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