00001 // @(#)root/mathmore:$Id: GSLMultiFitFunctionAdapter.h 33180 2010-04-25 10:14:07Z moneta $ 00002 // Authors: L. Moneta, Dec 2006 00003 00004 /********************************************************************** 00005 * * 00006 * Copyright (c) 2004 ROOT Foundation, CERN/PH-SFT * 00007 * * 00008 * This library is free software; you can redistribute it and/or * 00009 * modify it under the terms of the GNU General Public License * 00010 * as published by the Free Software Foundation; either version 2 * 00011 * of the License, or (at your option) any later version. * 00012 * * 00013 * This library is distributed in the hope that it will be useful, * 00014 * but WITHOUT ANY WARRANTY; without even the implied warranty of * 00015 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * 00016 * General Public License for more details. * 00017 * * 00018 * You should have received a copy of the GNU General Public License * 00019 * along with this library (see file COPYING); if not, write * 00020 * to the Free Software Foundation, Inc., 59 Temple Place, Suite * 00021 * 330, Boston, MA 02111-1307 USA, or contact the author. * 00022 * * 00023 **********************************************************************/ 00024 00025 // Header file for class GSLMultiMinFunctionAdapter 00026 // 00027 // Generic adapter for gsl_multimin_function signature 00028 // usable for any c++ class which defines operator( ) 00029 // 00030 // Created by: Lorenzo Moneta at Fri Nov 12 16:58:51 2004 00031 // 00032 // Last update: Fri Nov 12 16:58:51 2004 00033 // 00034 #ifndef ROOT_Math_GSLMultiFitFunctionAdapter 00035 #define ROOT_Math_GSLMultiFitFunctionAdapter 00036 00037 #include "gsl/gsl_vector.h" 00038 #include "gsl/gsl_matrix.h" 00039 00040 #include <cassert> 00041 00042 #include <iostream> 00043 00044 namespace ROOT { 00045 namespace Math { 00046 00047 00048 00049 /** 00050 Class for adapting a C++ functor class to C function pointers used by GSL MultiFit 00051 Algorithm 00052 The templated C++ function class must implement: 00053 00054 <em> double operator( const double * x)</em> 00055 and if the derivatives are required: 00056 <em> void Gradient( const double * x, double * g)</em> 00057 and 00058 <em> void FdF( const double * x, double &f, double * g)</em> 00059 00060 This class defines static methods with will be used to fill the 00061 \a gsl_multimin_function and 00062 \a gsl_multimin_function_fdf structs used by GSL. 00063 See for examples the 00064 <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Providing-a-function-to-minimize.html#Providing-a-function-to-minimize">GSL online manual</A> 00065 00066 @ingroup MultiMin 00067 */ 00068 00069 00070 template<class FuncVector> 00071 class GSLMultiFitFunctionAdapter { 00072 00073 public: 00074 00075 static int F( const gsl_vector * x, void * p, gsl_vector * f ) { 00076 // p is a pointer to an iterator of functions 00077 unsigned int n = f->size; 00078 // need to copy iterator otherwise next time the function is called it wont work 00079 FuncVector & funcVec = *( reinterpret_cast< FuncVector *> (p) ); 00080 if (n == 0) return -1; 00081 for (unsigned int i = 0; i < n ; ++i) { 00082 gsl_vector_set(f, i, (funcVec[i])(x->data) ); 00083 } 00084 return 0; 00085 } 00086 00087 00088 static int Df( const gsl_vector * x, void * p, gsl_matrix * h) { 00089 00090 // p is a pointer to an iterator of functions 00091 unsigned int n = h->size1; 00092 unsigned int npar = h->size2; 00093 if (n == 0) return -1; 00094 if (npar == 0) return -2; 00095 FuncVector & funcVec = *( reinterpret_cast< FuncVector *> (p) ); 00096 for (unsigned int i = 0; i < n ; ++i) { 00097 double * g = (h->data)+i*npar; //pointer to start of i-th row 00098 assert ( npar == (funcVec[i]).NDim() ); 00099 (funcVec[i]).Gradient(x->data, g); 00100 } 00101 return 0; 00102 } 00103 00104 /// evaluate derivative and function at the same time 00105 static int FDf( const gsl_vector * x, void * p, gsl_vector * f, gsl_matrix * h) { 00106 // should be implemented in the function 00107 // p is a pointer to an iterator of functions 00108 unsigned int n = h->size1; 00109 unsigned int npar = h->size2; 00110 if (n == 0) return -1; 00111 if (npar == 0) return -2; 00112 FuncVector & funcVec = *( reinterpret_cast< FuncVector *> (p) ); 00113 assert ( f->size == n); 00114 for (unsigned int i = 0; i < n ; ++i) { 00115 assert ( npar == (funcVec[i]).NDim() ); 00116 double fval = 0; 00117 double * g = (h->data)+i*npar; //pointer to start of i-th row 00118 (funcVec[i]).FdF(x->data, fval, g); 00119 gsl_vector_set(f, i, fval ); 00120 } 00121 return 0; 00122 } 00123 00124 }; 00125 00126 00127 } // namespace Math 00128 } // namespace ROOT 00129 00130 00131 #endif /* ROOT_Math_GSLMultiMinFunctionAdapter */