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