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