ADMB Documentation  11.1.1913
 All Classes Files Functions Variables Typedefs Friends Defines
f1b2locl.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: f1b2locl.cpp 1709 2014-02-28 21:48:21Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #include <df1b2fun.h>
00012 #include "f1b2locl.h"
00013 class newadkludge;
00014 extern newadkludge * newadkl;
00015 
00020 void local_init_df1b2variable::dot_calcs(local_dep_df1b2variable& v,int j)
00021 {
00022   if (adpool_stack_pointer<1)
00023   {
00024     cerr << "this can't happen" << endl;
00025     ad_exit(1);
00026   }
00027   int global_nvar=adpool_nvar_stack[adpool_stack_pointer-1];
00028   double * xd=p->get_u_dot();
00029   double * ud=v.get_u_dot();
00030   double * zd=(v.p)->get_u_dot();
00031   for (int i=0;i<global_nvar;i++)
00032   {
00033     //for (int j=0;j<nvar;j++)
00034     {
00035       *zd+=ud[j] * *xd;
00036     }
00037     zd++;
00038     xd++;
00039   }
00040 }
00041 
00046 local_init_df1b2variable::local_init_df1b2variable
00047   (void) : df1b2variable(newadkl)
00048 {
00049   p=0;
00050 }
00051 
00056 local_init_df1b2variable::local_init_df1b2variable
00057   (double _xu, double * _xdot) : df1b2variable()
00058 {
00059   *get_u()=_xu;
00060   p=0;
00061   xdot=_xdot;
00062   set_dot(num_vars-1);
00063 }
00064 
00069 local_init_df1b2variable::local_init_df1b2variable
00070   (const df1b2variable & _x) : df1b2variable(newadkl)
00071 {
00072   ADUNCONST(df1b2variable,x)
00073   p=&_x;
00074   xu=*(x.get_u());
00075 }
00076 
00081 local_init_df1b2vector::local_init_df1b2vector(const df1b2vector & _x)
00082 {
00083   ADUNCONST(df1b2vector,x)
00084   p=&_x;
00085   int mmin=x.indexmin();
00086   int mmax=x.indexmax();
00087   //int ind_index = x(mmin).get_ind_index();
00088   df1b2variable::noallocate=1;
00089   df1b2vector::allocate(mmin,mmax);
00090   df1b2variable::noallocate=0;
00091 }
00092 
00093  // #define USE_BARD_PEN
00094  // class newadkludge;
00095  // extern newadkludge * newadkl=0;
00096  //
00097  //
00098   typedef local_init_var  * PLOCAL_INIT_VAR;
00099  //
00100  // class laplace_approximation_calculator;
00101  // laplace_approximation_calculator * funnel_init_var::lapprox=0;
00102  // df1b2variable * funnel_init_var::funnel_constraints_penalty=0;
00103   int local_init_var::num_vars=0;
00104  // //int funnel_init_var::num_all_vars=0;
00105   int local_init_var::num_inactive_vars=0;
00106   int local_init_var::num_active_parameters=0;
00107  // //funnel_init_var ** funnel_init_var::all_list=new PFUNNEL_INIT_VAR[200];
00108   local_init_var ** local_init_var::list=new PLOCAL_INIT_VAR[200];
00109   local_init_var ** local_init_var::inactive_list=new PLOCAL_INIT_VAR[200];
00110   init_df1b2vector * local_init_var::py=0;
00111   imatrix * local_init_var::plist=0;
00112  //
00113  // void  xxx(init_df1b2vector & tmp,int x){;}
00114  //
00115 
00120   void local_init_var::add_to_list(void)
00121   {
00122     index=num_vars;
00123     list[num_vars++]=this;
00124     //all_list[num_all_vars++]=this;
00125   }
00126  //
00127  // void funnel_init_var::delete_from_list(void)
00128  // {
00129  //   if (index!=num_vars-1)
00130  //   {
00131  //     cerr << "can only delete last member" << endl;
00132  //     ad_exit(1);
00133  //   }
00134  //   num_vars--;
00135  //   index=-1;
00136  // }
00137  //
00138 
00143   void local_init_var::add_to_inactive_list(void)
00144   {
00145     index=-1;
00146     inactive_list[num_inactive_vars++]=this;
00147   }
00148 
00153   void local_init_var::allocate(void)
00154   {
00155     //cout << "In allocate" << endl;
00156   }
00157 
00162   void local_init_var::allocate_all(void)
00163   {
00164     num_active_parameters=local_init_var::nvarcalc_all();
00165     if (py)
00166     {
00167       if (py->indexmax() != num_active_parameters)
00168       {
00169         delete py;
00170         py=0;
00171       }
00172     }
00173 
00174     df1b2variable::save_adpool_pointer();
00175     adpool * tmppool=df1b2variable::pool;
00176     if (!localf1b2gradlist)
00177     {
00178       localf1b2gradlist = new df1b2_gradlist(400000U,20000U,800000U,40000U,
00179         200000U,10000U,adstring("lf1b2list1"));
00180       if (!localf1b2gradlist)
00181       {
00182         cerr << "Error allocating memory for local df1b2gradlist" << endl;
00183         ad_exit(1);
00184       }
00185     }
00186     globalf1b2gradlist=f1b2gradlist;
00187     f1b2gradlist=localf1b2gradlist;
00188 
00189     if (tmppool)
00190     {
00191       //cout << tmppool << endl;
00192       // check if current pool is the right size
00193       if (tmppool->nvar != num_active_parameters)
00194       {
00195         // check sizes of other pools
00196         int found_pool_flag=0;
00197         for (int i=0;i<df1b2variable::adpool_counter;i++)
00198         {
00199           if (df1b2variable::adpool_vector[i]->nvar == num_active_parameters)
00200           {
00201             adpool * tmp = df1b2variable::pool;
00202             df1b2variable::pool=df1b2variable::adpool_vector[i];
00203             if (!tmp->on_adpool_vector())
00204             {
00205               df1b2variable::adpool_vector[df1b2variable::adpool_counter]=tmp;
00206               df1b2variable::nvar_vector[df1b2variable::adpool_counter]=
00207                 tmp->nvar;
00208               //df1b2variable::adpool_counter++;
00209               df1b2variable::increment_adpool_counter();
00210               tmp->on_adpool_vector()=1;
00211             }
00212             found_pool_flag=1;
00213             break;
00214           }
00215         }
00216         if (!found_pool_flag)
00217         {
00218           if (df1b2variable::adpool_counter>=df1b2variable::adpool_vectorsize)
00219           {
00220              cerr << "Need to increase adpool_vectorsize" << endl;
00221              ad_exit(1);
00222           }
00223           if (!df1b2variable::pool->on_adpool_vector())
00224           {
00225             df1b2variable::adpool_vector[df1b2variable::adpool_counter]=
00226               df1b2variable::pool;
00227             df1b2variable::nvar_vector[df1b2variable::adpool_counter]=
00228               df1b2variable::pool->nvar;
00229             //df1b2variable::adpool_counter++;
00230             df1b2variable::increment_adpool_counter();
00231             df1b2variable::pool->on_adpool_vector()=1;
00232           }
00233           df1b2variable::pool=new adpool();
00234           if (!df1b2variable::pool)
00235           {
00236             cerr << "Memory allocation error" << endl;
00237             ad_exit(1);
00238           }
00239 
00240           df1b2variable::nvar=num_active_parameters;
00241           df1b2variable::set_blocksize();
00242 
00243           if (!df1b2variable::pool->on_adpool_vector())
00244           {
00245             df1b2variable::adpool_vector[df1b2variable::adpool_counter]=
00246               df1b2variable::pool;
00247             df1b2variable::nvar_vector[df1b2variable::adpool_counter]=
00248               df1b2variable::pool->nvar;
00249             //df1b2variable::adpool_counter++;
00250             df1b2variable::increment_adpool_counter();
00251             df1b2variable::pool->on_adpool_vector()=1;
00252           }
00253         }
00254       }
00255     }
00256     else
00257     {
00258       df1b2variable::pool=new adpool();
00259       if (!df1b2variable::pool)
00260       {
00261         cerr << "Memory allocation error" << endl;
00262         ad_exit(1);
00263       }
00264       df1b2variable::nvar=num_active_parameters;
00265       df1b2variable::set_blocksize();
00266       df1b2variable::adpool_vector[df1b2variable::adpool_counter]=
00267         df1b2variable::pool;
00268       df1b2variable::nvar_vector[df1b2variable::adpool_counter]=
00269         df1b2variable::pool->nvar;
00270       //df1b2variable::adpool_counter++;
00271       df1b2variable::increment_adpool_counter();
00272     }
00273 
00274       df1b2variable::nvar=num_active_parameters;
00275       df1b2variable::set_blocksize();
00276 
00277       if (!py)
00278       {
00279         py = new init_df1b2vector(1,num_active_parameters);
00280       }
00281       if (!py)
00282       {
00283         cerr << "memory allocation error" << endl;
00284         ad_exit(1);
00285       }
00286       //init_df1b2vector& tmp = *py;
00287 
00288     if (plist)
00289     {
00290       if (plist->indexmax() != num_active_parameters)
00291       {
00292         delete plist;
00293         plist=0;
00294       }
00295     }
00296     if (!plist)
00297     {
00298       plist = new imatrix(1,num_active_parameters,1,2);
00299     }
00300     if (!plist)
00301     {
00302       cerr << "memory allocation error" << endl;
00303       ad_exit(1);
00304     }
00305 
00306     int ii=1;
00307     int i;
00308       for(i=0;i<num_vars;i++)
00309       {
00310         list[i]->xinit(*py,ii);
00311       }
00312 
00313     ii=1;
00314     for(i=0;i<num_vars;i++)
00315     {
00316       list[i]->set_index(*plist,ii);
00317     }
00318 
00319     for(i=0;i<num_inactive_vars;i++)
00320     {
00321       inactive_list[i]->allocate();
00322     }
00323 
00324     local_init_var::reset(*py);
00325   }
00326 
00327  // funnel_init_df1b2variable::funnel_init_df1b2variable
00328  //   (const df1b2_init_number & _x) : df1b2variable(newadkl)
00329  //   //(df1b2_init_number & x) : df1b2variable()
00330  // {
00331  //   ADUNCONST(df1b2_init_number,x)
00332  //   type=0;
00333  //   pointer=0;
00334  //   ind_index=x.get_ind_index();
00335  //   if (ind_index<0)
00336  //   {
00337  //     add_to_inactive_list();
00338  //   }
00339  //   else
00340  //   {
00341  //     add_to_list();
00342  //     lapprox->used_flags(ind_index)+=1;
00343  //   }
00344  //   //cout << "ind_index = " << ind_index << endl;
00345  //   xu=*(x.get_u());
00346  // }
00347  //
00348  //
00349  // funnel_init_df1b2variable::funnel_init_df1b2variable
00350  //   (const random_effects_bounded_vector_info & _u)
00351  //   : df1b2variable(newadkl)
00352  // {
00353  //   ADUNCONST(random_effects_bounded_vector_info,u)
00354  //   df1b2variable& x = (*(u.pv)).df1b2vector::operator () (u.i);
00355  //
00356  //   type=1;
00357  //   pointer=u.pv;
00358  //   ind_index = x.get_ind_index();
00359  //   if (ind_index<0)
00360  //   {
00361  //     add_to_inactive_list();
00362  //   }
00363  //   else
00364  //   {
00365  //     add_to_list();
00366  //     lapprox->used_flags(ind_index)+=1;
00367  //   }
00368  //   xu=*(x.get_u());
00369  //
00370  // }
00371  //
00372 
00377   void local_init_df1b2variable::allocate(const df1b2variable& x)
00378   {
00379     cerr << "Haven't defined htis yet" << endl;
00380     ad_exit(1);
00381   }
00382  //
00383  // funnel_init_df1b2variable::funnel_init_df1b2variable
00384  //   (void) : df1b2variable(newadkl)
00385  // {
00386  //   type=0;
00387  //   pointer=0;
00388  //   ind_index = -1;
00389  //   if (ind_index<0)
00390  //   {
00391  //     add_to_inactive_list();
00392  //   }
00393  //   else
00394  //   {
00395  //     add_to_list();
00396  //   }
00397  // }
00398  //
00399 
00404   void local_init_df1b2variable::
00405     preallocate(const df1b2variable & _x)
00406   {
00407     ADUNCONST(df1b2variable,x)
00408     type=0;
00409     pointer=0;
00410     ind_index = x.get_ind_index();
00411     if (ind_index<0)
00412     {
00413       add_to_inactive_list();
00414     }
00415     else
00416     {
00417       add_to_list();
00418     }
00419     xu=*(x.get_u());
00420   }
00421 
00422  //
00423  // funnel_init_df1b2variable::funnel_init_df1b2variable
00424  //   (const df1b2variable & _x) : df1b2variable(newadkl)
00425  // {
00426  //   ADUNCONST(df1b2variable,x)
00427  //   type=0;
00428  //   pointer=0;
00429  //   ind_index = x.get_ind_index();
00430  //   if (ind_index<0)
00431  //   {
00432  //     add_to_inactive_list();
00433  //   }
00434  //   else
00435  //   {
00436  //     add_to_list();
00437  //   }
00438  //   xu=*(x.get_u());
00439  // }
00440  //
00441 
00446   void local_init_df1b2variable::allocate(void)
00447   {
00448     df1b2variable::allocate();
00449     *(get_u())=xu;
00450     if (index>=0)
00451       get_u_dot()[index]=1.0;
00452   }
00453  //
00454  // funnel_dependent_df1b2variable::funnel_dependent_df1b2variable
00455  //   (const df1b2variable& x)
00456  // {
00457  //   df1b2variable::operator = (x);
00458  //   if (!df1b2_gradlist::no_derivatives)
00459  //   {
00460  //     df1b2variable * tmp = (df1b2variable *) (this);
00461  //     //set_dependent_variable(*tmp);
00462  //   }
00463  //   df1b2_gradlist::set_no_derivatives();
00464  // }
00465  //
00466 
00471   void local_init_df1b2variable::set_value(const init_df1b2vector& _x,
00472     const int& _ii)
00473   {
00474     //df1b2variable pen=0.0;
00475     ADUNCONST(init_df1b2vector,x)
00476     ADUNCONST(int,ii)
00477     df1b2variable::operator = (x(ii++));
00478   }
00479 
00484   void local_init_df1b2variable::set_value(const init_df1b2vector& _x,
00485     const int& _ii,const df1b2variable& _pen)
00486   {
00487     ADUNCONST(init_df1b2vector,x)
00488     ADUNCONST(int,ii)
00489     ADUNCONST(df1b2variable,pen)
00490     if (!pointer)
00491     {
00492       df1b2variable::operator = (x(ii++));
00493     }
00494     else
00495     {
00496       switch (type)
00497       {
00498         case 1:   // vector
00499         {
00500           df1b2_init_bounded_vector & b = *(df1b2_init_bounded_vector*)pointer;
00501           //laplace_approximation_calculator * l =lapprox;
00502           //int uf=-1;
00503           if (ind_index>0)
00504           {
00505            // uf=l->used_flags(ind_index);
00506           }
00507           else
00508           {
00509             if (!initial_params::straight_through_flag)
00510             {
00511               df1b2variable::operator = (boundp(x(ii++),b.getminb(),b.getmaxb(),
00512                 pen));
00513             }
00514             else
00515             {
00516               df1b2variable::operator = (x(ii++));
00517               *get_u()=boundp(*get_u(),b.getminb(),b.getmaxb());
00518               double diff=b.getmaxb()-b.getminb();
00519               df1b2variable ss=((*this)-b.getminb())/diff;
00520   #           ifdef USE_BARD_PEN
00521                 const double l4=log(4.0);
00522                 double wght=.000001/diff;
00523                 pen-=wght*(log(ss+double(1.e-40))+log((double(1.0)-ss)
00524                    +double(1.e-40))+l4);
00525   #           else
00526                XXXX
00527   #           endif
00528             }
00529           }
00530           break;
00531         }
00532         case 2:  // matrix
00533         default:
00534         {
00535           cerr << "the bounded matrix case in "
00536             " void funnel_init_df1b2variable::xinit  has not bee implemented"
00537             << endl;
00538             ad_exit(1);
00539         }
00540       }
00541     }
00542   }
00543 
00548   int local_init_var::nvarcalc_all(void)
00549   {
00550     int n=0;
00551     for (int i=0;i<num_vars;i++)
00552     {
00553       n+=list[i]->nvar_calc();
00554     }
00555     return n;
00556   }
00557 
00562   void local_init_var::set_dot_all(void)
00563   {
00564     int ii=0;
00565     for (int i=0;i<num_vars;i++)
00566     {
00567       list[i]->set_dot(ii);
00568     }
00569   }
00570 
00575   void local_init_var::dot_calcs_all(local_dep_df1b2variable& v)
00576   {
00577     for (int i=0;i<num_vars;i++)
00578     {
00579       list[i]->dot_calcs(v,i);
00580     }
00581   }
00582 
00587   void local_init_df1b2variable::xinit(init_df1b2vector& y,int& ii)
00588   {
00589     y(ii)= xu;
00590     ii++;
00591   }
00592 
00597   void local_init_df1b2variable::xinit(dvector& y,int& ii)
00598   {
00599     if (!pointer)
00600     {
00601       y(ii)= xu;
00602     }
00603     else
00604     {
00605       switch (type)
00606       {
00607         case 1:   // vector
00608         {
00609           df1b2_init_bounded_vector & b = *(df1b2_init_bounded_vector*)pointer;
00610           y(ii)=boundpin(xu,b.getminb(),b.getmaxb());
00611           break;
00612         }
00613         case 2:  // matrix
00614         default:
00615         {
00616           cerr << "the bounded matrix case in "
00617             " void local_init_df1b2variable::xinit  has not bee implemented"
00618             << endl;
00619             ad_exit(1);
00620         }
00621       }
00622     }
00623     ii++;
00624   }
00625 
00626  // /*
00627  // void funnel_init_df1b2variable::xinit(dvector& y,int& ii)
00628  // {
00629  //   y(ii)= xu;
00630  //   ii++;
00631  // }
00632  // */
00633  //
00634 
00639   void local_init_df1b2variable::set_index(imatrix& y,int& ii)
00640   {
00641     //cout << "FUCK " << ind_index << " " << ii << endl;
00642     y(ii,1)= ind_index;
00643     y(ii,2)= ii;
00644     ii++;
00645   }
00646 
00651   void local_init_var::reset(init_df1b2vector& x)
00652   {
00653     int ii=1;
00654     //df1b2variable pen=0.0;
00655     for (int i=0;i<num_vars;i++)
00656     {
00657       //list[i]->set_value(x,ii,pen);
00658       list[i]->set_value(x,ii);
00659     }
00660   }
00661  //
00662  //
00663  // funnel_init_df1b2vector::funnel_init_df1b2vector(
00664  //   const df1b2_init_vector& _x)
00665  // {
00666  //   ADUNCONST(df1b2_init_vector,x)
00667  //   //type=0;
00668  //   //pointer=0;
00669  //   p=&_x;
00670  //   int mmin=p->indexmin();
00671  //   int mmax=p->indexmax();
00672  //   int ind_index = (*p)(mmin).get_ind_index();
00673  //   if (ind_index<0)
00674  //   {
00675  //     add_to_inactive_list();
00676  //   }
00677  //   else
00678  //   {
00679  //     add_to_list();
00680  //   }
00681  //   df1b2variable::noallocate=1;
00682  //   df1b2vector::allocate(mmin,mmax);
00683  //   df1b2variable::noallocate=0;
00684  // }
00685  //
00686  //
00687 
00692   int local_init_df1b2vector::nvar_calc(void)
00693   {
00694     return p->indexmax()-p->indexmin()+1;
00695   }
00696 
00697 
00702   void local_init_df1b2vector::xinit(init_df1b2vector& y,int& ii)
00703   {
00704     df1b2vector * vp = (df1b2vector*)(p);
00705     int mmin=vp->indexmin();
00706     int mmax=vp->indexmax();
00707     int i;
00708     for (i=mmin;i<=mmax;i++)
00709     {
00710       y(ii)= value((*vp)(i));
00711       ii++;
00712     }
00713   }
00714 
00719   void local_init_df1b2vector::set_index(imatrix& y,int& ii)
00720   {
00721     int mmin=p->indexmin();
00722     int mmax=p->indexmax();
00723     int i;
00724     for (i=mmin;i<=mmax;i++)
00725     {
00726       //y(ii,1)= ( *(df1b2_init_vector *)(p) )(i).get_ind_index();
00727       y(ii,2)= ii;
00728       ii++;
00729     }
00730   }
00731 
00736   void local_init_df1b2vector::set_value(const init_df1b2vector& _x,
00737     const int& _ii)
00738   {
00739     ADUNCONST(int,ii)
00740     ADUNCONST(init_df1b2vector,x)
00741     int mmin=p->indexmin();
00742     int mmax=p->indexmax();
00743     int i;
00744     for (i=mmin;i<=mmax;i++)
00745     {
00746       (*this)(i) = (x(ii++));
00747     }
00748   }
00749 
00754   void local_init_df1b2vector::set_value(const init_df1b2vector& _x,
00755     const int& _ii,const df1b2variable& pen)
00756   {
00757     ADUNCONST(int,ii)
00758     ADUNCONST(init_df1b2vector,x)
00759     int mmin=p->indexmin();
00760     int mmax=p->indexmax();
00761     int i;
00762     for (i=mmin;i<=mmax;i++)
00763     {
00764       (*this)(i) = (x(ii++));
00765     }
00766   }