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