ADMB Documentation  11.1.2249
 All Classes Files Functions Variables Typedefs Friends Defines
getbigs.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: getbigs.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  */
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, usize = 0;
00038     dmatrix minv(1,us,1,us);
00039     dmatrix Dux(1,us,1,nvar);
00040     dmatrix uhat_prime(1,us,1,nvar);
00041     uhat_prime.initialize();
00042     int Bnvar=nvar+us;
00043     // get l_uu and l_xu for covariance calculations
00044     if (lapprox->hesstype !=2)
00045     {
00046       tmpstring = ad_comm::adprogram_name + ".luu";
00047       uistream ifs1((char*)(tmpstring));
00048       ifs1 >> usize >> xsize;
00049       if (!ifs1)
00050       {
00051         cerr << "Error reading from file " << tmpstring << endl;
00052         ad_exit(1);
00053       }
00054       // check xsize and usize
00055       if (xsize !=nvar ||usize !=us)
00056       {
00057         cerr << "size error in file " << tmpstring << endl;
00058         ad_exit(1);
00059       }
00060       ifs1 >> minv;
00061       ifs1 >> Dux;
00062       if (!ifs1)
00063       {
00064         cerr << "Error reading from file " << tmpstring << endl;
00065         ad_exit(1);
00066       }
00067       uhat_prime=-minv*Dux;
00068     }
00069     else
00070     {
00071       if (nvar !=lapprox->xsize)
00072       {
00073         cerr << "error in sizes in mod_sd" << endl;
00074         ad_exit(1);
00075       }
00076       if (us !=lapprox->usize)
00077       {
00078         cerr << "error in sizes in mod_sd" << endl;
00079         ad_exit(1);
00080       }
00081       xsize=lapprox->xsize;
00082       usize=lapprox->usize;
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     int i;
00176 
00177     for (i=1;i<=xsize;i++)
00178     {
00179       for (int j=1;j<=xsize;j++)
00180       {
00181         BS(i,j)=S(i,j);
00182       }
00183     }
00184 
00185     for (i=xsize+1;i<=Bnvar;i++)
00186     {
00187       for (int j=1;j<=xsize;j++)
00188       {
00189         BS(i,j)=Sux(i-xsize,j);
00190         BS(j,i)=BS(i,j);
00191       }
00192     }
00193 
00194     for (i=xsize+1;i<=Bnvar;i++)
00195     {
00196       for (int j=xsize+1;j<=Bnvar;j++)
00197       {
00198         BS(i,j)=Suu(i-xsize,j-xsize);
00199       }
00200     }
00201  //
00202  //   int us=nvar1-nvar;
00203  //   int xsize,usize;
00204  //   // get l_uu and l_xu for covariance calculations
00205  //   tmpstring = ad_comm::adprogram_name + ".luu";
00206  //   uistream ifs1((char*)(tmpstring));
00207  //   ifs1 >> usize >> xsize;
00208  //   if (!ifs1)
00209  //   {
00210  //     cerr << "Error reading from file " << tmpstring << endl;
00211  //     ad_exit(1);
00212  //   }
00213  //   // check xsize and usize
00214  //   if (xsize !=nvar ||usize !=us)
00215  //   {
00216  //     cerr << "size error in file " << tmpstring << endl;
00217  //     ad_exit(1);
00218  //   }
00219  //   dmatrix minv(1,usize,1,usize);
00220  //   dmatrix Dux(1,usize,1,xsize);
00221  //   int Bnvar=xsize+usize;
00222  //   ifs1 >> minv;
00223  //   ifs1 >> Dux;
00224  //   if (!ifs1)
00225  //   {
00226  //     cerr << "Error reading from file " << tmpstring << endl;
00227  //     ad_exit(1);
00228  //   }
00229  //   dmatrix S(1,nvar,1,nvar);
00230  //   read_covariance_matrix(S,nvar);
00231  //   dmatrix uhat_prime=minv*Dux;
00232  //   dmatrix Sux=uhat_prime*S;
00233  //   dmatrix Suu=Sux*trans(uhat_prime);
00234  //   Suu+=minv;
00235  //   //Suu=minv;
00236  //   minv.deallocate();
00237  //   BS.initialize();
00238  //   // random effects are never bounded?
00239  //
00240  //   int i;
00241  //
00242  //   for (i=1;i<=xsize;i++)
00243  //   {
00244  //     for (int j=1;j<=xsize;j++)
00245  //     {
00246  //       BS(i,j)=S(i,j);
00247  //     }
00248  //   }
00249  //
00250  //   for (i=xsize+1;i<=Bnvar;i++)
00251  //   {
00252  //     for (int j=1;j<=xsize;j++)
00253  //     {
00254  //       BS(i,j)=Sux(i-xsize,j);
00255  //       BS(j,i)=BS(i,j);
00256  //     }
00257  //   }
00258  //
00259  //   for (i=xsize+1;i<=Bnvar;i++)
00260  //   {
00261  //     for (int j=xsize+1;j<=Bnvar;j++)
00262  //     {
00263  //       BS(i,j)=Suu(i-xsize,j-xsize);
00264  //     }
00265  //   }
00266 }