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