ADMB Documentation  11.1x.2711
 All Classes Files Functions Variables Typedefs Friends Defines
df13fun.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df13fun.cpp 1711 2014-02-28 22:46:44Z johnoel $
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(void)
00024   {
00025     v[0]=0;
00026     v[1]=0;
00027     v[2]=0;
00028     v[3]=0;
00029   }
00030 
00035   df1_three_variable::df1_three_variable(const df1_three_variable& x)
00036   {
00037     v[0]=x.v[0];
00038     v[1]=x.v[1];
00039     v[2]=x.v[2];
00040     v[3]=x.v[3];
00041   }
00042 
00043   df1_three_variable::df1_three_variable(void){}
00044 
00049  df1_three_vector::df1_three_vector(const df1_three_vector& m2)
00050  {
00051    index_min=m2.index_min;
00052    index_max=m2.index_max;
00053    shape=m2.shape;
00054    if (shape)
00055    {
00056      (shape->ncopies)++;
00057    }
00058    v = m2.v;
00059  }
00060 
00064 df1_three_vector::~df1_three_vector()
00065 {
00066   deallocate();
00067 }
00068 
00073  void df1_three_vector::deallocate(void)
00074  {
00075    if(shape)
00076    {
00077      if (shape->ncopies)
00078      {
00079        (shape->ncopies)--;
00080      }
00081      else
00082      {
00083        v = (df1_three_variable*) (shape->trueptr);
00084        delete [] v;
00085        v = NULL;
00086        delete shape;
00087        shape=0;
00088      }
00089    }
00090  }
00091 
00096  dvector value(const df1_three_vector& v)
00097  {
00098    int mmin=v.indexmin();
00099    int mmax=v.indexmax();
00100    dvector cv(mmin,mmax);
00101    for (int i=mmin;i<=mmax;i++)
00102    {
00103      cv(i)=value(v(i));
00104    }
00105    return cv;
00106  }
00107 
00112   void df1_three_vector::initialize(void)
00113   {
00114     int mmin=indexmin();
00115     int mmax=indexmax();
00116     for (int i=mmin;i<=mmax;i++)
00117     {
00118       (*this)(i)=0.0;
00119     }
00120   }
00121 
00126   df1_three_vector::df1_three_vector(void)
00127   {
00128     allocate();
00129   }
00130 
00135   df1_three_vector::df1_three_vector(int min,int max)
00136   {
00137     allocate(min,max);
00138   }
00139 
00144   void df1_three_vector::allocate(int min,int max)
00145   {
00146     index_min=min;
00147     index_max=max;
00148     v=new df1_three_variable[max-min+1];
00149     if (v==0)
00150     {
00151       cerr << "error allocating memory in df1_three_vector" << endl;
00152       ad_exit(1);
00153     }
00154     if ( (shape=new vector_shapex(min,max,v)) == NULL)
00155     {
00156       cerr << "Error trying to allocate memory for df1_three_vector"
00157            << endl;;
00158       ad_exit(1);
00159     }
00160     v-=min;
00161   }
00162 
00167   void df1_three_vector::allocate(void)
00168   {
00169     index_min=0;
00170     index_max=-1;
00171     v=0;
00172     shape=0;
00173   }
00174 
00179  dmatrix value(const df1_three_matrix& v)
00180  {
00181    int rmin=v.indexmin();
00182    int rmax=v.indexmax();
00183    dmatrix cm(rmin,rmax);
00184    for (int i=rmin;i<=rmax;i++)
00185    {
00186      int cmin=v(i).indexmin();
00187      int cmax=v(i).indexmax();
00188      cm(i).allocate(cmin,cmax);
00189      for (int j=cmin;j<=cmax;j++)
00190      {
00191        cm(i,j)=value(v(i,j));
00192      }
00193    }
00194    return cm;
00195  }
00196 
00201  df1_three_matrix::df1_three_matrix(const df1_three_matrix& m2)
00202  {
00203    index_min=m2.index_min;
00204    index_max=m2.index_max;
00205    shape=m2.shape;
00206    if (shape)
00207    {
00208      (shape->ncopies)++;
00209    }
00210    v = m2.v;
00211  }
00212 
00217  df1_three_matrix::~df1_three_matrix()
00218  {
00219    deallocate();
00220  }
00221 
00226  void df1_three_matrix::deallocate(void)
00227  {
00228    if (shape)
00229    {
00230      if (shape->ncopies)
00231      {
00232        (shape->ncopies)--;
00233      }
00234      else
00235      {
00236        v = (df1_three_vector*) (shape->get_pointer());
00237        delete [] v;
00238        v=0;
00239        delete shape;
00240        shape=0;
00241      }
00242    }
00243  }
00244 
00249   void df1_three_matrix::initialize(void)
00250   {
00251     int mmin=indexmin();
00252     int mmax=indexmax();
00253     for (int i=mmin;i<=mmax;i++)
00254     {
00255       (*this)(i).initialize();
00256     }
00257   }
00258 
00263   df1_three_matrix::df1_three_matrix(int rmin,int rmax,int cmin,int cmax)
00264   {
00265     index_min=rmin;
00266     index_max=rmax;
00267     v=new df1_three_vector[rmax-rmin+1];
00268     if (v==0)
00269     {
00270       cerr << "error allocating memory in df1_three_matrix" << endl;
00271       ad_exit(1);
00272     }
00273     if ( (shape=new mat_shapex(v)) == NULL)
00274     {
00275       cerr << "Error trying to allocate memory for df1_three_vector"
00276            << endl;;
00277     }
00278     v-=rmin;
00279 
00280     for (int i=rmin;i<=rmax;i++)
00281     {
00282       v[i].allocate(cmin,cmax);
00283     }
00284   }
00285 
00290 df1_three_variable& df1_three_variable::operator-=(const df1_three_variable& v)
00291   {
00292     *get_u() -= *v.get_u();
00293     *get_u_x() -= *v.get_u_x();
00294     *get_u_y() -= *v.get_u_y();
00295     *get_u_z() -= *v.get_u_z();
00296     return *this;
00297   }
00298 
00303 df1_three_variable operator-(const df1_three_variable& v)
00304 {
00305   df1_three_variable z;
00306 
00307   *z.get_u() = -(*v.get_u());
00308   *z.get_u_x() = -(*v.get_u_x());
00309   *z.get_u_y() = -(*v.get_u_y());
00310   *z.get_u_z() = -(*v.get_u_z());
00311 
00312   return z;
00313 }
00314 
00319 df1_three_variable& df1_three_variable::operator+=(const df1_three_variable& v)
00320   {
00321     *get_u() += *v.get_u();
00322     *get_u_x() += *v.get_u_x();
00323     *get_u_y() += *v.get_u_y();
00324     *get_u_z() += *v.get_u_z();
00325     return *this;
00326   }
00327 
00332   df1_three_variable& df1_three_variable::operator *= (double v)
00333   {
00334    /*
00335     df1_three_variable x=*this * v;
00336     *this=x;
00337     return *this;
00338    */
00339     *get_u()*=v;
00340     *get_u_x() = *get_u_x()*v;
00341     *get_u_y() = *get_u_y()*v;
00342     *get_u_z() = *get_u_z()*v;
00343     return *this;
00344   }
00345 
00350 df1_three_variable& df1_three_variable::operator*=(const df1_three_variable& y)
00351   {
00352    /*
00353     df1_three_variable x=*this * v;
00354     *this=x;
00355     return *this;
00356     */
00357     double tmp=value(y);
00358     *get_u_x() = *get_u_x()*tmp+ *get_u() * *y.get_u_x();
00359     *get_u_y() = *get_u_y()*tmp+ *get_u() * *y.get_u_y();
00360     *get_u_z() = *get_u_z()*tmp+ *get_u() * *y.get_u_z();
00361     // need to do this last
00362     *get_u()*=tmp;
00363     return *this;
00364   }
00365 
00370   df1_three_variable& df1_three_variable::operator /= (double y)
00371   {
00372    /*
00373     df1_three_variable x=*this * (1.0/y);
00374     *this=x;
00375     return *this;
00376    */
00377     double tmp=1.0/y;
00378     *get_u()*=tmp;
00379     *get_u_x() = *get_u_x()*tmp;
00380     *get_u_y() = *get_u_y()*tmp;
00381     *get_u_z() = *get_u_z()*tmp;
00382     return *this;
00383   }
00384 
00389   df1_three_variable& df1_three_variable::my_diveq (const df1_three_variable& y)
00390   {
00391     double tmp=1.0/value(y);
00392     *get_u()*=tmp;
00393     return *this;
00394   }
00395 
00400 df1_three_variable& df1_three_variable::operator/=(const df1_three_variable& y)
00401   {
00402    /*
00403     df1_three_variable x=*this * inv(y);
00404     *this=x;
00405     return *this;
00406    */
00407    // properly optimized code
00408     double tmp=1.0/value(y);
00409     *get_u()*=tmp;
00410     *get_u_x() = *get_u_x()*tmp- *get_u()*tmp* *y.get_u_x();
00411     *get_u_y() = *get_u_y()*tmp- *get_u()*tmp* *y.get_u_y();
00412     *get_u_z() = *get_u_z()*tmp- *get_u()*tmp* *y.get_u_z();
00413     return *this;
00414   }
00415 
00420   df1_three_variable& df1_three_variable::operator += (double v)
00421   {
00422     *get_u() += v;
00423 
00424     return *this;
00425   }
00426 
00431   df1_three_variable& df1_three_variable::operator -= (double v)
00432   {
00433     *get_u() -= v;
00434 
00435     return *this;
00436   }
00437 
00442 void set_derivatives(df1_three_variable& z,const df1_three_variable& x,
00443   double u, double zp)
00444 {
00445     //*z.get_u() = u;
00446     *z.get_u_x() = zp* *x.get_u_x();
00447     *z.get_u_y() = zp* *x.get_u_y();
00448     *z.get_u_z() = zp* *x.get_u_z();
00449 }
00450 
00455 void set_derivatives( df1_three_variable& z, const df1_three_variable& x,
00456   const df1_three_variable& y, double u,
00457   double f_u,double f_v)
00458 {
00459     *z.get_u() = u;
00460 
00461     *z.get_u_x() = f_u* *x.get_u_x()
00462                  + f_v* *y.get_u_x();
00463 
00464     *z.get_u_y() = f_u* *x.get_u_y()
00465                  + f_v* *y.get_u_y();
00466 
00467     *z.get_u_z() = f_u* *x.get_u_z()
00468                  + f_v* *y.get_u_z();
00469 }
00470 
00475   df1_three_variable sqrt(const df1_three_variable& x)
00476   {
00477     df1_three_variable z;
00478     double u=::sqrt(*x.get_u());
00479     *z.get_u()=u;
00480     //double xinv=1.0/(*x.get_u());
00481     double zp=0.5/u;
00482 
00483     set_derivatives(z,x,u,zp);
00484 
00485     return z;
00486   }
00487 
00492   df1_three_variable atan(const df1_three_variable& x)
00493   {
00494     df1_three_variable z;
00495     double cx=value(x);
00496     double d=1.0/(1+square(cx));
00497     //double d2=square(d);
00498     double u=::atan(cx);
00499     *z.get_u()=u;
00500     double zp=d;
00501 
00502     set_derivatives(z,x,u,zp);
00503     return z;
00504   }
00505 
00510   df1_three_variable square(const df1_three_variable& x)
00511   {
00512     df1_three_variable z;
00513     double u=value(x);
00514     *z.get_u()=u*u;
00515     double zp=2.0*u;
00516 
00517     set_derivatives(z,x,u,zp);
00518     return z;
00519   }
00520 
00525   df1_three_variable tan(const df1_three_variable& x)
00526   {
00527     df1_three_variable z;
00528     double u=::tan(*x.get_u());
00529     *z.get_u()=u;
00530     double v=1.0/::cos(*x.get_u());
00531     //double w=::sin(*x.get_u());
00532     double v2=v*v;
00533     double zp=v2;
00534 
00535     set_derivatives(z,x,u,zp);
00536     return z;
00537   }
00538 
00543   df1_three_variable sin(const df1_three_variable& x)
00544   {
00545     df1_three_variable z;
00546     double u=::sin(*x.get_u());
00547     *z.get_u()=u;
00548     double zp=::cos(*x.get_u());
00549 
00550     set_derivatives(z,x,u,zp);
00551     return z;
00552   }
00553 
00558   df1_three_variable fabs(const df1_three_variable& x)
00559   {
00560     df1_three_variable z;
00561     if (value(x)>=0.0)
00562       z=x;
00563     else
00564       z=-x;
00565     return z;
00566   }
00567 
00572   df1_three_variable log(const df1_three_variable& x)
00573   {
00574     df1_three_variable z;
00575     double u=::log(*x.get_u());
00576     *z.get_u()=u;
00577     double zp=1/(*x.get_u());
00578 
00579     set_derivatives(z,x,u,zp);
00580     return z;
00581   }
00582 
00587   df1_three_variable exp(const df1_three_variable& x)
00588   {
00589     df1_three_variable z;
00590     double u=::exp(*x.get_u());
00591     *z.get_u()=u;
00592     double zp=u;
00593 
00594     set_derivatives(z,x,u,zp);
00595     return z;
00596   }
00597 
00602   df1_three_variable inv(const df1_three_variable& x)
00603   {
00604     df1_three_variable z;
00605     double xinv=1.0/(*x.get_u());
00606     *z.get_u()=xinv;
00607     double zp=-xinv*xinv;
00608     set_derivatives(z,x,xinv,zp);
00609 
00610     return z;
00611   }
00612 
00617 df1_three_variable& df1_three_variable::operator=(const df1_three_variable& x)
00618   {
00619     *get_u() = *x.get_u();
00620     *get_u_x() = *x.get_u_x();
00621     *get_u_y() = *x.get_u_y();
00622     *get_u_z() = *x.get_u_z();
00623     return *this;
00624   }
00625 
00630   df1_three_variable& df1_three_variable::operator = (double x)
00631   {
00632     *get_u() = x;
00633     *get_u_x() =0.0;
00634     *get_u_y() =0.0;
00635     *get_u_z() =0.0;
00636     return *this;
00637   }
00638 
00643   df1_three_variable operator * (const df1_three_variable& x,
00644     const df1_three_variable& y)
00645   {
00646     df1_three_variable z;
00647     double u= *x.get_u() * *y.get_u();
00648     *z.get_u() = u;
00649     double f_u=*y.get_u();
00650     double f_v=*x.get_u();
00651     set_derivatives(z,x,y,u, f_u, f_v);
00652     return z;
00653   }
00654 
00659   df1_three_variable operator * (double x,
00660     const df1_three_variable& y)
00661   {
00662     df1_three_variable z;
00663     *z.get_u() = x *  *y.get_u();
00664     *z.get_u_x() = x * *y.get_u_x();
00665     *z.get_u_y() = x * *y.get_u_y();
00666     *z.get_u_z() = x * *y.get_u_z();
00667     return z;
00668   }
00669 
00674   df1_three_variable operator * (const df1_three_variable& y,
00675     double x)
00676   {
00677     df1_three_variable z;
00678     *z.get_u() = x *  *y.get_u();
00679     *z.get_u_x() = x * *y.get_u_x();
00680     *z.get_u_y() = x * *y.get_u_y();
00681     *z.get_u_z() = x * *y.get_u_z();
00682     return z;
00683   }
00684 
00689   df1_three_variable operator / (const df1_three_variable& x,
00690     double y)
00691   {
00692     double u=1/y;
00693     return x*u;
00694   }
00695 
00700   df1_three_variable operator / (const df1_three_variable& x,
00701     const df1_three_variable& y)
00702   {
00703     df1_three_variable u=inv(y);
00704     return x*u;
00705   }
00706 
00711   df1_three_variable operator / (const double x,
00712     const df1_three_variable& y)
00713   {
00714     df1_three_variable u=inv(y);
00715     return x*u;
00716   }
00717 
00722 df1_three_variable pow(const df1_three_variable& x,const df1_three_variable& y)
00723   {
00724     return exp(y*log(x));
00725   }
00726 
00731   df1_three_variable operator + (const double x,const df1_three_variable& y)
00732   {
00733     df1_three_variable z;
00734     *z.get_u() =  x + *y.get_u();
00735     *z.get_u_x() = *y.get_u_x();
00736     *z.get_u_y() = *y.get_u_y();
00737     *z.get_u_z() = *y.get_u_z();
00738     return z;
00739   }
00740 
00745   df1_three_variable operator + (const df1_three_variable& x,const double y)
00746   {
00747     df1_three_variable z;
00748     *z.get_u() =  *x.get_u() + y;
00749     *z.get_u_x() = *x.get_u_x();
00750     *z.get_u_y() = *x.get_u_y();
00751     *z.get_u_z() = *x.get_u_z();
00752     return z;
00753   }
00754 
00759   df1_three_variable operator + (const df1_three_variable& x,
00760     const df1_three_variable& y)
00761   {
00762     df1_three_variable z;
00763     *z.get_u() = *x.get_u() + *y.get_u();
00764     *z.get_u_x() = *x.get_u_x() + *y.get_u_x();
00765     *z.get_u_y() = *x.get_u_y()+*y.get_u_y();
00766     *z.get_u_z() = *x.get_u_z()+*y.get_u_z();
00767     return z;
00768   }
00769 
00774   df1_three_variable operator - (const df1_three_variable& x,
00775     const df1_three_variable& y)
00776   {
00777     df1_three_variable z;
00778     *z.get_u() = *x.get_u() - *y.get_u();
00779     *z.get_u_x() = *x.get_u_x()  - *y.get_u_x();
00780     *z.get_u_y() = *x.get_u_y() - *y.get_u_y();
00781     *z.get_u_z() = *x.get_u_z() - *y.get_u_z();
00782     return z;
00783   }
00784 
00789   df1_three_variable operator - (double x,
00790     const df1_three_variable& y)
00791   {
00792     df1_three_variable z;
00793     *z.get_u() = x - *y.get_u();
00794     *z.get_u_x() = - *y.get_u_x();
00795     *z.get_u_y() = - *y.get_u_y();
00796     *z.get_u_z() = - *y.get_u_z();
00797     return z;
00798   }
00799 
00804   df1_three_variable operator - (const df1_three_variable& x,
00805     double y)
00806   {
00807     df1_three_variable z;
00808     *z.get_u() = *x.get_u()-y;
00809     *z.get_u_x() = *x.get_u_x();
00810     *z.get_u_y() = *x.get_u_y();
00811     *z.get_u_z() = *x.get_u_z();
00812     return z;
00813   }
00814 
00815   df1_three_variable operator - (const df1_three_variable& x,
00816     const df1_three_variable& y);
00817   df1_three_variable operator / (const df1_three_variable& x,
00818     const df1_three_variable& y);
00819   df1_three_variable operator * (const df1_three_variable& x,
00820     const df1_three_variable& y);
00821 
00826   init_df1_three_variable::~init_df1_three_variable()
00827   {
00828     deallocate();
00829   }
00830 
00835   void init_df1_three_variable::deallocate(void)
00836   {
00837     num_ind_var=0;
00838   }
00839 
00844   init_df1_three_variable::init_df1_three_variable(const prevariable& _v)
00845   {
00846     ADUNCONST(prevariable,v)
00847     if (num_ind_var>2)
00848     {
00849       cerr << "can only have 2 independent_variables in df1_three_variable"
00850        " function" << endl;
00851       ad_exit(1);
00852     }
00853     ind_var[num_ind_var++]=&v;
00854     *get_u() =  value(v);
00855     switch(num_ind_var)
00856     {
00857     case 1:
00858       *get_u_x() = 1.0;
00859       *get_u_y() = 0.0;
00860       *get_u_z() = 0.0;
00861       break;
00862     case 2:
00863       *get_u_x() = 0.0;
00864       *get_u_y() = 1.0;
00865       *get_u_z() = 0.0;
00866       break;
00867     case 3:
00868       *get_u_x() = 0.0;
00869       *get_u_y() = 0.0;
00870       *get_u_z() = 1.0;
00871       break;
00872     default:
00873       cerr << "illegal num_ind_var value of " << num_ind_var
00874            << " in  df1_three_variable function" << endl;
00875       ad_exit(1);
00876     }
00877   }
00878 
00883   init_df1_three_variable::init_df1_three_variable(double v)
00884   {
00885     *get_u() =  v;
00886     *get_u_x() = 0.0;
00887     *get_u_y() = 0.0;
00888     *get_u_z() = 0.0;
00889   }
00890 
00895 df1_three_matrix choleski_decomp(const df1_three_matrix& MM)
00896 {
00897   // kludge to deal with constantness
00898   df1_three_matrix & M= (df1_three_matrix &) MM;
00899   int rmin=M.indexmin();
00900   int cmin=M(rmin).indexmin();
00901   int rmax=M.indexmax();
00902   int cmax=M(rmin).indexmax();
00903   if (rmin !=1 || cmin !=1)
00904   {
00905     cerr << "minimum row and column inidices must equal 1 in "
00906       "df1b2matrix choleski_decomp(const df1_three_atrix& MM)"
00907          << endl;
00908     ad_exit(1);
00909   }
00910   if (rmax !=cmax)
00911   {
00912     cerr << "Error in df1b2matrix choleski_decomp(const df1_three_matrix& MM)"
00913       " Matrix not square" << endl;
00914     ad_exit(1);
00915   }
00916 
00917   int n=rmax-rmin+1;
00918   df1_three_matrix L(1,n,1,n);
00919 #ifndef SAFE_INITIALIZE
00920     L.initialize();
00921 #endif
00922 
00923   int i,j,k;
00924   df1_three_variable tmp;
00925 
00926     if (value(M(1,1))<=0)
00927     {
00928       cerr << "Error matrix not positive definite in choleski_decomp"
00929         <<endl;
00930       ad_exit(1);
00931     }
00932 
00933   L(1,1)=sqrt(M(1,1));
00934   for (i=2;i<=n;i++)
00935   {
00936     L(i,1)=M(i,1)/L(1,1);
00937   }
00938 
00939   for (i=2;i<=n;i++)
00940   {
00941     for (j=2;j<=i-1;j++)
00942     {
00943       tmp=M(i,j);
00944       for (k=1;k<=j-1;k++)
00945       {
00946         tmp-=L(i,k)*L(j,k);
00947       }
00948       L(i,j)=tmp/L(j,j);
00949     }
00950     tmp=M(i,i);
00951     for (k=1;k<=i-1;k++)
00952     {
00953       tmp-=L(i,k)*L(i,k);
00954     }
00955 
00956     if (value(tmp)<=0)
00957     {
00958       cerr << "Error matrix not positive definite in choleski_decomp"
00959         <<endl;
00960       ad_exit(1);
00961     }
00962 
00963     L(i,i)=sqrt(tmp);
00964   }
00965 
00966   return L;
00967 }
00968 
00973 dvariable& dvariable::operator = (const df1_three_variable& v)
00974 {
00975   const prevariable * px=df1_three_variable::ind_var[0];
00976   const prevariable * py=df1_three_variable::ind_var[1];
00977   const prevariable * pz=df1_three_variable::ind_var[2];
00978   double  dfx= *v.get_u_x();
00979   double  dfy= *v.get_u_y();
00980   double  dfz= *v.get_u_z();
00981   value(*this)=*v.get_u();
00982 
00983   gradient_structure::GRAD_STACK1->set_gradient_stack(default_evaluation3ind,
00984     &(value(*this)),&(value(*px)),dfx,&(value(*py)),dfy,&(value(*pz)),
00985     dfz);
00986 
00987   return *this;
00988 }