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