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