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