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