ADMB Documentation  11.6rc.3396
 All Classes Files Functions Variables Typedefs Friends Defines
df13fun.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id$
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2009-2012 ADMB Foundation
00006  *
00007  */
00013 #include "df13fun.h"
00014 
00015   prevariable * df1_three_variable::ind_var[3];
00016 
00017   int df1_three_variable::num_ind_var=0;
00018 
00023 void df1_three_variable::initialize()
00024 {
00025   v[0] = 0;
00026   v[1] = 0;
00027   v[2] = 0;
00028   v[3] = 0;
00029 }
00033 df1_three_variable::df1_three_variable()
00034 {
00035   initialize();
00036 }
00040 df1_three_variable::df1_three_variable(const df1_three_variable& x)
00041 {
00042   v[0] = x.v[0];
00043   v[1] = x.v[1];
00044   v[2] = x.v[2];
00045   v[3] = x.v[3];
00046 }
00050 df1_three_vector::df1_three_vector(const df1_three_vector& m2)
00051 {
00052   index_min=m2.index_min;
00053   index_max=m2.index_max;
00054   shape=m2.shape;
00055   if (shape)
00056   {
00057     (shape->ncopies)++;
00058   }
00059   v = m2.v;
00060 }
00061 
00065 df1_three_vector::~df1_three_vector()
00066 {
00067   deallocate();
00068 }
00069 
00074  void df1_three_vector::deallocate(void)
00075  {
00076    if(shape)
00077    {
00078      if (shape->ncopies)
00079      {
00080        (shape->ncopies)--;
00081      }
00082      else
00083      {
00084        v = (df1_three_variable*) (shape->trueptr);
00085        delete [] v;
00086        v = NULL;
00087        delete shape;
00088        shape=0;
00089      }
00090    }
00091  }
00092 
00097  dvector value(const df1_three_vector& v)
00098  {
00099    int mmin=v.indexmin();
00100    int mmax=v.indexmax();
00101    dvector cv(mmin,mmax);
00102    for (int i=mmin;i<=mmax;i++)
00103    {
00104      cv(i)=value(v(i));
00105    }
00106    return cv;
00107  }
00108 
00113   void df1_three_vector::initialize(void)
00114   {
00115     int mmin=indexmin();
00116     int mmax=indexmax();
00117     for (int i=mmin;i<=mmax;i++)
00118     {
00119       (*this)(i)=0.0;
00120     }
00121   }
00122 
00127   df1_three_vector::df1_three_vector(void)
00128   {
00129     allocate();
00130   }
00131 
00136   df1_three_vector::df1_three_vector(int min,int max)
00137   {
00138     allocate(min,max);
00139   }
00140 
00145   void df1_three_vector::allocate(int min,int max)
00146   {
00147     index_min=min;
00148     index_max=max;
00149     v=new df1_three_variable[max-min+1];
00150     if (v==0)
00151     {
00152       cerr << "error allocating memory in df1_three_vector" << endl;
00153       ad_exit(1);
00154     }
00155     if ( (shape=new vector_shapex(min,max,v)) == NULL)
00156     {
00157       cerr << "Error trying to allocate memory for df1_three_vector"
00158            << endl;;
00159       ad_exit(1);
00160     }
00161     v-=min;
00162   }
00163 
00168   void df1_three_vector::allocate(void)
00169   {
00170     index_min=0;
00171     index_max=-1;
00172     v=0;
00173     shape=0;
00174   }
00175 
00180  dmatrix value(const df1_three_matrix& v)
00181  {
00182    int rmin=v.indexmin();
00183    int rmax=v.indexmax();
00184    dmatrix cm(rmin,rmax);
00185    for (int i=rmin;i<=rmax;i++)
00186    {
00187      int cmin=v(i).indexmin();
00188      int cmax=v(i).indexmax();
00189      cm(i).allocate(cmin,cmax);
00190      for (int j=cmin;j<=cmax;j++)
00191      {
00192        cm(i,j)=value(v(i,j));
00193      }
00194    }
00195    return cm;
00196  }
00197 
00202  df1_three_matrix::df1_three_matrix(const df1_three_matrix& m2)
00203  {
00204    index_min=m2.index_min;
00205    index_max=m2.index_max;
00206    shape=m2.shape;
00207    if (shape)
00208    {
00209      (shape->ncopies)++;
00210    }
00211    v = m2.v;
00212  }
00213 
00218  df1_three_matrix::~df1_three_matrix()
00219  {
00220    deallocate();
00221  }
00222 
00227  void df1_three_matrix::deallocate(void)
00228  {
00229    if (shape)
00230    {
00231      if (shape->ncopies)
00232      {
00233        (shape->ncopies)--;
00234      }
00235      else
00236      {
00237        v = (df1_three_vector*) (shape->get_pointer());
00238        delete [] v;
00239        v=0;
00240        delete shape;
00241        shape=0;
00242      }
00243    }
00244  }
00245 
00250   void df1_three_matrix::initialize(void)
00251   {
00252     int mmin=indexmin();
00253     int mmax=indexmax();
00254     for (int i=mmin;i<=mmax;i++)
00255     {
00256       (*this)(i).initialize();
00257     }
00258   }
00259 
00264   df1_three_matrix::df1_three_matrix(int rmin,int rmax,int cmin,int cmax)
00265   {
00266     index_min=rmin;
00267     index_max=rmax;
00268     v=new df1_three_vector[rmax-rmin+1];
00269     if (v==0)
00270     {
00271       cerr << "error allocating memory in df1_three_matrix" << endl;
00272       ad_exit(1);
00273     }
00274     if ( (shape=new mat_shapex(v)) == NULL)
00275     {
00276       cerr << "Error trying to allocate memory for df1_three_vector"
00277            << endl;;
00278     }
00279     v-=rmin;
00280 
00281     for (int i=rmin;i<=rmax;i++)
00282     {
00283       v[i].allocate(cmin,cmax);
00284     }
00285   }
00286 
00291 df1_three_variable& df1_three_variable::operator-=(const df1_three_variable& _v)
00292 {
00293   *get_u() -= *_v.get_u();
00294   *get_u_x() -= *_v.get_u_x();
00295   *get_u_y() -= *_v.get_u_y();
00296   *get_u_z() -= *_v.get_u_z();
00297 
00298   return *this;
00299 }
00300 
00305 df1_three_variable operator-(const df1_three_variable& v)
00306 {
00307   df1_three_variable z;
00308 
00309   *z.get_u() = -(*v.get_u());
00310   *z.get_u_x() = -(*v.get_u_x());
00311   *z.get_u_y() = -(*v.get_u_y());
00312   *z.get_u_z() = -(*v.get_u_z());
00313 
00314   return z;
00315 }
00316 
00321 df1_three_variable& df1_three_variable::operator+=(const df1_three_variable& _v)
00322 {
00323   *get_u() += *_v.get_u();
00324   *get_u_x() += *_v.get_u_x();
00325   *get_u_y() += *_v.get_u_y();
00326   *get_u_z() += *_v.get_u_z();
00327 
00328   return *this;
00329 }
00330 
00335 df1_three_variable& df1_three_variable::operator*=(double _v)
00336 {
00337 /*
00338 df1_three_variable x=*this * v;
00339 *this=x;
00340 return *this;
00341 */
00342   *get_u() *= _v;
00343   *get_u_x() = *get_u_x() * _v;
00344   *get_u_y() = *get_u_y() * _v;
00345   *get_u_z() = *get_u_z() * _v;
00346 
00347   return *this;
00348 }
00349 
00354 df1_three_variable& df1_three_variable::operator*=(const df1_three_variable& y)
00355   {
00356    /*
00357     df1_three_variable x=*this * v;
00358     *this=x;
00359     return *this;
00360     */
00361     double tmp=value(y);
00362     *get_u_x() = *get_u_x()*tmp+ *get_u() * *y.get_u_x();
00363     *get_u_y() = *get_u_y()*tmp+ *get_u() * *y.get_u_y();
00364     *get_u_z() = *get_u_z()*tmp+ *get_u() * *y.get_u_z();
00365     // need to do this last
00366     *get_u()*=tmp;
00367     return *this;
00368   }
00369 
00374   df1_three_variable& df1_three_variable::operator /= (double y)
00375   {
00376    /*
00377     df1_three_variable x=*this * (1.0/y);
00378     *this=x;
00379     return *this;
00380    */
00381     double tmp=1.0/y;
00382     *get_u()*=tmp;
00383     *get_u_x() = *get_u_x()*tmp;
00384     *get_u_y() = *get_u_y()*tmp;
00385     *get_u_z() = *get_u_z()*tmp;
00386     return *this;
00387   }
00388 
00393   df1_three_variable& df1_three_variable::my_diveq (const df1_three_variable& y)
00394   {
00395     double tmp=1.0/value(y);
00396     *get_u()*=tmp;
00397     return *this;
00398   }
00399 
00404 df1_three_variable& df1_three_variable::operator/=(const df1_three_variable& y)
00405   {
00406    /*
00407     df1_three_variable x=*this * inv(y);
00408     *this=x;
00409     return *this;
00410    */
00411    // properly optimized code
00412     double tmp=1.0/value(y);
00413     *get_u()*=tmp;
00414     *get_u_x() = *get_u_x()*tmp- *get_u()*tmp* *y.get_u_x();
00415     *get_u_y() = *get_u_y()*tmp- *get_u()*tmp* *y.get_u_y();
00416     *get_u_z() = *get_u_z()*tmp- *get_u()*tmp* *y.get_u_z();
00417     return *this;
00418   }
00419 
00424 df1_three_variable& df1_three_variable::operator+=(double _v)
00425 {
00426   *get_u() += _v;
00427 
00428   return *this;
00429 }
00430 
00435 df1_three_variable& df1_three_variable::operator-=(double _v)
00436 {
00437   *get_u() -= _v;
00438 
00439   return *this;
00440 }
00441 
00446 void set_derivatives(df1_three_variable& z,const df1_three_variable& x,
00447   double u, double zp)
00448 {
00449     //*z.get_u() = u;
00450     *z.get_u_x() = zp* *x.get_u_x();
00451     *z.get_u_y() = zp* *x.get_u_y();
00452     *z.get_u_z() = zp* *x.get_u_z();
00453 }
00454 
00459 void set_derivatives( df1_three_variable& z, const df1_three_variable& x,
00460   const df1_three_variable& y, double u,
00461   double f_u,double f_v)
00462 {
00463     *z.get_u() = u;
00464 
00465     *z.get_u_x() = f_u* *x.get_u_x()
00466                  + f_v* *y.get_u_x();
00467 
00468     *z.get_u_y() = f_u* *x.get_u_y()
00469                  + f_v* *y.get_u_y();
00470 
00471     *z.get_u_z() = f_u* *x.get_u_z()
00472                  + f_v* *y.get_u_z();
00473 }
00474 
00479   df1_three_variable sqrt(const df1_three_variable& x)
00480   {
00481     df1_three_variable z;
00482     double u=::sqrt(*x.get_u());
00483     *z.get_u()=u;
00484     //double xinv=1.0/(*x.get_u());
00485     double zp=0.5/u;
00486 
00487     set_derivatives(z,x,u,zp);
00488 
00489     return z;
00490   }
00491 
00496   df1_three_variable atan(const df1_three_variable& x)
00497   {
00498     df1_three_variable z;
00499     double cx=value(x);
00500     double d=1.0/(1+square(cx));
00501     //double d2=square(d);
00502     double u=::atan(cx);
00503     *z.get_u()=u;
00504     double zp=d;
00505 
00506     set_derivatives(z,x,u,zp);
00507     return z;
00508   }
00509 
00514   df1_three_variable square(const df1_three_variable& x)
00515   {
00516     df1_three_variable z;
00517     double u=value(x);
00518     *z.get_u()=u*u;
00519     double zp=2.0*u;
00520 
00521     set_derivatives(z,x,u,zp);
00522     return z;
00523   }
00524 
00529   df1_three_variable tan(const df1_three_variable& x)
00530   {
00531     df1_three_variable z;
00532     double u=::tan(*x.get_u());
00533     *z.get_u()=u;
00534     double v=1.0/::cos(*x.get_u());
00535     //double w=::sin(*x.get_u());
00536     double v2=v*v;
00537     double zp=v2;
00538 
00539     set_derivatives(z,x,u,zp);
00540     return z;
00541   }
00542 
00547   df1_three_variable sin(const df1_three_variable& x)
00548   {
00549     df1_three_variable z;
00550     double u=::sin(*x.get_u());
00551     *z.get_u()=u;
00552     double zp=::cos(*x.get_u());
00553 
00554     set_derivatives(z,x,u,zp);
00555     return z;
00556   }
00557 
00562   df1_three_variable fabs(const df1_three_variable& x)
00563   {
00564     df1_three_variable z;
00565     if (value(x)>=0.0)
00566       z=x;
00567     else
00568       z=-x;
00569     return z;
00570   }
00571 
00576   df1_three_variable log(const df1_three_variable& x)
00577   {
00578     df1_three_variable z;
00579     double u=::log(*x.get_u());
00580     *z.get_u()=u;
00581     double zp=1/(*x.get_u());
00582 
00583     set_derivatives(z,x,u,zp);
00584     return z;
00585   }
00586 
00591   df1_three_variable exp(const df1_three_variable& x)
00592   {
00593     df1_three_variable z;
00594     double u=::exp(*x.get_u());
00595     *z.get_u()=u;
00596     double zp=u;
00597 
00598     set_derivatives(z,x,u,zp);
00599     return z;
00600   }
00601 
00606   df1_three_variable inv(const df1_three_variable& x)
00607   {
00608     df1_three_variable z;
00609     double xinv=1.0/(*x.get_u());
00610     *z.get_u()=xinv;
00611     double zp=-xinv*xinv;
00612     set_derivatives(z,x,xinv,zp);
00613 
00614     return z;
00615   }
00616 
00621 df1_three_variable& df1_three_variable::operator=(const df1_three_variable& x)
00622   {
00623     *get_u() = *x.get_u();
00624     *get_u_x() = *x.get_u_x();
00625     *get_u_y() = *x.get_u_y();
00626     *get_u_z() = *x.get_u_z();
00627     return *this;
00628   }
00629 
00634   df1_three_variable& df1_three_variable::operator = (double x)
00635   {
00636     *get_u() = x;
00637     *get_u_x() =0.0;
00638     *get_u_y() =0.0;
00639     *get_u_z() =0.0;
00640     return *this;
00641   }
00642 
00647   df1_three_variable operator * (const df1_three_variable& x,
00648     const df1_three_variable& y)
00649   {
00650     df1_three_variable z;
00651     double u= *x.get_u() * *y.get_u();
00652     *z.get_u() = u;
00653     double f_u=*y.get_u();
00654     double f_v=*x.get_u();
00655     set_derivatives(z,x,y,u, f_u, f_v);
00656     return z;
00657   }
00658 
00663   df1_three_variable operator * (double x,
00664     const df1_three_variable& y)
00665   {
00666     df1_three_variable z;
00667     *z.get_u() = x *  *y.get_u();
00668     *z.get_u_x() = x * *y.get_u_x();
00669     *z.get_u_y() = x * *y.get_u_y();
00670     *z.get_u_z() = x * *y.get_u_z();
00671     return z;
00672   }
00673 
00678   df1_three_variable operator * (const df1_three_variable& y,
00679     double x)
00680   {
00681     df1_three_variable z;
00682     *z.get_u() = x *  *y.get_u();
00683     *z.get_u_x() = x * *y.get_u_x();
00684     *z.get_u_y() = x * *y.get_u_y();
00685     *z.get_u_z() = x * *y.get_u_z();
00686     return z;
00687   }
00688 
00693   df1_three_variable operator / (const df1_three_variable& x,
00694     double y)
00695   {
00696     double u=1/y;
00697     return x*u;
00698   }
00699 
00704   df1_three_variable operator / (const df1_three_variable& x,
00705     const df1_three_variable& y)
00706   {
00707     df1_three_variable u=inv(y);
00708     return x*u;
00709   }
00710 
00715   df1_three_variable operator / (const double x,
00716     const df1_three_variable& y)
00717   {
00718     df1_three_variable u=inv(y);
00719     return x*u;
00720   }
00721 
00726 df1_three_variable pow(const df1_three_variable& x,const df1_three_variable& y)
00727   {
00728     return exp(y*log(x));
00729   }
00730 
00735   df1_three_variable operator + (const double x,const df1_three_variable& y)
00736   {
00737     df1_three_variable z;
00738     *z.get_u() =  x + *y.get_u();
00739     *z.get_u_x() = *y.get_u_x();
00740     *z.get_u_y() = *y.get_u_y();
00741     *z.get_u_z() = *y.get_u_z();
00742     return z;
00743   }
00744 
00749   df1_three_variable operator + (const df1_three_variable& x,const double y)
00750   {
00751     df1_three_variable z;
00752     *z.get_u() =  *x.get_u() + y;
00753     *z.get_u_x() = *x.get_u_x();
00754     *z.get_u_y() = *x.get_u_y();
00755     *z.get_u_z() = *x.get_u_z();
00756     return z;
00757   }
00758 
00763   df1_three_variable operator + (const df1_three_variable& x,
00764     const df1_three_variable& y)
00765   {
00766     df1_three_variable z;
00767     *z.get_u() = *x.get_u() + *y.get_u();
00768     *z.get_u_x() = *x.get_u_x() + *y.get_u_x();
00769     *z.get_u_y() = *x.get_u_y()+*y.get_u_y();
00770     *z.get_u_z() = *x.get_u_z()+*y.get_u_z();
00771     return z;
00772   }
00773 
00778   df1_three_variable operator - (const df1_three_variable& x,
00779     const df1_three_variable& y)
00780   {
00781     df1_three_variable z;
00782     *z.get_u() = *x.get_u() - *y.get_u();
00783     *z.get_u_x() = *x.get_u_x()  - *y.get_u_x();
00784     *z.get_u_y() = *x.get_u_y() - *y.get_u_y();
00785     *z.get_u_z() = *x.get_u_z() - *y.get_u_z();
00786     return z;
00787   }
00788 
00793   df1_three_variable operator - (double x,
00794     const df1_three_variable& y)
00795   {
00796     df1_three_variable z;
00797     *z.get_u() = x - *y.get_u();
00798     *z.get_u_x() = - *y.get_u_x();
00799     *z.get_u_y() = - *y.get_u_y();
00800     *z.get_u_z() = - *y.get_u_z();
00801     return z;
00802   }
00803 
00808   df1_three_variable operator - (const df1_three_variable& x,
00809     double y)
00810   {
00811     df1_three_variable z;
00812     *z.get_u() = *x.get_u()-y;
00813     *z.get_u_x() = *x.get_u_x();
00814     *z.get_u_y() = *x.get_u_y();
00815     *z.get_u_z() = *x.get_u_z();
00816     return z;
00817   }
00818 
00819   df1_three_variable operator - (const df1_three_variable& x,
00820     const df1_three_variable& y);
00821   df1_three_variable operator / (const df1_three_variable& x,
00822     const df1_three_variable& y);
00823   df1_three_variable operator * (const df1_three_variable& x,
00824     const df1_three_variable& y);
00825 
00830   init_df1_three_variable::~init_df1_three_variable()
00831   {
00832     deallocate();
00833   }
00834 
00839   void init_df1_three_variable::deallocate(void)
00840   {
00841     num_ind_var=0;
00842   }
00843 
00848 init_df1_three_variable::init_df1_three_variable(const prevariable& _v)
00849 {
00850   if (num_ind_var > 2)
00851   {
00852     cerr << "can only have 2 independent_variables in df1_three_variable"
00853        " function" << endl;
00854     ad_exit(1);
00855   }
00856   else
00857   {
00858     ADUNCONST(prevariable,v)
00859     ind_var[num_ind_var++]=&v;
00860     *get_u() =  value(v);
00861     switch(num_ind_var)
00862     {
00863     case 1:
00864       *get_u_x() = 1.0;
00865       *get_u_y() = 0.0;
00866       *get_u_z() = 0.0;
00867       break;
00868     case 2:
00869       *get_u_x() = 0.0;
00870       *get_u_y() = 1.0;
00871       *get_u_z() = 0.0;
00872       break;
00873     case 3:
00874       *get_u_x() = 0.0;
00875       *get_u_y() = 0.0;
00876       *get_u_z() = 1.0;
00877       break;
00878     default:
00879       cerr << "illegal num_ind_var value of " << num_ind_var
00880            << " in  df1_three_variable function" << endl;
00881       ad_exit(1);
00882     }
00883   }
00884 }
00885 
00890   init_df1_three_variable::init_df1_three_variable(double v)
00891   {
00892     *get_u() =  v;
00893     *get_u_x() = 0.0;
00894     *get_u_y() = 0.0;
00895     *get_u_z() = 0.0;
00896   }
00897 
00902 df1_three_matrix choleski_decomp(const df1_three_matrix& MM)
00903 {
00904   // kludge to deal with constantness
00905   df1_three_matrix & M= (df1_three_matrix &) MM;
00906   int rmin=M.indexmin();
00907   int cmin=M(rmin).indexmin();
00908   int rmax=M.indexmax();
00909   int cmax=M(rmin).indexmax();
00910   if (rmin !=1 || cmin !=1)
00911   {
00912     cerr << "minimum row and column inidices must equal 1 in "
00913       "df1b2matrix choleski_decomp(const df1_three_atrix& MM)"
00914          << endl;
00915     ad_exit(1);
00916   }
00917   if (rmax !=cmax)
00918   {
00919     cerr << "Error in df1b2matrix choleski_decomp(const df1_three_matrix& MM)"
00920       " Matrix not square" << endl;
00921     ad_exit(1);
00922   }
00923 
00924   int n=rmax-rmin+1;
00925   df1_three_matrix L(1,n,1,n);
00926 #ifndef SAFE_INITIALIZE
00927     L.initialize();
00928 #endif
00929 
00930   int i,j,k;
00931   df1_three_variable tmp;
00932 
00933     if (value(M(1,1))<=0)
00934     {
00935       cerr << "Error matrix not positive definite in choleski_decomp"
00936         <<endl;
00937       ad_exit(1);
00938     }
00939 
00940   L(1,1)=sqrt(M(1,1));
00941   for (i=2;i<=n;i++)
00942   {
00943     L(i,1)=M(i,1)/L(1,1);
00944   }
00945 
00946   for (i=2;i<=n;i++)
00947   {
00948     for (j=2;j<=i-1;j++)
00949     {
00950       tmp=M(i,j);
00951       for (k=1;k<=j-1;k++)
00952       {
00953         tmp-=L(i,k)*L(j,k);
00954       }
00955       L(i,j)=tmp/L(j,j);
00956     }
00957     tmp=M(i,i);
00958     for (k=1;k<=i-1;k++)
00959     {
00960       tmp-=L(i,k)*L(i,k);
00961     }
00962 
00963     if (value(tmp)<=0)
00964     {
00965       cerr << "Error matrix not positive definite in choleski_decomp"
00966         <<endl;
00967       ad_exit(1);
00968     }
00969 
00970     L(i,i)=sqrt(tmp);
00971   }
00972 
00973   return L;
00974 }
00975 
00980 dvariable& dvariable::operator = (const df1_three_variable& v)
00981 {
00982   const prevariable * px=df1_three_variable::ind_var[0];
00983   const prevariable * py=df1_three_variable::ind_var[1];
00984   const prevariable * pz=df1_three_variable::ind_var[2];
00985   double  dfx= *v.get_u_x();
00986   double  dfy= *v.get_u_y();
00987   double  dfz= *v.get_u_z();
00988   value(*this)=*v.get_u();
00989 
00990   gradient_structure::GRAD_STACK1->set_gradient_stack(default_evaluation3ind,
00991     &(value(*this)),&(value(*px)),dfx,&(value(*py)),dfy,&(value(*pz)),
00992     dfz);
00993 
00994   return *this;
00995 }