ADMB Documentation  11.2.2828
 All Classes Files Functions Variables Typedefs Friends Defines
df12fun.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: df12fun.cpp 2786 2014-12-11 19:16:42Z johnoel $
00003  *
00004  * Author: David Fournier
00005  * Copyright (c) 2008-2012 Regents of the University of California
00006  */
00011 //class df1b2matrix;
00012 
00013 #include "df12fun.h"
00014 
00015 prevariable* df1_two_variable::ind_var[2];
00016 
00017 int df1_two_variable::num_ind_var=0;
00018 
00022 df1_two_variable::df1_two_variable(const df1_two_variable& x)
00023 {
00024   v[0] = x.v[0];
00025   v[1] = x.v[1];
00026   v[2] = x.v[2];
00027 }
00031 df1_two_vector::df1_two_vector(const df1_two_vector& m2)
00032 {
00033   index_min = m2.index_min;
00034   index_max = m2.index_max;
00035   shape  =m2.shape;
00036   if (shape)
00037   {
00038     (shape->ncopies)++;
00039   }
00040   v = m2.v;
00041 }
00045 df1_two_vector::~df1_two_vector()
00046 {
00047   deallocate();
00048 }
00049 
00054  void df1_two_vector::deallocate(void)
00055  {
00056    if(shape)
00057    {
00058      if (shape->ncopies)
00059      {
00060        (shape->ncopies)--;
00061      }
00062      else
00063      {
00064        v = (df1_two_variable*) (shape->trueptr);
00065        delete [] v;
00066        v = NULL;
00067        delete shape;
00068        shape=0;
00069      }
00070    }
00071  }
00072 
00077  dvector value(const df1_two_vector& v)
00078  {
00079    int mmin=v.indexmin();
00080    int mmax=v.indexmax();
00081    dvector cv(mmin,mmax);
00082    for (int i=mmin;i<=mmax;i++)
00083    {
00084      cv(i)=value(v(i));
00085    }
00086    return cv;
00087  }
00088 
00093   void df1_two_vector::initialize(void)
00094   {
00095     int mmin=indexmin();
00096     int mmax=indexmax();
00097     for (int i=mmin;i<=mmax;i++)
00098     {
00099       (*this)(i)=0.0;
00100     }
00101   }
00102 
00107   df1_two_vector::df1_two_vector(void)
00108   {
00109     allocate();
00110   }
00111 
00116   df1_two_vector::df1_two_vector(int min,int max)
00117   {
00118     allocate(min,max);
00119   }
00120 
00125   void df1_two_vector::allocate(int min,int max)
00126   {
00127     index_min=min;
00128     index_max=max;
00129     v=new df1_two_variable[max-min+1];
00130     if (v==0)
00131     {
00132       cerr << "error allocating memory in df1_two_vector" << endl;
00133       ad_exit(1);
00134     }
00135     if ( (shape=new vector_shapex(min,max,v)) == NULL)
00136     {
00137       cerr << "Error trying to allocate memory for df1_two_vector"
00138            << endl;;
00139       ad_exit(1);
00140     }
00141     v-=min;
00142   }
00143 
00148   void df1_two_vector::allocate(void)
00149   {
00150     index_min=0;
00151     index_max=-1;
00152     v=0;
00153     shape=0;
00154   }
00155 
00160  dmatrix value(const df1_two_matrix& v)
00161  {
00162    int rmin=v.indexmin();
00163    int rmax=v.indexmax();
00164    dmatrix cm(rmin,rmax);
00165    for (int i=rmin;i<=rmax;i++)
00166    {
00167      int cmin=v(i).indexmin();
00168      int cmax=v(i).indexmax();
00169      cm(i).allocate(cmin,cmax);
00170      for (int j=cmin;j<=cmax;j++)
00171      {
00172        cm(i,j)=value(v(i,j));
00173      }
00174    }
00175    return cm;
00176  }
00177 
00182  df1_two_matrix::df1_two_matrix(const df1_two_matrix& m2)
00183  {
00184    index_min=m2.index_min;
00185    index_max=m2.index_max;
00186    shape=m2.shape;
00187    if (shape)
00188    {
00189      (shape->ncopies)++;
00190    }
00191    v = m2.v;
00192  }
00193 
00198  df1_two_matrix::~df1_two_matrix()
00199  {
00200    deallocate();
00201  }
00202 
00207  void df1_two_matrix::deallocate(void)
00208  {
00209    if (shape)
00210    {
00211      if (shape->ncopies)
00212      {
00213        (shape->ncopies)--;
00214      }
00215      else
00216      {
00217        v = (df1_two_vector*) (shape->get_pointer());
00218        delete [] v;
00219        v=0;
00220        delete shape;
00221        shape=0;
00222      }
00223    }
00224  }
00225 
00230   void df1_two_matrix::initialize(void)
00231   {
00232     int mmin=indexmin();
00233     int mmax=indexmax();
00234     for (int i=mmin;i<=mmax;i++)
00235     {
00236       (*this)(i).initialize();
00237     }
00238   }
00239 
00244   df1_two_matrix::df1_two_matrix(int rmin,int rmax,int cmin,int cmax)
00245   {
00246     index_min=rmin;
00247     index_max=rmax;
00248     v=new df1_two_vector[rmax-rmin+1];
00249     if (v==0)
00250     {
00251       cerr << "error allocating memory in df1_two_matrix" << endl;
00252       ad_exit(1);
00253     }
00254     if ( (shape=new mat_shapex(v)) == NULL)
00255     {
00256       cerr << "Error trying to allocate memory for df1_two_vector"
00257            << endl;;
00258     }
00259     v-=rmin;
00260 
00261     for (int i=rmin;i<=rmax;i++)
00262     {
00263       v[i].allocate(cmin,cmax);
00264     }
00265   }
00266 
00271   df1_two_variable& df1_two_variable::operator -= (const df1_two_variable& v)
00272   {
00273     *get_u() -= *v.get_u();
00274     *get_u_x() -= *v.get_u_x();
00275     *get_u_y() -= *v.get_u_y();
00276     return *this;
00277   }
00278 
00283   df1_two_variable& df1_two_variable::operator /= (const df1_two_variable& y)
00284   {
00285    /*
00286     df1_three_variable x=*this * inv(y);
00287     *this=x;
00288     return *this;
00289    */
00290    // properly optimized code
00291     double tmp=1.0/value(y);
00292     *get_u()*=tmp;
00293     *get_u_x() = *get_u_x()*tmp- *get_u()*tmp* *y.get_u_x();
00294     *get_u_y() = *get_u_y()*tmp- *get_u()*tmp* *y.get_u_y();
00295     return *this;
00296   }
00297 
00302 df1_two_variable operator-(const df1_two_variable& v)
00303 {
00304   df1_two_variable z;
00305 
00306   *z.get_u() = -(*v.get_u());
00307   *z.get_u_x() = -(*v.get_u_x());
00308   *z.get_u_y() = -(*v.get_u_y());
00309 
00310   return z;
00311 }
00312 
00317   df1_two_variable& df1_two_variable::operator += (const df1_two_variable& v)
00318   {
00319     *get_u() += *v.get_u();
00320     *get_u_x() += *v.get_u_x();
00321     *get_u_y() += *v.get_u_y();
00322     return *this;
00323   }
00324 
00325   df1_two_variable& df1_two_variable::operator *= (double v)
00326   {
00327    /*
00328     df1_three_variable x=*this * v;
00329     *this=x;
00330     return *this;
00331    */
00332     *get_u()*=v;
00333     *get_u_x() = *get_u_x()*v;
00334     *get_u_y() = *get_u_y()*v;
00335     return *this;
00336   }
00337 
00338 
00343   df1_two_variable& df1_two_variable::operator *= (const df1_two_variable& v)
00344   {
00345     df1_two_variable x=*this * v;
00346     *this=x;
00347     return *this;
00348   }
00349 
00354   df1_two_variable& df1_two_variable::operator += (double v)
00355   {
00356     *get_u() += v;
00357 
00358     return *this;
00359   }
00360 
00365   df1_two_variable& df1_two_variable::operator -= (double v)
00366   {
00367     *get_u() -= v;
00368 
00369     return *this;
00370   }
00371 
00376 void set_derivatives( df1_two_variable& z,const df1_two_variable& x,double u,
00377   double zp)
00378 {
00379     //*z.get_u() = u;
00380     *z.get_u_x() = zp* *x.get_u_x();
00381     *z.get_u_y() = zp* *x.get_u_y();
00382 }
00383 
00388 void set_derivatives( df1_two_variable& z, const df1_two_variable& x,
00389   const df1_two_variable& y, double u,
00390   double f_u,double f_v)
00391 {
00392     *z.get_u() = u;
00393 
00394     *z.get_u_x() = f_u* *x.get_u_x()
00395                  + f_v* *y.get_u_x();
00396 
00397     *z.get_u_y() = f_u* *x.get_u_y()
00398                  + f_v* *y.get_u_y();
00399 }
00400 
00405   df1_two_variable sqrt(const df1_two_variable& x)
00406   {
00407     df1_two_variable z;
00408     double u=::sqrt(*x.get_u());
00409     *z.get_u()=u;
00410     //double xinv=1.0/(*x.get_u());
00411     double zp=0.5/u;
00412 
00413 
00414     set_derivatives(z,x,u,zp);
00415 
00416     return z;
00417   }
00418 
00423   df1_two_variable atan(const df1_two_variable& x)
00424   {
00425     df1_two_variable z;
00426     double cx=value(x);
00427     double d=1.0/(1+square(cx));
00428     //double d2=square(d);
00429     double u=::atan(cx);
00430     *z.get_u()=u;
00431     double zp=d;
00432 
00433     set_derivatives(z,x,u,zp);
00434     return z;
00435   }
00436 
00441   df1_two_variable square(const df1_two_variable& x)
00442   {
00443     df1_two_variable z;
00444     double u=value(x);
00445     *z.get_u()=u*u;
00446     double zp=2.0*u;
00447 
00448     set_derivatives(z,x,u,zp);
00449     return z;
00450   }
00451 
00456   df1_two_variable tan(const df1_two_variable& x)
00457   {
00458     df1_two_variable z;
00459     double u=::tan(*x.get_u());
00460     *z.get_u()=u;
00461     double v=1.0/::cos(*x.get_u());
00462     //double w=::sin(*x.get_u());
00463     double v2=v*v;
00464     double zp=v2;
00465 
00466     set_derivatives(z,x,u,zp);
00467     return z;
00468   }
00469 
00474   df1_two_variable sin(const df1_two_variable& x)
00475   {
00476     df1_two_variable z;
00477     double u=::sin(*x.get_u());
00478     *z.get_u()=u;
00479     double zp=::cos(*x.get_u());
00480 
00481     set_derivatives(z,x,u,zp);
00482     return z;
00483   }
00484 
00489   df1_two_variable log(const df1_two_variable& x)
00490   {
00491     df1_two_variable z;
00492     double u=::log(*x.get_u());
00493     *z.get_u()=u;
00494     double zp=1/(*x.get_u());
00495 
00496     set_derivatives(z,x,u,zp);
00497     return z;
00498   }
00499 
00504   df1_two_variable exp(const df1_two_variable& x)
00505   {
00506     df1_two_variable z;
00507     double u=::exp(*x.get_u());
00508     *z.get_u()=u;
00509     double zp=u;
00510 
00511     set_derivatives(z,x,u,zp);
00512     return z;
00513   }
00514 
00519   df1_two_variable pow(const df1_two_variable& x,const df1_two_variable& y)
00520   {
00521     return exp(y*log(x));
00522   }
00523 
00528   df1_two_variable pow(double x,const df1_two_variable& y)
00529   {
00530     return exp(y*log(x));
00531   }
00532 
00537   df1_two_variable inv(const df1_two_variable& x)
00538   {
00539     df1_two_variable z;
00540     double xinv=1.0/(*x.get_u());
00541     *z.get_u()=xinv;
00542     double zp=-xinv*xinv;
00543     set_derivatives(z,x,xinv,zp);
00544 
00545     return z;
00546   }
00547 
00552   df1_two_variable& df1_two_variable::operator = (const df1_two_variable& x)
00553   {
00554     *get_u() = *x.get_u();
00555     *get_u_x() = *x.get_u_x();
00556     *get_u_y() = *x.get_u_y();
00557     return *this;
00558   }
00559 
00564   df1_two_variable& df1_two_variable::operator = (double x)
00565   {
00566     *get_u() = x;
00567     *get_u_x() =0.0;
00568     *get_u_y() =0.0;
00569     return *this;
00570   }
00571 
00576   df1_two_variable operator * (const df1_two_variable& x,
00577     const df1_two_variable& y)
00578   {
00579     df1_two_variable z;
00580     double u= *x.get_u() * *y.get_u();
00581     *z.get_u() = u;
00582     double f_u=*y.get_u();
00583     double f_v=*x.get_u();
00584     set_derivatives(z,x,y,u, f_u, f_v);
00585     return z;
00586   }
00587 
00592   df1_two_variable operator * (double x,
00593     const df1_two_variable& y)
00594   {
00595     df1_two_variable z;
00596     *z.get_u() = x *  *y.get_u();
00597     *z.get_u_x() = x * *y.get_u_x();
00598     *z.get_u_y() = x * *y.get_u_y();
00599     return z;
00600   }
00601 
00606   df1_two_variable operator * (const df1_two_variable& y,
00607     double x)
00608   {
00609     df1_two_variable z;
00610     *z.get_u() = x *  *y.get_u();
00611     *z.get_u_x() = x * *y.get_u_x();
00612     *z.get_u_y() = x * *y.get_u_y();
00613     return z;
00614   }
00615 
00620   df1_two_variable operator / (const df1_two_variable& x,
00621     double y)
00622   {
00623     double u=1/y;
00624     return x*u;
00625   }
00626 
00631   df1_two_variable operator / (const df1_two_variable& x,
00632     const df1_two_variable& y)
00633   {
00634     df1_two_variable u=inv(y);
00635     return x*u;
00636   }
00637 
00642   df1_two_variable operator / (const double x,
00643     const df1_two_variable& y)
00644   {
00645     df1_two_variable u=inv(y);
00646     return x*u;
00647   }
00648 
00653   df1_two_variable operator + (const double x,const df1_two_variable& y)
00654   {
00655     df1_two_variable z;
00656     *z.get_u() =  x + *y.get_u();
00657     *z.get_u_x() = *y.get_u_x();
00658     *z.get_u_y() = *y.get_u_y();
00659     return z;
00660   }
00661 
00666   df1_two_variable operator + (const df1_two_variable& x,const double y)
00667   {
00668     df1_two_variable z;
00669     *z.get_u() =  *x.get_u() + y;
00670     *z.get_u_x() = *x.get_u_x();
00671     *z.get_u_y() = *x.get_u_y();
00672     return z;
00673   }
00674 
00679   df1_two_variable operator + (const df1_two_variable& x,
00680     const df1_two_variable& y)
00681   {
00682     df1_two_variable z;
00683     *z.get_u() = *x.get_u() + *y.get_u();
00684     *z.get_u_x() = *x.get_u_x() + *y.get_u_x();
00685     *z.get_u_y() = *x.get_u_y()+*y.get_u_y();
00686     return z;
00687   }
00688 
00693   df1_two_variable operator - (const df1_two_variable& x,
00694     const df1_two_variable& y)
00695   {
00696     df1_two_variable z;
00697     *z.get_u() = *x.get_u() - *y.get_u();
00698     *z.get_u_x() = *x.get_u_x()  - *y.get_u_x();
00699     *z.get_u_y() = *x.get_u_y() - *y.get_u_y();
00700     return z;
00701   }
00702 
00707   df1_two_variable operator - (double x,
00708     const df1_two_variable& y)
00709   {
00710     df1_two_variable z;
00711     *z.get_u() = x - *y.get_u();
00712     *z.get_u_x() = - *y.get_u_x();
00713     *z.get_u_y() = - *y.get_u_y();
00714     return z;
00715   }
00716 
00721   df1_two_variable operator - (const df1_two_variable& x,
00722     double y)
00723   {
00724     df1_two_variable z;
00725     *z.get_u() = *x.get_u()-y;
00726     *z.get_u_x() = *x.get_u_x();
00727     *z.get_u_y() = *x.get_u_y();
00728     return z;
00729   }
00730 
00735   init_df1_two_variable::~init_df1_two_variable()
00736   {
00737     deallocate();
00738   }
00739 
00744   void init_df1_two_variable::deallocate(void)
00745   {
00746     num_ind_var=0;
00747   }
00748 
00753 init_df1_two_variable::init_df1_two_variable(const prevariable& _v)
00754 {
00755   if (num_ind_var > 1)
00756   {
00757     cerr << "can only have 2 independent_variables in df1_two_variable"
00758        " function" << endl;
00759     ad_exit(1);
00760   }
00761   else
00762   {
00763     ADUNCONST(prevariable,v)
00764     ind_var[num_ind_var++]=&v;
00765     *get_u() =  value(v);
00766     switch(num_ind_var)
00767     {
00768     case 1:
00769       *get_u_x() = 1.0;
00770       *get_u_y() = 0.0;
00771       break;
00772     case 2:
00773       *get_u_x() = 0.0;
00774       *get_u_y() = 1.0;
00775       break;
00776     default:
00777       cerr << "illegal num_ind_var value of " << num_ind_var
00778            << " in  df1_two_variable function" << endl;
00779       ad_exit(1);
00780     }
00781   }
00782 }
00783 
00788   init_df1_two_variable::init_df1_two_variable(double v)
00789   {
00790     *get_u() =  v;
00791     *get_u_x() = 0.0;
00792     *get_u_y() = 0.0;
00793   }
00794 
00799   df1_two_variable::df1_two_variable(void)
00800   {
00801   }
00802 
00807 df1_two_matrix choleski_decomp(const df1_two_matrix& MM)
00808 {
00809   // kludge to deal with constantness
00810   df1_two_matrix & M= (df1_two_matrix &) MM;
00811   int rmin=M.indexmin();
00812   int cmin=M(rmin).indexmin();
00813   int rmax=M.indexmax();
00814   int cmax=M(rmin).indexmax();
00815   if (rmin !=1 || cmin !=1)
00816   {
00817     cerr << "minimum row and column inidices must equal 1 in "
00818       "df1b2matrix choleski_decomp(const df1_two_atrix& MM)"
00819          << endl;
00820     ad_exit(1);
00821   }
00822   if (rmax !=cmax)
00823   {
00824     cerr << "Error in df1b2matrix choleski_decomp(const df1_two_matrix& MM)"
00825       " Matrix not square" << endl;
00826     ad_exit(1);
00827   }
00828 
00829   int n=rmax-rmin+1;
00830   df1_two_matrix L(1,n,1,n);
00831 #ifndef SAFE_INITIALIZE
00832     L.initialize();
00833 #endif
00834 
00835   int i,j,k;
00836   df1_two_variable tmp;
00837 
00838     if (value(M(1,1))<=0)
00839     {
00840       cerr << "Error matrix not positive definite in choleski_decomp"
00841         <<endl;
00842       ad_exit(1);
00843     }
00844 
00845   L(1,1)=sqrt(M(1,1));
00846   for (i=2;i<=n;i++)
00847   {
00848     L(i,1)=M(i,1)/L(1,1);
00849   }
00850 
00851   for (i=2;i<=n;i++)
00852   {
00853     for (j=2;j<=i-1;j++)
00854     {
00855       tmp=M(i,j);
00856       for (k=1;k<=j-1;k++)
00857       {
00858         tmp-=L(i,k)*L(j,k);
00859       }
00860       L(i,j)=tmp/L(j,j);
00861     }
00862     tmp=M(i,i);
00863     for (k=1;k<=i-1;k++)
00864     {
00865       tmp-=L(i,k)*L(i,k);
00866     }
00867 
00868     if (value(tmp)<=0)
00869     {
00870       cerr << "Error matrix not positive definite in choleski_decomp"
00871         <<endl;
00872       ad_exit(1);
00873     }
00874 
00875     L(i,i)=sqrt(tmp);
00876   }
00877 
00878   return L;
00879 }
00880 
00881   df1_two_variable fabs(const df1_two_variable& x)
00882   {
00883     df1_two_variable z;
00884     if (value(x)>=0.0)
00885       z=x;
00886     else
00887       z=-x;
00888     return z;
00889   }