ADMB Documentation  11.1.2270
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2f15.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2f15.cpp 1919 2014-04-22 22:02:01Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #include <df1b2fun.h>
00012 
00017 void df1b2_init_vector::allocate(int mmin,int mmax,const adstring&s)
00018 {
00019   allocate(mmin,mmax,1,s);
00020 }
00021 
00026 void df1b2_init_vector::allocate(int mmin,int mmax,int ps,const adstring&s)
00027 {
00028   set_phase_start(ps);
00029   df1b2vector::allocate(mmin,mmax);
00030 }
00031 
00036 void initial_df1b2params::add_to_list(void)
00037 {
00038   if (num_initial_df1b2params>=1000)
00039   {
00040     cerr << " This version of ADMB only supports 1000 initial df1b2parameter"
00041       " objects" << endl;
00042     ad_exit(1);
00043   }
00044   varsptr[num_initial_df1b2params++]= this;
00045 }
00046 
00051 void initial_df1b2params::reset(const init_df1b2vector& x,
00052   const df1b2variable& _pen)
00053 {
00054   ADUNCONST(df1b2variable,pen)
00055   //df1b2variable::nvar=global_nvar;
00056   pen=0.0;
00057   df1b2variable pen1=0.0;
00058   int ii=1;
00059   for (int i=0;i<num_initial_df1b2params;i++)
00060   {
00061     if (withinbound(0,(varsptr[i])->phase_start,current_phase))
00062     {
00063       (varsptr[i])->set_value(x,ii,pen1);
00064     }
00065   }
00066   pen+=pen1;
00067 }
00068 
00069 void set_value(const df1b2vector& _x,const init_df1b2vector& _v, const int& _ii,
00070   double fmin,double fmax,const df1b2variable& fpen,double s);
00075  double initial_df1b2params::get_scalefactor(void)
00076  {
00077    return scalefactor;
00078  }
00079 
00085 void initial_df1b2params::set_scalefactor(const double sf)
00086 {
00087   scalefactor=sf;
00088 }
00089 
00090 int initial_df1b2params::set_index(void)
00091 {
00092   int ii=1;
00093   for (int i=0;i<num_initial_df1b2params;i++)
00094   {
00095     if (withinbound(0,(varsptr[i])->phase_start,current_phase))
00096     {
00097       (varsptr[i])->set_index(1,ii);
00098     }
00099     else
00100     {
00101       (varsptr[i])->set_index(0,ii);
00102     }
00103   }
00104   return ii-1;
00105 }
00106 
00111 void df1b2_init_vector::set_value(const init_df1b2vector& _x,
00112   const int& _ii,const df1b2variable& pen)
00113 {
00114   ADUNCONST(init_df1b2vector,x)
00115   ADUNCONST(int,ii)
00116 
00117   int mmin=indexmin();
00118   int mmax=indexmax();
00119   if (scalefactor==0.0)
00120   {
00121     for (int i=mmin;i<=mmax;i++)
00122     {
00123       (*this)(i) = (x(ii++));
00124     }
00125   }
00126   else
00127   {
00128     for (int i=mmin;i<=mmax;i++)
00129     {
00130       (*this)(i) = (x(ii++)/scalefactor);
00131     }
00132   }
00133 }
00134 
00139 void df1b2_init_vector::set_index(int aflag,const int& _ii)
00140 {
00141   ADUNCONST(int,ii)
00142 
00143   int mmin=indexmin();
00144   int mmax=indexmax();
00145   if (aflag)
00146   {
00147     for (int i=mmin;i<=mmax;i++)
00148     {
00149       (*this)(i).get_ind_index() = ii++;
00150     }
00151   }
00152   else
00153   {
00154     for (int i=mmin;i<=mmax;i++)
00155     {
00156       (*this)(i).get_ind_index() = -1;
00157     }
00158   }
00159 }
00160 
00165 void df1b2_init_number::set_index(int aflag,const int& _ii)
00166 {
00167   ADUNCONST(int,ii)
00168   if (aflag)
00169   {
00170     get_ind_index()=ii++;
00171   }
00172   else
00173   {
00174     get_ind_index()=-1;
00175   }
00176 }
00177 
00182 void df1b2_init_matrix::set_index(int aflag,const int& _ii)
00183 {
00184   ADUNCONST(int,ii)
00185 
00186   int rmin=indexmin();
00187   int rmax=indexmax();
00188   if (aflag)
00189   {
00190     for (int i=rmin;i<=rmax;i++)
00191     {
00192       int cmin=(*this)(i).indexmin();
00193       int cmax=(*this)(i).indexmax();
00194       {
00195         for (int j=cmin;j<=cmax;j++)
00196         {
00197           (*this)(i,j).get_ind_index() = ii++;
00198         }
00199       }
00200     }
00201   }
00202   else
00203   {
00204     for (int i=rmin;i<=rmax;i++)
00205     {
00206       int cmin=(*this)(i).indexmin();
00207       int cmax=(*this)(i).indexmax();
00208       {
00209         for (int j=cmin;j<=cmax;j++)
00210         {
00211           (*this)(i,j).get_ind_index() = -1;
00212         }
00213       }
00214     }
00215   }
00216 }
00217 
00222 void df1b2_init_matrix::set_value(const init_df1b2vector& _x,
00223   const int& _ii,const df1b2variable& pen)
00224 {
00225   ADUNCONST(init_df1b2vector,x)
00226   ADUNCONST(int,ii)
00227 
00228   int rmin=indexmin();
00229   int rmax=indexmax();
00230   for (int i=rmin;i<=rmax;i++)
00231   {
00232     int cmin=(*this)(i).indexmin();
00233     int cmax=(*this)(i).indexmax();
00234     {
00235       for (int j=cmin;j<=cmax;j++)
00236       {
00237         (*this)(i,j) = (x(ii++));
00238       }
00239     }
00240   }
00241 }
00242 
00247 void df1b2_init_number::set_value(const init_df1b2vector& _x,const int& _ii,
00248     const df1b2variable& pen)
00249 {
00250   ADUNCONST(init_df1b2vector,x)
00251   ADUNCONST(int,ii)
00252   if (scalefactor==0.0)
00253     operator = (x(ii++));
00254   else
00255     operator = (x(ii++)/scalefactor);
00256 }
00257 
00262 void df1b2_init_number::allocate(const adstring& s)
00263 {
00264   if (!df1b2variable::noallocate)
00265   {
00266     df1b2variable::allocate();
00267   }
00268   set_phase_start(1);
00269 }
00270 
00275 void df1b2_init_number::allocate(int _n,const adstring& s)
00276 {
00277   df1b2variable::allocate();
00278   set_phase_start(_n);
00279 }
00280 
00285 void df1b2variable::allocate(const char *s)
00286 {
00287   df1b2variable::allocate();
00288 }
00289 
00294 df1b2_init_number::df1b2_init_number() :  df1b2variable(do_naught_kludge_a)
00295 {
00296   //cout << "Here" << endl;
00297 }
00298 
00303 void df1b2_init_number::operator = (const df1b2variable& _x)
00304 {
00305   df1b2variable::operator = (_x);
00306 }
00307 
00312 void df1b2_init_bounded_number::allocate(double _minb,double _maxb,
00313  int _n,const char * s)
00314 {
00315   minb=_minb;
00316   maxb=_maxb;
00317   df1b2_init_number::allocate(s);
00318 
00319   set_phase_start(_n);
00320 }
00321 
00326 void df1b2_init_bounded_number::allocate(double _minb,
00327   double _maxb,const char * s)
00328 {
00329   minb=_minb;
00330   maxb=_maxb;
00331   df1b2_init_number::allocate(s);
00332   set_phase_start(1);
00333 }
00334 
00339 void df1b2_init_bounded_number::set_value(const init_df1b2vector& _x,
00340   const int& _ii,const df1b2variable& pen)
00341 {
00342   ADUNCONST(init_df1b2vector,x)
00343   ADUNCONST(int,ii)
00344   if (scalefactor==0.0)
00345     ::set_value(*this,x,ii,minb,maxb,pen);
00346   else
00347     ::set_value(*this,x,ii,minb,maxb,pen,scalefactor);
00348 }
00349 
00354 void set_value(const df1b2variable& _u,const init_df1b2vector& _x,
00355   const int& _ii,double fmin,double fmax,const df1b2variable& _fpen)
00356 {
00357   ADUNCONST(init_df1b2vector,x)
00358   ADUNCONST(int,ii)
00359   ADUNCONST(df1b2variable,u)
00360   ADUNCONST(df1b2variable,fpen)
00361   if (!initial_params::straight_through_flag)
00362   {
00363     u=boundp(x(ii++),fmin,fmax,fpen);
00364   }
00365   else
00366   {
00367     u=x(ii);
00368     value(u)=boundp(value(x(ii++)),fmin,fmax);
00369     double diff=fmax-fmin;
00370     //t=fmin + diff*ss;
00371     df1b2variable ss=(u-fmin)/diff;
00372 #   ifdef USE_BARD_PEN
00373       const double l4=log(4.0);
00374       double pen=.000001/diff;
00375       fpen-=pen*(log(ss+1.e-40)+log((1.0-ss)+1.e-40)+l4);
00376 #   else
00377           XXXX
00378 #   endif
00379   }
00380 }
00381 
00382 void set_value(const df1b2variable& _u,const init_df1b2vector& _x,
00383   const int& _ii,double fmin,double fmax,const df1b2variable& _fpen,
00384   double s)
00385 {
00386   ADUNCONST(init_df1b2vector,x)
00387   ADUNCONST(int,ii)
00388   ADUNCONST(df1b2variable,u)
00389   ADUNCONST(df1b2variable,fpen)
00390   if (!initial_params::straight_through_flag)
00391   {
00392     u=boundp(x(ii++),fmin,fmax,fpen,s);
00393   }
00394   else
00395   {
00396     u=x(ii);
00397     value(u)=boundp(value(x(ii++)),fmin,fmax,s);
00398     double diff=fmax-fmin;
00399     //t=fmin + diff*ss;
00400     df1b2variable ss=(u-fmin)/diff;
00401 #   ifdef USE_BARD_PEN
00402       const double l4=log(4.0);
00403       double pen=.000001/diff;
00404       fpen-=pen*(log(ss+1.e-40)+log((1.0-ss)+1.e-40)+l4);
00405 #   else
00406           XXXX
00407 #   endif
00408   }
00409 }
00410 
00418 df1b2variable boundp(const df1b2variable& x, double fmin, double fmax,
00419   const df1b2variable& _fpen)
00420 {
00421   ADUNCONST(df1b2variable,fpen)
00422   df1b2variable t;
00423   //df1b2variable y;
00424   //y=x;
00425   double diff=fmax-fmin;
00426   const double l4=log(4.0);
00427   df1b2variable ss=0.49999999999999999*sin(x*1.57079632679489661)+0.50;
00428   t=fmin + diff*ss;
00429 
00430 #ifdef USE_BARD_PEN
00431   double pen=.000001/diff;
00432   fpen-=pen*(log(ss+1.e-40)+log((1.0-ss)+1.e-40)+l4);
00433 #else
00434   if (x < -.9999)
00435   {
00436     fpen+=cube(-0.9999-x);
00437     if (x < -1.)
00438     {
00439       fpen+=1.e+6*cube(-1.0-x);
00440       if (x < -1.02)
00441       {
00442         fpen+=1.e+10*cube(-1.02-x);
00443       }
00444     }
00445   }
00446   if (x > 0.9999)
00447   {
00448     fpen+=cube(x-0.9999);
00449     if (x > 1.)
00450     {
00451       fpen+=1.e+6*cube(x-1.);
00452       if (x > 1.02)
00453       {
00454         fpen+=1.e+10*cube(x-1.02);
00455       }
00456     }
00457   }
00458 #endif
00459   return(t);
00460 }
00461 
00469 df1b2variable boundp(const df1b2variable& _x, double fmin, double fmax,
00470   const df1b2variable& _fpen,double s)
00471 {
00472   ADUNCONST(df1b2variable,fpen)
00473   df1b2variable t;
00474   df1b2variable x=_x/s;
00475   //df1b2variable y;
00476   //y=x;
00477   double diff=fmax-fmin;
00478   const double l4=log(4.0);
00479 
00480   // ss is underlying varialbe on [0,1] and t lives in [fmin,fmax]
00481   df1b2variable ss=0.49999999999999999*sin(x*1.57079632679489661)+0.50;
00482   t=fmin + diff*ss;
00483 
00484   // Add penalty
00485 #ifdef USE_BARD_PEN
00486   double pen=.000001/diff;
00487   fpen-=pen*(log(ss+1.e-40)+log((1.0-ss)+1.e-40)+l4);
00488 #else
00489   if (x < -.9999)
00490   {
00491     fpen+=cube(-0.9999-x);
00492     if (x < -1.)
00493     {
00494       fpen+=1.e+6*cube(-1.0-x);
00495       if (x < -1.02)
00496       {
00497         fpen+=1.e+10*cube(-1.02-x);
00498       }
00499     }
00500   }
00501   if (x > 0.9999)
00502   {
00503     fpen+=cube(x-0.9999);
00504     if (x > 1.)
00505     {
00506       fpen+=1.e+6*cube(x-1.);
00507       if (x > 1.02)
00508       {
00509         fpen+=1.e+10*cube(x-1.02);
00510       }
00511     }
00512   }
00513 #endif
00514   return(t);
00515 }
00516 
00521 df1b2variable boundp(const df1b2variable& x, double fmin, double fmax)
00522 {
00523   df1b2variable t;
00524   //df1b2variable y;
00525   //y=x;
00526   double diff=fmax-fmin;
00527   df1b2variable ss=0.49999999999999999*sin(x*1.57079632679489661)+0.50;
00528   t=fmin + diff*ss;
00529 
00530   return(t);
00531 }
00532 
00537 void set_value_mc(const df1b2vector& _x,const init_df1b2vector& _v,
00538   const int& _ii, double fmin,double fmax)
00539 {
00540   ADUNCONST(int,ii)
00541   ADUNCONST(df1b2vector,x)
00542   ADUNCONST(init_df1b2vector,v)
00543   int min=x.indexmin();
00544   int max=x.indexmax();
00545   for (int i=min;i<=max;i++)
00546   {
00547     //x(i)=boundp(v(ii++),fmin,fmax,fpen);
00548     const double pinv=1./PI;
00549     df1b2variable y=atan(v(ii++))*pinv+0.5;
00550     x(i)=fmin+(fmax-fmin)*y;
00551   }
00552 }
00553 
00558 void set_value(const df1b2vector& _x,const init_df1b2vector& _v,
00559   const int& _ii, double fmin,double fmax,const df1b2variable& fpen)
00560 {
00561   ADUNCONST(int,ii)
00562   ADUNCONST(df1b2vector,x)
00563   ADUNCONST(init_df1b2vector,v)
00564   int min=x.indexmin();
00565   int max=x.indexmax();
00566   for (int i=min;i<=max;i++)
00567   {
00568     x(i)=boundp(v(ii++),fmin,fmax,fpen);
00569     //cout << setprecision(15) << fpen << " " << fmin << " " << fmax
00570      //  << " " << v(ii-1) << " " << x(i) << endl;
00571   }
00572 }
00573 void set_value(const df1b2vector& _x,const init_df1b2vector& _v,
00574   const int& _ii, double fmin,double fmax,const df1b2variable& fpen,double s)
00575 {
00576   ADUNCONST(int,ii)
00577   ADUNCONST(df1b2vector,x)
00578   ADUNCONST(init_df1b2vector,v)
00579   int min=x.indexmin();
00580   int max=x.indexmax();
00581   for (int i=min;i<=max;i++)
00582   {
00583     x(i)=boundp(v(ii++),fmin,fmax,fpen,s);
00584     //cout << setprecision(15) << fpen << " " << fmin << " " << fmax
00585      //  << " " << v(ii-1) << " " << x(i) << endl;
00586   }
00587 }
00588 
00593 void df1b2_init_bounded_vector::set_value(const init_df1b2vector& _x,
00594   const int& ii,const df1b2variable& pen)
00595 {
00596   init_df1b2vector& x=(init_df1b2vector&) _x;
00597 
00598   if (scalefactor==0.0)
00599   {
00600     if (initial_params::mc_phase)
00601     {
00602       ::set_value_mc(*this,x,ii,minb,maxb);
00603     }
00604     else
00605     {
00606       ::set_value(*this,x,ii,minb,maxb,pen);
00607     }
00608   }
00609   else
00610   {
00611     if (initial_params::mc_phase)
00612     {
00613       // should this have scalefactor ??
00614       //::set_value_mc(*this,x,ii,minb,maxb,scalefactor);
00615       ::set_value_mc(*this,x,ii,minb,maxb);
00616     }
00617     else
00618     {
00619       ::set_value(*this,x,ii,minb,maxb,pen,scalefactor);
00620       //::set_value(*this,x,ii,minb,maxb,pen);
00621     }
00622   }
00623 }
00624 
00629 re_df1b2_init_bounded_vector::re_df1b2_init_bounded_vector(void)
00630 {
00631   initial_df1b2params::have_bounded_random_effects=1;
00632 }
00633 
00638 void re_df1b2_init_bounded_vector::set_value(const dvector& _x,
00639   const int& _ii)
00640 {
00641   df1b2_init_bounded_vector::set_value(_x,_ii);
00642 }
00643 
00648 void re_df1b2_init_bounded_vector::set_value(const init_df1b2vector& _x,
00649   const int& _ii,const df1b2variable& _pen)
00650 {
00651   ADUNCONST(int,ii)
00652   ADUNCONST(df1b2variable,pen)
00653   init_df1b2vector& x=(init_df1b2vector&) _x;
00654   int mmin=indexmin();
00655   int mmax=indexmax();
00656   for (int i=mmin;i<=mmax;i++)
00657   {
00658     df1b2variable& y=df1b2vector::operator()(i);
00659     if (!initial_params::straight_through_flag)
00660     {
00661       y = (boundp(x(ii++),getminb(),getmaxb(),pen));
00662     }
00663     else
00664     {
00665       y = (x(ii++));
00666       *y.get_u()=boundp(*y.get_u(),getminb(),getmaxb());
00667       double diff=getmaxb()-getminb();
00668       df1b2variable ss=(y-getminb())/diff;
00669 #     ifdef USE_BARD_PEN
00670         const double l4=log(4.0);
00671         double wght=.000001/diff;
00672         pen-=wght*(log(ss+1.e-40)+log((1.0-ss)+1.e-40)+l4);
00673 #     else
00674        XXXX
00675 #     endif
00676     }
00677   }
00678 }
00679 
00684 void df1b2_init_bounded_dev_vector::set_value(const init_df1b2vector& x,
00685   const int& ii,const df1b2variable& _pen)
00686 {
00687   ADUNCONST(df1b2variable,pen);
00688   df1b2_init_bounded_vector::set_value(x,ii,pen);
00689   df1b2variable m=mean(*this);
00690   pen+=1000.0*square(m);
00691 }
00692 
00697 void df1b2_init_bounded_vector::allocate(int mmin,int mmax,
00698   double _minb,double _maxb)
00699 {
00700   minb=_minb;
00701   maxb=_maxb;
00702   df1b2_init_vector::allocate(mmin,mmax);
00703 }
00704 
00709 void df1b2_init_bounded_vector::allocate(int mmin,int mmax,
00710   double _minb,double _maxb,const char * s)
00711 {
00712   allocate(mmin,mmax,_minb,_maxb);
00713 }
00714 
00719 void df1b2_init_bounded_vector::allocate(int mmin,int mmax,
00720   double _minb,double _maxb,int n,const char * s)
00721 {
00722   allocate(mmin,mmax,_minb,_maxb);
00723   set_phase_start(n);
00724 }
00725 
00730 int active(const df1b2_init_vector& x)
00731 {
00732    if (x.current_phase >= x.phase_start)
00733      return 1;
00734    else
00735      return 0;
00736 }
00737 
00742 int active(const df1b2_init_number& x)
00743 {
00744    if (x.current_phase >= x.phase_start)
00745      return 1;
00746    else
00747      return 0;
00748 }
00749 
00754 int active(const df1b2_init_matrix& x)
00755 {
00756    if (x.current_phase >= x.phase_start)
00757      return 1;
00758    else
00759      return 0;
00760 }
00761 
00766 void random_effects_bounded_vector_info::set_value
00767   (const init_df1b2vector& _x,const int& _ii,const df1b2variable& _pen)
00768 {
00769   ADUNCONST(int,ii)
00770   ADUNCONST(init_df1b2vector,x)
00771   ADUNCONST(df1b2variable,pen);
00772   df1b2variable& y=pv->df1b2vector::operator()(i);
00773   if (!initial_params::straight_through_flag)
00774   {
00775     // df1b2variable& tmp = boundp(x(ii++),b.getminb(),b.getmaxb(),pen);
00776     // df1b2variable::operator = (tmp);
00777     //df1b2variable::operator =
00778     y = (boundp(x(ii++),pv->getminb(),pv->getmaxb(),pen));
00779   }
00780   else
00781   {
00782     y = (x(ii++));
00783     *y.get_u()=boundp(*y.get_u(),pv->getminb(),pv->getmaxb());
00784     double diff=pv->getmaxb()-pv->getminb();
00785     df1b2variable ss=(y-pv->getminb())/diff;
00786 #   ifdef USE_BARD_PEN
00787       const double l4=log(4.0);
00788       double wght=.000001/diff;
00789       pen-=wght*(log(ss+1.e-40)+log((1.0-ss)+1.e-40)+l4);
00790 #   else
00791      XXXX
00792 #   endif
00793   }
00794 }