ADMB Documentation  11.1.1890
 All Classes Files Functions Variables Typedefs Friends Defines
test_trust.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: test_trust.cpp 1709 2014-02-28 21:48:21Z 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 <fvar.hpp>
00013 #  include <admodel.h>
00014 #  include <df1b2fun.h>
00015 #  include <adrndeff.h>
00016 double evaluate_function(const dvector& x,function_minimizer * pfmin);
00017 void get_second_ders(int xs,int us,const init_df1b2vector y,dmatrix& Hess,
00018   dmatrix& Dux, df1b2_gradlist * f1b2gradlist,function_minimizer * pfmin,
00019   laplace_approximation_calculator* lap);
00020 double calculate_laplace_approximation(const dvector& x,const dvector& u0,
00021   const dmatrix& Hess,const dvector& _xadjoint,const dvector& _uadjoint,
00022   const dmatrix& _Hessadjoint,function_minimizer * pmin);
00023 
00024 double calculate_importance_sample(const dvector& x,const dvector& u0,
00025   const dmatrix& Hess,const dvector& _xadjoint,const dvector& _uadjoint,
00026   const dmatrix& _Hessadjoint,function_minimizer * pmin);
00027 
00028 dmatrix choleski_decomp_positive(const dmatrix& M,double b);
00029 
00030 //dvector laplace_approximation_calculator::default_calculations
00031  // (const dvector& _x,const double& _f,function_minimizer * pfmin)
00032 
00037 dvector laplace_approximation_calculator::test_trust_region_method
00038   (const dvector& _x,const double& _f,function_minimizer * pfmin)
00039 {
00040   // for use when there is no separability
00041   ADUNCONST(dvector,x)
00042   //ADUNCONST(double,f)
00043   //int i,j;
00044   int i;
00045 
00046   initial_params::set_inactive_only_random_effects();
00047   gradient_structure::set_NO_DERIVATIVES();
00048   initial_params::reset(x);    // get current x values into the model
00049   gradient_structure::set_YES_DERIVATIVES();
00050 
00051  /*
00052   initial_params::set_active_only_random_effects();
00053   int lmn_flag=0;
00054   if (ad_comm::time_flag)
00055   {
00056     if (ad_comm::ptm1)
00057     {
00058       ad_comm::ptm1->get_elapsed_time_and_reset();
00059     }
00060     if (ad_comm::ptm)
00061     {
00062       ad_comm::ptm->get_elapsed_time_and_reset();
00063     }
00064   }
00065   if (ad_comm::time_flag)
00066   {
00067     if (ad_comm::ptm)
00068     {
00069       double time=ad_comm::ptm->get_elapsed_time();
00070       if (ad_comm::global_logfile)
00071       {
00072         (*ad_comm::global_logfile) << " Time pos 0 "
00073           << time << endl;
00074       }
00075     }
00076   }
00077 
00078   double maxg;
00079   double maxg_save;
00080   dvector uhat_old(1,usize);
00081   double f_from_1=0.0;
00082   if (!inner_lmnflag)
00083   {
00084     if (!ADqd_flag)
00085     {
00086       uhat=get_uhat_quasi_newton(x,pfmin);
00087       maxg=fabs(fmc1.gmax);
00088       f_from_1=fmc1.bestf;
00089     }
00090     else
00091       uhat=get_uhat_quasi_newton_qd(x,pfmin);
00092   }
00093   else
00094   {
00095     uhat=get_uhat_lm_newton(x,pfmin);
00096   }
00097 
00098   if (ad_comm::time_flag)
00099   {
00100     if (ad_comm::ptm)
00101     {
00102       double time=ad_comm::ptm->get_elapsed_time_and_reset();
00103       if (ad_comm::global_logfile)
00104       {
00105         (*ad_comm::global_logfile) << " Time in inner optimization "
00106           << time << endl;
00107       }
00108     }
00109   }
00110  */
00111   for (i=1;i<=xsize;i++)
00112   {
00113     y(i)=x(i);
00114   }
00115   for (i=1;i<=usize;i++)
00116   {
00117     y(i+xsize)=uhat(i);
00118   }
00119   int n=xsize+usize;
00120   dvector xx(1,n);
00121   for (i=1;i<=n;i++)
00122   {
00123     xx(i)=value(y(i));
00124   }
00125   double bestf=do_one_feval(xx,pfmin);
00126   double Delta=10;
00127   double lambda;
00128   cout << "input Delta" << endl;
00129   cin >> Delta;
00130   cout << "input lambda" << endl;
00131   cin >> lambda;
00132   int outer_iter=0;
00133   int max_iteration=10;
00134   do
00135   {
00136     outer_iter++;
00137     //int ierr=0;
00138     //int niters=0;
00139     dvector g(1,n);
00140     dmatrix H(1,n,1,n);
00141 
00142     // test newton raphson
00143     get_complete_hessian(H,g,pfmin);
00144 
00145     double tol=1.e-6;
00146     //int itmax=1000;
00147     //int itol=1;
00148     //int iter=0;
00149     //double err=0;
00150     //lambda=1;
00151 
00152     //cout << "input Delta" << endl;
00153     //cin >> Delta;
00154     //cout << "input lambda" << endl;
00155     //cin >> lambda;
00156 
00157     int i;
00158 
00159     for (i=1;i<=n;i++)
00160     {
00161       H(i,i)+=lambda;
00162     }
00163 
00164     cout << "initial fun value - " << double(0.5)*xx*(H*xx)+xx*g;
00165     double truef=do_one_feval(xx,pfmin);
00166     cout << "real function value = " << truef << endl;
00167     double estdiff;
00168     double truediff;
00169     int iflag=0;
00170     int inner_iter=0;
00171     //int oldbest=(int)bestf;
00172     int maxfn=15;
00173     dvector xret=lincg(xx,g,H,tol,Delta,pfmin,truef,estdiff,
00174      truediff,bestf,iflag,inner_iter,maxfn);
00175     cout << " norm(g) = " << norm(g)
00176          << "  norm(H*xret+g) = " << norm(H*xret+g)
00177          << "  norm(H*xret-g) = " << norm(H*xret-g)
00178          << " inner_iter = " << inner_iter << endl;
00179     if (iflag==1)
00180     {
00181       Delta/=5.0;
00182     }
00183     for (i=1;i<=n;i++)
00184     {
00185       y(i)+=xret(i);
00186     }
00187     for (i=1;i<=n;i++)
00188     {
00189       xx(i)=value(y(i));
00190     }
00191 
00192     f1b2gradlist->reset();
00193     f1b2gradlist->list.initialize();
00194     f1b2gradlist->list2.initialize();
00195     f1b2gradlist->list3.initialize();
00196     f1b2gradlist->nlist.initialize();
00197     f1b2gradlist->nlist2.initialize();
00198     f1b2gradlist->nlist3.initialize();
00199   }
00200   while(outer_iter<=max_iteration);
00201   return xx;
00202 }
00203 
00208 void laplace_approximation_calculator::get_complete_hessian
00209   (dmatrix& H,dvector& g,function_minimizer * pfmin)
00210 {
00211   int i,j,ip;
00212 
00213   if (ad_comm::time_flag)
00214   {
00215     if (ad_comm::ptm)
00216     {
00217         (*ad_comm::global_logfile) << " Starting Newton-Raphson "
00218           <<  endl;
00219     }
00220   }
00221 
00222   for (ip=1;ip<=num_der_blocks;ip++)
00223   {
00224     df1b2variable::minder=minder(ip);
00225     df1b2variable::maxder=maxder(ip);
00226 
00227     set_u_dot(ip);
00228 
00229     // do we need to reallocate memory for df1b2variables?
00230     check_for_need_to_reallocate(ip);
00231     df1b2_gradlist::set_yes_derivatives();
00232     //cout << re_objective_function_value::pobjfun << endl;
00233     //cout << re_objective_function_value::pobjfun->ptr << endl;
00234     (*re_objective_function_value::pobjfun)=double(0.0);
00235 
00236 #if defined(USE_DDOUBLE)
00237 #undef double
00238     df1b2variable pen=double(0.0);
00239     df1b2variable zz=double(0.0);
00240 #define double dd_real
00241 #else
00242     df1b2variable pen=0.0;
00243     df1b2variable zz=0.0;
00244 #endif
00245     initial_df1b2params::reset(y,pen);
00246      //cout << setprecision(15) << y << endl;
00247     if (initial_df1b2params::separable_flag)
00248     {
00249       H.initialize();
00250       grad.initialize();
00251     }
00252 
00253     double time1 = 0;
00254     if (ad_comm::time_flag)
00255     {
00256       if (ad_comm::ptm)
00257       {
00258         time1=ad_comm::ptm->get_elapsed_time();
00259       }
00260     }
00261 
00262     pfmin->user_function();
00263 
00264     if (ad_comm::time_flag)
00265     {
00266       if (ad_comm::ptm)
00267       {
00268         if (ad_comm::global_logfile)
00269         {
00270           double time=ad_comm::ptm->get_elapsed_time();
00271           (*ad_comm::global_logfile) << "       Time in user_function() "
00272             <<  ip << "  " << time-time1 << endl;
00273         }
00274       }
00275     }
00276 
00277     re_objective_function_value::fun_without_pen
00278       =value(*re_objective_function_value::pobjfun);
00279 
00280     (*re_objective_function_value::pobjfun)+=pen;
00281     (*re_objective_function_value::pobjfun)+=zz;
00282 
00283     //cout << setprecision(15) << *re_objective_function_value::pobjfun << endl;
00284     //cout << setprecision(15) << pen << endl;
00285     if (!initial_df1b2params::separable_flag)
00286     {
00287       set_dependent_variable(*re_objective_function_value::pobjfun);
00288       df1b2_gradlist::set_no_derivatives();
00289       df1b2variable::passnumber=1;
00290       df1b2_gradcalc1();
00291       int mind=y(1).minder;
00292       int jmin=max(mind,1);
00293       int jmax=min(y(1).maxder,xsize+usize);
00294       for (i=1;i<=xsize+usize;i++)
00295         for (j=jmin;j<=jmax;j++)
00296           H(i,j)=y(i).u_bar[j-mind];
00297 
00298     // ****************************************************************
00299     // ****************************************************************
00300     // ****************************************************************
00301     // temporary shit
00302      /*
00303       for (i=1;i<=usize;i++)
00304       {
00305         for (j=jmin;j<=jmax;j++)
00306         {
00307           //Hess(i,j-xsize)=y(i+xsize).u_bar[j-mind];
00308           y(i+xsize).get_u_bar_tilde()[j-mind]=Hessadjoint(i,j-xsize);
00309         }
00310       }
00311       cout << Hess << endl;
00312       df1b2variable::passnumber=2;
00313       df1b2_gradcalc1();
00314 
00315       df1b2variable::passnumber=3;
00316       df1b2_gradcalc1();
00317      */
00318     // ****************************************************************
00319     // ****************************************************************
00320     // ****************************************************************
00321     // ****************************************************************
00322       for (j=jmin;j<=jmax;j++)
00323         g(j)= re_objective_function_value::pobjfun->u_dot[j-mind];
00324     }
00325     if (ip<num_der_blocks)
00326       f1b2gradlist->reset();
00327   }
00328 
00329 
00330   // just to match master pvm routine
00331   if (ad_comm::time_flag)
00332   {
00333     if (ad_comm::ptm)
00334     {
00335       /*double time=*/ad_comm::ptm->get_elapsed_time();
00336     }
00337   }
00338 }
00339 
00340 
00341 
00342 #define EPS 1.0e-14
00343 
00344 
00345 
00346 /*
00347 int main(int argc,char * argv[])
00348 {
00349   int n=4;
00350   dmatrix H(1,n,1,n);
00351 
00352   double lambda;
00353   lambda=1.0;
00354 
00355   dvector g(1,n);
00356   dvector x(1,n);
00357   double tol=1.e-6;
00358   int itmax=1000;
00359   int itol=1;
00360   int iter=0;
00361   double err=0;
00362   double Delta=20;
00363   H.initialize();
00364   lambda=1;
00365 
00366   g=1.0;
00367   x=0.0;
00368   H(1,1)=1;
00369   H(2,2)=2;
00370   H(3,3)=3;
00371   H(4,4)=3;
00372   //cout << "input lambda" << endl;
00373   //cin >> lambda;
00374 
00375   for (int i=1;i<=n;i++)
00376   {
00377     H(i,i)+=lambda;
00378   }
00379   cout << "initial fun value - " << 0.5*x*(H*x)+x*g;
00380   double truef=do_one_feval(x_k,pfmin);
00381   cout << "real function value = " << truef << endl;
00382   dvector xret=lincg(x,g,H,tol,Delta,truef,estdiff);
00383 
00384   cout << H*xret << "   " <<  g << endl;
00385   exit(0);
00386 
00387 }
00388 */
00389 
00390 /*
00391 dvector laplace_approximation_calculator::lincg(dvector& x,
00392   dvector& c, dmatrix& H,double tol,double Delta,function_minimizer * pfmin,
00393   double& truef,double& estdiff,double& truediff,double& bestf,
00394   int& iflag,int& inner_iter,int maxfn)
00395 {
00396   iflag=0;
00397   int n=c.indexmax();
00398   dvector g_k(1,n);
00399   dvector p_k(1,n);
00400   dvector x_k(1,n);
00401   dvector xbest_k(1,n);
00402   double alpha_k;
00403   double beta_k;
00404   double n2gk=0.0;
00405   x_k=0.0;
00406   g_k=c;
00407   //g_k=H*x+c;
00408   p_k=-g_k;
00409   n2gk=norm2(g_k);
00410   int improve_flag=0;
00411   inner_iter=0;
00412   xbest_k=0;
00413 
00414   do
00415   {
00416     if (++inner_iter > maxfn)
00417       return xbest_k;
00418 
00419     dvector v=H*p_k;
00420     double tmp2=p_k*v;
00421     cout << tmp2 << endl;
00422 
00423     if  (tmp2 <=0.0)
00424     {
00425       cout << "matrix not positive definite " << tmp2 << endl;
00426      // find point at boubndary
00427 
00428      double a=norm2(p_k);
00429      double b=2.0*x_k*p_k;
00430      double cc=norm2(x_k)-square(Delta);
00431      double d=sqrt(b*b-4*a*cc);
00432      double r1=(-b+d)/(2.0*a);
00433      double r2=(-b-d)/(2.0*a);
00434 
00435      cout << " the roots are " << r1 << "  " << r2 << endl;
00436      if (r1>0)
00437      {
00438        x_k=x_k+r1*p_k;
00439      }
00440      else
00441      {
00442        x_k=x_k+r2*p_k;
00443      }
00444      cout << "function value = "  << 0.5*(x_k*(H*x_k)) + x_k*c <<
00445        " norm(x_k) =  " << norm(x_k) << endl;
00446      break;
00447     }
00448 
00449     alpha_k=n2gk/(tmp2);
00450     x_k+=alpha_k*p_k;
00451 
00452     double nxk=norm(x_k);
00453     if (nxk>Delta)
00454     {
00455       // we have goneoutside the trust region
00456       cout << " we have gone outside the trust region " << endl;
00457       // find point at boubndary
00458 
00459       double a=norm2(p_k);
00460       double b=2.0*x_k*p_k;
00461       double cc=norm2(x_k)-square(Delta);
00462       double d=b*b-4*a*cc;
00463       if (d<0)
00464         cout <<" d = " << d << endl;
00465       d=sqrt(d);
00466       double r1=(-b+d)/(2.0*a);
00467       double r2=(-b-d)/(2.0*a);
00468 
00469       cout << " the roots are " << r1 << "  " << r2 << endl;
00470       dvector y_1=x_k+r1*p_k;
00471       dvector y_2=x_k+r2*p_k;
00472       if (fabs(r1)<fabs(r2))
00473       {
00474         x_k=x_k+r1*p_k;
00475       }
00476       else
00477       {
00478         x_k=x_k+r2*p_k;
00479       }
00480       cout << norm(y_1) << " " << norm(y_2) << endl;
00481       cout << "r1 function value = "  << 0.5*(y_1*(H*y_1)) + y_1*c <<
00482           "  r2 functiomn value = "  << 0.5*(y_2*(H*y_2)) + y_2*c <<
00483               " normp_k = " << norm(p_k)  << sqrt(n2gk) << endl;
00484       break;
00485     }
00486 
00487     g_k+= alpha_k*v;
00488     double tmp=n2gk;
00489     n2gk=norm2(g_k);
00490     beta_k=n2gk/tmp;
00491     p_k=-g_k+beta_k*p_k;
00492     cout << "estimated function value = "
00493          << truef + 0.5*(x_k*(H*x_k)) + x_k*c << endl;
00494     double tf=do_one_feval(x+x_k,pfmin);
00495     cout << "real function value = " << tf << endl;
00496     if (tf<=bestf)
00497     {
00498      improve_flag=1;
00499      xbest_k=x_k;
00500      cout << "norm(H*xbest_k+c)" << norm(H*xbest_k+c)  << endl;
00501      bestf=tf;
00502     }
00503     else
00504     {
00505       if (improve_flag==0) iflag=1;
00506       return xbest_k;
00507     }
00508   }
00509   while (n2gk>tol);
00510   estdiff=0.5*(x_k*(H*x_k)) + x_k*c;
00511   cout << "estimated function value = "
00512     << truef + estdiff << endl;
00513   double tf=do_one_feval(x+x_k,pfmin);
00514   truediff=tf-truef;
00515   truef=tf;
00516   cout << "real function value = " << tf << endl;
00517   return xbest_k;
00518 }
00519 */
00520 
00525 double laplace_approximation_calculator::do_one_feval
00526   (const dvector& x,function_minimizer * pfmin)
00527 {
00528   double f=0.0;
00529   //double fb=1.e+100;
00530   dvector g(1,usize);
00531   dvector ub(1,usize);
00532   initial_params::set_active_random_effects();
00533   int nvar=initial_params::nvarcalc();
00534   int nvar1=initial_params::nvarcalc_all();
00535   cout << nvar << " " << nvar1 << endl;
00536   gradient_structure::set_NO_DERIVATIVES();
00537   dvariable vf=0.0;
00538   dvariable pen=initial_params::reset(dvar_vector(x));
00539   *objective_function_value::pobjfun=0.0;
00540   pfmin->AD_uf_inner();
00541   vf+=*objective_function_value::pobjfun;
00542   vf+=pen;
00543   f=value(vf);
00544   return f;
00545 }
00546 
00547 /*
00548 dvector laplace_approximation_calculator::lincg(dvector& x,
00549   dvector& c, dmatrix& H,double tol,double Delta,function_minimizer * pfmin,
00550   double& truef,double& estdiff,double& truediff,double& bestf,
00551   int& iflag,int& inner_iter,int maxfn)
00552 {
00553   iflag=0;
00554   int n=c.indexmax();
00555   dvector g_k(1,n);
00556   dvector p_k(1,n);
00557   dvector x_k(1,n);
00558   dvector xbest_k(1,n);
00559   double alpha_k;
00560   double beta_k;
00561   double n2gk=0.0;
00562   x_k=0.0;
00563   g_k=c;
00564   p_k=-g_k;
00565   n2gk=norm2(g_k);
00566   int improve_flag=0;
00567   inner_iter=0;
00568   xbest_k=0;
00569   cout << norm2(H*p_k+c) << endl;
00570 
00571   do
00572   {
00573     dvector v=H*p_k;
00574     double tmp2=p_k*v;
00575     alpha_k=n2gk/(tmp2);
00576     x_k+=alpha_k*p_k;
00577     double nxk=norm(x_k);
00578     g_k+= alpha_k*v;
00579     double tmp=n2gk;
00580     n2gk=norm2(g_k);
00581     beta_k=n2gk/tmp;
00582     p_k=-g_k+beta_k*p_k;
00583     cout << norm2(H*p_k+c) << endl;
00584   }
00585   while(1);
00586 }
00587 */
00588 
00589 /*
00590 dvector laplace_approximation_calculator::lincg(dvector& xinit,
00591   dvector& c, dmatrix& H,double tol,double Delta,function_minimizer * pfmin,
00592   double& truef,double& estdiff,double& truediff,double& bestf,
00593   int& iflag,int& inner_iter,int maxfn)
00594 {
00595   iflag=0;
00596   int n=c.indexmax();
00597   dmatrix g(0,20,1,n);
00598   dmatrix q(0,20,1,n);
00599   dmatrix p(0,20,1,n);
00600   dmatrix x(0,20,1,n);
00601   dmatrix r(0,20,1,n);
00602   dvector alpha(0,20);
00603   dvector beta(0,20);
00604   int k=0;
00605   x(0)=0.0;
00606   g(0)=H*x(0)+c;
00607   p(0)=-g(0);
00608   cout << norm(H*x(k)+c) << endl;
00609 
00610   do
00611   {
00612     alpha(k)=norm2(g(k))/(p(k)*H*p(k));
00613     x(k+1)=x(k)+alpha(k)*p(k);
00614     g(k+1)=g(k)+alpha(k)*H*p(k);
00615     beta(k)=norm2(g(k+1))/norm2(g(k));
00616     p(k+1)=-g(k+1)+beta(k)*p(k);
00617     k++;
00618     cout << norm(H*x(k)+c) << endl;
00619   }
00620   while(1);
00621 }
00622 */
00623 
00628 dvector laplace_approximation_calculator::lincg(dvector& xinit,
00629   dvector& c, dmatrix& H1,double tol,double Delta,function_minimizer * pfmin,
00630   double& truef,double& estdiff,double& truediff,double& bestf,
00631   int& iflag,int& inner_iter,int maxfn)
00632 {
00633   iflag=0;
00634   int n=c.indexmax();
00635   dmatrix r(1,20,1,n);
00636   dmatrix p(1,20,1,n);
00637   dmatrix x(1,20,1,n);
00638   dvector alpha(1,20);
00639   dvector beta(1,20);
00640   dmatrix H=-H1;
00641   int k=1;
00642   x(1)=0.0;
00643   r(1)=-c;
00644   p(1)=r(1);
00645   cout << norm(H*x(k)+c) << endl;
00646   cout << norm(H*x(k)-c) << endl;
00647 
00648   do
00649   {
00650     dvector w=H*p(k);
00651     alpha(k)=(r(k)*H*r(k))/norm2(w);
00652     r(k+1)=r(k)-alpha(k)*(H*p(k));
00653     beta(k)=(r(k+1)*H*r(k+1))/(r(k)*H*r(k));
00654     p(k+1)=r(k)+beta(k)*p(k);
00655     x(k+1)=x(k)+alpha(k)*p(k);
00656     k++;
00657     cout << norm(H*x(k)+c) << endl;
00658     cout << norm(H*x(k)-c) << endl;
00659   }
00660   while(1);
00661   return 0;
00662 }
00663 
00664 
00665 #endif //#if defined(USE_LAPLACE)