// @(#)root/mathmore:$Id$
// Author: L. Moneta Wed Sep  6 09:52:26 2006

/**********************************************************************
 *                                                                    *
 * Copyright (c) 2006  LCG ROOT Math Team, CERN/PH-SFT                *
 *                                                                    *
 *                                                                    *
 **********************************************************************/

// implementation file for class WrappedTF1 and WrappedMultiTF1

#include "Math/WrappedTF1.h"
#include "Math/WrappedMultiTF1.h"
#include "TClass.h"   // needed to copy the TF1 pointer

#include <cmath>


namespace ROOT {

namespace Math {

// static value for epsilon used in derivative calculations
double WrappedTF1::fgEps      = 0.001;
double WrappedMultiTF1::fgEps = 0.001;


WrappedTF1::WrappedTF1 ( TF1 & f  )  :
   fLinear(false),
   fPolynomial(false),
   fFunc(&f),
   fX () 
   //fParams(f.GetParameters(),f.GetParameters()+f.GetNpar())
{
   // constructor from a TF1 function pointer.

   // init the pointers for CINT
   //if (fFunc->GetMethodCall() )  fFunc->InitArgs(fX, &fParams.front() );
   if (fFunc->GetMethodCall() )  fFunc->InitArgs(fX, 0 );
   // distinguish case of polynomial functions and linear functions
   if (fFunc->GetNumber() >= 300 && fFunc->GetNumber() < 310) {
      fLinear = true;
      fPolynomial = true;
   }
   // check that in case function is linear the linear terms are not zero
   if (fFunc->IsLinear() ) {
      int ip = 0; 
      fLinear = true;
      while (fLinear && ip < fFunc->GetNpar() )  { 
         fLinear &= (fFunc->GetLinearPart(ip) != 0) ;
         ip++;
      }
   }
}

WrappedTF1::WrappedTF1(const WrappedTF1 & rhs) :
   BaseFunc(),
   BaseGradFunc(),
   IGrad(),
   fLinear(rhs.fLinear),
   fPolynomial(rhs.fPolynomial),
   fFunc(rhs.fFunc),
   fX()
   //fParams(rhs.fParams)
{
   // copy constructor
   fFunc->InitArgs(fX, 0 );
}

WrappedTF1 & WrappedTF1::operator = (const WrappedTF1 & rhs) {
   // assignment operator
   if (this == &rhs) return *this;  // time saving self-test
   fLinear = rhs.fLinear;
   fPolynomial = rhs.fPolynomial; 
   fFunc = rhs.fFunc; 
   fFunc->InitArgs(fX, 0 );
   //fParams = rhs.fParams;
   return *this;
}

void  WrappedTF1::ParameterGradient(double x, const double * par, double * grad ) const {
   // evaluate the derivative of the function with respect to the parameters
   if (!fLinear) {
      // need to set parameter values
      fFunc->SetParameters( par );
      // no need to call InitArgs (it is called in TF1::GradientPar)
      fFunc->GradientPar(&x,grad,fgEps);
   }
   else {
      unsigned int np = NPar();
      for (unsigned int i = 0; i < np; ++i)
         grad[i] = DoParameterDerivative(x, par, i);
   }
}

double WrappedTF1::DoDerivative( double  x  ) const {
   // return the function derivatives w.r.t. x

   // parameter are passed as non-const in Derivative
   //double * p =  (fParams.size() > 0) ? const_cast<double *>( &fParams.front()) : 0;
   return  fFunc->Derivative(x,(double *) 0,fgEps);
}

double WrappedTF1::DoParameterDerivative(double x, const double * p, unsigned int ipar ) const {
   // evaluate the derivative of the function with respect to the parameters
   //IMPORTANT NOTE: TF1::GradientPar returns 0 for fixed parameters to avoid computing useless derivatives
   //  BUT the TLinearFitter wants to have the derivatives also for fixed parameters.
   //  so in case of fLinear (or fPolynomial) a non-zero value will be returned for fixed parameters

   if (! fLinear ) {
      fFunc->SetParameters( p );
      return fFunc->GradientPar(ipar, &x,fgEps);
   }
   else if (fPolynomial) {
      // case of polynomial function (no parameter dependency)
      return std::pow(x, static_cast<int>(ipar) );
   }
   else {
      // case of general linear function (built in TFormula with ++ )
      const TFormula * df = dynamic_cast<const TFormula*>( fFunc->GetLinearPart(ipar) );
      assert(df != 0);
      fX[0] = x;
      // hack since TFormula::EvalPar is not const
      return (const_cast<TFormula*> ( df) )->Eval( x ) ; // derivatives should not depend on parameters since func is linear 
   }
}

void WrappedTF1::SetDerivPrecision(double eps) { fgEps = eps; }

double WrappedTF1::GetDerivPrecision( ) { return fgEps; }



// impelmentations for WrappedMultiTF1


WrappedMultiTF1::WrappedMultiTF1 (TF1 & f, unsigned int dim  )  : 
   fLinear(false), 
   fPolynomial(false), 
   fOwnFunc(false),
   fFunc(&f),
   fDim(dim)
   //fParams(f.GetParameters(),f.GetParameters()+f.GetNpar())
{ 
   // constructor of WrappedMultiTF1 
   // pass a dimension if dimension specified in TF1 does not correspond to real dimension
   // for example in case of multi-dimensional TF1 objects defined as TF1 (i.e. for functions with dims > 3 )
   if (fDim == 0) fDim = fFunc->GetNdim();

   // check that in case function is linear the linear terms are not zero
   // function is linear when is a TFormula created with "++"
   // hyperplane are not yet existing in TFormula
   if (fFunc->IsLinear() ) {
      int ip = 0; 
      fLinear = true;
      while (fLinear && ip < fFunc->GetNpar() )  { 
         fLinear &= (fFunc->GetLinearPart(ip) != 0) ; 
         ip++;
      }
   }
   // distinguish case of polynomial functions and linear functions
   if (fDim == 1 && fFunc->GetNumber() >= 300 && fFunc->GetNumber() < 310) {
      fLinear = true;
      fPolynomial = true;
   }
}


WrappedMultiTF1::WrappedMultiTF1(const WrappedMultiTF1 & rhs) :
   BaseFunc(),
   BaseParamFunc(),
   fLinear(rhs.fLinear), 
   fPolynomial(rhs.fPolynomial), 
   fOwnFunc(rhs.fOwnFunc),
   fFunc(rhs.fFunc),
   fDim(rhs.fDim)
   //fParams(rhs.fParams) 
{
   // copy constructor
   if (fOwnFunc) SetAndCopyFunction(rhs.fFunc);
}


WrappedMultiTF1 & WrappedMultiTF1::operator= (const WrappedMultiTF1 & rhs) {
   // Assignment operator
   if (this == &rhs) return *this;  // time saving self-test
   fLinear = rhs.fLinear;  
   fPolynomial = rhs.fPolynomial;  
   fOwnFunc = rhs.fOwnFunc;
   fDim = rhs.fDim;
   //fParams = rhs.fParams;
   return *this;
}


void  WrappedMultiTF1::ParameterGradient(const double * x, const double * par, double * grad ) const {
   // evaluate the gradient of the function with respect to the parameters
   //IMPORTANT NOTE: TF1::GradientPar returns 0 for fixed parameters to avoid computing useless derivatives
   //  BUT the TLinearFitter wants to have the derivatives also for fixed parameters.
   //  so in case of fLinear (or fPolynomial) a non-zero value will be returned for fixed parameters

   if (!fLinear) {
      // need to set parameter values
      fFunc->SetParameters( par );
      // no need to call InitArgs (it is called in TF1::GradientPar)
      fFunc->GradientPar(x,grad,fgEps);
   }
   else {  // case of linear functions
      unsigned int np = NPar();
      for (unsigned int i = 0; i < np; ++i)
         grad[i] = DoParameterDerivative(x, par, i);
   }
}

double WrappedMultiTF1::DoParameterDerivative(const double * x, const double * p, unsigned int ipar ) const {
   // evaluate the derivative of the function with respect to parameter ipar
   // see note above concerning the fixed parameters
   if (! fLinear ) {
      fFunc->SetParameters( p );
      return fFunc->GradientPar(ipar, x,fgEps);
   }
   if (fPolynomial) {
      // case of polynomial function (no parameter dependency)  (case for dim = 1)
      assert (fDim == 1);
      if (ipar == 0) return 1.0;
      return std::pow(x[0], static_cast<int>(ipar) );  
   }
   else {
      // case of general linear function (built in TFormula with ++ )
      const TFormula * df = dynamic_cast<const TFormula*>( fFunc->GetLinearPart(ipar) );
      assert(df != 0);
      return (const_cast<TFormula*> ( df) )->EvalPar( x ) ; // derivatives should not depend on parameters since
                                                            // function  is linear
   }
}

void WrappedMultiTF1::SetDerivPrecision(double eps) { fgEps = eps; }

double WrappedMultiTF1::GetDerivPrecision( ) { return fgEps; }

void WrappedMultiTF1::SetAndCopyFunction(const TF1 * f) {
   const TF1 * funcToCopy = (f) ? f : fFunc; 
   TF1 * fnew = (TF1*) funcToCopy->IsA()->New(); 
   funcToCopy->Copy(*fnew); 
   fFunc = fnew; 
   fOwnFunc = true; 
}


} // end namespace Fit

} // end namespace ROOT


 WrappedTF1.cxx:1
 WrappedTF1.cxx:2
 WrappedTF1.cxx:3
 WrappedTF1.cxx:4
 WrappedTF1.cxx:5
 WrappedTF1.cxx:6
 WrappedTF1.cxx:7
 WrappedTF1.cxx:8
 WrappedTF1.cxx:9
 WrappedTF1.cxx:10
 WrappedTF1.cxx:11
 WrappedTF1.cxx:12
 WrappedTF1.cxx:13
 WrappedTF1.cxx:14
 WrappedTF1.cxx:15
 WrappedTF1.cxx:16
 WrappedTF1.cxx:17
 WrappedTF1.cxx:18
 WrappedTF1.cxx:19
 WrappedTF1.cxx:20
 WrappedTF1.cxx:21
 WrappedTF1.cxx:22
 WrappedTF1.cxx:23
 WrappedTF1.cxx:24
 WrappedTF1.cxx:25
 WrappedTF1.cxx:26
 WrappedTF1.cxx:27
 WrappedTF1.cxx:28
 WrappedTF1.cxx:29
 WrappedTF1.cxx:30
 WrappedTF1.cxx:31
 WrappedTF1.cxx:32
 WrappedTF1.cxx:33
 WrappedTF1.cxx:34
 WrappedTF1.cxx:35
 WrappedTF1.cxx:36
 WrappedTF1.cxx:37
 WrappedTF1.cxx:38
 WrappedTF1.cxx:39
 WrappedTF1.cxx:40
 WrappedTF1.cxx:41
 WrappedTF1.cxx:42
 WrappedTF1.cxx:43
 WrappedTF1.cxx:44
 WrappedTF1.cxx:45
 WrappedTF1.cxx:46
 WrappedTF1.cxx:47
 WrappedTF1.cxx:48
 WrappedTF1.cxx:49
 WrappedTF1.cxx:50
 WrappedTF1.cxx:51
 WrappedTF1.cxx:52
 WrappedTF1.cxx:53
 WrappedTF1.cxx:54
 WrappedTF1.cxx:55
 WrappedTF1.cxx:56
 WrappedTF1.cxx:57
 WrappedTF1.cxx:58
 WrappedTF1.cxx:59
 WrappedTF1.cxx:60
 WrappedTF1.cxx:61
 WrappedTF1.cxx:62
 WrappedTF1.cxx:63
 WrappedTF1.cxx:64
 WrappedTF1.cxx:65
 WrappedTF1.cxx:66
 WrappedTF1.cxx:67
 WrappedTF1.cxx:68
 WrappedTF1.cxx:69
 WrappedTF1.cxx:70
 WrappedTF1.cxx:71
 WrappedTF1.cxx:72
 WrappedTF1.cxx:73
 WrappedTF1.cxx:74
 WrappedTF1.cxx:75
 WrappedTF1.cxx:76
 WrappedTF1.cxx:77
 WrappedTF1.cxx:78
 WrappedTF1.cxx:79
 WrappedTF1.cxx:80
 WrappedTF1.cxx:81
 WrappedTF1.cxx:82
 WrappedTF1.cxx:83
 WrappedTF1.cxx:84
 WrappedTF1.cxx:85
 WrappedTF1.cxx:86
 WrappedTF1.cxx:87
 WrappedTF1.cxx:88
 WrappedTF1.cxx:89
 WrappedTF1.cxx:90
 WrappedTF1.cxx:91
 WrappedTF1.cxx:92
 WrappedTF1.cxx:93
 WrappedTF1.cxx:94
 WrappedTF1.cxx:95
 WrappedTF1.cxx:96
 WrappedTF1.cxx:97
 WrappedTF1.cxx:98
 WrappedTF1.cxx:99
 WrappedTF1.cxx:100
 WrappedTF1.cxx:101
 WrappedTF1.cxx:102
 WrappedTF1.cxx:103
 WrappedTF1.cxx:104
 WrappedTF1.cxx:105
 WrappedTF1.cxx:106
 WrappedTF1.cxx:107
 WrappedTF1.cxx:108
 WrappedTF1.cxx:109
 WrappedTF1.cxx:110
 WrappedTF1.cxx:111
 WrappedTF1.cxx:112
 WrappedTF1.cxx:113
 WrappedTF1.cxx:114
 WrappedTF1.cxx:115
 WrappedTF1.cxx:116
 WrappedTF1.cxx:117
 WrappedTF1.cxx:118
 WrappedTF1.cxx:119
 WrappedTF1.cxx:120
 WrappedTF1.cxx:121
 WrappedTF1.cxx:122
 WrappedTF1.cxx:123
 WrappedTF1.cxx:124
 WrappedTF1.cxx:125
 WrappedTF1.cxx:126
 WrappedTF1.cxx:127
 WrappedTF1.cxx:128
 WrappedTF1.cxx:129
 WrappedTF1.cxx:130
 WrappedTF1.cxx:131
 WrappedTF1.cxx:132
 WrappedTF1.cxx:133
 WrappedTF1.cxx:134
 WrappedTF1.cxx:135
 WrappedTF1.cxx:136
 WrappedTF1.cxx:137
 WrappedTF1.cxx:138
 WrappedTF1.cxx:139
 WrappedTF1.cxx:140
 WrappedTF1.cxx:141
 WrappedTF1.cxx:142
 WrappedTF1.cxx:143
 WrappedTF1.cxx:144
 WrappedTF1.cxx:145
 WrappedTF1.cxx:146
 WrappedTF1.cxx:147
 WrappedTF1.cxx:148
 WrappedTF1.cxx:149
 WrappedTF1.cxx:150
 WrappedTF1.cxx:151
 WrappedTF1.cxx:152
 WrappedTF1.cxx:153
 WrappedTF1.cxx:154
 WrappedTF1.cxx:155
 WrappedTF1.cxx:156
 WrappedTF1.cxx:157
 WrappedTF1.cxx:158
 WrappedTF1.cxx:159
 WrappedTF1.cxx:160
 WrappedTF1.cxx:161
 WrappedTF1.cxx:162
 WrappedTF1.cxx:163
 WrappedTF1.cxx:164
 WrappedTF1.cxx:165
 WrappedTF1.cxx:166
 WrappedTF1.cxx:167
 WrappedTF1.cxx:168
 WrappedTF1.cxx:169
 WrappedTF1.cxx:170
 WrappedTF1.cxx:171
 WrappedTF1.cxx:172
 WrappedTF1.cxx:173
 WrappedTF1.cxx:174
 WrappedTF1.cxx:175
 WrappedTF1.cxx:176
 WrappedTF1.cxx:177
 WrappedTF1.cxx:178
 WrappedTF1.cxx:179
 WrappedTF1.cxx:180
 WrappedTF1.cxx:181
 WrappedTF1.cxx:182
 WrappedTF1.cxx:183
 WrappedTF1.cxx:184
 WrappedTF1.cxx:185
 WrappedTF1.cxx:186
 WrappedTF1.cxx:187
 WrappedTF1.cxx:188
 WrappedTF1.cxx:189
 WrappedTF1.cxx:190
 WrappedTF1.cxx:191
 WrappedTF1.cxx:192
 WrappedTF1.cxx:193
 WrappedTF1.cxx:194
 WrappedTF1.cxx:195
 WrappedTF1.cxx:196
 WrappedTF1.cxx:197
 WrappedTF1.cxx:198
 WrappedTF1.cxx:199
 WrappedTF1.cxx:200
 WrappedTF1.cxx:201
 WrappedTF1.cxx:202
 WrappedTF1.cxx:203
 WrappedTF1.cxx:204
 WrappedTF1.cxx:205
 WrappedTF1.cxx:206
 WrappedTF1.cxx:207
 WrappedTF1.cxx:208
 WrappedTF1.cxx:209
 WrappedTF1.cxx:210
 WrappedTF1.cxx:211
 WrappedTF1.cxx:212
 WrappedTF1.cxx:213
 WrappedTF1.cxx:214
 WrappedTF1.cxx:215
 WrappedTF1.cxx:216
 WrappedTF1.cxx:217
 WrappedTF1.cxx:218
 WrappedTF1.cxx:219
 WrappedTF1.cxx:220
 WrappedTF1.cxx:221
 WrappedTF1.cxx:222
 WrappedTF1.cxx:223
 WrappedTF1.cxx:224
 WrappedTF1.cxx:225
 WrappedTF1.cxx:226
 WrappedTF1.cxx:227
 WrappedTF1.cxx:228
 WrappedTF1.cxx:229
 WrappedTF1.cxx:230
 WrappedTF1.cxx:231
 WrappedTF1.cxx:232
 WrappedTF1.cxx:233
 WrappedTF1.cxx:234
 WrappedTF1.cxx:235
 WrappedTF1.cxx:236
 WrappedTF1.cxx:237
 WrappedTF1.cxx:238
 WrappedTF1.cxx:239
 WrappedTF1.cxx:240
 WrappedTF1.cxx:241
 WrappedTF1.cxx:242
 WrappedTF1.cxx:243
 WrappedTF1.cxx:244
 WrappedTF1.cxx:245
 WrappedTF1.cxx:246
 WrappedTF1.cxx:247
 WrappedTF1.cxx:248
 WrappedTF1.cxx:249
 WrappedTF1.cxx:250
 WrappedTF1.cxx:251
 WrappedTF1.cxx:252
 WrappedTF1.cxx:253
 WrappedTF1.cxx:254
 WrappedTF1.cxx:255