ADMB Documentation  11.5.3152
 All Classes Files Functions Variables Typedefs Friends Defines
fvar_m42.cpp
Go to the documentation of this file.
00001 
00007 #include <fvar.hpp>
00008 #ifdef __TURBOC__
00009   #pragma hdrstop
00010   #include <iomanip.h>
00011 #endif
00012 
00013 #ifdef __ZTC__
00014   #include <iomanip.hpp>
00015 #endif
00016 
00017 #define TINY 1.0e-20;
00018 void dfinvpret(void);
00019 
00020 // int min(int a,int b);
00021 void df_xldet(void);
00022 
00023 #if defined(max)
00024 #undef max
00025 #endif
00026 #if defined(min)
00027 #undef min
00028 #endif
00029 
00030 dvariable ln_det(const dvar_matrix& aa, int& sgn);
00031 
00032 dvariable ln_det(const dvar_matrix& a)
00033 {
00034   int sgn = 1;
00035   return ln_det(a, sgn);
00036 }
00042 dvariable ln_det(const dvar_matrix& aa, int& sgn)
00043 {
00044   int errflag=0;
00045   int i,j,k,n;
00046   n=aa.colsize();
00047   int lb=aa.colmin();
00048   int ub=aa.colmax();
00049   ivector indx(lb,ub);
00050   if (lb!=aa.rowmin()||ub!=aa.colmax())
00051   {
00052     cerr << "Error matrix not square in det()"<<endl;
00053     ad_exit(1);
00054   }
00055   int One=1;
00056   indx.fill_seqadd(lb,One);
00057   double ld;
00058   double big,dum,sum,temp;
00059   dvar_matrix_position dmp(aa,1);
00060   dmatrix bb=value(aa);
00061   dvector vv(lb,ub);
00062   dvector part_prod(lb,ub);
00063 
00064   ld=0.0;
00065   for (i=lb;i<=ub;i++)
00066   {
00067     big=0.0;
00068     for (j=lb;j<=ub;j++)
00069     {
00070       temp=fabs(bb.elem(i,j));
00071       if (temp > big)
00072       {
00073         big=temp;
00074       }
00075     }
00076     if (big == 0.0)
00077     {
00078       cerr << "Error in matrix inverse -- matrix singular in "
00079       "inv(dvar_matrix)\n";
00080       big=1.e+10;
00081       errflag=1;
00082     }
00083     vv[i]=1.0/big;
00084   }
00085 
00086   for (j=lb;j<=ub;j++)
00087   {
00088     for (i=lb;i<j;i++)
00089     {
00090       sum=bb.elem(i,j);
00091       for (k=lb;k<i;k++)
00092       {
00093         sum = sum - bb.elem(i,k)*bb.elem(k,j);
00094       }
00095       //a[i][j]=sum;
00096       bb(i,j)=sum;
00097     }
00098     int imax = j;
00099     big=0.0;
00100     for (i=j;i<=ub;i++)
00101     {
00102       sum=bb.elem(i,j);
00103       for (k=lb;k<j;k++)
00104       {
00105         sum = sum - bb(i,k)*bb(k,j);
00106       }
00107       bb(i,j)=sum;
00108       dum=vv.elem(i)*fabs(sum);
00109       if ( dum >= big)
00110       {
00111         big=dum;
00112         imax=i;
00113       }
00114     }
00115     if (j != imax)
00116     {
00117       for (k=lb;k<=ub;k++)
00118       {
00119         dum=bb.elem(imax,k);
00120         bb.elem(imax,k)=bb.elem(j,k);
00121         bb.elem(j,k)=dum;
00122       }
00123       //d = -1.*d;
00124       sgn=-1*sgn;
00125       vv.elem(imax)=vv.elem(j);
00126 
00127       //if (j<ub)
00128       {
00129         int itemp=indx.elem(imax);
00130         indx.elem(imax)=indx.elem(j);
00131         indx.elem(j)=itemp;
00132       }
00133       //cout << "indx= " <<indx<<endl;
00134     }
00135 
00136     if (bb.elem(j,j) == 0.0)
00137     {
00138       bb(j,j)=TINY;
00139     }
00140 
00141     if (j != n)
00142     {
00143       dum=1.0/bb(j,j);
00144       for (i=j+1;i<=ub;i++)
00145       {
00146         bb.elem(i,j) *= dum;
00147       }
00148     }
00149   }
00150   if (bb(1,1)>0)
00151     part_prod(1)=ld+log(bb(1,1));
00152   else
00153   {
00154     part_prod(1)=ld+log(-bb(1,1));
00155     sgn=-sgn;
00156   }
00157   for (j=lb+1;j<=ub;j++)
00158   {
00159     if (bb(j,j)>0)
00160       part_prod(j)=part_prod(j-1)+log(bb(j,j));
00161     else
00162     {
00163       part_prod(j)=part_prod(j-1)+log(-bb(j,j));
00164       sgn=-sgn;
00165     }
00166   }
00167   double ldet=part_prod(ub);
00168   dvariable rdet=nograd_assign(ldet);
00169   save_identifier_string("PLACE7");
00170   part_prod.save_dvector_value();
00171   part_prod.save_dvector_position();
00172   indx.save_ivector_value();
00173   indx.save_ivector_position();
00174   save_identifier_string("PLACE3");
00175   aa.save_dvar_matrix_position();
00176   rdet.save_prevariable_position();
00177   bb.save_dmatrix_value();
00178   save_identifier_string("PLACE2");
00179   bb.save_dmatrix_position();
00180   save_identifier_string("PLACE1");
00181   save_double_value(ld);
00182   save_identifier_string("PLACE0");
00183   gradient_structure::GRAD_STACK1->
00184       set_gradient_stack(df_xldet);
00185   if (errflag) sgn=-1;
00186   return rdet;
00187 }
00188 
00190 void df_xldet(void)
00191 {
00192   verify_identifier_string("PLACE0");
00193   /*double ld=*/restore_double_value();
00194   verify_identifier_string("PLACE1");
00195   dmatrix_position bpos=restore_dmatrix_position();
00196   verify_identifier_string("PLACE2");
00197   dmatrix b=restore_dmatrix_value(bpos);
00198   //dvar_matrix_position rdet_pos=restore_prevariable_position();
00199   double dfdet=restore_prevariable_derivative();
00200   dvar_matrix_position a_pos=restore_dvar_matrix_position();
00201   verify_identifier_string("PLACE3");
00202   ivector_position indx_pos=restore_ivector_position();
00203   ivector indx=restore_ivector_value(indx_pos);
00204   dvector_position part_prod_pos=restore_dvector_position();
00205   dvector part_prod=restore_dvector_value(part_prod_pos);
00206   verify_identifier_string("PLACE7");
00207   int lb=b.colmin();
00208   int ub=b.colmax();
00209   dmatrix dfb(lb,ub,lb,ub);
00210 
00211   dvector dfpart_prod(lb,ub);
00212 
00213   #ifndef SAFE_INITIALIZE
00214     dfb.initialize();
00215     dfpart_prod.initialize();
00216   #endif
00217 
00218 
00219   dfpart_prod(ub)=dfdet;
00220   int j;
00221   for (j=ub;j>=lb+1;j--)
00222   {
00223     if (b(j,j)>0)
00224     {
00225       // part_prod(j)=part_prod(j-1)+log(b(j,j));
00226       dfpart_prod(j-1)+=dfpart_prod(j);
00227       dfb(j,j)+=dfpart_prod(j)/b(j,j);
00228     }
00229     else
00230     {
00231       // part_prod(j)=part_prod(j-1)+log(-b(j,j));
00232       dfpart_prod(j-1)+=dfpart_prod(j);
00233       dfb(j,j)+=dfpart_prod(j)/b(j,j);
00234     }
00235     dfpart_prod(j)=0.;
00236   }
00237   //part_prod(1)=ld+log(b(lb,lb));
00238   dfb(lb,lb)+=dfpart_prod(lb)/b(lb,lb);
00239   dfpart_prod(lb)=0.;
00240 
00241   double dfsum=0.;
00242   for (j=ub;j>=lb;j--)
00243   {
00244     for (int i=ub;i>=lb;i--)
00245     {
00246       if (i<=j)
00247       {
00248         // b(i,j)=sum;
00249         dfsum+=dfb(i,j);
00250         dfb(i,j)=0.;
00251       }
00252       else
00253       {
00254         // b(i,j)=sum/b(j,j);
00255         dfsum+=dfb(i,j)/b(j,j);
00256         dfb(j,j)-=dfb(i,j)*b(i,j)/b(j,j);
00257         dfb(i,j)=0.;
00258       }
00259 
00260       for (int k=min(i-1,j-1);k>=lb;k--)
00261       {
00262         // sum-=b(i,k)*b(k,j);
00263         dfb(i,k)-=dfsum*b(k,j);
00264         dfb(k,j)-=dfsum*b(i,k);
00265       }
00266       // sum=value(a(indx(i),j);
00267       save_dmatrix_derivatives(a_pos,dfsum,indx(i),j); // like this
00268       dfsum=0.;
00269     }
00270   }
00271 }
00272 
00273 #undef TINY