ADMB Documentation  11.1.2547
 All Classes Files Functions Variables Typedefs Friends Defines
mod_sd.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: mod_sd.cpp 2491 2014-10-22 08:43:24Z 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 #  include <df1b2fun.h>
00014 #include <admodel.h>
00015 
00016 dmatrix * GAUSS_varcovariance_matrix = NULL;
00017 
00018 void  set_gauss_covariance_matrix(const dll_data_matrix& m)
00019 {
00020   GAUSS_varcovariance_matrix = &((dmatrix&)(m) );
00021 }
00022 
00023 void  set_gauss_covariance_matrix(const dmatrix& m)
00024 {
00025   GAUSS_varcovariance_matrix = &((dmatrix&)(m) );
00026 }
00027 
00028 void  set_covariance_matrix(const dll_data_matrix& m)
00029 {
00030   GAUSS_varcovariance_matrix = &((dmatrix&)(m) );
00031 }
00032 
00033 void  set_covariance_matrix(const dmatrix& m)
00034 {
00035   GAUSS_varcovariance_matrix = &((dmatrix&)(m) );
00036 }
00037 
00038 void function_minimizer::sd_routine(void)
00039 {
00040   int nvar=initial_params::nvarcalc(); // get the number of active parameters
00041   dvector x(1,nvar);
00042   initial_params::xinit(x); // get the number of active parameters
00043 
00044   initial_params::restore_start_phase();
00045   int nvar1=initial_params::nvarcalc(); // get the number of active parameters
00046   int num_sdrep_types=stddev_params::num_stddev_params +
00047     initial_params::num_active_calc();
00048 
00049   param_labels.allocate(1,num_sdrep_types);
00050   param_size.allocate(1,num_sdrep_types);
00051 
00052   int ii=1;
00053   unsigned int max_name_length = 0;
00054   int i;
00055   for (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 (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 (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 (i=1;i<=ndvar;i++)
00156           {
00157             int j;
00158             for (j=1;j<=nvar;j++)
00159             {
00160               tmpsub(j)=(tv(i)*S(j))*scale(j);
00161             }
00162             ofs << tmpsub << "  ";
00163             tmpsub=tv(i)*S;
00164             for (j=1;j<=i;j++)
00165             {
00166               tmp(nvar+j)=tmpsub*tv(j);
00167               ofs << tmp(nvar+j) << " ";
00168             }
00169             diag(i+nvar)=tmp(i+nvar);
00170 
00171             if (diag(i+nvar)<=0.0)
00172             {
00173               cerr << "Estimated covariance matrix may not"
00174                " be positive definite" << endl;
00175               cerr << sort(eigenvalues(S)) << endl;
00176             }
00177             ofs << endl;
00178           }
00179         }
00180       }
00181       else  // have random effects
00182       {
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       }
00263     }
00264     // *******************************************************
00265     #else
00266     // *******************************************************
00267     {
00268       for (i=1;i<=nvar;i++)
00269       {
00270         for (int j=1;j<=i;j++)
00271         {
00272           tmp(j)=S(i,j)*scale(i)*scale(j);
00273           ofs << tmp(j) << " ";
00274         }
00275         ofs << endl;
00276         diag(i)=tmp(i);
00277       }
00278       dvector tv(1,nvar);
00279       adstring tmpstring="admodel.dep";
00280       if (ad_comm::wd_flag)
00281          tmpstring = ad_comm::adprogram_name + ".dep";
00282       cifstream cif((char*)tmpstring);
00283       int tmp_nvar,tmp_ndvar;
00284       cif >> tmp_nvar >> tmp_ndvar;
00285       dvector tmpsub(1,nvar);
00286       for (i=1;i<=ndvar;i++)
00287       {
00288         cif >> tv;  // v(i)
00289         for (int j=1;j<=nvar;j++)
00290         {
00291           tmpsub(j)=(tv*S(j))*scale(j);
00292         }
00293         ofs << tmpsub << "  ";
00294         tmpsub=tv*S;
00295         cif.seekg(0,ios::beg);
00296         cif >> tmp_nvar >> tmp_ndvar;
00297         for (j=1;j<=i;j++)
00298         {
00299           cif >> v;
00300           tmp(nvar+j)=tmpsub*v;
00301           ofs << tmp(nvar+j) << " ";
00302         }
00303         diag(i+nvar)=tmp(i+nvar);
00304 
00305         if (diag(i+nvar)<=0.0)
00306         {
00307           cerr << "Estimated covariance matrix may not be positive definite"
00308              << endl;
00309         }
00310         ofs << endl;
00311       }
00312     }
00313     // *******************************************************
00314     // *******************************************************
00315    #endif
00316   }
00317 
00318 
00319   {
00320     cifstream cif("admodel.tmp");
00321     //ofstream ofs("admodel.cor");
00322     ofstream ofs((char*)(ad_comm::adprogram_name + adstring(".cor")));
00323     ofstream ofsd((char*)(ad_comm::adprogram_name + adstring(".std")));
00324 
00325     int offset=1;
00326     dvector param_values(1,nvar1+ndvar);
00327     initial_params::copy_all_values(param_values,offset);
00328     stddev_params::copy_all_values(param_values,offset);
00329 
00330     int i;
00331     for (i=1;i<=nvar1;i++)
00332     {
00333       if (diag(i)<=0.0)
00334       {
00335         cerr << "Estimated covariance matrix may not be positive definite"
00336        << endl;
00337         exit(1);
00338       }
00339       else
00340       {
00341         diag(i)=sqrt(diag(i));
00342       }
00343     }
00344     for (i=nvar1+1;i<=nvar1+ndvar;i++)
00345     {
00346       if (diag(i)<0.0)
00347       {
00348         cerr << "Estimated covariance matrix may not be positive definite"
00349        << endl;
00350         exit(1);
00351       }
00352       else if (diag(i)==0.0)
00353       {
00354         diag(i)=0.0;
00355       }
00356       else
00357       {
00358         diag(i)=sqrt(diag(i));
00359       }
00360     }
00361 
00362     {
00363       dvector dd=diag(nvar1+1,nvar1+ndvar);
00364       dd.shift(1);
00365       int ii=0;
00366       stddev_params::get_all_sd_values(dd,ii);
00367     }
00368 
00369     int lc=1;
00370     int ic=1;
00371     ofs << " The logarithm of the determinant of the hessian = " << lndet
00372         << endl;
00373     ofs << " index  ";
00374     ofsd << " index  ";
00375     ofs << " name  ";
00376     ofsd << " name ";
00377     unsigned int inmax = max_name_length > 5 ? max_name_length - 5 : 0;
00378     for (unsigned int in = 1;in <= inmax; in++)
00379     {
00380       ofs << " ";
00381       ofsd << " ";
00382     }
00383     ofs << "  value      std.dev   ";
00384     ofsd << "  value      std.dev";
00385     for (i=1;i<=nvar+ndvar;i++)
00386     {
00387       ofs << " " << setw(4) << i << "   ";
00388     }
00389     ofs << endl;
00390     ofsd << endl;
00391 
00392     if (GAUSS_varcovariance_matrix) (*GAUSS_varcovariance_matrix).initialize();
00393 
00394     for (i=1;i<=nvar1+ndvar;i++)
00395     {
00396       int j;
00397       for (j=1;j<=i;j++)
00398       {
00399         cif >> tmp(j);
00400       }
00401       for (j=1;j<=i;j++)
00402       {
00403         if (diag(i)==0.0 || diag(j)==0.0)
00404         {
00405           tmp(j)=0.0;
00406         }
00407         else
00408         {
00409           if (i!=j)
00410           {
00411             tmp(j)/=(diag(i)*diag(j));
00412           }
00413           else
00414           {
00415             tmp(j)=1;
00416           }
00417         }
00418       }
00419       ofs << "  " << setw(4) << i << "   ";
00420       ofsd << "  " << setw(4) << i << "   ";
00421       ofs << param_labels[lc];
00422       ofsd << param_labels[lc];
00423 // get the std dev of profiles likelihood variables into the right slots
00424       if (start_stdlabels <= lc && end_stdlabels >= lc)
00425       {
00426         for (int ix=0;ix<likeprof_params::num_likeprof_params;ix++)
00427         {
00428           if (param_labels[lc]==likeprof_params::likeprofptr[ix]->label())
00429           {
00430             likeprof_params::likeprofptr[ix]->get_sigma()=diag(i);
00431           }
00432         }
00433       }
00434 
00435       unsigned int inmax = max_name_length + 1 > param_labels[lc].size()
00436                            ? max_name_length + 1 - param_labels[lc].size()
00437                            : 0;
00438       for (unsigned int in = 1; in <= inmax; in++)
00439       {
00440         ofs << " ";
00441         ofsd << " ";
00442       }
00443       if (++ic> param_size[lc])
00444       {
00445         lc++;
00446         ic=1;
00447       }
00448       ofs << setscientific() << setw(11) << setprecision(4) << param_values(i)
00449           << " ";
00450       ofs << setscientific() << setw(10) << setprecision(4) << diag(i) << " ";
00451 
00452       if (GAUSS_varcovariance_matrix)
00453       {
00454         if (GAUSS_varcovariance_matrix->indexmax()>=i)
00455           (*GAUSS_varcovariance_matrix) (i,1)=diag(i);
00456       }
00457 
00458       ofsd << setscientific() << setw(11) << setprecision(4) << param_values(i)
00459            << " ";
00460       ofsd << setscientific() << setw(10) << setprecision(4) << diag(i);
00461       for (j=1;j<=i;j++)
00462       {
00463         ofs << " " << setfixed() << setprecision(4) << setw(7)
00464             << tmp(j);
00465         if (GAUSS_varcovariance_matrix)
00466         {
00467           if (GAUSS_varcovariance_matrix->indexmax()>=i  &&
00468             (*GAUSS_varcovariance_matrix)(i).indexmax()>j)
00469           {
00470             (*GAUSS_varcovariance_matrix) (i,j+1)=tmp(j);
00471           }
00472         }
00473       }
00474       ofs << endl;
00475       ofsd << endl;
00476     }
00477   }
00478 #if defined(_MSC_VER)
00479   if (system("del admodel.tmp") == -1)
00480 #else
00481   if (unlink("admodel.tmp") == -1)
00482 #endif
00483   {
00484     char msg[40] = {"Error trying to delete temporary file "};
00485     cerr << msg << "admodel.tmp" << endl;
00486   }
00487 }