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