ADMB Documentation  11.2.2853
 All Classes Files Functions Variables Typedefs Friends Defines
mod_sd.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: mod_sd.cpp 2691 2014-11-19 18:53:17Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00007 #ifndef _MSC_VER
00008   #include <unistd.h>
00009 #endif
00010 #if !defined(DOS386)
00011   #define DOS386
00012 #endif
00013 
00014 #include <df1b2fun.h>
00015 #include <admodel.h>
00016 
00017 dmatrix * GAUSS_varcovariance_matrix = NULL;
00018 
00019 void  set_gauss_covariance_matrix(const dll_data_matrix& m)
00020 {
00021   GAUSS_varcovariance_matrix = &((dmatrix&)(m) );
00022 }
00023 
00024 void  set_gauss_covariance_matrix(const dmatrix& m)
00025 {
00026   GAUSS_varcovariance_matrix = &((dmatrix&)(m) );
00027 }
00028 
00029 void  set_covariance_matrix(const dll_data_matrix& m)
00030 {
00031   GAUSS_varcovariance_matrix = &((dmatrix&)(m) );
00032 }
00033 
00034 void  set_covariance_matrix(const dmatrix& m)
00035 {
00036   GAUSS_varcovariance_matrix = &((dmatrix&)(m) );
00037 }
00038 
00039 void function_minimizer::sd_routine(void)
00040 {
00041   int nvar=initial_params::nvarcalc(); // get the number of active parameters
00042   dvector x(1,nvar);
00043   initial_params::xinit(x); // get the number of active parameters
00044 
00045   initial_params::restore_start_phase();
00046   int nvar1=initial_params::nvarcalc(); // get the number of active parameters
00047   int num_sdrep_types=stddev_params::num_stddev_params +
00048     initial_params::num_active_calc();
00049 
00050   param_labels.allocate(1,num_sdrep_types);
00051   param_size.allocate(1,num_sdrep_types);
00052 
00053   int ii=1;
00054   size_t max_name_length = 0;
00055   for (int i=0;i<initial_params::num_initial_params;i++)
00056   {
00057     //if ((initial_params::varsptr[i])->phase_start
00058      // <= initial_params::current_phase)
00059     if (withinbound(0,(initial_params::varsptr[i])->phase_start,
00060       initial_params::current_phase))
00061     {
00062       param_labels[ii]=
00063        (initial_params::varsptr[i])->label();
00064       param_size[ii]=
00065        (initial_params::varsptr[i])->size_count();
00066       if (max_name_length<param_labels[ii].size())
00067       {
00068         max_name_length=param_labels[ii].size();
00069       }
00070       ii++;
00071     }
00072   }
00073 
00074   int start_stdlabels=ii;
00075   for (int i=0;i< stddev_params::num_stddev_params;i++)
00076   {
00077     param_labels[ii]=
00078       stddev_params::stddevptr[i]->label();
00079     param_size[ii]=
00080       stddev_params::stddevptr[i]->size_count();
00081     if (max_name_length<param_labels[ii].size())
00082     {
00083       max_name_length=param_labels[ii].size();
00084     }
00085     ii++;
00086   }
00087   int end_stdlabels=ii-1;
00088 
00089   int ndvar=stddev_params::num_stddev_calc();
00090   dvector scale(1,nvar1);   // need to get scale from somewhere
00091   dvector v(1,nvar);  // need to read in v from model.rep
00092   dmatrix S(1,nvar,1,nvar);
00093   {
00094     uistream cif("admodel.cov");
00095     int tmp_nvar = 0;
00096     cif >> tmp_nvar;
00097     if (nvar !=tmp_nvar)
00098     {
00099       cerr << "Incorrect number of independent variables in file"
00100         " model.cov" << endl;
00101       exit(1);
00102     }
00103     cif >> S;
00104     if (!cif)
00105     {
00106       cerr << "error reading covariance matrix from model.cov" << endl;
00107       exit(1);
00108     }
00109   }
00110   int sgn;
00111   initial_params::stddev_scale(scale,x);
00112   double lndet=-ln_det(S,sgn)-2.0*sum(log(scale));
00113   initial_params::set_active_random_effects();
00114   //int nvar1=initial_params::nvarcalc();
00115   dvector diag(1,nvar1+ndvar);
00116   dvector tmp(1,nvar1+ndvar);
00117 
00118   {
00119     ofstream ofs("admodel.tmp");
00120 
00121     #if defined(__GNU__) || defined(DOS386)  || defined(__GNUDOS__)
00122     // *******************************************************
00123     // *******************************************************
00124     {
00125       if (nvar==nvar1)  // no random effects
00126       {
00127         for (int i=1;i<=nvar;i++)
00128         {
00129           for (int j=1;j<=i;j++)
00130           {
00131             tmp(j)=S(i,j)*scale(i)*scale(j);
00132             ofs << tmp(j) << " ";
00133           }
00134           ofs << endl;
00135           diag(i)=tmp(i);
00136         }
00137         dmatrix tv(1,ndvar,1,nvar1);
00138         adstring tmpstring="admodel.dep";
00139         if (ad_comm::wd_flag)
00140            tmpstring = ad_comm::adprogram_name + ".dep";
00141         cifstream cif((char*)tmpstring);
00142 
00143         int tmp_nvar = 0, tmp_ndvar = 0;
00144         cif >> tmp_nvar >> tmp_ndvar;
00145         if (tmp_nvar!=nvar1)
00146         {
00147           cerr << " tmp_nvar != nvar1 in file " << tmpstring
00148                  << endl;
00149           ad_exit(1);
00150         }
00151         if (ndvar>0)
00152         {
00153           cif >> tv;
00154           dvector tmpsub(1,nvar);
00155           for (int i=1;i<=ndvar;i++)
00156           {
00157             for (int j=1;j<=nvar;j++)
00158             {
00159               tmpsub(j)=(tv(i)*S(j))*scale(j);
00160             }
00161             ofs << tmpsub << "  ";
00162             tmpsub=tv(i)*S;
00163             for (int j=1;j<=i;j++)
00164             {
00165               tmp(nvar+j)=tmpsub*tv(j);
00166               ofs << tmp(nvar+j) << " ";
00167             }
00168             diag(i+nvar)=tmp(i+nvar);
00169 
00170             if (diag(i+nvar)<=0.0)
00171             {
00172               cerr << "Estimated covariance matrix may not"
00173                " be positive definite" << endl;
00174               cerr << sort(eigenvalues(S)) << endl;
00175             }
00176             ofs << endl;
00177           }
00178         }
00179       }
00180       else  // have random effects
00181       {
00182         dmatrix tv(1,ndvar,1,nvar1);
00183         adstring tmpstring="admodel.dep";
00184         if (ad_comm::wd_flag)
00185            tmpstring = ad_comm::adprogram_name + ".dep";
00186         cifstream cif((char*)tmpstring);
00187 
00188         int tmp_nvar = 0, tmp_ndvar = 0;
00189         cif >> tmp_nvar >> tmp_ndvar;
00190         if (tmp_nvar!=nvar1)
00191         {
00192           cerr << " tmp_nvar != nvar1 in file " << tmpstring
00193                  << endl;
00194           ad_exit(1);
00195         }
00196 
00197         dmatrix BS(1,nvar1,1,nvar1);
00198         BS.initialize();
00199         get_bigS(ndvar,nvar1,nvar,S,BS,scale);
00200 
00201         {
00202           tmpstring = ad_comm::adprogram_name + ".bgs";
00203           uostream uos((char*)(tmpstring));
00204           if (!uos)
00205           {
00206             cerr << "error opening file " << tmpstring << endl;
00207             ad_exit(1);
00208           }
00209           uos << nvar1;
00210           uos << BS;
00211           if (!uos)
00212           {
00213             cerr << "error writing to file " << tmpstring << endl;
00214             ad_exit(1);
00215           }
00216         }
00217 
00218         for (int i=1;i<=nvar1;i++)
00219         {
00220           for (int j=1;j<=i;j++)
00221           {
00222             tmp(j)=BS(i,j)*scale(i)*scale(j);
00223             ofs << tmp(j) << " ";
00224           }
00225           ofs << endl;
00226           diag(i)=tmp(i);
00227         }
00228 
00229         if (ndvar>0)
00230         {
00231           cif >> tv;
00232           dvector tmpsub(1,nvar1);
00233           for (int i=1;i<=ndvar;i++)
00234           {
00235             for (int j=1;j<=nvar1;j++)
00236             {
00237               tmpsub(j)=(tv(i)*BS(j))*scale(j);
00238             }
00239             ofs << tmpsub << "  ";
00240             tmpsub=tv(i)*BS;
00241             for (int j=1;j<=i;j++)
00242             {
00243               tmp(nvar1+j)=tmpsub*tv(j);
00244               ofs << tmp(nvar1+j) << " ";
00245             }
00246             diag(i+nvar1)=tmp(i+nvar1);
00247 
00248             if (diag(i+nvar1)<=0.0)
00249             {
00250               if (norm(tv(i))>1.e-100)
00251               {
00252                 cerr << "Estimated covariance matrix may not"
00253                  " be positive definite" << endl;
00254                 cerr << sort(eigenvalues(BS)) << endl;
00255               }
00256             }
00257             ofs << endl;
00258           }
00259         }
00260       }
00261     }
00262     // *******************************************************
00263     #else
00264     // *******************************************************
00265     {
00266       for (int i=1;i<=nvar;i++)
00267       {
00268         for (int j=1;j<=i;j++)
00269         {
00270           tmp(j)=S(i,j)*scale(i)*scale(j);
00271           ofs << tmp(j) << " ";
00272         }
00273         ofs << endl;
00274         diag(i)=tmp(i);
00275       }
00276       dvector tv(1,nvar);
00277       adstring tmpstring="admodel.dep";
00278       if (ad_comm::wd_flag)
00279          tmpstring = ad_comm::adprogram_name + ".dep";
00280       cifstream cif((char*)tmpstring);
00281       int tmp_nvar,tmp_ndvar;
00282       cif >> tmp_nvar >> tmp_ndvar;
00283       dvector tmpsub(1,nvar);
00284       for (int i=1;i<=ndvar;i++)
00285       {
00286         cif >> tv;  // v(i)
00287         for (int j=1;j<=nvar;j++)
00288         {
00289           tmpsub(j)=(tv*S(j))*scale(j);
00290         }
00291         ofs << tmpsub << "  ";
00292         tmpsub=tv*S;
00293         cif.seekg(0,ios::beg);
00294         cif >> tmp_nvar >> tmp_ndvar;
00295         for (int j=1;j<=i;j++)
00296         {
00297           cif >> v;
00298           tmp(nvar+j)=tmpsub*v;
00299           ofs << tmp(nvar+j) << " ";
00300         }
00301         diag(i+nvar)=tmp(i+nvar);
00302 
00303         if (diag(i+nvar)<=0.0)
00304         {
00305           cerr << "Estimated covariance matrix may not be positive definite"
00306              << endl;
00307         }
00308         ofs << endl;
00309       }
00310     }
00311     // *******************************************************
00312     // *******************************************************
00313    #endif
00314   }
00315 
00316 
00317   {
00318     cifstream cif("admodel.tmp");
00319     //ofstream ofs("admodel.cor");
00320     ofstream ofs((char*)(ad_comm::adprogram_name + adstring(".cor")));
00321     ofstream ofsd((char*)(ad_comm::adprogram_name + adstring(".std")));
00322 
00323     int offset=1;
00324     dvector param_values(1,nvar1+ndvar);
00325     initial_params::copy_all_values(param_values,offset);
00326     stddev_params::copy_all_values(param_values,offset);
00327 
00328     for (int i=1;i<=nvar1;i++)
00329     {
00330       if (diag(i)<=0.0)
00331       {
00332         cerr << "Estimated covariance matrix may not be positive definite"
00333        << endl;
00334         exit(1);
00335       }
00336       else
00337       {
00338         diag(i)=sqrt(diag(i));
00339       }
00340     }
00341     for (int i=nvar1+1;i<=nvar1+ndvar;i++)
00342     {
00343       if (diag(i)<0.0)
00344       {
00345         cerr << "Estimated covariance matrix may not be positive definite"
00346        << endl;
00347         exit(1);
00348       }
00349       else if (diag(i)==0.0)
00350       {
00351         diag(i)=0.0;
00352       }
00353       else
00354       {
00355         diag(i)=sqrt(diag(i));
00356       }
00357     }
00358 
00359     {
00360       dvector dd=diag(nvar1+1,nvar1+ndvar);
00361       dd.shift(1);
00362       ii=0;
00363       stddev_params::get_all_sd_values(dd,ii);
00364     }
00365 
00366     int lc=1;
00367     int ic=1;
00368     ofs << " The logarithm of the determinant of the hessian = " << lndet
00369         << endl;
00370     ofs << " index  ";
00371     ofsd << " index  ";
00372     ofs << " name  ";
00373     ofsd << " name ";
00374     size_t inmax = max_name_length > 5 ? max_name_length - 5 : 0;
00375     for (size_t in = 1;in <= inmax; in++)
00376     {
00377       ofs << " ";
00378       ofsd << " ";
00379     }
00380     ofs << "  value      std.dev   ";
00381     ofsd << "  value      std.dev";
00382     for (int i=1;i<=nvar+ndvar;i++)
00383     {
00384       ofs << " " << setw(4) << i << "   ";
00385     }
00386     ofs << endl;
00387     ofsd << endl;
00388 
00389     if (GAUSS_varcovariance_matrix) (*GAUSS_varcovariance_matrix).initialize();
00390 
00391     for (int i=1;i<=nvar1+ndvar;i++)
00392     {
00393       for (int j=1;j<=i;j++)
00394       {
00395         cif >> tmp(j);
00396       }
00397       for (int j=1;j<=i;j++)
00398       {
00399         if (diag(i)==0.0 || diag(j)==0.0)
00400         {
00401           tmp(j)=0.0;
00402         }
00403         else
00404         {
00405           if (i!=j)
00406           {
00407             tmp(j)/=(diag(i)*diag(j));
00408           }
00409           else
00410           {
00411             tmp(j)=1;
00412           }
00413         }
00414       }
00415       ofs << "  " << setw(4) << i << "   ";
00416       ofsd << "  " << setw(4) << i << "   ";
00417       ofs << param_labels[lc];
00418       ofsd << param_labels[lc];
00419 // get the std dev of profiles likelihood variables into the right slots
00420       if (start_stdlabels <= lc && end_stdlabels >= lc)
00421       {
00422         for (int ix=0;ix<likeprof_params::num_likeprof_params;ix++)
00423         {
00424           if (param_labels[lc]==likeprof_params::likeprofptr[ix]->label())
00425           {
00426             likeprof_params::likeprofptr[ix]->get_sigma()=diag(i);
00427           }
00428         }
00429       }
00430 
00431       inmax = max_name_length + 1 > param_labels[lc].size()
00432               ? max_name_length + 1 - param_labels[lc].size() : 0;
00433       for (size_t in = 1; in <= inmax; in++)
00434       {
00435         ofs << " ";
00436         ofsd << " ";
00437       }
00438       if (++ic> param_size[lc])
00439       {
00440         lc++;
00441         ic=1;
00442       }
00443       ofs << setscientific() << setw(11) << setprecision(4) << param_values(i)
00444           << " ";
00445       ofs << setscientific() << setw(10) << setprecision(4) << diag(i) << " ";
00446 
00447       if (GAUSS_varcovariance_matrix)
00448       {
00449         if (GAUSS_varcovariance_matrix->indexmax()>=i)
00450           (*GAUSS_varcovariance_matrix) (i,1)=diag(i);
00451       }
00452 
00453       ofsd << setscientific() << setw(11) << setprecision(4) << param_values(i)
00454            << " ";
00455       ofsd << setscientific() << setw(10) << setprecision(4) << diag(i);
00456       for (int j=1;j<=i;j++)
00457       {
00458         ofs << " " << setfixed() << setprecision(4) << setw(7)
00459             << tmp(j);
00460         if (GAUSS_varcovariance_matrix)
00461         {
00462           if (GAUSS_varcovariance_matrix->indexmax()>=i  &&
00463             (*GAUSS_varcovariance_matrix)(i).indexmax()>j)
00464           {
00465             (*GAUSS_varcovariance_matrix) (i,j+1)=tmp(j);
00466           }
00467         }
00468       }
00469       ofs << endl;
00470       ofsd << endl;
00471     }
00472   }
00473 #if defined(_MSC_VER)
00474   if (system("del admodel.tmp") == -1)
00475 #else
00476   if (unlink("admodel.tmp") == -1)
00477 #endif
00478   {
00479     char msg[40] = {"Error trying to delete temporary file "};
00480     cerr << msg << "admodel.tmp" << endl;
00481   }
00482 }