ADMB Documentation  11.1.1913
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2lp9.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2lp9.cpp 1667 2014-02-24 23:53:59Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #if defined(USE_LAPLACE)
00012 #  include <admodel.h>
00013 #  include <df1b2fun.h>
00014 #  include <adrndeff.h>
00015         //int fcount =0;
00016   static int no_stuff=0;
00017               static void crap(void)
00018               {
00019               }
00020               static void crap(double ff,dvector& uuu,dvector& gg)
00021               {
00022                 //cout << setprecision(10) << setw(19) << ff << " "
00023                  //    << setw(19) << uuu   << "  "  << setw(19) << gg << endl;
00024               }
00025 
00026 typedef fmm * pfmm;
00027 
00032 dvector laplace_approximation_calculator::get_uhat_quasi_newton_block_diagonal
00033   (const dvector& x,function_minimizer * pfmin)
00034 {
00035   if (separable_function_difference)
00036   {
00037     delete separable_function_difference;
00038     separable_function_difference=0;
00039   }
00040   separable_function_difference = new dvector(1,num_separable_calls);
00041 
00042   fmm ** pfmc1 = new pfmm[num_separable_calls];
00043   pfmc1--;
00044   int i;
00045   ivector ishape(1,num_separable_calls);
00046   dvector gmax(1,num_separable_calls);
00047 
00048   for (i=1;i<=num_separable_calls;i++)
00049   {
00050     int m=(*derindex)(i).indexmax();
00051     ishape(i)=m;
00052     if (m>0)
00053     {
00054     pfmc1[i] = new fmm(m);
00055     pfmc1[i]->iprint=0;
00056     pfmc1[i]->crit=inner_crit;
00057     pfmc1[i]->ireturn=0;
00058     pfmc1[i]->itn=0;
00059     pfmc1[i]->ifn=0;
00060     pfmc1[i]->ialph=0;
00061     pfmc1[i]->ihang=0;
00062     pfmc1[i]->ihflag=0;
00063     pfmc1[i]->maxfn=100;
00064     pfmc1[i]->gmax=1.e+100;
00065     pfmc1[i]->use_control_c=0;
00066     }
00067     else
00068     {
00069       pfmc1[i]= (fmm *)(0);
00070     }
00071   }
00072   dmatrix gg(1,num_separable_calls,1,ishape);
00073   dmatrix ggb(1,num_separable_calls,1,ishape);
00074   dmatrix uu(1,num_separable_calls,1,ishape);
00075   dmatrix uub(1,num_separable_calls,1,ishape);
00076   dvector ff(1,num_separable_calls);
00077   dvector ffb(1,num_separable_calls);
00078   ivector icon(1,num_separable_calls);
00079   icon.initialize();
00080   ffb=1.e+100;
00081 
00082   double f=0.0;
00083   double fb=1.e+100;
00084   dvector g(1,usize);
00085   dvector ub(1,usize);
00086   independent_variables u(1,usize);
00087   gradcalc(0,g);
00088   fmc1.itn=0;
00089   fmc1.ifn=0;
00090   fmc1.ireturn=0;
00091   initial_params::xinit(u);    // get the initial values into the
00092   fmc1.ialph=0;
00093   fmc1.ihang=0;
00094   fmc1.ihflag=0;
00095 
00096   if (init_switch)
00097   {
00098     u.initialize();
00099   }
00100 
00101   for (int ii=1;ii<=2;ii++)
00102   {
00103     // get the initial u into the uu's
00104     for (i=1;i<=num_separable_calls;i++)
00105     {
00106       int m=(*derindex)(i).indexmax();
00107       for (int j=1;j<=m;j++)
00108       {
00109         uu(i,j)=u((*derindex)(i)(j));
00110       }
00111     }
00112     fmc1.dfn=1.e-2;
00113     dvariable pen=0.0;
00114     int converged=0;
00115     int initrun_flag=1;
00116     int loop_counter=0;
00117     int loop_flag=0;
00118 
00119     while (converged==0)
00120     {
00121       if (loop_flag) loop_counter++;
00122       if (loop_counter>18)
00123       {
00124         cout << loop_counter;
00125       }
00126       if (!initrun_flag)
00127       {
00128         converged=1;
00129       }
00130       for (int i=1;i<=num_separable_calls;i++)
00131       {
00132         if (ishape(i)>0) //check to see if there are any active randoem effects
00133         {               // in this function call
00134           if (!icon(i))
00135           {
00136             independent_variables& uuu=*(independent_variables*)(&(uu(i)));
00137             if (i==19)
00138               crap(ff[i],uuu,gg[i]);
00139             (pfmc1[i])->fmin(ff[i],uuu,gg(i));
00140             if (i==19)
00141               crap();
00142             gmax(i)=fabs(pfmc1[i]->gmax);
00143             if (!initrun_flag)
00144             {
00145               if (gmax(i)<1.e-4  || pfmc1[i]->ireturn<=0)
00146               {
00147                 icon(i)=1;
00148               }
00149               else
00150               {
00151                 converged=0;
00152               }
00153             }
00154           }
00155         }
00156       }
00157       initrun_flag=0;
00158       for (int i2=1;i2<=num_separable_calls;i2++)
00159       {
00160         int m=(*derindex)(i2).indexmax();
00161         for (int j=1;j<=m;j++)
00162         {
00163           u((*derindex)(i2)(j))=uu(i2,j);
00164         }
00165       }
00166       // put the
00167       //if (fmc1.ireturn>0)
00168       {
00169         dvariable vf=0.0;
00170         pen=initial_params::reset(dvar_vector(u));
00171         *objective_function_value::pobjfun=0.0;
00172 
00173         //num_separable_calls=0;
00174 
00175         pmin->inner_opt_flag=1;
00176         pfmin->AD_uf_inner();
00177         pmin->inner_opt_flag=0;
00178 
00179         if (saddlepointflag)
00180         {
00181           *objective_function_value::pobjfun*=-1.0;
00182         }
00183         if ( no_stuff==0 && quadratic_prior::get_num_quadratic_prior()>0)
00184         {
00185           quadratic_prior::get_M_calculations();
00186         }
00187         vf+=*objective_function_value::pobjfun;
00188 
00189         objective_function_value::fun_without_pen=value(vf);
00190         vf+=pen;
00191 
00192         gradcalc(usize,g);
00193         for (int i=1;i<=num_separable_calls;i++)
00194         {
00195           int m=(*derindex)(i).indexmax();
00196           for (int j=1;j<=m;j++)
00197           {
00198             gg(i,j)=g((*derindex)(i)(j));
00199           }
00200         }
00201         {
00202           ofstream ofs("l:/temp1.dat");
00203           ofs << g.indexmax() << " " << setprecision(15) << g << endl;
00204         }
00205         if (saddlepointflag==2)
00206         {
00207           ff[1]=-(*separable_function_difference)(1);
00208           for (int i=2;i<=num_separable_calls;i++)
00209           {
00210             ff[i]=-(*separable_function_difference)(i);
00211             //ff[i]=-(*separable_function_difference)(i)
00212              // +(*separable_function_difference)(i-1);
00213 
00214             if (ff[i] < ffb[i])
00215             {
00216               ffb[i]=ff[i];
00217               uub[i]=uu[i];
00218               ggb[i]=gg[i];
00219             }
00220           }
00221         }
00222         else
00223         {
00224           ff[1]=(*separable_function_difference)(1);
00225           for (int i=2;i<=num_separable_calls;i++)
00226           {
00227             ff[i]=(*separable_function_difference)(i);
00228             //ff[i]=(*separable_function_difference)(i)
00229              // -(*separable_function_difference)(i-1);
00230 
00231             if (ff[i] < ffb[i])
00232             {
00233               ffb[i]=ff[i];
00234               uub[i]=uu[i];
00235               ggb[i]=gg[i];
00236             }
00237           }
00238         }
00239         f=0.0;
00240         for (int i2=1;i2<=num_separable_calls;i2++)
00241         {
00242           f+=ff[i2];
00243         }
00244         if (f<fb)
00245         {
00246           fb=f;
00247           ub=u;
00248         }
00249       }
00250       u=ub;
00251     }
00252     double tmax=max(gmax);
00253     cout <<  " inner maxg = " << tmax << endl;
00254 
00255     if (tmax< 1.e-4) break;
00256   }
00257   fmc1.ireturn=0;
00258   fmc1.fbest=fb;
00259   gradient_structure::set_NO_DERIVATIVES();
00260   //num_separable_calls=0;
00261   pmin->inner_opt_flag=1;
00262   pfmin->AD_uf_inner();
00263   pmin->inner_opt_flag=0;
00264   if ( no_stuff==0 && quadratic_prior::get_num_quadratic_prior()>0)
00265   {
00266     quadratic_prior::get_M_calculations();
00267   }
00268   gradient_structure::set_YES_DERIVATIVES();
00269   for (i=1;i<=num_separable_calls;i++)
00270   {
00271     if (pfmc1[i])
00272     {
00273       delete pfmc1[i];
00274     }
00275   }
00276   pfmc1++;
00277   delete [] pfmc1;
00278   pfmc1 = 0;
00279   return u;
00280 }
00281 #endif  // #if defined(USE_LAPLACE)