ADMB Documentation  11.5.3260
 All Classes Files Functions Variables Typedefs Friends Defines
df11fun.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id$
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00007 
00008 char df12fun_notice[50]="copyright (c) 2006 otter research ltd";
00009 
00010 #include "df11fun.h"
00011 
00012 prevariable* df1_one_variable::ind_var[1] = { NULL };
00013 int df1_one_variable::num_ind_var = 0;
00014 
00018 df1_one_variable::df1_one_variable()
00019 {
00020   v[0] = 0;
00021   v[1] = 0;
00022 }
00026 df1_one_variable::df1_one_variable(const df1_one_variable& x)
00027 {
00028   v[0] = x.v[0];
00029   v[1] = x.v[1];
00030 }
00034 df1_one_vector::df1_one_vector(const df1_one_vector& m2)
00035 {
00036   index_min = m2.index_min;
00037   index_max = m2.index_max;
00038   shape = m2.shape;
00039   if (shape)
00040   {
00041     (shape->ncopies)++;
00042   }
00043   v = m2.v;
00044 }
00048 df1_one_vector::~df1_one_vector()
00049 {
00050   deallocate();
00051 }
00052 
00053  void df1_one_vector::deallocate(void)
00054  {
00055    if(shape)
00056    {
00057      if (shape->ncopies)
00058      {
00059        (shape->ncopies)--;
00060      }
00061      else
00062      {
00063        v = (df1_one_variable*) (shape->trueptr);
00064        delete [] v;
00065        v = NULL;
00066        delete shape;
00067        shape=0;
00068      }
00069    }
00070  }
00071 
00072  dvector value(const df1_one_vector& v)
00073  {
00074    int mmin=v.indexmin();
00075    int mmax=v.indexmax();
00076    dvector cv(mmin,mmax);
00077    for (int i=mmin;i<=mmax;i++)
00078    {
00079      cv(i)=value(v(i));
00080    }
00081    return cv;
00082  }
00083 
00084   void df1_one_vector::initialize(void)
00085   {
00086     int mmin=indexmin();
00087     int mmax=indexmax();
00088     for (int i=mmin;i<=mmax;i++)
00089     {
00090       (*this)(i)=0.0;
00091     }
00092   }
00093 
00094   df1_one_vector::df1_one_vector(void)
00095   {
00096     allocate();
00097   }
00098 
00099   df1_one_vector::df1_one_vector(int min,int max)
00100   {
00101     allocate(min,max);
00102   }
00103 
00104   void df1_one_vector::allocate(int min,int max)
00105   {
00106     index_min=min;
00107     index_max=max;
00108     v=new df1_one_variable[max-min+1];
00109     if (v==0)
00110     {
00111       cerr << "error allocating memory in df1_one_vector" << endl;
00112       ad_exit(1);
00113     }
00114     if ( (shape=new vector_shapex(min,max,v)) == NULL)
00115     {
00116       cerr << "Error trying to allocate memory for df1_one_vector"
00117            << endl;;
00118       ad_exit(1);
00119     }
00120     v-=min;
00121   }
00122 
00123   void df1_one_vector::allocate(void)
00124   {
00125     index_min=0;
00126     index_max=-1;
00127     v=0;
00128     shape=0;
00129   }
00130 
00131  dmatrix value(const df1_one_matrix& v)
00132  {
00133    int rmin=v.indexmin();
00134    int rmax=v.indexmax();
00135    dmatrix cm(rmin,rmax);
00136    for (int i=rmin;i<=rmax;i++)
00137    {
00138      int cmin=v(i).indexmin();
00139      int cmax=v(i).indexmax();
00140      cm(i).allocate(cmin,cmax);
00141      for (int j=cmin;j<=cmax;j++)
00142      {
00143        cm(i,j)=value(v(i,j));
00144      }
00145    }
00146    return cm;
00147  }
00148 
00149  df1_one_matrix::df1_one_matrix(const df1_one_matrix& m2)
00150  {
00151    index_min=m2.index_min;
00152    index_max=m2.index_max;
00153    shape=m2.shape;
00154    if (shape)
00155    {
00156      (shape->ncopies)++;
00157    }
00158    v = m2.v;
00159  }
00160 
00161  df1_one_matrix::~df1_one_matrix()
00162  {
00163    deallocate();
00164  }
00165 
00166  void df1_one_matrix::deallocate(void)
00167  {
00168    if (shape)
00169    {
00170      if (shape->ncopies)
00171      {
00172        (shape->ncopies)--;
00173      }
00174      else
00175      {
00176        v = (df1_one_vector*) (shape->get_pointer());
00177        delete [] v;
00178        v=0;
00179        delete shape;
00180        shape=0;
00181      }
00182    }
00183  }
00184 
00185 
00186   void df1_one_matrix::initialize(void)
00187   {
00188     int mmin=indexmin();
00189     int mmax=indexmax();
00190     for (int i=mmin;i<=mmax;i++)
00191     {
00192       (*this)(i).initialize();
00193     }
00194   }
00195 
00196 
00197   df1_one_matrix::df1_one_matrix(int rmin,int rmax,int cmin,int cmax)
00198   {
00199     index_min=rmin;
00200     index_max=rmax;
00201     v=new df1_one_vector[rmax-rmin+1];
00202     if (v==0)
00203     {
00204       cerr << "error allocating memory in df1_one_matrix" << endl;
00205       ad_exit(1);
00206     }
00207     if ( (shape=new mat_shapex(v)) == NULL)
00208     {
00209       cerr << "Error trying to allocate memory for df1_one_vector"
00210            << endl;;
00211     }
00212     v-=rmin;
00213 
00214     for (int i=rmin;i<=rmax;i++)
00215     {
00216       v[i].allocate(cmin,cmax);
00217     }
00218   }
00219 
00220 df1_one_variable& df1_one_variable::operator-=(const df1_one_variable& _v)
00221 {
00222   *get_u() -= *_v.get_u();
00223   *get_u_x() -= *_v.get_u_x();
00224 
00225   return *this;
00226 }
00227 df1_one_variable operator-(const df1_one_variable& v)
00228 {
00229   df1_one_variable z;
00230 
00231   *z.get_u() = -(*v.get_u());
00232   *z.get_u_x() = -(*v.get_u_x());
00233 
00234   return z;
00235 }
00236 df1_one_variable& df1_one_variable::operator+=(const df1_one_variable& _v)
00237 {
00238   *get_u() += *_v.get_u();
00239   *get_u_x() += *_v.get_u_x();
00240 
00241   return *this;
00242 }
00243 df1_one_variable& df1_one_variable::operator*=(const df1_one_variable& _v)
00244 {
00245   df1_one_variable x = *this * _v;
00246   *this = x;
00247 
00248   return *this;
00249 }
00250 df1_one_variable& df1_one_variable::operator+=(double _v)
00251 {
00252   *get_u() += _v;
00253 
00254   return *this;
00255 }
00256 df1_one_variable& df1_one_variable::operator-=(double _v)
00257 {
00258   *get_u() -= _v;
00259 
00260   return *this;
00261 }
00262 
00263 void set_derivatives( df1_one_variable& z,const df1_one_variable& x,double u,
00264   double zp)
00265 {
00266     //*z.get_u() = u;
00267     *z.get_u_x() = zp* *x.get_u_x();
00268 }
00269 
00270 void set_derivatives( df1_one_variable& z, const df1_one_variable& x,
00271   const df1_one_variable& y, double u,
00272   double f_u,double f_v)
00273 {
00274     *z.get_u() = u;
00275 
00276     *z.get_u_x() = f_u* *x.get_u_x()
00277                  + f_v* *y.get_u_x();
00278 }
00279 
00280   df1_one_variable sqrt(const df1_one_variable& x)
00281   {
00282     df1_one_variable z;
00283     double u=::sqrt(*x.get_u());
00284     *z.get_u()=u;
00285     //double xinv=1.0/(*x.get_u());
00286     double zp=0.5/u;
00287 
00288     set_derivatives(z,x,u,zp);
00289 
00290     return z;
00291   }
00292 
00293 
00294 
00295   df1_one_variable atan(const df1_one_variable& x)
00296   {
00297     df1_one_variable z;
00298     double cx=value(x);
00299     double d=1.0/(1+square(cx));
00300     //double d2=square(d);
00301     double u=::atan(cx);
00302     *z.get_u()=u;
00303     double zp=d;
00304 
00305     set_derivatives(z,x,u,zp);
00306     return z;
00307   }
00308 
00309   df1_one_variable square(const df1_one_variable& x)
00310   {
00311     df1_one_variable z;
00312     double u=value(x);
00313     *z.get_u()=u*u;
00314     double zp=2.0*u;
00315 
00316     set_derivatives(z,x,u,zp);
00317     return z;
00318   }
00319 
00320   df1_one_variable tan(const df1_one_variable& x)
00321   {
00322     df1_one_variable z;
00323     double u=::tan(*x.get_u());
00324     *z.get_u()=u;
00325     double v=1.0/::cos(*x.get_u());
00326     //double w=::sin(*x.get_u());
00327     double v2=v*v;
00328     double zp=v2;
00329 
00330     set_derivatives(z,x,u,zp);
00331     return z;
00332   }
00333 
00334   df1_one_variable sin(const df1_one_variable& x)
00335   {
00336     df1_one_variable z;
00337     double u=::sin(*x.get_u());
00338     *z.get_u()=u;
00339     double zp=::cos(*x.get_u());
00340 
00341     set_derivatives(z,x,u,zp);
00342     return z;
00343   }
00344 
00345   df1_one_variable log(const df1_one_variable& x)
00346   {
00347     df1_one_variable z;
00348     double u=::log(*x.get_u());
00349     *z.get_u()=u;
00350     double zp=1/(*x.get_u());
00351 
00352     set_derivatives(z,x,u,zp);
00353     return z;
00354   }
00355 
00356   df1_one_variable exp(const df1_one_variable& x)
00357   {
00358     df1_one_variable z;
00359     double u=::exp(*x.get_u());
00360     *z.get_u()=u;
00361     double zp=u;
00362 
00363     set_derivatives(z,x,u,zp);
00364     return z;
00365   }
00366 
00367   df1_one_variable inv(const df1_one_variable& x)
00368   {
00369     df1_one_variable z;
00370     double xinv=1.0/(*x.get_u());
00371     *z.get_u()=xinv;
00372     double zp=-xinv*xinv;
00373     set_derivatives(z,x,xinv,zp);
00374 
00375     return z;
00376   }
00377 
00382  dvector first_derivatives(const df1_one_vector& v)
00383  {
00384    int mmin=v.indexmin();
00385    int mmax=v.indexmax();
00386    dvector cv(mmin,mmax);
00387    for (int i=mmin;i<=mmax;i++)
00388    {
00389      cv(i)=*(v(i).get_udot());
00390    }
00391    return cv;
00392  }
00393 
00398  dmatrix first_derivatives(const df1_one_matrix& v)
00399  {
00400    int rmin=v.indexmin();
00401    int rmax=v.indexmax();
00402    dmatrix cm(rmin,rmax);
00403    for (int i=rmin;i<=rmax;i++)
00404    {
00405      int cmin=v(i).indexmin();
00406      int cmax=v(i).indexmax();
00407      cm(i).allocate(cmin,cmax);
00408      cm(i)=first_derivatives(v(i));
00409    }
00410    return cm;
00411  }
00412 
00413   df1_one_variable& df1_one_variable::operator = (const df1_one_variable& x)
00414   {
00415     *get_u() = *x.get_u();
00416     *get_u_x() = *x.get_u_x();
00417     return *this;
00418   }
00419 
00420   df1_one_variable& df1_one_variable::operator = (double x)
00421   {
00422     *get_u() = x;
00423     *get_u_x() =0.0;
00424     return *this;
00425   }
00426 
00427 
00428 df1_one_variable operator*(
00429   const df1_one_variable& x,
00430   const df1_one_variable& y)
00431 {
00432   df1_one_variable z;
00433   double u= *x.get_u() * *y.get_u();
00434   *z.get_u() = u;
00435   double f_u=*y.get_u();
00436   double f_v=*x.get_u();
00437   set_derivatives(z,x,y,u, f_u, f_v);
00438   return z;
00439 }
00440 
00441   df1_one_variable operator * (double x,
00442     const df1_one_variable& y)
00443   {
00444     df1_one_variable z;
00445     *z.get_u() = x *  *y.get_u();
00446     *z.get_u_x() = x * *y.get_u_x();
00447     return z;
00448   }
00449 
00450 
00451   df1_one_variable operator * (const df1_one_variable& y,
00452     double x)
00453   {
00454     df1_one_variable z;
00455     *z.get_u() = x *  *y.get_u();
00456     *z.get_u_x() = x * *y.get_u_x();
00457     return z;
00458   }
00459 
00460 
00461 
00462   df1_one_variable operator / (const df1_one_variable& x,
00463     double y)
00464   {
00465     double u=1/y;
00466     return x*u;
00467   }
00468 
00469   df1_one_variable operator / (const df1_one_variable& x,
00470     const df1_one_variable& y)
00471   {
00472     df1_one_variable u=inv(y);
00473     return x*u;
00474   }
00475 
00476   df1_one_variable operator / (const double x,
00477     const df1_one_variable& y)
00478   {
00479     df1_one_variable u=inv(y);
00480     return x*u;
00481   }
00482 
00483 
00484   df1_one_variable operator + (const double x,const df1_one_variable& y)
00485   {
00486     df1_one_variable z;
00487     *z.get_u() =  x + *y.get_u();
00488     *z.get_u_x() = *y.get_u_x();
00489     return z;
00490   }
00491 
00492   df1_one_variable operator + (const df1_one_variable& x,const double y)
00493   {
00494     df1_one_variable z;
00495     *z.get_u() =  *x.get_u() + y;
00496     *z.get_u_x() = *x.get_u_x();
00497     return z;
00498   }
00499 
00500 
00501 
00502   df1_one_variable operator + (const df1_one_variable& x,
00503     const df1_one_variable& y)
00504   {
00505     df1_one_variable z;
00506     *z.get_u() = *x.get_u() + *y.get_u();
00507     *z.get_u_x() = *x.get_u_x() + *y.get_u_x();
00508     return z;
00509   }
00510 
00511   df1_one_variable operator - (const df1_one_variable& x,
00512     const df1_one_variable& y)
00513   {
00514     df1_one_variable z;
00515     *z.get_u() = *x.get_u() - *y.get_u();
00516     *z.get_u_x() = *x.get_u_x()  - *y.get_u_x();
00517     return z;
00518   }
00519 
00520   df1_one_variable operator - (double x,
00521     const df1_one_variable& y)
00522   {
00523     df1_one_variable z;
00524     *z.get_u() = x - *y.get_u();
00525     *z.get_u_x() = - *y.get_u_x();
00526     return z;
00527   }
00528 
00529   df1_one_variable operator - (const df1_one_variable& x,
00530     double y)
00531   {
00532     df1_one_variable z;
00533     *z.get_u() = *x.get_u()-y;
00534     *z.get_u_x() = *x.get_u_x();
00535     return z;
00536   }
00537 
00541 init_df1_one_variable::~init_df1_one_variable()
00542 {
00543   deallocate();
00544 }
00548 void init_df1_one_variable::deallocate(void)
00549 {
00550   num_ind_var=0;
00551 }
00552 
00553 init_df1_one_variable::init_df1_one_variable(const prevariable& _v)
00554 {
00555   ADUNCONST(prevariable,v)
00556   if (num_ind_var > 0)
00557   {
00558     cerr << "can only have 1 independent_variables in df1_one_variable"
00559        " function" << endl;
00560     ad_exit(1);
00561   }
00562   else
00563   {
00564     ind_var[num_ind_var++]=&v;
00565     *get_u() =  value(v);
00566     switch(num_ind_var)
00567     {
00568     case 1:
00569       *get_u_x() = 1.0;
00570       break;
00571     default:
00572       cerr << "illegal num_ind_var value of " << num_ind_var
00573            << " in  df1_one_variable function" << endl;
00574       ad_exit(1);
00575     }
00576   }
00577 }
00578 
00579   init_df1_one_variable::init_df1_one_variable(double v)
00580   {
00581     *get_u() =  v;
00582     *get_u_x() = 0.0;
00583   }
00584 
00585 df1_one_matrix choleski_decomp(const df1_one_matrix& MM)
00586 {
00587   // kludge to deal with constantness
00588   df1_one_matrix & M= (df1_one_matrix &) MM;
00589   int rmin=M.indexmin();
00590   int cmin=M(rmin).indexmin();
00591   int rmax=M.indexmax();
00592   int cmax=M(rmin).indexmax();
00593   if (rmin !=1 || cmin !=1)
00594   {
00595     cerr << "minimum row and column inidices must equal 1 in "
00596       "df1b2matrix choleski_decomp(const df1_one_atrix& MM)"
00597          << endl;
00598     ad_exit(1);
00599   }
00600   if (rmax !=cmax)
00601   {
00602     cerr << "Error in df1b2matrix choleski_decomp(const df1_one_matrix& MM)"
00603       " Matrix not square" << endl;
00604     ad_exit(1);
00605   }
00606 
00607   int n=rmax-rmin+1;
00608   df1_one_matrix L(1,n,1,n);
00609 #ifndef SAFE_INITIALIZE
00610     L.initialize();
00611 #endif
00612 
00613   int i,j,k;
00614   df1_one_variable tmp;
00615 
00616     if (value(M(1,1))<=0)
00617     {
00618       cerr << "Error matrix not positive definite in choleski_decomp"
00619         <<endl;
00620       ad_exit(1);
00621     }
00622 
00623   L(1,1)=sqrt(M(1,1));
00624   for (i=2;i<=n;i++)
00625   {
00626     L(i,1)=M(i,1)/L(1,1);
00627   }
00628 
00629   for (i=2;i<=n;i++)
00630   {
00631     for (j=2;j<=i-1;j++)
00632     {
00633       tmp=M(i,j);
00634       for (k=1;k<=j-1;k++)
00635       {
00636         tmp-=L(i,k)*L(j,k);
00637       }
00638       L(i,j)=tmp/L(j,j);
00639     }
00640     tmp=M(i,i);
00641     for (k=1;k<=i-1;k++)
00642     {
00643       tmp-=L(i,k)*L(i,k);
00644     }
00645 
00646     if (value(tmp)<=0)
00647     {
00648       cerr << "Error matrix not positive definite in choleski_decomp"
00649         <<endl;
00650       ad_exit(1);
00651     }
00652 
00653     L(i,i)=sqrt(tmp);
00654   }
00655 
00656   return L;
00657 }