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