ADMB Documentation  11.2.2814
 All Classes Files Functions Variables Typedefs Friends Defines
getbigs.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: getbigs.cpp 2787 2014-12-11 19:33:05Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #  include <df1b2fun.h>
00012 #include <admodel.h>
00013 
00018 void function_minimizer::get_bigS(int ndvar,int nvar1,int nvar,
00019   dmatrix& S,dmatrix& BS,dvector& scale)
00020 {
00021   dmatrix tv(1,ndvar,1,nvar1);
00022   adstring tmpstring="admodel.dep";
00023   if (ad_comm::wd_flag)
00024      tmpstring = ad_comm::adprogram_name + ".dep";
00025   cifstream cif((char*)tmpstring);
00026 
00027   int tmp_nvar = 0, tmp_ndvar = 0;
00028   cif >> tmp_nvar >> tmp_ndvar;
00029   if (tmp_nvar!=nvar1)
00030   {
00031     cerr << " tmp_nvar != nvar1 in file " << tmpstring
00032            << endl;
00033     ad_exit(1);
00034   }
00035 
00036     int us=nvar1-nvar;
00037     int xsize = 0;
00038     dmatrix minv(1,us,1,us);
00039     dmatrix uhat_prime(1,us,1,nvar);
00040     uhat_prime.initialize();
00041     int Bnvar=nvar+us;
00042     // get l_uu and l_xu for covariance calculations
00043     if (lapprox->hesstype !=2)
00044     {
00045       dmatrix Dux(1,us,1,nvar);
00046       int usize = 0;
00047       tmpstring = ad_comm::adprogram_name + ".luu";
00048       uistream ifs1((char*)(tmpstring));
00049       ifs1 >> usize >> xsize;
00050       if (!ifs1)
00051       {
00052         cerr << "Error reading from file " << tmpstring << endl;
00053         ad_exit(1);
00054       }
00055       // check xsize and usize
00056       if (xsize !=nvar ||usize !=us)
00057       {
00058         cerr << "size error in file " << tmpstring << endl;
00059         ad_exit(1);
00060       }
00061       ifs1 >> minv;
00062       ifs1 >> Dux;
00063       if (!ifs1)
00064       {
00065         cerr << "Error reading from file " << tmpstring << endl;
00066         ad_exit(1);
00067       }
00068       uhat_prime=-minv*Dux;
00069     }
00070     else
00071     {
00072       if (nvar !=lapprox->xsize)
00073       {
00074         cerr << "error in sizes in mod_sd" << endl;
00075         ad_exit(1);
00076       }
00077       if (us !=lapprox->usize)
00078       {
00079         cerr << "error in sizes in mod_sd" << endl;
00080         ad_exit(1);
00081       }
00082       xsize=lapprox->xsize;
00083       // calculate uhat_prime from the block diagnal matrix
00084       d3_array & H=*(lapprox->block_diagonal_hessian);
00085       d3_array & Dux=*(lapprox->block_diagonal_Dux);
00086       int mmin=H.indexmin();
00087       int mmax=H.indexmax();
00088       for (int i=mmin;i<=mmax;i++)
00089       {
00090         if (allocated(H(i)))
00091         {
00092           dmatrix tmp=-inv(H(i))*Dux(i);
00093           int rmin=H(i).indexmin();
00094           int rmax=H(i).indexmax();
00095           int tmpmin=Dux(i).indexmin();
00096           int cmin=Dux(i,tmpmin).indexmin();
00097           int cmax=Dux(i,tmpmin).indexmax();
00098 
00099           for (int j=rmin;j<=rmax;j++)
00100           {
00101             for (int k=cmin;k<=cmax;k++)
00102             {
00103               int i1=(*lapprox->block_diagonal_re_list)(i,j)-nvar;
00104               uhat_prime(i1,(*lapprox->block_diagonal_fe_list)(i,k))=
00105                 tmp(j,k);
00106             }
00107           }
00108         }
00109       }
00110       // rescale uhat_prime to be der wrt x
00111       {
00112         int rmin=uhat_prime.indexmin();
00113         int rmax=uhat_prime.indexmax();
00114         int cmin=uhat_prime(rmin).indexmin();
00115         int cmax=uhat_prime(rmin).indexmax();
00116         for (int i=rmin;i<=rmax;i++)
00117         {
00118           for (int j=cmin;j<=cmax;j++)
00119           {
00120             uhat_prime(i,j)*=scale(j);
00121           }
00122         }
00123       }
00124     }
00125     dmatrix Sux=uhat_prime*S;
00126     dmatrix Suu=Sux*trans(uhat_prime);
00127     if (lapprox->hesstype !=2)
00128     {
00129       if (lapprox->saddlepointflag==2)
00130       {
00131         Suu-=minv;
00132       }
00133       else
00134       {
00135         Suu+=minv;
00136       }
00137     }
00138     else
00139     {
00140       d3_array & H=*(lapprox->block_diagonal_hessian);
00141       int mmin=H.indexmin();
00142       int mmax=H.indexmax();
00143       for (int i=mmin;i<=mmax;i++)
00144       {
00145         if (allocated(H(i)))
00146         {
00147           dmatrix tmp=inv(H(i));
00148           int rmin=H(i).indexmin();
00149           int rmax=H(i).indexmax();
00150 
00151           for (int j=rmin;j<=rmax;j++)
00152           {
00153             for (int k=rmin;k<=rmax;k++)
00154             {
00155               int j1=(*lapprox->block_diagonal_re_list)(i,j)-nvar;
00156               int k1=(*lapprox->block_diagonal_re_list)(i,k)-nvar;
00157               if (lapprox->saddlepointflag==2)
00158               {
00159                 Suu(j1,k1)-=tmp(j,k);
00160               }
00161               else
00162               {
00163                 Suu(j1,k1)+=tmp(j,k);
00164               }
00165             }
00166           }
00167         }
00168       }
00169     }
00170     minv.deallocate();
00171     BS.initialize();
00172     // random effects are never bounded?
00173     scale(xsize+1,Bnvar)=1.0;
00174 
00175     for (int i=1;i<=xsize;i++)
00176     {
00177       for (int j=1;j<=xsize;j++)
00178       {
00179         BS(i,j)=S(i,j);
00180       }
00181     }
00182 
00183     for (int i=xsize+1;i<=Bnvar;i++)
00184     {
00185       for (int j=1;j<=xsize;j++)
00186       {
00187         BS(i,j)=Sux(i-xsize,j);
00188         BS(j,i)=BS(i,j);
00189       }
00190     }
00191 
00192     for (int i=xsize+1;i<=Bnvar;i++)
00193     {
00194       for (int j=xsize+1;j<=Bnvar;j++)
00195       {
00196         BS(i,j)=Suu(i-xsize,j-xsize);
00197       }
00198     }
00199  //
00200  //   int us=nvar1-nvar;
00201  //   int xsize,usize;
00202  //   // get l_uu and l_xu for covariance calculations
00203  //   tmpstring = ad_comm::adprogram_name + ".luu";
00204  //   uistream ifs1((char*)(tmpstring));
00205  //   ifs1 >> usize >> xsize;
00206  //   if (!ifs1)
00207  //   {
00208  //     cerr << "Error reading from file " << tmpstring << endl;
00209  //     ad_exit(1);
00210  //   }
00211  //   // check xsize and usize
00212  //   if (xsize !=nvar ||usize !=us)
00213  //   {
00214  //     cerr << "size error in file " << tmpstring << endl;
00215  //     ad_exit(1);
00216  //   }
00217  //   dmatrix minv(1,usize,1,usize);
00218  //   dmatrix Dux(1,usize,1,xsize);
00219  //   int Bnvar=xsize+usize;
00220  //   ifs1 >> minv;
00221  //   ifs1 >> Dux;
00222  //   if (!ifs1)
00223  //   {
00224  //     cerr << "Error reading from file " << tmpstring << endl;
00225  //     ad_exit(1);
00226  //   }
00227  //   dmatrix S(1,nvar,1,nvar);
00228  //   read_covariance_matrix(S,nvar);
00229  //   dmatrix uhat_prime=minv*Dux;
00230  //   dmatrix Sux=uhat_prime*S;
00231  //   dmatrix Suu=Sux*trans(uhat_prime);
00232  //   Suu+=minv;
00233  //   //Suu=minv;
00234  //   minv.deallocate();
00235  //   BS.initialize();
00236  //   // random effects are never bounded?
00237  //
00238  //   int i;
00239  //
00240  //   for (i=1;i<=xsize;i++)
00241  //   {
00242  //     for (int j=1;j<=xsize;j++)
00243  //     {
00244  //       BS(i,j)=S(i,j);
00245  //     }
00246  //   }
00247  //
00248  //   for (i=xsize+1;i<=Bnvar;i++)
00249  //   {
00250  //     for (int j=1;j<=xsize;j++)
00251  //     {
00252  //       BS(i,j)=Sux(i-xsize,j);
00253  //       BS(j,i)=BS(i,j);
00254  //     }
00255  //   }
00256  //
00257  //   for (i=xsize+1;i<=Bnvar;i++)
00258  //   {
00259  //     for (int j=xsize+1;j<=Bnvar;j++)
00260  //     {
00261  //       BS(i,j)=Suu(i-xsize,j-xsize);
00262  //     }
00263  //   }
00264 }