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