ADMB Documentation  11.1x.2738
 All Classes Files Functions Variables Typedefs Friends Defines
df1b2lp9.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df1b2lp9.cpp 2575 2014-11-06 19:12:44Z 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   ivector ishape(1,num_separable_calls);
00044   dvector gmax(1,num_separable_calls);
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             if (i==19)
00136               crap(ff[i],uuu,gg[i]);
00137             (pfmc1[i])->fmin(ff[i],uuu,gg(i));
00138             if (i==19)
00139               crap();
00140             gmax(i)=fabs(pfmc1[i]->gmax);
00141             if (!initrun_flag)
00142             {
00143               if (gmax(i)<1.e-4  || pfmc1[i]->ireturn<=0)
00144               {
00145                 icon(i)=1;
00146               }
00147               else
00148               {
00149                 converged=0;
00150               }
00151             }
00152           }
00153         }
00154       }
00155       initrun_flag=0;
00156       for (int i2=1;i2<=num_separable_calls;i2++)
00157       {
00158         int m=(*derindex)(i2).indexmax();
00159         for (int j=1;j<=m;j++)
00160         {
00161           u((*derindex)(i2)(j))=uu(i2,j);
00162         }
00163       }
00164       // put the
00165       //if (fmc1.ireturn>0)
00166       {
00167         dvariable vf=0.0;
00168         pen=initial_params::reset(dvar_vector(u));
00169         *objective_function_value::pobjfun=0.0;
00170 
00171         //num_separable_calls=0;
00172 
00173         pmin->inner_opt_flag=1;
00174         pfmin->AD_uf_inner();
00175         pmin->inner_opt_flag=0;
00176 
00177         if (saddlepointflag)
00178         {
00179           *objective_function_value::pobjfun*=-1.0;
00180         }
00181         if ( no_stuff==0 && quadratic_prior::get_num_quadratic_prior()>0)
00182         {
00183           quadratic_prior::get_M_calculations();
00184         }
00185         vf+=*objective_function_value::pobjfun;
00186 
00187         objective_function_value::fun_without_pen=value(vf);
00188         vf+=pen;
00189 
00190         gradcalc(usize,g);
00191         for (int i=1;i<=num_separable_calls;i++)
00192         {
00193           int m=(*derindex)(i).indexmax();
00194           for (int j=1;j<=m;j++)
00195           {
00196             gg(i,j)=g((*derindex)(i)(j));
00197           }
00198         }
00199         {
00200           ofstream ofs("l:/temp1.dat");
00201           ofs << g.indexmax() << " " << setprecision(15) << g << endl;
00202         }
00203         if (saddlepointflag==2)
00204         {
00205           ff[1]=-(*separable_function_difference)(1);
00206           for (int i=2;i<=num_separable_calls;i++)
00207           {
00208             ff[i]=-(*separable_function_difference)(i);
00209             //ff[i]=-(*separable_function_difference)(i)
00210              // +(*separable_function_difference)(i-1);
00211 
00212             if (ff[i] < ffb[i])
00213             {
00214               ffb[i]=ff[i];
00215               uub[i]=uu[i];
00216               ggb[i]=gg[i];
00217             }
00218           }
00219         }
00220         else
00221         {
00222           ff[1]=(*separable_function_difference)(1);
00223           for (int i=2;i<=num_separable_calls;i++)
00224           {
00225             ff[i]=(*separable_function_difference)(i);
00226             //ff[i]=(*separable_function_difference)(i)
00227              // -(*separable_function_difference)(i-1);
00228 
00229             if (ff[i] < ffb[i])
00230             {
00231               ffb[i]=ff[i];
00232               uub[i]=uu[i];
00233               ggb[i]=gg[i];
00234             }
00235           }
00236         }
00237         f=0.0;
00238         for (int i2=1;i2<=num_separable_calls;i2++)
00239         {
00240           f+=ff[i2];
00241         }
00242         if (f<fb)
00243         {
00244           fb=f;
00245           ub=u;
00246         }
00247       }
00248       u=ub;
00249     }
00250     double tmax=max(gmax);
00251     cout <<  " inner maxg = " << tmax << endl;
00252 
00253     if (tmax< 1.e-4) break;
00254   }
00255   fmc1.ireturn=0;
00256   fmc1.fbest=fb;
00257   gradient_structure::set_NO_DERIVATIVES();
00258   //num_separable_calls=0;
00259   pmin->inner_opt_flag=1;
00260   pfmin->AD_uf_inner();
00261   pmin->inner_opt_flag=0;
00262   if ( no_stuff==0 && quadratic_prior::get_num_quadratic_prior()>0)
00263   {
00264     quadratic_prior::get_M_calculations();
00265   }
00266   gradient_structure::set_YES_DERIVATIVES();
00267   for (int i=1;i<=num_separable_calls;i++)
00268   {
00269     if (pfmc1[i])
00270     {
00271       delete pfmc1[i];
00272     }
00273   }
00274   pfmc1++;
00275   delete [] pfmc1;
00276   pfmc1 = 0;
00277   return u;
00278 }