ADMB Documentation  11.1.1927
 All Classes Files Functions Variables Typedefs Friends Defines
fvar_m47.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: fvar_m47.cpp 1919 2014-04-22 22:02:01Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 #include "fvar.hpp"
00012 
00013 #ifdef __TURBOC__
00014 #pragma hdrstop
00015 #include <iostream.h>
00016 #endif
00017 
00018 #ifdef __ZTC__
00019 #include <iostream.hpp>
00020 #endif
00021 
00022 #ifdef __SUN__
00023 #include <iostream.h>
00024 #endif
00025 #ifdef __NDPX__
00026 #include <iostream.h>
00027 #endif
00028 void dfcholeski_decomp_positive(void);
00029 dvar_matrix choleski_decomp_positive(const dvar_matrix& MM, double eps,
00030   dvariable& _fpen);
00031 
00036 dvar_matrix positive_definite_matrix(const dvar_matrix& MM, double eps,
00037   dvariable& _fpen)
00038 {
00039   dvar_matrix ch_m=choleski_decomp_positive(MM,eps,_fpen);
00040   return ch_m*trans(ch_m);
00041 }
00042 
00047 dvar_matrix choleski_decomp_positive(const dvar_matrix& MM, double eps,
00048   dvariable& _fpen)
00049 {
00050   // kludge to deal with constantness
00051   dmatrix M=value(MM);
00052   double fpen;
00053   if (M.colsize() != M.rowsize())
00054   {
00055     cerr << "Error in chol_decomp. Matrix not square" << endl;
00056     ad_exit(1);
00057   }
00058   int rowsave=M.rowmin();
00059   int colsave=M.colmin();
00060   M.rowshift(1);
00061   M.colshift(1);
00062   int n=M.rowmax();
00063 
00064   dmatrix L(1,n,1,n);
00065   //dmatrix C(1,n,1,n);
00066   //imatrix B(1,n,1,n);
00067   //B.initialize();
00068 #ifndef SAFE_INITIALIZE
00069     L.initialize();
00070 #endif
00071 
00072   int i,j,k;
00073   double tmp;
00074   double ptmp;
00075   fpen=0.0;
00076   ptmp=posfun(M(1,1),eps,fpen);
00077   L(1,1)=sqrt(ptmp);
00078   for (i=2;i<=n;i++)
00079   {
00080     L(i,1)=M(i,1)/L(1,1);
00081   }
00082 
00083   for (i=2;i<=n;i++)
00084   {
00085     for (j=2;j<=i-1;j++)
00086     {
00087       tmp=M(i,j);
00088       for (k=1;k<=j-1;k++)
00089       {
00090         tmp-=L(i,k)*L(j,k);
00091       }
00092       L(i,j)=tmp/L(j,j);
00093     }
00094     tmp=M(i,i);
00095     for (k=1;k<=i-1;k++)
00096     {
00097       tmp-=L(i,k)*L(i,k);
00098     }
00099     ptmp=posfun(tmp,eps,fpen);
00100     L(i,i)=sqrt(ptmp);
00101   }
00102   L.rowshift(rowsave);
00103   L.colshift(colsave);
00104 
00105   value(_fpen)=fpen;
00106   dvar_matrix vc=nograd_assign(L);
00107   save_identifier_string("qs");
00108   _fpen.save_prevariable_position();
00109   save_double_value(eps);
00110   vc.save_dvar_matrix_position();
00111   MM.save_dvar_matrix_value();
00112   save_identifier_string("rl");
00113   MM.save_dvar_matrix_position();
00114   save_identifier_string("lo");
00115   gradient_structure::GRAD_STACK1->
00116       set_gradient_stack(dfcholeski_decomp_positive);
00117   return vc;
00118 }
00119 
00124 void dfcholeski_decomp_positive(void)
00125 {
00126   verify_identifier_string("lo");
00127   dvar_matrix_position MMpos=restore_dvar_matrix_position();
00128   verify_identifier_string("rl");
00129   dmatrix M=restore_dvar_matrix_value(MMpos);
00130   dvar_matrix_position vcpos=restore_dvar_matrix_position();
00131   double eps=restore_double_value();
00132   prevariable_position fpenpos=restore_prevariable_position();
00133   verify_identifier_string("qs");
00134   dmatrix dfL=restore_dvar_matrix_derivatives(vcpos);
00135   double dfpen=restore_prevariable_derivative(fpenpos);
00136 
00137   if (M.colsize() != M.rowsize())
00138   {
00139     cerr << "Error in chol_decomp. Matrix not square" << endl;
00140     ad_exit(1);
00141   }
00142   int rowsave=M.rowmin();
00143   int colsave=M.colmin();
00144   dfL.rowshift(1);
00145   dfL.colshift(1);
00146   M.rowshift(1);
00147   M.colshift(1);
00148   int n=M.rowmax();
00149 
00150   dmatrix L(1,n,1,n);
00151   dvector tmp(1,n);
00152   dvector ptmp(1,n);
00153   dmatrix tmp1(1,n,1,n);
00154   dmatrix dftmp1(1,n,1,n);
00155   dmatrix dfM(1,n,1,n);
00156   dvector dfptmp(1,n);
00157   dvector dftmp(1,n);
00158   tmp.initialize();
00159   tmp1.initialize();
00160   dftmp.initialize();
00161   dftmp1.initialize();
00162   dfM.initialize();
00163 #ifndef SAFE_INITIALIZE
00164     L.initialize();
00165 #endif
00166 
00167   double fpen;
00168   int i,j,k;
00169   ptmp(1)=posfun(M(1,1),eps,fpen);
00170   L(1,1)=sqrt(ptmp(1));
00171   for (i=2;i<=n;i++)
00172   {
00173     L(i,1)=M(i,1)/L(1,1);
00174   }
00175 
00176   for (i=2;i<=n;i++)
00177   {
00178     for (j=2;j<=i-1;j++)
00179     {
00180       tmp1(i,j)=M(i,j);
00181       for (k=1;k<=j-1;k++)
00182       {
00183         tmp1(i,j)-=L(i,k)*L(j,k);
00184       }
00185       L(i,j)=tmp1(i,j)/L(j,j);
00186     }
00187     tmp(i)=M(i,i);
00188     for (k=1;k<=i-1;k++)
00189     {
00190       tmp(i)-=L(i,k)*L(i,k);
00191     }
00192     double pen;
00193     ptmp(i)=posfun(tmp(i),eps,pen);
00194     L(i,i)=sqrt(ptmp(i));
00195   }
00196 
00197   dfptmp.initialize();
00198 
00199  //*******************************************************************8
00200  //*******************************************************************8
00201  //*******************************************************************8
00202   for (i=n;i>=2;i--)
00203   {
00204     //L(i,i)=sqrt(ptmp(i));
00205     dfptmp(i)+=dfL(i,i)/(2.0*L(i,i));
00206     dfL(i,i)=0.0;
00207     // ptmp(i)=posfun(tmp(i),eps);
00208     dftmp(i)=dfptmp(i)*dfposfun(tmp(i),eps);
00209     dftmp(i)+=dfpen*dfposfun1(tmp(i),eps);
00210     dfptmp(i)=0.0;
00211 
00212     for (k=i-1;k>=1;k--)
00213     {
00214       //tmp(i)-=L(i,k)*L(i,k);
00215       dfL(i,k)-=2.*dftmp(i)*L(i,k);
00216     }
00217     //tmp(i)=M(i,i);
00218     dfM(i,i)+=dftmp(i);
00219     dftmp(i)=0.0;
00220     for (j=i-1;j>=2;j--)
00221     {
00222       //L(i,j)=tmp1(i,j)/L(j,j);
00223       double linv=1./L(j,j);
00224       dftmp1(i,j)+=dfL(i,j)*linv;
00225       dfL(j,j)-=dfL(i,j)*tmp1(i,j)*linv*linv;
00226       dfL(i,j)=0.0;
00227       for (k=j-1;k>=1;k--)
00228       {
00229         //tmp1(i,j)-=L(i,k)*L(j,k);
00230         dfL(i,k)-=dftmp1(i,j)*L(j,k);
00231         dfL(j,k)-=dftmp1(i,j)*L(i,k);
00232       }
00233       //tmp1(i,j)=M(i,j);
00234       dfM(i,j)+=dftmp1(i,j);
00235       dftmp1(i,j)=0.0;
00236     }
00237   }
00238   double linv=1./L(1,1);
00239   for (i=n;i>=2;i--)
00240   {
00241     //L(i,1)=M(i,1)/L(1,1);
00242     dfM(i,1)+=dfL(i,1)*linv;
00243     dfL(1,1)-=dfL(i,1)*M(i,1)*linv*linv;
00244     dfL(i,1)=0.0;
00245   }
00246   //L(1,1)=sqrt(ptmp(1));
00247   dfptmp(1)+=dfL(1,1)/(2.*L(1,1));
00248   dfL(1,1)=0.0;
00249   //ptmp(1)=posfun(M(1,1),eps,fpen);
00250   dfM(1,1)+=dfptmp(1)*dfposfun(M(1,1),eps);
00251   dfM(1,1)+=dfpen*dfposfun1(M(1,1),eps);;
00252   dfptmp(1)=0.0;
00253 
00254 
00255  //*******************************************************************8
00256  //*******************************************************************8
00257  //*******************************************************************8
00258   dfM.rowshift(rowsave);
00259   dfM.colshift(colsave);
00260 
00261   dfM.save_dmatrix_derivatives(MMpos);
00262 }