ADMB Documentation  11.5.3297
 All Classes Files Functions Variables Typedefs Friends Defines
df12fun.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  */
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()
00023 {
00024   v[0] = 0;
00025   v[1] = 0;
00026   v[2] = 0;
00027 }
00028 
00032 df1_two_variable::df1_two_variable(const df1_two_variable& x)
00033 {
00034   v[0] = x.v[0];
00035   v[1] = x.v[1];
00036   v[2] = x.v[2];
00037 }
00041 df1_two_vector::df1_two_vector(const df1_two_vector& m2)
00042 {
00043   index_min = m2.index_min;
00044   index_max = m2.index_max;
00045   shape  =m2.shape;
00046   if (shape)
00047   {
00048     (shape->ncopies)++;
00049   }
00050   v = m2.v;
00051 }
00055 df1_two_vector::~df1_two_vector()
00056 {
00057   deallocate();
00058 }
00059 
00064  void df1_two_vector::deallocate(void)
00065  {
00066    if(shape)
00067    {
00068      if (shape->ncopies)
00069      {
00070        (shape->ncopies)--;
00071      }
00072      else
00073      {
00074        v = (df1_two_variable*) (shape->trueptr);
00075        delete [] v;
00076        v = NULL;
00077        delete shape;
00078        shape=0;
00079      }
00080    }
00081  }
00082 
00087  dvector value(const df1_two_vector& v)
00088  {
00089    int mmin=v.indexmin();
00090    int mmax=v.indexmax();
00091    dvector cv(mmin,mmax);
00092    for (int i=mmin;i<=mmax;i++)
00093    {
00094      cv(i)=value(v(i));
00095    }
00096    return cv;
00097  }
00098 
00103   void df1_two_vector::initialize(void)
00104   {
00105     int mmin=indexmin();
00106     int mmax=indexmax();
00107     for (int i=mmin;i<=mmax;i++)
00108     {
00109       (*this)(i)=0.0;
00110     }
00111   }
00112 
00117   df1_two_vector::df1_two_vector(void)
00118   {
00119     allocate();
00120   }
00121 
00126   df1_two_vector::df1_two_vector(int min,int max)
00127   {
00128     allocate(min,max);
00129   }
00130 
00135   void df1_two_vector::allocate(int min,int max)
00136   {
00137     index_min=min;
00138     index_max=max;
00139     v=new df1_two_variable[max-min+1];
00140     if (v==0)
00141     {
00142       cerr << "error allocating memory in df1_two_vector" << endl;
00143       ad_exit(1);
00144     }
00145     if ( (shape=new vector_shapex(min,max,v)) == NULL)
00146     {
00147       cerr << "Error trying to allocate memory for df1_two_vector"
00148            << endl;;
00149       ad_exit(1);
00150     }
00151     v-=min;
00152   }
00153 
00158   void df1_two_vector::allocate(void)
00159   {
00160     index_min=0;
00161     index_max=-1;
00162     v=0;
00163     shape=0;
00164   }
00165 
00170  dmatrix value(const df1_two_matrix& v)
00171  {
00172    int rmin=v.indexmin();
00173    int rmax=v.indexmax();
00174    dmatrix cm(rmin,rmax);
00175    for (int i=rmin;i<=rmax;i++)
00176    {
00177      int cmin=v(i).indexmin();
00178      int cmax=v(i).indexmax();
00179      cm(i).allocate(cmin,cmax);
00180      for (int j=cmin;j<=cmax;j++)
00181      {
00182        cm(i,j)=value(v(i,j));
00183      }
00184    }
00185    return cm;
00186  }
00187 
00192  df1_two_matrix::df1_two_matrix(const df1_two_matrix& m2)
00193  {
00194    index_min=m2.index_min;
00195    index_max=m2.index_max;
00196    shape=m2.shape;
00197    if (shape)
00198    {
00199      (shape->ncopies)++;
00200    }
00201    v = m2.v;
00202  }
00203 
00208  df1_two_matrix::~df1_two_matrix()
00209  {
00210    deallocate();
00211  }
00212 
00217  void df1_two_matrix::deallocate(void)
00218  {
00219    if (shape)
00220    {
00221      if (shape->ncopies)
00222      {
00223        (shape->ncopies)--;
00224      }
00225      else
00226      {
00227        v = (df1_two_vector*) (shape->get_pointer());
00228        delete [] v;
00229        v=0;
00230        delete shape;
00231        shape=0;
00232      }
00233    }
00234  }
00235 
00240   void df1_two_matrix::initialize(void)
00241   {
00242     int mmin=indexmin();
00243     int mmax=indexmax();
00244     for (int i=mmin;i<=mmax;i++)
00245     {
00246       (*this)(i).initialize();
00247     }
00248   }
00249 
00254   df1_two_matrix::df1_two_matrix(int rmin,int rmax,int cmin,int cmax)
00255   {
00256     index_min=rmin;
00257     index_max=rmax;
00258     v=new df1_two_vector[rmax-rmin+1];
00259     if (v==0)
00260     {
00261       cerr << "error allocating memory in df1_two_matrix" << endl;
00262       ad_exit(1);
00263     }
00264     if ( (shape=new mat_shapex(v)) == NULL)
00265     {
00266       cerr << "Error trying to allocate memory for df1_two_vector"
00267            << endl;;
00268     }
00269     v-=rmin;
00270 
00271     for (int i=rmin;i<=rmax;i++)
00272     {
00273       v[i].allocate(cmin,cmax);
00274     }
00275   }
00276 
00281 df1_two_variable& df1_two_variable::operator-=(const df1_two_variable& _v)
00282 {
00283   *get_u() -= *_v.get_u();
00284   *get_u_x() -= *_v.get_u_x();
00285   *get_u_y() -= *_v.get_u_y();
00286 
00287   return *this;
00288 }
00289 
00294   df1_two_variable& df1_two_variable::operator /= (const df1_two_variable& y)
00295   {
00296    /*
00297     df1_three_variable x=*this * inv(y);
00298     *this=x;
00299     return *this;
00300    */
00301    // properly optimized code
00302     double tmp=1.0/value(y);
00303     *get_u()*=tmp;
00304     *get_u_x() = *get_u_x()*tmp- *get_u()*tmp* *y.get_u_x();
00305     *get_u_y() = *get_u_y()*tmp- *get_u()*tmp* *y.get_u_y();
00306     return *this;
00307   }
00308 
00313 df1_two_variable operator-(const df1_two_variable& v)
00314 {
00315   df1_two_variable z;
00316 
00317   *z.get_u() = -(*v.get_u());
00318   *z.get_u_x() = -(*v.get_u_x());
00319   *z.get_u_y() = -(*v.get_u_y());
00320 
00321   return z;
00322 }
00323 
00328 df1_two_variable& df1_two_variable::operator += (const df1_two_variable& _v)
00329 {
00330   *get_u() += *_v.get_u();
00331   *get_u_x() += *_v.get_u_x();
00332   *get_u_y() += *_v.get_u_y();
00333   return *this;
00334 }
00335 
00339 df1_two_variable& df1_two_variable::operator*=(double _v)
00340 {
00341 /*
00342   df1_three_variable x=*this * v;
00343   *this=x;
00344   return *this;
00345 */
00346   *get_u() *= _v;
00347   *get_u_x() *= _v;
00348   *get_u_y() *= _v;
00349 
00350   return *this;
00351 }
00356   df1_two_variable& df1_two_variable::operator *= (const df1_two_variable& v)
00357   {
00358     df1_two_variable x=*this * v;
00359     *this=x;
00360     return *this;
00361   }
00362 
00367 df1_two_variable& df1_two_variable::operator+=(double _v)
00368 {
00369   *get_u() += _v;
00370 
00371   return *this;
00372 }
00373 
00378 df1_two_variable& df1_two_variable::operator-=(double _v)
00379 {
00380   *get_u() -= _v;
00381 
00382   return *this;
00383 }
00384 
00389 void set_derivatives( df1_two_variable& z,const df1_two_variable& x,double u,
00390   double zp)
00391 {
00392     //*z.get_u() = u;
00393     *z.get_u_x() = zp* *x.get_u_x();
00394     *z.get_u_y() = zp* *x.get_u_y();
00395 }
00396 
00401 void set_derivatives( df1_two_variable& z, const df1_two_variable& x,
00402   const df1_two_variable& y, double u,
00403   double f_u,double f_v)
00404 {
00405     *z.get_u() = u;
00406 
00407     *z.get_u_x() = f_u* *x.get_u_x()
00408                  + f_v* *y.get_u_x();
00409 
00410     *z.get_u_y() = f_u* *x.get_u_y()
00411                  + f_v* *y.get_u_y();
00412 }
00413 
00418   df1_two_variable sqrt(const df1_two_variable& x)
00419   {
00420     df1_two_variable z;
00421     double u=::sqrt(*x.get_u());
00422     *z.get_u()=u;
00423     //double xinv=1.0/(*x.get_u());
00424     double zp=0.5/u;
00425 
00426 
00427     set_derivatives(z,x,u,zp);
00428 
00429     return z;
00430   }
00431 
00436   df1_two_variable atan(const df1_two_variable& x)
00437   {
00438     df1_two_variable z;
00439     double cx=value(x);
00440     double d=1.0/(1+square(cx));
00441     //double d2=square(d);
00442     double u=::atan(cx);
00443     *z.get_u()=u;
00444     double zp=d;
00445 
00446     set_derivatives(z,x,u,zp);
00447     return z;
00448   }
00449 
00454   df1_two_variable square(const df1_two_variable& x)
00455   {
00456     df1_two_variable z;
00457     double u=value(x);
00458     *z.get_u()=u*u;
00459     double zp=2.0*u;
00460 
00461     set_derivatives(z,x,u,zp);
00462     return z;
00463   }
00464 
00469   df1_two_variable tan(const df1_two_variable& x)
00470   {
00471     df1_two_variable z;
00472     double u=::tan(*x.get_u());
00473     *z.get_u()=u;
00474     double v=1.0/::cos(*x.get_u());
00475     //double w=::sin(*x.get_u());
00476     double v2=v*v;
00477     double zp=v2;
00478 
00479     set_derivatives(z,x,u,zp);
00480     return z;
00481   }
00482 
00487   df1_two_variable sin(const df1_two_variable& x)
00488   {
00489     df1_two_variable z;
00490     double u=::sin(*x.get_u());
00491     *z.get_u()=u;
00492     double zp=::cos(*x.get_u());
00493 
00494     set_derivatives(z,x,u,zp);
00495     return z;
00496   }
00497 
00502   df1_two_variable log(const df1_two_variable& x)
00503   {
00504     df1_two_variable z;
00505     double u=::log(*x.get_u());
00506     *z.get_u()=u;
00507     double zp=1/(*x.get_u());
00508 
00509     set_derivatives(z,x,u,zp);
00510     return z;
00511   }
00512 
00517   df1_two_variable exp(const df1_two_variable& x)
00518   {
00519     df1_two_variable z;
00520     double u=::exp(*x.get_u());
00521     *z.get_u()=u;
00522     double zp=u;
00523 
00524     set_derivatives(z,x,u,zp);
00525     return z;
00526   }
00527 
00532   df1_two_variable pow(const df1_two_variable& x,const df1_two_variable& y)
00533   {
00534     return exp(y*log(x));
00535   }
00536 
00541   df1_two_variable pow(double x,const df1_two_variable& y)
00542   {
00543     return exp(y*log(x));
00544   }
00545 
00550   df1_two_variable inv(const df1_two_variable& x)
00551   {
00552     df1_two_variable z;
00553     double xinv=1.0/(*x.get_u());
00554     *z.get_u()=xinv;
00555     double zp=-xinv*xinv;
00556     set_derivatives(z,x,xinv,zp);
00557 
00558     return z;
00559   }
00560 
00565   df1_two_variable& df1_two_variable::operator = (const df1_two_variable& x)
00566   {
00567     *get_u() = *x.get_u();
00568     *get_u_x() = *x.get_u_x();
00569     *get_u_y() = *x.get_u_y();
00570     return *this;
00571   }
00572 
00577   df1_two_variable& df1_two_variable::operator = (double x)
00578   {
00579     *get_u() = x;
00580     *get_u_x() =0.0;
00581     *get_u_y() =0.0;
00582     return *this;
00583   }
00584 
00589   df1_two_variable operator * (const df1_two_variable& x,
00590     const df1_two_variable& y)
00591   {
00592     df1_two_variable z;
00593     double u= *x.get_u() * *y.get_u();
00594     *z.get_u() = u;
00595     double f_u=*y.get_u();
00596     double f_v=*x.get_u();
00597     set_derivatives(z,x,y,u, f_u, f_v);
00598     return z;
00599   }
00600 
00605   df1_two_variable operator * (double x,
00606     const df1_two_variable& y)
00607   {
00608     df1_two_variable z;
00609     *z.get_u() = x *  *y.get_u();
00610     *z.get_u_x() = x * *y.get_u_x();
00611     *z.get_u_y() = x * *y.get_u_y();
00612     return z;
00613   }
00614 
00619   df1_two_variable operator * (const df1_two_variable& y,
00620     double x)
00621   {
00622     df1_two_variable z;
00623     *z.get_u() = x *  *y.get_u();
00624     *z.get_u_x() = x * *y.get_u_x();
00625     *z.get_u_y() = x * *y.get_u_y();
00626     return z;
00627   }
00628 
00633   df1_two_variable operator / (const df1_two_variable& x,
00634     double y)
00635   {
00636     double u=1/y;
00637     return x*u;
00638   }
00639 
00644   df1_two_variable operator / (const df1_two_variable& x,
00645     const df1_two_variable& y)
00646   {
00647     df1_two_variable u=inv(y);
00648     return x*u;
00649   }
00650 
00655   df1_two_variable operator / (const double x,
00656     const df1_two_variable& y)
00657   {
00658     df1_two_variable u=inv(y);
00659     return x*u;
00660   }
00661 
00666   df1_two_variable operator + (const double x,const df1_two_variable& y)
00667   {
00668     df1_two_variable z;
00669     *z.get_u() =  x + *y.get_u();
00670     *z.get_u_x() = *y.get_u_x();
00671     *z.get_u_y() = *y.get_u_y();
00672     return z;
00673   }
00674 
00679   df1_two_variable operator + (const df1_two_variable& x,const double y)
00680   {
00681     df1_two_variable z;
00682     *z.get_u() =  *x.get_u() + y;
00683     *z.get_u_x() = *x.get_u_x();
00684     *z.get_u_y() = *x.get_u_y();
00685     return z;
00686   }
00687 
00692   df1_two_variable operator + (const df1_two_variable& x,
00693     const df1_two_variable& y)
00694   {
00695     df1_two_variable z;
00696     *z.get_u() = *x.get_u() + *y.get_u();
00697     *z.get_u_x() = *x.get_u_x() + *y.get_u_x();
00698     *z.get_u_y() = *x.get_u_y()+*y.get_u_y();
00699     return z;
00700   }
00701 
00706   df1_two_variable operator - (const df1_two_variable& x,
00707     const df1_two_variable& y)
00708   {
00709     df1_two_variable z;
00710     *z.get_u() = *x.get_u() - *y.get_u();
00711     *z.get_u_x() = *x.get_u_x()  - *y.get_u_x();
00712     *z.get_u_y() = *x.get_u_y() - *y.get_u_y();
00713     return z;
00714   }
00715 
00720   df1_two_variable operator - (double x,
00721     const df1_two_variable& y)
00722   {
00723     df1_two_variable z;
00724     *z.get_u() = x - *y.get_u();
00725     *z.get_u_x() = - *y.get_u_x();
00726     *z.get_u_y() = - *y.get_u_y();
00727     return z;
00728   }
00729 
00734   df1_two_variable operator - (const df1_two_variable& x,
00735     double y)
00736   {
00737     df1_two_variable z;
00738     *z.get_u() = *x.get_u()-y;
00739     *z.get_u_x() = *x.get_u_x();
00740     *z.get_u_y() = *x.get_u_y();
00741     return z;
00742   }
00743 
00748   init_df1_two_variable::~init_df1_two_variable()
00749   {
00750     deallocate();
00751   }
00752 
00757   void init_df1_two_variable::deallocate(void)
00758   {
00759     num_ind_var=0;
00760   }
00761 
00766 init_df1_two_variable::init_df1_two_variable(const prevariable& _v)
00767 {
00768   if (num_ind_var > 1)
00769   {
00770     cerr << "can only have 2 independent_variables in df1_two_variable"
00771        " function" << endl;
00772     ad_exit(1);
00773   }
00774   else
00775   {
00776     ADUNCONST(prevariable,v)
00777     ind_var[num_ind_var++]=&v;
00778     *get_u() =  value(v);
00779     switch(num_ind_var)
00780     {
00781     case 1:
00782       *get_u_x() = 1.0;
00783       *get_u_y() = 0.0;
00784       break;
00785     case 2:
00786       *get_u_x() = 0.0;
00787       *get_u_y() = 1.0;
00788       break;
00789     default:
00790       cerr << "illegal num_ind_var value of " << num_ind_var
00791            << " in  df1_two_variable function" << endl;
00792       ad_exit(1);
00793     }
00794   }
00795 }
00796 
00801   init_df1_two_variable::init_df1_two_variable(double v)
00802   {
00803     *get_u() =  v;
00804     *get_u_x() = 0.0;
00805     *get_u_y() = 0.0;
00806   }
00807 
00808 
00813 df1_two_matrix choleski_decomp(const df1_two_matrix& MM)
00814 {
00815   // kludge to deal with constantness
00816   df1_two_matrix & M= (df1_two_matrix &) MM;
00817   int rmin=M.indexmin();
00818   int cmin=M(rmin).indexmin();
00819   int rmax=M.indexmax();
00820   int cmax=M(rmin).indexmax();
00821   if (rmin !=1 || cmin !=1)
00822   {
00823     cerr << "minimum row and column inidices must equal 1 in "
00824       "df1b2matrix choleski_decomp(const df1_two_atrix& MM)"
00825          << endl;
00826     ad_exit(1);
00827   }
00828   if (rmax !=cmax)
00829   {
00830     cerr << "Error in df1b2matrix choleski_decomp(const df1_two_matrix& MM)"
00831       " Matrix not square" << endl;
00832     ad_exit(1);
00833   }
00834 
00835   int n=rmax-rmin+1;
00836   df1_two_matrix L(1,n,1,n);
00837 #ifndef SAFE_INITIALIZE
00838     L.initialize();
00839 #endif
00840 
00841   int i,j,k;
00842   df1_two_variable tmp;
00843 
00844     if (value(M(1,1))<=0)
00845     {
00846       cerr << "Error matrix not positive definite in choleski_decomp"
00847         <<endl;
00848       ad_exit(1);
00849     }
00850 
00851   L(1,1)=sqrt(M(1,1));
00852   for (i=2;i<=n;i++)
00853   {
00854     L(i,1)=M(i,1)/L(1,1);
00855   }
00856 
00857   for (i=2;i<=n;i++)
00858   {
00859     for (j=2;j<=i-1;j++)
00860     {
00861       tmp=M(i,j);
00862       for (k=1;k<=j-1;k++)
00863       {
00864         tmp-=L(i,k)*L(j,k);
00865       }
00866       L(i,j)=tmp/L(j,j);
00867     }
00868     tmp=M(i,i);
00869     for (k=1;k<=i-1;k++)
00870     {
00871       tmp-=L(i,k)*L(i,k);
00872     }
00873 
00874     if (value(tmp)<=0)
00875     {
00876       cerr << "Error matrix not positive definite in choleski_decomp"
00877         <<endl;
00878       ad_exit(1);
00879     }
00880 
00881     L(i,i)=sqrt(tmp);
00882   }
00883 
00884   return L;
00885 }
00886 
00887   df1_two_variable fabs(const df1_two_variable& x)
00888   {
00889     df1_two_variable z;
00890     if (value(x)>=0.0)
00891       z=x;
00892     else
00893       z=-x;
00894     return z;
00895   }