ADMB Documentation  11.1x.2711
 All Classes Files Functions Variables Typedefs Friends Defines
newmodm5.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: newmodm5.cpp 2581 2014-11-07 03:21:24Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00007 #include <admodel.h>
00008 #include <df1b2fun.h>
00009 #include <adrndeff.h>
00010 
00011 #ifdef ISZERO
00012   #undef ISZERO
00013 #endif
00014 #define ISZERO(d) ((d)==0.0)
00015 
00016 void function_minimizer::prof_minimize_re(int iprof, double sigma,
00017   double new_value, const double& _fprof,const int underflow_flag,
00018   double global_min, const double& _penalties, const double& _final_weight)
00019    {
00020      double& penalties=(double&) _penalties;
00021      double& fprof=(double&) _fprof;
00022      double& final_weight=(double&) _final_weight;
00023     if (!underflow_flag)
00024     {
00025      int max_profile_phases=3;
00026      int profile_phase=1;
00027      initial_params::current_phase = initial_params::max_number_phases;
00028      while (profile_phase <= max_profile_phases)
00029      {
00030       {
00031 // ****************************************************************
00032 // ****************************************************************
00033 // ****************************************************************
00034 // ****************************************************************
00035        initial_params::set_active_only_random_effects();
00036        //cout << nvar << endl;
00037        /*int unvar=*/initial_params::nvarcalc(); // get the number of active
00038        initial_params::restore_start_phase();
00039        initial_params::set_inactive_random_effects();
00040        int nvar=initial_params::nvarcalc(); // get the number of active
00041 
00042 // ****************************************************************
00043 // ****************************************************************
00044 // ****************************************************************
00045 // ****************************************************************
00046 
00047 
00048        dvector g(1,nvar);
00049        independent_variables x(1,nvar);
00050        initial_params::xinit(x);    // get the initial values into the
00051                                     // x vector
00052        fmm fmc(nvar);
00053        fmc.maxfn= maxfn;
00054        fmc.iprint= iprint;
00055        fmc.crit = crit;
00056        fmc.imax = imax;
00057        fmc.dfn= dfn;
00058        fmc.scroll_flag= scroll_flag;
00059        fmc.min_improve=min_improve;
00060        double f=0.0;
00061        gradient_structure::set_YES_DERIVATIVES();
00062        // set convergence criterion for this phase
00063        if (!(!convergence_criteria))
00064        {
00065          int ind=min(convergence_criteria.indexmax(),
00066            initial_params::current_phase);
00067          fmc.crit=convergence_criteria(ind);
00068        }
00069        if (!(!maximum_function_evaluations))
00070        {
00071          int ind=min(maximum_function_evaluations.indexmax(),
00072            initial_params::current_phase);
00073          fmc.maxfn=int(maximum_function_evaluations(ind));
00074        }
00075        int itnsave=0;
00076        //double weight=pow(50.0,profile_phase)/(sigma*sigma);
00077        double weight = pow(120.0,profile_phase);
00078        if (!ISZERO(sigma))
00079        {
00080          weight /= (sigma*sigma);
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 (lapprox)
00201         {
00202           if (lapprox->hesstype==2)
00203           {
00204             //lapprox->num_separable_calls=0;
00205             lapprox->separable_calls_counter=0;
00206           }
00207         }
00208       //**********************************************************
00209       //**********************************************************
00210       //**********************************************************
00211       //**********************************************************
00212         userfunction();
00213         dvariable tv=likeprof_params::likeprofptr[iprof]->variable();
00214         vf+=weight*square(new_value-tv);
00215         f+=value(vf);
00216         gradcalc(nvar,g1);
00217         g+=g1;
00218       }
00219     }
00220     initial_params::set_inactive_only_random_effects();
00221   }
00222 
00223 
00224 
00225 
00226 // ****************************************************************
00227 // ****************************************************************
00228 // ****************************************************************
00229 // ****************************************************************
00230 
00231 
00232 
00233        gradient_structure::set_NO_DERIVATIVES();
00234        iexit=fmc.iexit;
00235        ihflag=fmc.ihflag;
00236        ihang=fmc.ihang;
00237        maxfn_flag=fmc.maxfn_flag;
00238        quit_flag=fmc.quit_flag;
00239        fprof=value(initial_params::reset(dvar_vector(x)));
00240        *objective_function_value::pobjfun=0.0;
00241       //**********************************************************
00242       //**********************************************************
00243       //**********************************************************
00244         if (lapprox)
00245         {
00246           if (lapprox->hesstype==2)
00247           {
00248             //lapprox->num_separable_calls=0;
00249             lapprox->separable_calls_counter=0;
00250           }
00251         }
00252       //**********************************************************
00253       //**********************************************************
00254       //**********************************************************
00255       //**********************************************************
00256        userfunction();
00257        double tv=value(likeprof_params::likeprofptr[iprof]->variable());
00258        fprof+=value(*objective_function_value::pobjfun);
00259        penalties=weight*(new_value-tv)*(new_value-tv);
00260        fprof+=penalties;
00261        if (quit_flag=='Q') break;
00262        if (!quit_flag || quit_flag == 'N')
00263        {
00264          profile_phase++;
00265        }
00266       }
00267      }
00268     }
00269     else
00270     {
00271       fprof=global_min+20.0;
00272     }
00273    }