ADMB Documentation  11.4.2891
 All Classes Files Functions Variables Typedefs Friends Defines
montebds.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id$
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00007 #include <admodel.h>
00008 
00009 /*
00010 ofstream ofss("variance");
00011 double ssmin(double x, double y)
00012 {
00013   if (x<y) return x;
00014   return y;
00015 }
00016 double ssmax(double y, double x)
00017 {
00018   if (x>y) return x;
00019   return y;
00020 }
00021 */
00022 
00023 const double simbdsmax=1.e+90;
00024 
00025   void initial_params::set_all_simulation_bounds(const dmatrix& symbds)
00026   {
00027     int ii=1;
00028     for (int i=0;i<num_initial_params;i++)
00029     {
00030       //if ((varsptr[i])->phase_start <= current_phase)
00031       if (withinbound(0,(varsptr[i])->phase_start,current_phase))
00032         (varsptr[i])->set_simulation_bounds(symbds,ii);
00033     }
00034   }
00035 
00036 void param_init_number::set_simulation_bounds(const dmatrix& _symbds,
00037   const int& _ii)
00038   {
00039     int& ii=(int&) _ii;
00040     dmatrix& symbds=(dmatrix&) _symbds;
00041     symbds(1,ii)= -simbdsmax;
00042     symbds(2,ii)= simbdsmax;
00043     ii++;
00044   }
00045 
00046 void param_init_bounded_number::set_simulation_bounds(const dmatrix& _symbds,
00047   const int& _ii)
00048   {
00049     dmatrix& symbds=(dmatrix&) _symbds;
00050     int& ii=(int&) _ii;
00051     symbds(1,ii)= minb-value(*this);
00052     symbds(2,ii)= maxb-value(*this);
00053     ii++;
00054   }
00055 
00056 void param_init_vector::set_simulation_bounds(const dmatrix& _symbds,
00057   const int& _ii)
00058   {
00059     dmatrix& symbds=(dmatrix&) _symbds;
00060     int& ii=(int&) _ii;
00061     int mmin=indexmin();
00062     int mmax=indexmax();
00063     for (int i=mmin;i<=mmax;i++)
00064     {
00065       symbds(1,ii)= -simbdsmax;
00066       symbds(2,ii)= simbdsmax;
00067       ii++;
00068     }
00069   }
00070 
00071 void param_init_bounded_vector::set_simulation_bounds(const dmatrix& _symbds,
00072   const int& _ii)
00073   {
00074     dmatrix& symbds=(dmatrix&) _symbds;
00075     int& ii=(int&) _ii;
00076     int mmin=indexmin();
00077     int mmax=indexmax();
00078     for (int i=mmin;i<=mmax;i++)
00079     {
00080       symbds(1,ii)= minb-value((*this)(i));
00081       symbds(2,ii)= maxb-value((*this)(i));
00082       ii++;
00083     }
00084   }
00085 
00086 void param_init_matrix::set_simulation_bounds(const dmatrix& _symbds,
00087   const int& _ii)
00088   {
00089     dmatrix& symbds=(dmatrix&) _symbds;
00090     int& ii=(int&) _ii;
00091     int rmin=rowmin();
00092     int rmax=rowmax();
00093     for (int i=rmin;i<=rmax;i++)
00094     {
00095       int cmin=(*this)(i).indexmin();
00096       int cmax=(*this)(i).indexmax();
00097       for (int j=cmin;j<=cmax;j++)
00098       {
00099         symbds(1,ii)= -simbdsmax;
00100         symbds(2,ii)= simbdsmax;
00101         ii++;
00102       }
00103     }
00104   }
00105 
00106 void param_init_bounded_matrix::set_simulation_bounds(const dmatrix& _symbds,
00107   const int& _ii)
00108   {
00109     dmatrix& symbds=(dmatrix&) _symbds;
00110     int& ii=(int&) _ii;
00111     int rmin=rowmin();
00112     int rmax=rowmax();
00113     for (int i=rmin;i<=rmax;i++)
00114     {
00115       int cmin=(*this)(i).indexmin();
00116       int cmax=(*this)(i).indexmax();
00117       for (int j=cmin;j<=cmax;j++)
00118       {
00119         symbds(1,ii)= minb-value((*this)(i,j));
00120         symbds(2,ii)= maxb-value((*this)(i,j));
00121         ii++;
00122       }
00123     }
00124   }
00125 
00126 void param_init_number::add_value(const dvector& y, const dvector& ndev,
00127   const int& _ii, const double& s, const dvector& diag)
00128   {
00129     int& ii=(int&) _ii;
00130     (*this)+=diag(ii)*ndev(ii);
00131     ii++;
00132   }
00133 
00134 double new_value_mc(const double& _jac,double x,double min,double max,
00135   double eps, double sig)
00136 {
00137   double& jac=(double&) _jac;
00138   double y;
00139   double y1;
00140   double mult;
00141   double z;
00142   double z1;
00143   if (eps>0)
00144   {
00145     double d=max-x;
00146     double d2=.5*d;
00147     if (sig>d2)
00148     {
00149       mult=d2;
00150     }
00151     else
00152     {
00153       mult=sig;
00154     }
00155     y=mult*eps;
00156     y1=mult*(eps+1.e-6);
00157     if (y>0.9*d)
00158     {
00159       z=y-0.9*d;
00160       z1=y1-0.9*d;
00161       y=0.9*d+0.1*d*z/(1.+z);
00162       y1=0.9*d+0.1*d*z1/(1.+z1);
00163     }
00164     y+=x;
00165     y1+=x;
00166   }
00167   else
00168   {
00169     double d=x-min;
00170     double d2=.5*d;
00171     if (sig>d2)
00172     {
00173       mult=d2;
00174     }
00175     else
00176     {
00177       mult=sig;
00178     }
00179     y=-mult*eps;
00180     y1=-mult*(eps+1.e-6);
00181     if (y>0.9*d)
00182     {
00183       z=y-0.9*d;
00184       z1=y1-0.9*d;
00185       y=0.9*d+0.1*d*z/(1.+z);
00186       y1=0.9*d+0.1*d*z1/(1.+z1);
00187     }
00188     y=x-y;
00189     y1=x-y1;
00190   }
00191   jac=(y1-y)*1.e+6;
00192   return y;
00193 }
00194 
00195 void param_init_bounded_number::add_value(const dvector& _y,
00196   const dvector& ndev, const int& _ii, const double& s, const dvector& diag)
00197   {
00198     int& ii=(int&) _ii;
00199     *this = new_value_mc(_y(ii),value(*this),minb,maxb,ndev(ii),diag(ii));
00200     ii++;
00201   }
00202 
00203 void param_init_vector::add_value(const dvector& y, const dvector& ndev,
00204   const int& _ii, const double& s, const dvector& diag)
00205   {
00206     int& ii=(int&) _ii;
00207     int mmin=indexmin();
00208     int mmax=indexmax();
00209     for (int i=mmin;i<=mmax;i++)
00210     {
00211       (*this)(i)+=diag(ii)*ndev(ii);
00212       ii++;
00213     }
00214   }
00215 
00216 void param_init_bounded_vector::add_value(const dvector& y,
00217   const dvector& ndev, const int& _ii, const double& s, const dvector& diag)
00218   {
00219     int& ii=(int&) _ii;
00220     int mmin=indexmin();
00221     int mmax=indexmax();
00222     for (int i=mmin;i<=mmax;i++)
00223     {
00224       (*this)(i)= new_value_mc(y(ii),value((*this)(i)),minb,maxb,ndev(ii),
00225         diag(ii));
00226       ii++;
00227     }
00228   }
00229 
00230 /*
00231 void param_init_bounded_vector::add_value(const dvector& y,
00232   const dvector& ndev, const int& ii, const double& s)
00233   {
00234     int mmin=indexmin();
00235     int mmax=indexmax();
00236     dvector tmp(mmin,mmax);
00237     int jj=1;
00238     set_value_inv(tmp,jj);
00239     double hpi=.5*PI;
00240     for (int i=mmin;i<=mmax;i++)
00241     {
00242       tmp(i)=atan(y(ii)+ndev(ii))/hpi;
00243       s-= log(1./(hpi*(1.0+square(y(ii)+ndev(ii)))));
00244       ii++;
00245     }
00246     jj=1;
00247     dvariable pen=0.0;
00248     set_value(tmp,jj,pen);
00249   }
00250 */
00251 
00252 void param_init_matrix::add_value(const dvector& y, const dvector& ndev,
00253   const int& _ii, const double& s, const dvector& diag)
00254   {
00255     int& ii=(int&) _ii;
00256     int rmin=rowmin();
00257     int rmax=rowmax();
00258     for (int i=rmin;i<=rmax;i++)
00259     {
00260       int cmin=(*this)(i).indexmin();
00261       int cmax=(*this)(i).indexmax();
00262       for (int j=cmin;j<=cmax;j++)
00263       {
00264         (*this)(i,j)+=diag(ii)*ndev(ii);
00265         ii++;
00266       }
00267     }
00268   }
00269 
00270 void param_init_bounded_matrix::add_value(const dvector& y, const dvector& ndev,
00271   const int& _ii, const double& s, const dvector& diag)
00272   {
00273     int& ii=(int&) _ii;
00274     int rmin=rowmin();
00275     int rmax=rowmax();
00276     for (int i=rmin;i<=rmax;i++)
00277     {
00278       int cmin=(*this)(i).indexmin();
00279       int cmax=(*this)(i).indexmax();
00280       for (int j=cmin;j<=cmax;j++)
00281       {
00282         (*this)(i,j)= new_value_mc(y(ii),value((*this)(i,j)),minb,maxb,
00283           ndev(ii),diag(ii));
00284         ii++;
00285       }
00286     }
00287   }
00288 
00289 void param_init_d3array::add_value(const dvector& y, const dvector& ndev,
00290   const int& _ii, const double& s, const dvector& diag)
00291 {
00292   int& ii=(int&) _ii;
00293   int smin=indexmin();
00294   int smax=indexmax();
00295   for (int k=smin;k<=smax;k++)
00296   {
00297     int rmin=(*this)(k).indexmin();
00298     int rmax=(*this)(k).indexmax();
00299     for (int i=rmin;i<=rmax;i++)
00300     {
00301       int cmin=(*this)(k,i).indexmin();
00302       int cmax=(*this)(k,i).indexmax();
00303       for (int j=cmin;j<=cmax;j++)
00304       {
00305         (*this)(k,i,j)+=diag(ii)*ndev(ii);
00306         ii++;
00307       }
00308     }
00309   }
00310 }
00311 
00312 void param_init_d3array::add_value(const dvector& ndev, const int& _ii)
00313 {
00314   int& ii=(int&) _ii;
00315   int smin=indexmin();
00316   int smax=indexmax();
00317   for (int k=smin;k<=smax;k++)
00318   {
00319     int rmin=(*this)(k).indexmin();
00320     int rmax=(*this)(k).indexmax();
00321     for (int i=rmin;i<=rmax;i++)
00322     {
00323       int cmin=(*this)(k,i).indexmin();
00324       int cmax=(*this)(k,i).indexmax();
00325       for (int j=cmin;j<=cmax;j++)
00326       {
00327         (*this)(k,i,j)+=ndev(ii);
00328         ii++;
00329       }
00330     }
00331   }
00332 }
00333 
00334 void param_init_number::get_jacobian(const dvector& y, const dvector& ndev,
00335   const int& _ii)
00336   {
00337     int& ii=(int&) _ii;
00338     (*this)+=ndev(ii);
00339     ii++;
00340   }
00341 
00342 void param_init_bounded_number::get_jacobian(const dvector& y,
00343   const dvector& ndev, const int& _ii)
00344   {
00345     int& ii=(int&) _ii;
00346     (*this)+=ndev(ii);
00347     ii++;
00348   }
00349 
00350 void param_init_vector::get_jacobian(const dvector& y, const dvector& ndev,
00351   const int& _ii)
00352   {
00353     int& ii=(int&) _ii;
00354     int mmin=indexmin();
00355     int mmax=indexmax();
00356     for (int i=mmin;i<=mmax;i++)
00357     {
00358       (*this)(i)+=ndev(ii);
00359       ii++;
00360     }
00361   }
00362 
00363 void param_init_matrix::get_jacobian(const dvector& y, const dvector& ndev,
00364    const int& _ii)
00365   {
00366     int& ii=(int&) _ii;
00367     int rmin=rowmin();
00368     int rmax=rowmax();
00369     for (int i=rmin;i<=rmax;i++)
00370     {
00371       int cmin=(*this)(i).indexmin();
00372       int cmax=(*this)(i).indexmax();
00373       for (int j=cmin;j<=cmax;j++)
00374       {
00375         (*this)(i,j)+=ndev(ii);
00376         ii++;
00377       }
00378     }
00379   }
00380 
00381 void param_init_bounded_matrix::get_jacobian(const dvector& y,
00382   const dvector& ndev, const int& _ii)
00383   {
00384     int& ii=(int&) _ii;
00385     int rmin=rowmin();
00386     int rmax=rowmax();
00387     for (int i=rmin;i<=rmax;i++)
00388     {
00389       int cmin=(*this)(i).indexmin();
00390       int cmax=(*this)(i).indexmax();
00391       for (int j=cmin;j<=cmax;j++)
00392       {
00393         (*this)(i,j)+=ndev(ii);
00394         ii++;
00395       }
00396     }
00397   }
00398 
00399 void param_init_bounded_vector::get_jacobian(const dvector& _y,
00400   const dvector& _jac, const int& _ii)
00401   {
00402     dvector& y=(dvector&) _y;
00403     dvector& jac=(dvector&) _jac;
00404     int& ii=(int&) _ii;
00405     int mmin=indexmin();
00406     int mmax=indexmax();
00407     dvector tmp(mmin,mmax);
00408     int jj=1;
00409     set_value_inv(tmp,jj);
00410     double hpi=.5*PI;
00411     for (int i=mmin;i<=mmax;i++)
00412     {
00413       y(ii)=tan(hpi*tmp(i));
00414       jac(ii) = 1./(hpi*(1.0+square(y(ii))));
00415       ii++;
00416     }
00417   }
00418 
00419 void param_init_d3array::set_simulation_bounds(const dmatrix& _symbds,
00420   const int& _ii)
00421 {
00422   dmatrix& symbds=(dmatrix&) _symbds;
00423   int& ii=(int&) _ii;
00424   int smin=indexmin();
00425   int smax=indexmax();
00426   for (int k=smin;k<=smax;k++)
00427   {
00428     int rmin=(*this)(k).indexmin();
00429     int rmax=(*this)(k).indexmax();
00430     for (int i=rmin;i<=rmax;i++)
00431     {
00432       int cmin=(*this)(k,i).indexmin();
00433       int cmax=(*this)(k,i).indexmax();
00434       for (int j=cmin;j<=cmax;j++)
00435       {
00436         symbds(1,ii)= -simbdsmax;
00437         symbds(2,ii)= simbdsmax;
00438         ii++;
00439       }
00440     }
00441   }
00442 }