ADMB Documentation  11.1.2445
 All Classes Files Functions Variables Typedefs Friends Defines
model.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: model.cpp 2424 2014-09-29 20:53:15Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00007 #include <admodel.h>
00008 
00009 int initial_params::num_initial_params=0;
00010 
00011 #if !defined(BIG_INIT_PARAMS)
00012 const int initial_params::max_num_initial_params=4000;
00013   #if (__BORLANDC__  >= 0x0550)
00014 // this should be a resizeable array
00015 initial_params* initial_params::varsptr[4001];
00016   #else
00017 // this should be a resizeable array
00018 initial_params*
00019 initial_params::varsptr[initial_params::max_num_initial_params + 1];
00020   #endif
00021 #else
00022   const int initial_params::max_num_initial_params=250;
00023   adlist_ptr initial_params::varsptr(initial_params::max_num_initial_params);
00024 #endif
00025  int initial_params::max_number_phases=1;
00026  int initial_params::current_phase=1;
00027  int initial_params::restart_phase=0;
00028  int initial_params::sd_phase=0;
00029  int initial_params::mc_phase=0;
00030  int initial_params::mceval_phase=0;
00031  int AD_gaussflag=0;
00032  int ADqd_flag=0;
00033 
00034  double initial_params::get_scalefactor(void)
00035  {
00036    return scalefactor;
00037  }
00038  void initial_params::set_scalefactor(const double sf)
00039  {
00040    scalefactor=sf;
00041  }
00042 
00043  dvector function_minimizer::convergence_criteria;
00044  dvector function_minimizer::maximum_function_evaluations;
00045  int function_minimizer::sd_flag;
00046  adstring function_minimizer::user_data_file;
00047  adstring function_minimizer::user_par_file;
00048 
00049   int withinbound(int lb,int n,int ub)
00050   {
00051     if (lb<=n && n<=ub)
00052       return 1;
00053     else
00054      return 0;
00055   }
00056 
00060 initial_params::~initial_params()
00061 {
00062   num_initial_params--;
00063 #if defined(USE_SHARE_FLAGS)
00064   if (share_flags)
00065   {
00066     delete share_flags;
00067     share_flags = 0;
00068   }
00069 #endif
00070 }
00071 
00072 extern int* pointer_to_phase;
00073   initial_params::initial_params(void)
00074   {
00075 #  if defined(USE_SHARE_FLAGS)
00076      share_flags=0;
00077 #  endif
00078     phase_start=0;
00079     phase_save=-9999;
00080     initial_value=0;
00081     initial_value_flag=0;
00082     active_flag=0;
00083     scalefactor=0;
00084     pointer_to_phase=&initial_params::current_phase;
00085   }
00086 
00087   void initial_params::set_initial_value(double x)
00088   {
00089     initial_value=x;
00090     initial_value_flag=1;
00091   }
00092 
00093   double initial_params::get_initial_value(void)
00094   {
00095     return initial_value;
00096   }
00097 
00098   void initial_params::set_phase_start(int i)
00099   {
00100     phase_start=i;
00101     phase_save=i;
00102   }
00103 
00104   int initial_params::get_phase_start(void)
00105   {
00106     return phase_start;
00107   }
00108 
00109 
00110   void model_name_tag::allocate(const char * s)
00111   {
00112     name=s;
00113   }
00114 
00115   void initial_params::allocate(int _phase_start)
00116   {
00117 #if !defined(BIG_INIT_PARAMS)
00118     add_to_list();
00119 #else
00120     varsptr.add_to_list(this);
00121     num_initial_params++;
00122 #endif
00123     phase_start=_phase_start;
00124     phase_save=phase_start;
00125     if (max_number_phases<_phase_start) max_number_phases=_phase_start;
00126   }
00127 
00128   void initial_params::add_to_list(void)
00129   {
00130     if (num_initial_params>=initial_params::max_num_initial_params)
00131     {
00132       cerr << " This version of ADMB only supports "
00133           << initial_params::max_num_initial_params << " initial parameter"
00134         " objects" << endl;
00135       ad_exit(1);
00136     }
00137     varsptr[num_initial_params++]= this; // this is the list of
00138                                          // fundamental objects
00139   }
00140 
00141 int initial_params::correct_for_dev_objects(const dmatrix& H)
00142   {
00143     cout << H << endl << endl;
00144     int ii=1;
00145     for (int i=0;i<num_initial_params;i++)
00146     {
00147       if (withinbound(0,(varsptr[i])->phase_start,current_phase))
00148       {
00149         (varsptr[i])->dev_correction(H,ii);
00150       }
00151     }
00152     cout << H << endl << endl;
00153     return ii-1;
00154   }
00155 
00156 size_t initial_params::nvarcalc()
00157 {
00158   size_t nvar = 0;
00159   for (int i = 0; i < num_initial_params; i++)
00160   {
00161     //if ((varsptr[i])->phase_start <= current_phase)
00162 #if defined(USE_SHARE_FLAGS)
00163     if (varsptr[i]->share_flags != 0)
00164     {
00165       nvar += (varsptr[i])->shared_size_count();
00166     }
00167     else
00168     {
00169 #endif
00170       if (withinbound(0,(varsptr[i])->phase_start,current_phase))
00171       {
00172         nvar += (varsptr[i])->size_count();
00173       }
00174 #if defined(USE_SHARE_FLAGS)
00175     }
00176 #endif
00177   }
00178   return nvar;
00179 }
00180 
00181   int initial_params::num_active_calc()
00182   {
00183     int ntypes=0;
00184     for (int i=0;i<num_initial_params;i++)
00185     {
00186       //if ((varsptr[i])->phase_start <= current_phase)
00187       if (withinbound(0,(varsptr[i])->phase_start,current_phase))
00188       {
00189         ntypes++;
00190       }
00191     }
00192     return ntypes;
00193   }
00194 
00195   int initial_params::stddev_vscale(const dvar_vector& d,const dvar_vector& x)
00196   {
00197     int ii=1;
00198     for (int i=0;i<num_initial_params;i++)
00199     {
00200       if (withinbound(0,(varsptr[i])->phase_start,current_phase))
00201         (varsptr[i])->sd_vscale(d,x,ii);
00202     }
00203     return ii-1;
00204   }
00205 
00206 int initial_params::stddev_scale(const dvector& d, const dvector& x)
00207   {
00208     int ii=1;
00209     for (int i=0;i<num_initial_params;i++)
00210     {
00211       //if ((varsptr[i])->phase_start <= current_phase)
00212       if (withinbound(0,(varsptr[i])->phase_start,current_phase))
00213         (varsptr[i])->sd_scale(d,x,ii);
00214     }
00215     return ii-1;
00216   }
00217 
00218 int initial_params::stddev_curvscale(const dvector& d, const dvector& x)
00219   {
00220     int ii=1;
00221     for (int i=0;i<num_initial_params;i++)
00222     {
00223       //if ((varsptr[i])->phase_start <= current_phase)
00224       if (withinbound(0,(varsptr[i])->phase_start,current_phase))
00225         (varsptr[i])->curv_scale(d,x,ii);
00226     }
00227     return ii-1;
00228   }
00229 
00230   void initial_params::xinit(const dvector& x)
00231   {
00232     int ii=1;
00233     for (int i=0;i<num_initial_params;i++)
00234     {
00235       //if ((varsptr[i])->phase_start <= current_phase)
00236 #  if defined(USE_SHARE_FLAGS)
00237        if (varsptr[i]->share_flags !=0)
00238        {
00239           (varsptr[i])->shared_set_value_inv(x,ii);
00240        }
00241        else
00242        {
00243 #  endif
00244          if (withinbound(0,(varsptr[i])->phase_start,current_phase))
00245          {
00246            (varsptr[i])->set_value_inv(x,ii);
00247            (varsptr[i])->set_active_flag();
00248          }
00249 #  if defined(USE_SHARE_FLAGS)
00250         }
00251 #  endif
00252     }
00253   }
00254 
00255   void initial_params::set_active_only_random_effects(void)
00256   {
00257     for (int i=0;i<num_initial_params;i++)
00258     {
00259       (varsptr[i])->set_only_random_effects_active();
00260     }
00261   }
00262 
00263   void initial_params::set_inactive_only_random_effects(void)
00264   {
00265     for (int i=0;i<num_initial_params;i++)
00266     {
00267       (varsptr[i])->set_only_random_effects_inactive();
00268     }
00269   }
00270 
00271   void initial_params::set_active_random_effects(void)
00272   {
00273     for (int i=0;i<num_initial_params;i++)
00274     {
00275       (varsptr[i])->set_random_effects_active();
00276     }
00277   }
00278 
00279   void initial_params::restore_start_phase(void)
00280   {
00281     for (int i=0;i<num_initial_params;i++)
00282     {
00283       (varsptr[i])->restore_phase_start();
00284     }
00285   }
00286 
00287   void initial_params::restore_phase_start(void)
00288   {
00289     phase_start=phase_save;
00290   }
00291 
00292   void initial_params::set_inactive_random_effects(void)
00293   {
00294     for (int i=0;i<num_initial_params;i++)
00295     {
00296       (varsptr[i])->set_random_effects_inactive();
00297     }
00298   }
00299 
00300 void initial_params::xinit1(const dvector& _x, const dvector& g)
00301   {
00302     int ii=1;
00303     dvector& x=(dvector&) _x;
00304     for (int i=0;i<num_initial_params;i++)
00305     {
00306       //if ((varsptr[i])->phase_start <= current_phase)
00307       if (withinbound(0,(varsptr[i])->phase_start,current_phase))
00308       {
00309         (varsptr[i])->set_value_inv(x,ii);
00310         (varsptr[i])->set_active_flag();
00311       }
00312     }
00313     x=elem_prod(x,g);
00314   }
00315 
00316 dvariable initial_params::reset(const dvar_vector& x, const dvector& __pen)
00317   {
00318     dvector& _pen=(dvector&) __pen;
00319     int ii=1;
00320     dvariable pen=0.0;
00321     dvariable pen1;
00322     for (int i=0;i<num_initial_params;i++)
00323     {
00324       if (withinbound(0,(varsptr[i])->phase_start,current_phase))
00325       {
00326         (varsptr[i])->set_value(x,ii,pen1);
00327         _pen(ii-1)=value(pen1);
00328         pen+=pen1;
00329       }
00330     }
00331     return pen;
00332   }
00333 
00334 dvariable initial_params::reset1(const dvar_vector& _x, const dvector& g)
00335   {
00336     dvar_vector& x=(dvar_vector&) _x;
00337     int ii=1;
00338     dvariable pen=0.0;
00339     x=elem_div(x,g);
00340     for (int i=0;i<num_initial_params;i++)
00341     {
00342       //if ((varsptr[i])->phase_start <= current_phase)
00343       if (withinbound(0,(varsptr[i])->phase_start,current_phase))
00344         (varsptr[i])->set_value(x,ii,pen);
00345     }
00346     return pen;
00347   }
00348 
00349 dvariable initial_params::reset(const dvar_vector& x)
00350   {
00351     int ii=1;
00352     dvariable pen=0.0;
00353     for (int i=0;i<num_initial_params;i++)
00354     {
00355 #  if defined(USE_SHARE_FLAGS)
00356       if (varsptr[i]->share_flags !=0)
00357       {
00358          (varsptr[i])->shared_set_value(x,ii,pen);
00359       }
00360       else
00361       {
00362 #  endif
00363         //if ((varsptr[i])->phase_start <= current_phase)
00364         if (withinbound(0,(varsptr[i])->phase_start,current_phase))
00365         (varsptr[i])->set_value(x,ii,pen);
00366 #  if defined(USE_SHARE_FLAGS)
00367       }
00368 #  endif
00369     }
00370     return pen;
00371   }
00372 
00373 dvariable initial_params::reset(const dvector& x)
00374   {
00375     int ii=1;
00376     dvariable pen=0.0;
00377     for (int i=0;i<num_initial_params;i++)
00378     {
00379       //if ((varsptr[i])->phase_start <= current_phase)
00380       if (withinbound(0,(varsptr[i])->phase_start,current_phase))
00381         (varsptr[i])->set_value(x,ii,pen);
00382     }
00383     return pen;
00384   }
00385 
00386   void initial_params::save()
00387   {
00388     adstring tmp;
00389     if (current_phase==max_number_phases)
00390     {
00391       tmp="ar";
00392     }
00393     else if (current_phase>=10)
00394     {
00395       tmp=str(current_phase);
00396     }
00397     else
00398     {
00399       tmp="0" + str(current_phase);
00400     }
00401     {
00402       adstring tadstring=ad_comm::adprogram_name + adstring(".p") + tmp;
00403       ad_comm::global_savefile = new ofstream((char*)tadstring);
00404       if (ad_comm::global_savefile)
00405       {
00406         *(ad_comm::global_savefile) << setshowpoint();
00407         *(ad_comm::global_savefile) <<  "# Number of parameters = "
00408            << initial_params::nvarcalc() << " ";
00409         *(ad_comm::global_savefile) <<  " Objective function value = "
00410            << *objective_function_value::pobjfun;
00411         *(ad_comm::global_savefile) <<  "  Maximum gradient component = "
00412            << objective_function_value::gmax << endl;
00413         *(ad_comm::global_savefile) << setshowpoint();
00414 
00415         for (int i=0;i<num_initial_params;i++)
00416         {
00417            (varsptr[i])->save_value();
00418         }
00419         delete ad_comm::global_savefile;
00420         ad_comm::global_savefile=NULL;
00421       }
00422     }
00423     {
00424       adstring tadstring=ad_comm::adprogram_name + adstring(".b") + tmp;
00425       ad_comm::global_bsavefile = new uostream((char*)tadstring);
00426 
00427       if (ad_comm::global_bsavefile)
00428       {
00429         for (int i=0;i<num_initial_params;i++)
00430         {
00431            (varsptr[i])->bsave_value();
00432         }
00433         delete ad_comm::global_bsavefile;
00434         ad_comm::global_bsavefile=NULL;
00435       }
00436     }
00437   }
00438 
00439   void initial_params::set_active_flag(void)
00440   {
00441     active_flag=1;
00442   }
00443 
00444   void initial_params::set_inactive_flag(void)
00445   {
00446     active_flag=0;
00447   }
00448 
00449   int active(const initial_params& ip)
00450   {
00451     return ip.active_flag;
00452   }
00453 
00454 void param_init_number::set_value(const dvar_vector& x, const int& ii,
00455    const dvariable& pen)
00456   {
00457     if (!scalefactor)
00458       ::set_value(*this,x,ii);
00459     else
00460       ::set_value(*this,x,ii,scalefactor);
00461   }
00462 
00463 void param_init_number::set_value_inv(const dvector& x, const int& ii)
00464   {
00465     if (!scalefactor)
00466       ::set_value_inv(*this,(const dvector&)(x),ii);
00467     else
00468       ::set_value_inv(*this,(const dvector&)(x),ii,scalefactor);
00469   }
00470 
00471   param_init_number::param_init_number() : named_dvariable() ,
00472     initial_params()
00473   {
00474     //add_to_list();
00475   }
00476 
00477   int param_init_number::size_count(void)
00478   {
00479     return 1;
00480   }
00481 
00488   void param_init_bounded_number::allocate(const data_vector & v,
00489    const char * _s)
00490    {
00491     double lb = v(1);
00492     double ub = v(2);
00493     int   phz = int(v(3));
00494 
00495     allocate(lb,ub,phz,_s);
00496    }
00497 
00498   void param_init_bounded_number::allocate(double _minb,
00499     double _maxb,int _phase_start,const char * _s)
00500   {
00501     minb=_minb;
00502     maxb=_maxb;
00503     if (minb>maxb)
00504     {
00505       cerr << "Error allocating init_bounded_number " << endl
00506         << " minb = " << minb << " maxb = " << maxb << endl;
00507       ad_exit(1);
00508     }
00509     named_dvariable::allocate(_s);
00510     initial_params::allocate(_phase_start);
00511     if (ad_comm::global_bparfile)
00512     {
00513       *(ad_comm::global_bparfile) >> value(*this);
00514         if (!(*(ad_comm::global_bparfile)))
00515         {
00516           cerr << "error reading parameters from binary file "
00517                <<  endl;
00518           ad_exit(1);
00519         }
00520     }
00521     else if (ad_comm::global_parfile)
00522     {
00523       *(ad_comm::global_parfile) >> value(*this);
00524         if (!(*(ad_comm::global_parfile)))
00525         {
00526           cerr << "error reading parameters from file "
00527                << ad_comm::global_parfile->get_file_name() << endl;
00528           ad_exit(1);
00529         }
00530     }
00531     else
00532     {
00533       if ((!initial_value_flag) || initial_value <=minb
00534            || initial_value >= maxb)
00535       {
00536         prevariable::operator=((minb+maxb)/2.);
00537       }
00538       else
00539       {
00540         prevariable::operator=(initial_value);
00541       }
00542     }
00543   }
00544 void param_init_bounded_number::set_value(const dvar_vector& x,
00545   const int& ii, const dvariable& pen)
00546   {
00547     if (!scalefactor)
00548       ::set_value(*this,x,ii,minb,maxb,pen);
00549     else
00550       ::set_value(*this,x,ii,minb,maxb,pen,scalefactor);
00551   }
00552 
00553 void param_init_bounded_number::set_value_inv(const dvector& x, const int& ii)
00554   {
00555     if (!scalefactor)
00556       ::set_value_inv(*this,x,ii,minb,maxb);
00557     else
00558       ::set_value_inv(*this,x,ii,minb,maxb,scalefactor);
00559   }
00560   void param_init_bounded_number::allocate(double _minb,double _maxb,
00561     const char * _s)
00562   {
00563     allocate(_minb,_maxb,1,_s);
00564   }
00565 
00566   void check_datafile_pointer(void * p)
00567   {
00568     if (!p)
00569     {
00570       cerr << " Error trying to read in model data " << endl;
00571       cerr << " This is usual caused by a missing DAT file " << endl;
00572       ad_exit(1);
00573     }
00574   }
00575 
00576 data_number& data_number::operator=(const double& v)
00577   {
00578     val=v;
00579     return *this;
00580   }
00581 
00582   void data_number::allocate(const char * s)
00583   {
00584     check_datafile_pointer(ad_comm::global_datafile);
00585     model_name_tag::allocate(s);
00586     *(ad_comm::global_datafile) >> val;
00587   }
00588 
00589   void named_dvariable::allocate(const char * s)
00590   {
00591     model_name_tag::allocate(s);
00592   }
00593 
00594   void named_dvector::allocate(int mmin,int mmax,const char * s)
00595   {
00596     dvector::allocate(mmin,mmax);
00597     model_name_tag::allocate(s);
00598   }
00599 
00600   void named_dvector::allocate(const char * s)
00601   {
00602     dvector::allocate();
00603     model_name_tag::allocate(s);
00604   }
00605 
00606 
00607   void named_ivector::allocate(int mmin,int mmax,const char * s)
00608   {
00609     ivector::allocate(mmin,mmax);
00610     model_name_tag::allocate(s);
00611   }
00612 
00613   void named_dvar_vector::allocate(int mmin,int mmax,const char * s)
00614   {
00615     dvar_vector::allocate(mmin,mmax);
00616     model_name_tag::allocate(s);
00617   }
00618   void named_dvar_vector::allocate(const char * s)
00619   {
00620     dvar_vector::allocate();
00621     model_name_tag::allocate(s);
00622   }
00623 
00624   void param_init_number::allocate( int phasestart,const char * s)
00625   {
00626     named_dvariable::allocate(s);
00627     initial_params::allocate(phasestart);
00628     if (ad_comm::global_bparfile)
00629     {
00630       *(ad_comm::global_bparfile) >> value(*this);
00631         if (!(*(ad_comm::global_bparfile)))
00632         {
00633           cerr << "error reading parameters from binary file "
00634                << endl;
00635           ad_exit(1);
00636         }
00637     }
00638     else if (ad_comm::global_parfile)
00639     {
00640       *(ad_comm::global_parfile) >> value(*this);
00641         if (!(*(ad_comm::global_parfile)))
00642         {
00643           cerr << "error reading parameters from file "
00644                << ad_comm::global_parfile->get_file_name() << endl;
00645           ad_exit(1);
00646         }
00647     }
00648     else
00649     {
00650       prevariable::operator=(initial_value);
00651     }
00652   }
00653 
00654   void param_init_number::allocate(const char * _s)
00655   {
00656     allocate(1,_s);
00657   }
00658 
00659   void param_init_number::save_value(void)
00660   {
00661     #ifndef __ZTC__
00662       *(ad_comm::global_savefile) << label_class(this->label())
00663         << setprecision(12) << dvariable(*this) << endl;
00664     #else
00665       *(ad_comm::global_savefile) <<  label_class(this->label())
00666         << *this << endl;
00667     #endif
00668   }
00669 
00670   void param_init_number::bsave_value(void)
00671   {
00672     dvariable tmp=*this;
00673     *(ad_comm::global_bsavefile) << tmp;
00674   }
00675 
00676 void param_init_vector::set_value(const dvar_vector& x,
00677   const int& ii, const dvariable& pen)
00678   {
00679     if (!scalefactor)
00680       ::set_value(*this,x,ii);
00681     else
00682       ::set_value(*this,x,ii,scalefactor);
00683   }
00684 
00685   void param_init_vector::set_value_inv(const dvector& x, const int& ii)
00686   {
00687 #  if defined(USE_SHARE_FLAGS)
00688     if (share_flags)
00689     {
00690       int ndim=share_flags->get_shareflags()->dimension();
00691       if (ndim!=1)
00692       {
00693         cerr << "grouping flags dimension error" << endl;
00694         ad_exit(1);
00695       }
00696     }
00697     else
00698     {
00699 #  endif
00700     if (!scalefactor)
00701       ::set_value_inv(*this,x,ii);
00702     else
00703       ::set_value_inv(*this,x,ii,scalefactor);
00704 #  if defined(USE_SHARE_FLAGS)
00705     }
00706 #  endif
00707   }
00708 
00709   int param_init_vector::size_count(void)
00710   {
00711     return ::size_count(*this);
00712   }
00713 
00714   param_init_vector::param_init_vector(void) : named_dvar_vector() ,
00715     initial_params()
00716   {
00717     //add_to_list();
00718   }
00719 
00720   void param_init_vector::save_value(void)
00721   {
00722     if (!(!(*this)))
00723       *(ad_comm::global_savefile) << label_class(this->label())
00724         << dvar_vector(*this) << endl;
00725   }
00726 
00727   void param_init_vector::bsave_value(void)
00728   {
00729     if (!(!(*this)))
00730     {
00731       dvar_vector& tmp=*this;
00732       *(ad_comm::global_bsavefile) << tmp;
00733     }
00734   }
00735 
00736  /*
00737   void param_init_vector::allocate(int imin,int imax,
00738      const ivector& ishare,const char * s)
00739   {
00740     share_flags=new index_type(ishare);
00741     allocate(imin,imax,1,s);
00742   }
00743  */
00744 
00745   void param_init_vector::allocate(int imin,int imax,int phase_start,
00746      const char * s)
00747   {
00748     named_dvar_vector::allocate(imin,imax,s);
00749     if (!(!(*this)))
00750     {
00751       initial_params::allocate(phase_start);
00752       if (ad_comm::global_bparfile)
00753       {
00754         *(ad_comm::global_bparfile) >> dvar_vector(*this);
00755         if (!(*(ad_comm::global_bparfile)))
00756         {
00757           cerr << "error reading parameters from binary file "
00758                << endl;
00759           ad_exit(1);
00760         }
00761       }
00762       else if (ad_comm::global_parfile)
00763       {
00764         *(ad_comm::global_parfile) >> dvar_vector(*this);
00765         if (!(*(ad_comm::global_parfile)))
00766         {
00767           cerr << "error reading parameters from file "
00768                << ad_comm::global_parfile->get_file_name() << endl;
00769           ad_exit(1);
00770         }
00771       }
00772       else
00773       {
00774         dvar_vector::operator=(initial_value);
00775       }
00776     }
00777     else
00778     {
00779       initial_params::allocate(-1);
00780     }
00781   }
00782 
00783   void param_init_matrix::allocate(int rmin,int rmax,int cmin,int cmax,
00784     const char * s)
00785   {
00786     allocate(rmin,rmax,cmin,cmax,1,s);
00787   }
00788 
00789  /*
00790   void param_init_matrix::allocate(int rmin,int rmax,int cmin,int cmax,
00791     const imatrix& jshare,const char * s)
00792   {
00793     share_flags=new index_type(jshare);
00794     allocate(rmin,rmax,cmin,cmax,1,s);
00795   }
00796  */
00797 
00798   void param_init_matrix::allocate(int rmin,int rmax,int cmin,int cmax,
00799     int _phase_start, const char * s)
00800   {
00801     named_dvar_matrix::allocate(rmin,rmax,cmin,cmax,s);
00802     if (!(!(*this)))
00803     {
00804       initial_params::allocate(_phase_start);
00805       if (ad_comm::global_bparfile)
00806       {
00807         *(ad_comm::global_bparfile) >> dvar_matrix(*this);
00808         if (!(*(ad_comm::global_bparfile)))
00809         {
00810           cerr << "error reading parameters from binary file "
00811                << endl;
00812           ad_exit(1);
00813         }
00814       }
00815       else if (ad_comm::global_parfile)
00816       {
00817         *(ad_comm::global_parfile) >> dvar_matrix(*this);
00818         if (!(*(ad_comm::global_parfile)))
00819         {
00820           cerr << "error reading parameters from file "
00821                << ad_comm::global_parfile->get_file_name() << endl;
00822           ad_exit(1);
00823         }
00824       }
00825       else
00826       {
00827         dvar_matrix::operator=(initial_value);
00828       }
00829     }
00830     else
00831     {
00832       initial_params::allocate(-1);
00833     }
00834   }
00835 
00836   void param_init_vector::allocate(int imin,int imax,const char * s)
00837   {
00838     allocate(imin,imax,1,s);
00839   }
00840 
00841 void param_init_bounded_vector::set_value(const dvar_vector& x,
00842   const int& ii, const dvariable& pen)
00843   {
00844     if (!(!(*this)))
00845     {
00846       if (initial_params::mc_phase)
00847       {
00848         set_value_mc(*this,x,ii,minb,maxb);
00849       }
00850       else
00851       {
00852         if (!scalefactor)
00853           ::set_value(*this,x,ii,minb,maxb,pen);
00854         else
00855           ::set_value(*this,x,ii,minb,maxb,pen,scalefactor);
00856       }
00857     }
00858   }
00859 
00860   void param_init_bounded_vector::set_value_inv(const dvector& x, const int& ii)
00861   {
00862     if (!(!(*this)))
00863     {
00864       if (initial_params::mc_phase)
00865       {
00866         set_value_inv_mc(*this,x,ii,minb,maxb);
00867       }
00868       else
00869       {
00870         if (!scalefactor)
00871           ::set_value_inv(*this,x,ii,minb,maxb);
00872         else
00873           ::set_value_inv(*this,x,ii,minb,maxb,scalefactor);
00874       }
00875     }
00876   }
00877 
00878   int param_init_bounded_vector::size_count(void)
00879   {
00880     return ::size_count(*this);
00881   }
00882 
00883   param_init_bounded_vector::param_init_bounded_vector(void) :
00884     named_dvar_vector() , initial_params()
00885   {
00886     //add_to_list();
00887   }
00888 
00889   param_init_bounded_number::param_init_bounded_number(void) :
00890     param_init_number() {}
00891   //{
00892    // add_to_list();
00893   //}
00894 
00895   void param_init_bounded_vector::allocate(int imin,int imax,
00896     double _minb,double _maxb,int _phase_start,const char * s)
00897   {
00898     minb=_minb;
00899     maxb=_maxb;
00900     named_dvar_vector::allocate(imin,imax,s);
00901     initial_params::allocate(_phase_start);
00902     if (!(!(*this)))
00903     {
00904       if (ad_comm::global_bparfile)
00905       {
00906         *(ad_comm::global_bparfile) >> dvar_vector(*this);
00907         if (!(*(ad_comm::global_bparfile)))
00908         {
00909           cerr << "error reading parameters from binary file "
00910                << endl;
00911           ad_exit(1);
00912         }
00913       }
00914       else if (ad_comm::global_parfile)
00915       {
00916         *(ad_comm::global_parfile) >> dvar_vector(*this);
00917         if (!(*(ad_comm::global_parfile)))
00918         {
00919           cerr << "error reading parameters from file "
00920                << ad_comm::global_parfile->get_file_name() << endl;
00921           ad_exit(1);
00922         }
00923       }
00924       else
00925       {
00926         if ((!initial_value_flag) || initial_value <=minb
00927              || initial_value >= maxb)
00928         {
00929           initial_value=(minb+maxb)/2.;
00930         }
00931         dvar_vector::operator=(initial_value);
00932       }
00933     }
00934     else
00935     {
00936       initial_params::allocate(-1);
00937     }
00938   }
00939 
00940   void param_init_bounded_vector::allocate(int imin,int imax,
00941     double _minb,double _maxb,const char * s)
00942   {
00943     allocate(imin,imax,_minb,_maxb,1,s);
00944   }
00945 
00946   void param_init_bounded_vector::save_value(void)
00947   {
00948     if (!(!(*this)))
00949       *(ad_comm::global_savefile) << label_class(this->label())
00950                     << dvar_vector(*this) << endl;
00951   }
00952 
00953   void param_init_bounded_vector::bsave_value(void)
00954   {
00955     if (!(!(*this)))
00956     {
00957       dvar_vector& tmp=*this;
00958       *(ad_comm::global_bsavefile) << tmp;
00959     }
00960   }
00961 
00962 void param_init_matrix::set_value(const dvar_vector& x, const int& ii,
00963   const dvariable& pen)
00964   {
00965     if (!scalefactor)
00966       ::set_value(*this,x,ii);
00967     else
00968       ::set_value(*this,x,ii,scalefactor);
00969   }
00970 
00971 
00972 void param_init_matrix::set_value_inv(const dvector& x, const int& ii)
00973   {
00974     if (!scalefactor)
00975       ::set_value_inv(*this,x,ii);
00976     else
00977       ::set_value_inv(*this,x,ii,scalefactor);
00978   }
00979 
00980   param_init_matrix::param_init_matrix() : named_dvar_matrix() ,
00981     initial_params()
00982   {
00983     //add_to_list();
00984   }
00985 
00986   int param_init_matrix::size_count(void)
00987   {
00988     return ::size_count(*this);
00989   }
00990 
00991   void param_init_matrix::save_value(void)
00992   {
00993     if (!(!(*this)))
00994       *(ad_comm::global_savefile) << label_class(this->label())
00995         << dvar_matrix(*this) << endl;
00996   }
00997 
00998   void param_init_matrix::bsave_value(void)
00999   {
01000     if (!(!(*this)))
01001     {
01002       dvar_matrix& tmp=*this;
01003       *(ad_comm::global_bsavefile) << tmp;
01004     }
01005   }
01006 
01007 
01008   void data_matrix::allocate(int rmin,int rmax,int cmin,int cmax,
01009     const char * s)
01010   {
01011     named_dmatrix::allocate(rmin,rmax,cmin,cmax,s);
01012     if (!(!(*this)))
01013     {
01014       check_datafile_pointer(ad_comm::global_datafile);
01015       *(ad_comm::global_datafile) >> *this;
01016     }
01017   }
01018 
01019 void data_matrix::allocate(int rmin,int rmax, const ivector& cmin,
01020   const ivector& cmax, const char * s)
01021   {
01022     named_dmatrix::allocate(rmin,rmax,cmin,cmax,s);
01023     if (!(!(*this)))
01024     {
01025       check_datafile_pointer(ad_comm::global_datafile);
01026       *(ad_comm::global_datafile) >> *this;
01027     }
01028   }
01029 
01030 void data_matrix::allocate(int rmin,int rmax,int cmin, const ivector& cmax,
01031   const char * s)
01032   {
01033     named_dmatrix::allocate(rmin,rmax,cmin,cmax,s);
01034     if (!(!(*this)))
01035     {
01036       check_datafile_pointer(ad_comm::global_datafile);
01037       *(ad_comm::global_datafile) >> *this;
01038     }
01039   }
01040 
01041 void data_matrix::allocate(int rmin,int rmax, const ivector& cmin,int cmax,
01042   const char * s)
01043   {
01044     named_dmatrix::allocate(rmin,rmax,cmin,cmax,s);
01045     if (!(!(*this)))
01046     {
01047       check_datafile_pointer(ad_comm::global_datafile);
01048       *(ad_comm::global_datafile) >> *this;
01049     }
01050   }
01051 
01052   void data_3array::allocate(int hsl,int hsu,int rmin,int rmax,
01053     int cmin,int cmax,const char * s)
01054   {
01055     named_d3_array::allocate(hsl,hsu,rmin,rmax,cmin,cmax,s);
01056     if (!(!(*this)))
01057     {
01058       check_datafile_pointer(ad_comm::global_datafile);
01059       *(ad_comm::global_datafile) >> d3_array(*this);
01060     }
01061   }
01062 
01063 void data_3array::allocate(int hsl, int hsu, const index_type& rmin,
01064   const index_type& rmax, const index_type& cmin,
01065   const index_type& cmax, const char * s)
01066   {
01067     named_d3_array::allocate(hsl,hsu,rmin,rmax,cmin,cmax,s);
01068     if (!(!(*this)))
01069     {
01070       check_datafile_pointer(ad_comm::global_datafile);
01071       *(ad_comm::global_datafile) >> d3_array(*this);
01072     }
01073   }
01074 
01075 void data_3array::allocate(int hsl, int hsu, const ivector& rmin,
01076   int rmax, int cmin, int cmax, const char * s)
01077   {
01078     named_d3_array::allocate(hsl,hsu,rmin,rmax,cmin,cmax,s);
01079     if (!(!(*this)))
01080     {
01081       check_datafile_pointer(ad_comm::global_datafile);
01082       *(ad_comm::global_datafile) >> d3_array(*this);
01083     }
01084   }
01085 
01086  void data_3array::allocate(int hsl, int hsu, int rmin,
01087    const ivector& rmax, int cmin, int cmax, const char * s)
01088   {
01089     named_d3_array::allocate(hsl,hsu,rmin,rmax,cmin,cmax,s);
01090     if (!(!(*this)))
01091     {
01092       check_datafile_pointer(ad_comm::global_datafile);
01093       *(ad_comm::global_datafile) >> d3_array(*this);
01094     }
01095   }
01096 
01097 void data_3array::allocate(int hsl,int hsu, const ivector& rmin,
01098   const ivector& rmax,int cmin,int cmax,const char * s)
01099   {
01100     named_d3_array::allocate(hsl,hsu,rmin,rmax,cmin,cmax,s);
01101     if (!(!(*this)))
01102     {
01103       check_datafile_pointer(ad_comm::global_datafile);
01104       *(ad_comm::global_datafile) >> d3_array(*this);
01105     }
01106   }
01107 
01108 void data_3array::allocate(int hsl, int hsu, int rmin, int rmax,
01109   const ivector& cmin, int cmax, const char * s)
01110   {
01111     named_d3_array::allocate(hsl,hsu,rmin,rmax,cmin,cmax,s);
01112     if (!(!(*this)))
01113     {
01114       check_datafile_pointer(ad_comm::global_datafile);
01115       *(ad_comm::global_datafile) >> d3_array(*this);
01116     }
01117   }
01118 
01119 void data_3array::allocate(int hsl, int hsu, int rmin, int rmax,
01120   int cmin, const ivector& cmax, const char * s)
01121   {
01122     named_d3_array::allocate(hsl,hsu,rmin,rmax,cmin,cmax,s);
01123     if (!(!(*this)))
01124     {
01125       check_datafile_pointer(ad_comm::global_datafile);
01126       *(ad_comm::global_datafile) >> d3_array(*this);
01127     }
01128   }
01129 
01130  void data_3array::allocate(int hsl,int hsu,int rmin,int rmax,
01131    const ivector& cmin, const ivector& cmax,const char * s)
01132   {
01133     named_d3_array::allocate(hsl,hsu,rmin,rmax,cmin,cmax,s);
01134     if (!(!(*this)))
01135     {
01136       check_datafile_pointer(ad_comm::global_datafile);
01137       *(ad_comm::global_datafile) >> d3_array(*this);
01138     }
01139   }
01140 
01141   void data_imatrix::allocate(int rmin,int rmax,int cmin,int cmax,
01142    const char * s)
01143   {
01144     named_imatrix::allocate(rmin,rmax,cmin,cmax,s);
01145     if (!(!(*this)))
01146     {
01147       check_datafile_pointer(ad_comm::global_datafile);
01148       *(ad_comm::global_datafile) >> imatrix(*this);
01149     }
01150   }
01151 
01152 void data_imatrix::allocate(int rmin,int rmax, const index_type& cmin,
01153   const index_type& cmax, const char * s)
01154   {
01155     named_imatrix::allocate(rmin,rmax,cmin,cmax,s);
01156     if (!(!(*this)))
01157     {
01158       check_datafile_pointer(ad_comm::global_datafile);
01159       *(ad_comm::global_datafile) >> imatrix(*this);
01160     }
01161   }
01162 
01163   void data_int::allocate(int n,const char * s)
01164   {
01165     val=n;
01166     model_name_tag::allocate(s);
01167   }
01168 
01169   void data_int::allocate(const char * s)
01170   {
01171     check_datafile_pointer(ad_comm::global_datafile);
01172     *(ad_comm::global_datafile) >> val;
01173     model_name_tag::allocate(s);
01174   }
01175 
01176 void named_imatrix::allocate(int rmin,int rmax,int cmin,int cmax,
01177   const char * s)
01178 {
01179   imatrix::allocate(rmin,rmax,cmin,cmax);
01180   model_name_tag::allocate(s);
01181 }
01182 
01183 void named_imatrix::allocate(int rmin,int rmax, const index_type& cmin,
01184   const index_type& cmax, const char * s)
01185 {
01186   imatrix::allocate(rmin,rmax,cmin,cmax);
01187   model_name_tag::allocate(s);
01188 }
01189 
01190   void data_vector::allocate(int imin,int imax,const char * s)
01191   {
01192     named_dvector::allocate(imin,imax,s);
01193     if (!(!(*this)))
01194     {
01195       check_datafile_pointer(ad_comm::global_datafile);
01196       *(ad_comm::global_datafile) >> dvector(*this);
01197     }
01198   }
01199 
01200 
01201   void data_ivector::allocate(int imin,int imax,const char * s)
01202   {
01203     named_ivector::allocate(imin,imax,s);
01204     if (!(!(*this)))
01205     {
01206       check_datafile_pointer(ad_comm::global_datafile);
01207       *(ad_comm::global_datafile) >> ivector(*this);
01208     }
01209   }
01210 
01211 
01212 adstring initial_params::get_reportfile_name(void)
01213 {
01214   adstring tmp;
01215   if (current_phase==max_number_phases)
01216     tmp="rep";
01217   else if (current_phase>=10)
01218     tmp="r" + str(current_phase);
01219   else
01220     tmp="r0" + str(current_phase);
01221   tmp = "." + tmp;
01222   return tmp;
01223 }
01224 
01225 void initial_params::set_only_random_effects_active(void)
01226 {
01227   //phase_save=phase_start;
01228   phase_start=-1;
01229 }
01230 
01231 void initial_params::set_only_random_effects_inactive(void)
01232 {
01233   phase_start=phase_save;
01234 }
01235 void initial_params::set_random_effects_active(void) {;}
01236 
01237 void initial_params::set_random_effects_inactive(void) {;}
01238 
01239 pinitial_params& adlist_ptr::operator[](int i)
01240 {
01241   return (pinitial_params&)ptr[i];
01242 }
01246 adlist_ptr::adlist_ptr(int init_size)
01247 {
01248   current = 0;
01249   ptr = new ptovoid[init_size];
01250   if (ptr == 0)
01251   {
01252     cerr << "Error: allocating memory in adlist_ptr" << endl;
01253   }
01254   current_size = init_size;
01255 }
01256 void adlist_ptr::initialize()
01257 {
01258   for (int i = 0; i < current_size; i++)
01259   {
01260     ptr[i] = 0;
01261   }
01262   //reset current index to beginning
01263   current = 0;
01264 }
01268 void adlist_ptr::resize(void)
01269 {
01270   current_size *= 2;
01271   ptovoid* tmp = new ptovoid[current_size];
01272   if (tmp == 0)
01273   {
01274     cerr << "Error: allocating memory in adlist_ptr" << endl;
01275   }
01276   for (int i=0;i<current;i++)
01277   {
01278     tmp[i] = ptr[i];
01279   }
01280   delete [] ptr;
01281   ptr = tmp;
01282   tmp = 0;
01283 }
01287 void adlist_ptr::add_to_list(void* p)
01288 {
01289   if (current > current_size)
01290   {
01291     cerr << "This can't happen in adlist_ptr" << endl;
01292     exit(1);
01293   }
01294   if (current == current_size)
01295   {
01296     resize();
01297   }
01298   ptr[current++] = p;
01299 }
01303 adlist_ptr::~adlist_ptr()
01304 {
01305   current = 0;
01306   current_size = -1;
01307   if (ptr)
01308   {
01309     delete [] ptr;
01310     ptr = 0;
01311   }
01312 }