ADMB Documentation  11.6rc.3396
 All Classes Files Functions Variables Typedefs Friends Defines
newmodm5.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id$
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)
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     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       lapprox->get_hessian_components_banded_lme(this);
00171     }
00172 
00173     if (negdirections==0)
00174     {
00175       dvector g1(1,nvar);
00176       while (fmc.ireturn>=0)
00177       {
00178         fmc.fmin(f,x,g);
00179         double diff =
00180         new_value-value(likeprof_params::likeprofptr[iprof]->variable());
00181         if (fmc.itn>itnsave && diff < pow(.1,iprof)*sigma)
00182         {
00183           fmc.ifn=fmc.imax;
00184         }
00185         if (fmc.ireturn>0)
00186         {
00187           g=(*lapprox)(x,f,this);
00188         }
00189         dvariable vf=0.0;
00190         vf=initial_params::reset(dvar_vector(x));
00191         *objective_function_value::pobjfun=0.0;
00192       //**********************************************************
00193       //**********************************************************
00194       //**********************************************************
00195         if (lapprox)
00196         {
00197           if (lapprox->hesstype==2)
00198           {
00199             //lapprox->num_separable_calls=0;
00200             lapprox->separable_calls_counter=0;
00201           }
00202         }
00203       //**********************************************************
00204       //**********************************************************
00205       //**********************************************************
00206       //**********************************************************
00207         userfunction();
00208         dvariable tv=likeprof_params::likeprofptr[iprof]->variable();
00209         vf+=weight*square(new_value-tv);
00210         f+=value(vf);
00211         gradcalc(nvar,g1);
00212         g+=g1;
00213       }
00214     }
00215     initial_params::set_inactive_only_random_effects();
00216   }
00217 
00218 
00219 
00220 
00221 // ****************************************************************
00222 // ****************************************************************
00223 // ****************************************************************
00224 // ****************************************************************
00225 
00226 
00227 
00228        gradient_structure::set_NO_DERIVATIVES();
00229        iexit=fmc.iexit;
00230        ihflag=fmc.ihflag;
00231        ihang=fmc.ihang;
00232        maxfn_flag=fmc.maxfn_flag;
00233        quit_flag=fmc.quit_flag;
00234        fprof=value(initial_params::reset(dvar_vector(x)));
00235        *objective_function_value::pobjfun=0.0;
00236       //**********************************************************
00237       //**********************************************************
00238       //**********************************************************
00239         if (lapprox)
00240         {
00241           if (lapprox->hesstype==2)
00242           {
00243             //lapprox->num_separable_calls=0;
00244             lapprox->separable_calls_counter=0;
00245           }
00246         }
00247       //**********************************************************
00248       //**********************************************************
00249       //**********************************************************
00250       //**********************************************************
00251        userfunction();
00252        double tv=value(likeprof_params::likeprofptr[iprof]->variable());
00253        fprof+=value(*objective_function_value::pobjfun);
00254        penalties=weight*(new_value-tv)*(new_value-tv);
00255        fprof+=penalties;
00256        if (quit_flag=='Q') break;
00257        if (!quit_flag || quit_flag == 'N')
00258        {
00259          profile_phase++;
00260        }
00261       }
00262      }
00263     }
00264     else
00265     {
00266       fprof=global_min+20.0;
00267     }
00268    }