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