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