casa  $Rev:20696$
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
AutoDiff.h
Go to the documentation of this file.
00001 //# AutoDiff.h: An automatic differentiating class for functions
00002 //# Copyright (C) 1995,1998,1999,2001,2002
00003 //# Associated Universities, Inc. Washington DC, USA.
00004 //#
00005 //# This library is free software; you can redistribute it and/or modify it
00006 //# under the terms of the GNU Library General Public License as published by
00007 //# the Free Software Foundation; either version 2 of the License, or (at your
00008 //# option) any later version.
00009 //#
00010 //# This library is distributed in the hope that it will be useful, but WITHOUT
00011 //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00012 //# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
00013 //# License for more details.
00014 //#
00015 //# You should have received a copy of the GNU Library General Public License
00016 //# along with this library; if not, write to the Free Software Foundation,
00017 //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
00018 //#
00019 //# Correspondence concerning AIPS++ should be addressed as follows:
00020 //#        Internet email: aips2-request@nrao.edu.
00021 //#        Postal address: AIPS++ Project Office
00022 //#                        National Radio Astronomy Observatory
00023 //#                        520 Edgemont Road
00024 //#                        Charlottesville, VA 22903-2475 USA
00025 //#
00026 //#
00027 //# $Id: AutoDiff.h 21024 2011-03-01 11:46:18Z gervandiepen $
00028 
00029 #ifndef SCIMATH_AUTODIFF_H
00030 #define SCIMATH_AUTODIFF_H
00031 
00032 //# Includes
00033 #include <casa/aips.h>
00034 #include <scimath/Mathematics/AutoDiffRep.h>
00035 #include <casa/Containers/ObjectPool.h>
00036 #include <casa/OS/Mutex.h>
00037 
00038 namespace casa { //# NAMESPACE CASA - BEGIN
00039 
00040 //# Forward declarations
00041 template <class T> class Vector;
00042 template <class T> class AutoDiff;
00043 
00044 // <summary>
00045 // Class that computes partial derivatives by automatic differentiation.
00046 // </summary>
00047 //
00048 // <use visibility=export>
00049 //
00050 // <reviewed reviewer="UNKNOWN" date="before2004/08/25" tests="tAutoDiff.cc" demos="dAutoDiff.cc">
00051 // </reviewed>
00052 //
00053 // <prerequisite>
00054 // <li> 
00055 // </prerequisite>
00056 //
00057 // <etymology>
00058 // Class that computes partial derivatives by automatic differentiation, thus
00059 // AutoDiff.
00060 // </etymology>
00061 //
00062 // <synopsis>
00063 // Class that computes partial derivatives by automatic differentiation.
00064 // It does this by storing the value of a function and the values of its first 
00065 // derivatives with respect to its independent parameters.  When a mathematical
00066 // operation is applied to an AutoDiff object, the derivative values of the 
00067 // resulting new object are computed according to chain rules 
00068 // of differentiation.
00069 //
00070 // Suppose we have a function f(x0,x1,...,xn) and its differential is
00071 // <srcblock> 
00072 // df = (df/dx0)*dx0 + (df/dx1)*dx1 + ... + (df/dxn)*dxn
00073 // </srcblock> 
00074 // We can build a class that has the value of the function, 
00075 // f(x0,x1,...,xn), and the values of the derivatives, (df/dx0), (df/dx1), 
00076 // ..., (df/dxn) at (x0,x1,...,xn), as class members. 
00077 //
00078 // Now if we have another function, g(x0,x1,...,xn) and its differential is
00079 // dg = (dg/dx0)*dx0 + (dg/dx1)*dx1 + ... + (dg/dxn)*dxn,
00080 // since 
00081 // <srcblock> 
00082 // d(f+g) = df + dg, 
00083 // d(f*g) = g*df + f*dg, 
00084 // d(f/g) = df/g - fdg/g^2,
00085 // dsin(f) = cos(f)df, 
00086 // ...,
00087 // </srcblock> 
00088 // we can calculate
00089 // <srcblock>  
00090 // d(f+g), d(f*g), ..., 
00091 // </srcblock>  based on our information on
00092 // <srcblock>  
00093 // df/dx0, df/dx1, ..., dg/dx0, dg/dx1, ..., dg/dxn. 
00094 // </srcblock>  
00095 // All we need to do is to define the operators and derivatives of common 
00096 // mathematical functions.
00097 //
00098 // To be able to use the class as an automatic differentiator of a function
00099 // object, we need a templated function object, i.e. an object with:
00100 // <ul>
00101 // <li> a <src> template <class T> T operator()(const T)</src>
00102 // <li> or multiple variable input like: 
00103 //              <src> template <class T> T operator()(const Vector<T> &)</src>
00104 // <li> all variables and constants used in the calculation of the function
00105 //      value should have been typed with T
00106 // </ul>
00107 // A simple example of such a function object could be:
00108 // <srcblock>
00109 //      template <class T> f {
00110 //      public:
00111 //        T operator()(const T &x, const T &a, const T &b) {
00112 //              return a*b*x; }
00113 //      };
00114 //      // Instantiate the following versions:
00115 //      template class f<Double>;
00116 //      template class f<AutoDiff<Double> >;
00117 // </srcblock>
00118 // A call with values will produce the function value:
00119 // <srcblock>
00120 //      cout << f(7.0, 2.0, 3.0) << endl;
00121 //      // will produce the value at x=7 for a=2; b=3:
00122 //      42
00123 //      // But a call indicating that we want derivatives to a and b:
00124 //      cout << f(AutoDiff<Double>(7.0), AutoDiff<Double>(2.0, 2, 0),
00125 //                AutoDiff<Double>(3.0, 2, 1)) << endl;
00126 //      // will produce the value at x=7 for a=2; b=3:
00127 //      // and the partial derivatives wrt a and b at x=7:
00128 //      (42, [21, 14])
00129 //      // The following will calculate the derivate wrt x:
00130 //      cout << f(AutoDiff<Double>(7.0, 1, 0), AutoDiff<Double>(2.0),
00131 //                AutoDiff<Double>(3.0)) << endl;
00132 //      (42,[6])
00133 // </srcblock>
00134 // In actual practice, there are a few rules to obey for the structure of
00135 // the function object if you want to use the function object and its
00136 // derivatives in least squares fitting procedures in the Fitting
00137 // module. The major one is to view the function object having 'fixed' and
00138 // 'variable' parameters. I.e., rather than viewing the function as
00139 // depending on parameters <em>a, b, x</em> (<src>f(a,b,x)</src>), the 
00140 // function is considered to be <src>f(x; a,b)</src>, where <em>a, b</em> 
00141 // are 'fixed' parameters, and <em>x</em> a variable parameter.
00142 // Fixed parameters should be contained in a 
00143 // <linkto class=FunctionParam>FunctionParam</linkto> container object;
00144 // while the variable parameter(s) are given in the function
00145 // <src>operator()</src>. See <linkto class=Function>Function</linkto> class
00146 // for details.
00147 //
00148 // A Gaussian spectral profile would in general have the center frequency,
00149 // the width and the amplitude as fixed parameters, and the frequency as
00150 // a variable. Given a spectrum, you would solve for the fixed parameters,
00151 // given spectrum values. However, in other cases the role of the
00152 // parameters could be reversed. An example could be a whole stack of
00153 // observed (in the laboratory) spectra at different temperatures at
00154 // one frequency. In that case the width would be the variable parameter,
00155 // and the frequency one of the fixed (and to be solved for)parameters. 
00156 //
00157 // Since the calculation of the derivatives is done with simple overloading, 
00158 // the calculation of second (and higher) derivatives is easy. It should be
00159 // noted that higher deivatives are inefficient in the current incarnation
00160 // (there is no knowledge e.g. about symmetry in the Jacobian). However,
00161 // it is a very good way to get the correct answers of the derivatives. In
00162 // practice actual production code will be better off with specialization
00163 // of the <src>f<AutoDiff<> ></src> implementation.
00164 //
00165 // The <src>AutoDiff</src> class is the class the user communicates with.
00166 // Alias classes (<linkto class=AutoDiffA>AutoDiffA</linkto> and
00167 // <linkto class=AutoDiffA>AutoDiffX</linkto>) exists
00168 // to make it possible to have different incarnations of a templated
00169 // method (e.g. a generic one and a specialized one). See the
00170 // <src>dAutoDiff</src> demo for an example of its use.
00171 //
00172 // All operators and functions are declared in <linkto file=AutoDiffMath.h> 
00173 // AutoDiffMath</linkto>. The output operator in 
00174 // <linkto file=AutoDiffIO.h>AutoDiffIO</linkto>. The actual structure of the
00175 // data block used by <src>AutoDiff</src> is described in 
00176 // <linkto class=AutoDiffRep>AutoDiffRep</linkto>.
00177 // </synopsis>
00178 //
00179 // <example>
00180 // <srcblock>
00181 // // First a simple example.
00182 // // We have a function of the form f(x,y,z); and want to know the
00183 // // value of the function for x=10; y=20; z=30; and for
00184 // // the derivatives at those point.
00185 // // Specify the values; and indicate 3 derivatives:
00186 //      AutoDiff<Double> x(10.0, 3, 0); 
00187 //      AutoDiff<Double> y(20.0, 3, 1); 
00188 //      AutoDiff<Double> z(30.0, 3, 2);
00189 // // The result will be:
00190 //      AutoDiff<Double> result = x*y + sin(z);
00191 //      cout << result.value() << endl;
00192 // // 199.012
00193 //      cout << result.derivatives() << endl;
00194 // // [20, 10, 0.154251]
00195 // // Note: sin(30) = -0.988; cos(30) =  0.154251;
00196 // </srcblock>
00197 //
00198 // See for an extensive example the demo program dAutoDiff. It is
00199 // based on the example given above, and shows also the use of second
00200 // derivatives (which is just using <src>AutoDiff<AutoDiff<Double> ></src>
00201 // as template argument). 
00202 // <srcblock>
00203 //      // The function, with fixed parameters a,b:
00204 //      template <class T> class f {
00205 //      public:
00206 //        T operator()(const T& x) { return a_p*a_p*a_p*b_p*b_p*x; }
00207 //        void set(const T& a, const T& b) { a_p = a; b_p = b; }
00208 //      private:
00209 //        T a_p;
00210 //        T b_p;
00211 //      };
00212 //      // Call it with different template arguments:
00213 //        Double a0(2), b0(3), x0(7);
00214 //        f<Double> f0; f0.set(a0, b0);
00215 //        cout << "Value:     " << f0(x0) << endl;
00216 //      
00217 //        AutoDiff<Double> a1(2,2,0), b1(3,2,1), x1(7);
00218 //        f<AutoDiff<Double> > f1; f1.set(a1, b1);
00219 //        cout << "Diff a,b:   " << f1(x1) << endl;
00220 //      
00221 //        AutoDiff<Double> a2(2), b2(3), x2(7,1,0);
00222 //        f<AutoDiff<Double> > f2; f2.set(a2, b2);
00223 //        cout << "Diff x:     " << f2(x2) << endl;
00224 //      
00225 //        AutoDiff<AutoDiff<Double> > a3(AutoDiff<Double>(2,2,0),2,0),
00226 //          b3(AutoDiff<Double>(3,2,1),2,1), x3(AutoDiff<Double>(7),2);
00227 //        f<AutoDiff<AutoDiff<Double> > > f3; f3.set(a3, b3);
00228 //        cout << "Diff2 a,b:  " << f3(x3) << endl;
00229 //      
00230 //        AutoDiff<AutoDiff<Double> > a4(AutoDiff<Double>(2),1),
00231 //          b4(AutoDiff<Double>(3),1),
00232 //          x4(AutoDiff<Double>(7,1,0),1,0);
00233 //        f<AutoDiff<AutoDiff<Double> > > f4; f4.set(a4, b4);
00234 //        cout << "Diff2 x:    " << f4(x4) << endl;
00235 //
00236 //    // Result will be:
00237 //    // Value:     504
00238 //    // Diff a,b:  (504, [756, 336])
00239 //    // Diff x:    (504, [72])
00240 //    // Diff2 a,b: ((504, [756, 336]), [(756, [756, 504]), (336, [504, 112])])
00241 //    // Diff2 x:   ((504, [72]), [(72, [0])])
00242 //
00243 //    // It needed the template instantiations definitions:
00244 //      template class f<Double>;
00245 //      template class f<AutoDiff<Double> >;
00246 //      template class f<AutoDiff<AutoDiff<Double> > >;
00247 // </srcblock>
00248 // </example>
00249 //
00250 // <motivation>
00251 // The creation of the class was motivated by least-squares non-linear fit where
00252 // partial derivatives of a fitted function are needed. It would be tedious
00253 // to create functionals for all partial derivatives of a function.
00254 // </motivation>
00255 //
00256 // <templating arg=T>
00257 //  <li> any class that has the standard mathematical and comparisons
00258 //      defined
00259 // </templating>
00260 //
00261 // <todo asof="2001/06/07">
00262 // <li> Nothing I know
00263 // </todo>
00264 
00265 template <class T> class AutoDiff {
00266  public:
00267   //# Typedefs
00268   typedef T                     value_type;
00269   typedef value_type&           reference;
00270   typedef const value_type&     const_reference;
00271   typedef value_type*           iterator;
00272   typedef const value_type*     const_iterator;
00273 
00274   //# Constructors
00275   // Construct a constant with a value of zero.  Zero derivatives.
00276   AutoDiff();
00277 
00278   // Construct a constant with a value of v.  Zero derivatives.
00279   AutoDiff(const T &v);
00280 
00281   // A function f(x0,x1,...,xn,...) with a value of v.  The 
00282   // total number of derivatives is ndiffs, the nth derivative is one, and all 
00283   // others are zero. 
00284   AutoDiff(const T &v, const uInt ndiffs, const uInt n); 
00285 
00286   // A function f(x0,x1,...,xn,...) with a value of v.  The 
00287   // total number of derivatives is ndiffs.
00288   // All derivatives are zero. 
00289   AutoDiff(const T &v, const uInt ndiffs); 
00290 
00291   // Construct one from another
00292   AutoDiff(const AutoDiff<T> &other);
00293 
00294   // Construct a function f(x0,x1,...,xn) of a value v and a vector of 
00295   // derivatives derivs(0) = df/dx0, derivs(1) = df/dx1, ...
00296   AutoDiff(const T &v, const Vector<T> &derivs);
00297 
00298   ~AutoDiff();
00299 
00300   // Assignment operator.  Assign a constant to variable.  All derivatives
00301   // are zero.
00302   AutoDiff<T> &operator=(const T &v);
00303 
00304   // Assign one to another.
00305   AutoDiff<T> &operator=(const AutoDiff<T> &other);
00306 
00307   // Assignment operators
00308   // <group>
00309   void operator*=(const AutoDiff<T> &other);
00310   void operator/=(const AutoDiff<T> &other);
00311   void operator+=(const AutoDiff<T> &other);
00312   void operator-=(const AutoDiff<T> &other);
00313   void operator*=(const T other);
00314   void operator/=(const T other);
00315   void operator+=(const T other);
00316   void operator-=(const T other);
00317   // </group>
00318 
00319   // Returns the pointer to the structure of value and derivatives.
00320   // <group>
00321   AutoDiffRep<T> *theRep() { return rep_p; }
00322   const AutoDiffRep<T> *theRep() const { return rep_p; }
00323   // </group>
00324 
00325   // Returns the value of the function
00326   // <group>
00327   T &value() { return rep_p->val_p; }
00328   const T &value() const { return rep_p->val_p; }
00329   // </group>
00330   
00331   // Returns a vector of the derivatives of an AutoDiff
00332   // <group>
00333   Vector<T> derivatives() const;
00334   void derivatives(Vector<T> &res) const;
00335   // </group>
00336   
00337   // Returns a specific derivative. The second set does not check for
00338   // a valid which; the first set does through Vector addressing.
00339   // <group>
00340   T &derivative(uInt which) { return rep_p->grad_p(which); }
00341   const T &derivative(uInt which) const { return rep_p->grad_p(which); }
00342   T &deriv(uInt which) { return rep_p->grad_p[which]; }
00343   const T &deriv(uInt which) const { return rep_p->grad_p[which]; }
00344   // </group>
00345   
00346   // Return total number of derivatives
00347   uInt nDerivatives() const { return rep_p->nd_p; }
00348   
00349   // Is it a constant, i.e., with zero derivatives?
00350   Bool isConstant() const { return rep_p->nd_p == 0; }
00351 
00352   // Indicate that we are going to use a temporary value for the last time
00353   // (e.g. at the of a function returning by value). This way superfluous
00354   // copying can be circumvented.
00355   const AutoDiff<T> &ref() { rep_p->nocopy_p = True; return *this; }
00356   
00357  private:
00358   //# Data
00359   // Pool of data blocks
00360   static ObjectPool<AutoDiffRep<T>, uInt> theirPool;
00361   // Mutex for thread-safe access to theirPool.
00362   static Mutex theirMutex;
00363   // Value representation
00364   AutoDiffRep<T> *rep_p;
00365 
00366   //# Methods
00367   // Release a struct of value and derivative data
00368   void release() {
00369     if (!rep_p->nocopy_p) {
00370       ScopedMutexLock locker(theirMutex);
00371       theirPool.release(rep_p, rep_p->nd_p);
00372     } else {
00373       rep_p->nocopy_p = False;
00374     }
00375   }
00376 
00377 };
00378 
00379 
00380 } //# NAMESPACE CASA - END
00381 
00382 #ifndef CASACORE_NO_AUTO_TEMPLATES
00383 #include <scimath/Mathematics/AutoDiff.tcc>
00384 #endif //# CASACORE_NO_AUTO_TEMPLATES
00385 #endif