ADMB Documentation  11.1.1932
 All Classes Files Functions Variables Typedefs Friends Defines
newmodm5.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: newmodm5.cpp 1919 2014-04-22 22:02:01Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00007 #if defined(USE_LAPLACE)
00008 #  include <admodel.h>
00009 #  include <df1b2fun.h>
00010 #  include <adrndeff.h>
00011 
00012 void function_minimizer::prof_minimize_re(int iprof, double sigma,
00013   double new_value, const double& _fprof,const int underflow_flag,
00014   double global_min, const double& _penalties, const double& _final_weight)
00015    {
00016      double& penalties=(double&) _penalties;
00017      double& fprof=(double&) _fprof;
00018      double& final_weight=(double&) _final_weight;
00019     if (!underflow_flag)
00020     {
00021      int max_profile_phases=3;
00022      int profile_phase=1;
00023      initial_params::current_phase = initial_params::max_number_phases;
00024      while (profile_phase <= max_profile_phases)
00025      {
00026       {
00027 // ****************************************************************
00028 // ****************************************************************
00029 // ****************************************************************
00030 // ****************************************************************
00031        initial_params::set_active_only_random_effects();
00032        //cout << nvar << endl;
00033        /*int unvar=*/initial_params::nvarcalc(); // get the number of active
00034        initial_params::restore_start_phase();
00035        initial_params::set_inactive_random_effects();
00036        int nvar=initial_params::nvarcalc(); // get the number of active
00037 
00038 // ****************************************************************
00039 // ****************************************************************
00040 // ****************************************************************
00041 // ****************************************************************
00042 
00043 
00044        dvector g(1,nvar);
00045        independent_variables x(1,nvar);
00046        initial_params::xinit(x);    // get the initial values into the
00047                                     // x vector
00048        fmm fmc(nvar);
00049        fmc.maxfn= maxfn;
00050        fmc.iprint= iprint;
00051        fmc.crit = crit;
00052        fmc.imax = imax;
00053        fmc.dfn= dfn;
00054        fmc.scroll_flag= scroll_flag;
00055        fmc.min_improve=min_improve;
00056        double f=0.0;
00057        gradient_structure::set_YES_DERIVATIVES();
00058        // set convergence criterion for this phase
00059        if (!(!convergence_criteria))
00060        {
00061          int ind=min(convergence_criteria.indexmax(),
00062            initial_params::current_phase);
00063          fmc.crit=convergence_criteria(ind);
00064        }
00065        if (!(!maximum_function_evaluations))
00066        {
00067          int ind=min(maximum_function_evaluations.indexmax(),
00068            initial_params::current_phase);
00069          fmc.maxfn=int(maximum_function_evaluations(ind));
00070        }
00071        int itnsave=0;
00072        //double weight=pow(50.0,profile_phase)/(sigma*sigma);
00073        double weight;
00074        if (sigma)
00075        {
00076          weight=pow(120.0,profile_phase)/(sigma*sigma);
00077        }
00078        else
00079        {
00080          weight=pow(120.0,profile_phase);
00081        }
00082        final_weight=weight;
00083 
00084 // ****************************************************************
00085 // ****************************************************************
00086 // ****************************************************************
00087 // ****************************************************************
00088   {
00089     // calculate the number of random effects unvar
00090     // this turns on random effects variables and turns off
00091     // everything else
00092     //cout << nvar << endl;
00093     initial_params::set_active_only_random_effects();
00094     //cout << nvar << endl;
00095     int unvar=initial_params::nvarcalc(); // get the number of active
00096     //df1b2_gradlist::set_no_derivatives();
00097 
00098     if (funnel_init_var::py)
00099     {
00100       delete funnel_init_var::py;
00101       funnel_init_var::py=0;
00102     }
00103     if (lapprox)
00104     {
00105       delete lapprox;
00106       lapprox=0;
00107       df1b2variable::pool->deallocate();
00108 
00109       for (int i=1;i<df1b2variable::adpool_counter;i++)
00110       {
00111         delete df1b2variable::adpool_vector[i];
00112         df1b2variable::adpool_vector[i]=0;
00113         df1b2variable::nvar_vector[i]=0;
00114       }
00115       df1b2variable::adpool_counter=0;
00116     }
00117     lapprox=new laplace_approximation_calculator(nvar,unvar,1,nvar+unvar,
00118       this);
00119     if (lapprox==0)
00120     {
00121       cerr << "Error allocating memory for lapprox" << endl;
00122       ad_exit(1);
00123     }
00124     initial_df1b2params::current_phase=initial_params::current_phase;
00125 
00126     initial_df1b2params::save_varsptr();
00127     allocate();
00128     initial_df1b2params::restore_varsptr();
00129 
00130     df1b2_gradlist::set_no_derivatives();
00131     int nvar=initial_params::nvarcalc_all();
00132     dvector y(1,nvar);
00133     initial_params::xinit_all(y);
00134     initial_df1b2params::reset_all(y);
00135 
00136     gradient_structure::set_NO_DERIVATIVES();
00137     //vmon_begin();
00138     // see what kind of hessian we are dealing with
00139     int on=0;
00140     int nopt=0;
00141     //  DF Nov 27 11
00142     initial_params::set_inactive_only_random_effects();
00143     nvar=initial_params::nvarcalc(); // get the number of active
00144 
00145     if (lapprox->have_users_hesstype==0)
00146     {
00147       if (initial_df1b2params::separable_flag)
00148       {
00149         if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-newht",nopt))>-1)
00150         {
00151           lapprox->check_hessian_type2(this);
00152         }
00153         else
00154         {
00155           lapprox->check_hessian_type(this);
00156         }
00157       }
00158       cout << "Hesstype = " << lapprox->hesstype << endl;
00159     }
00160 
00161    /*
00162     cout << "NEED to remove !!!! " << endl;
00163     initial_df1b2params::separable_flag=0;
00164     lapprox->hesstype=1;
00165    */
00166 
00167     // linear mixed effects optimization
00168     if (laplace_approximation_calculator::variance_components_vector)
00169     {
00170       if (!lapprox)
00171       {
00172         cerr << "this can't happen" << endl;
00173         ad_exit(1);
00174       }
00175       lapprox->get_hessian_components_banded_lme(this);
00176     }
00177 
00178     if (negdirections==0)
00179     {
00180       dvector g1(1,nvar);
00181       while (fmc.ireturn>=0)
00182       {
00183         fmc.fmin(f,x,g);
00184         double diff =
00185         new_value-value(likeprof_params::likeprofptr[iprof]->variable());
00186         if (fmc.itn>itnsave && diff < pow(.1,iprof)*sigma)
00187         {
00188           fmc.ifn=fmc.imax;
00189         }
00190         if (fmc.ireturn>0)
00191         {
00192           g=(*lapprox)(x,f,this);
00193         }
00194         dvariable vf=0.0;
00195         vf=initial_params::reset(dvar_vector(x));
00196         *objective_function_value::pobjfun=0.0;
00197       //**********************************************************
00198       //**********************************************************
00199       //**********************************************************
00200       #if defined(USE_LAPLACE)
00201         if (lapprox)
00202         {
00203           if (lapprox->hesstype==2)
00204           {
00205             //lapprox->num_separable_calls=0;
00206             lapprox->separable_calls_counter=0;
00207           }
00208         }
00209       #endif
00210       //**********************************************************
00211       //**********************************************************
00212       //**********************************************************
00213       //**********************************************************
00214         userfunction();
00215         dvariable tv=likeprof_params::likeprofptr[iprof]->variable();
00216         vf+=weight*square(new_value-tv);
00217         f+=value(vf);
00218         gradcalc(nvar,g1);
00219         g+=g1;
00220       }
00221     }
00222     initial_params::set_inactive_only_random_effects();
00223   }
00224 
00225 
00226 
00227 
00228 // ****************************************************************
00229 // ****************************************************************
00230 // ****************************************************************
00231 // ****************************************************************
00232 
00233 
00234 
00235        gradient_structure::set_NO_DERIVATIVES();
00236        iexit=fmc.iexit;
00237        ihflag=fmc.ihflag;
00238        ihang=fmc.ihang;
00239        maxfn_flag=fmc.maxfn_flag;
00240        quit_flag=fmc.quit_flag;
00241        fprof=value(initial_params::reset(dvar_vector(x)));
00242        *objective_function_value::pobjfun=0.0;
00243       //**********************************************************
00244       //**********************************************************
00245       //**********************************************************
00246       #if defined(USE_LAPLACE)
00247         if (lapprox)
00248         {
00249           if (lapprox->hesstype==2)
00250           {
00251             //lapprox->num_separable_calls=0;
00252             lapprox->separable_calls_counter=0;
00253           }
00254         }
00255       #endif
00256       //**********************************************************
00257       //**********************************************************
00258       //**********************************************************
00259       //**********************************************************
00260        userfunction();
00261        double tv=value(likeprof_params::likeprofptr[iprof]->variable());
00262        fprof+=value(*objective_function_value::pobjfun);
00263        penalties=weight*(new_value-tv)*(new_value-tv);
00264        fprof+=penalties;
00265        if (quit_flag=='Q') break;
00266        if (!quit_flag || quit_flag == 'N')
00267        {
00268          profile_phase++;
00269        }
00270       }
00271      }
00272     }
00273     else
00274     {
00275       fprof=global_min+20.0;
00276     }
00277    }
00278 #endif //#if defined(USE_LAPLACE)