library: libFumili #include "TFumili.h"

# TFumili

class description - source file - inheritance tree (.pdf)

## class TFumili : public TVirtualFitter

| Inheritance Chart:
`TObject`
| <- | `TNamed`
| <- | `TVirtualFitter`
| <- | `TFumili` |

        public:

                            TFumili(Int_t maxpar = 25)
                            TFumili(const TFumili&)
                    virtual ~TFumili()
                       void BuildArrays()
           virtual Double_t Chisquare(Int_t npar, Double_t* params) const
             static TClass* Class()
               virtual void Clear(Option_t* opt = "")
                       void DeleteArrays()
                       void Derivatives(Double_t*, Double_t*)
                      Int_t Eval(Int_t& npar, Double_t* grad, Double_t& fval, Double_t* par, Int_t flag)
                   Double_t EvalTFN(Double_t*, Double_t*)
              virtual Int_t ExecuteCommand(const char* command, Double_t* args, Int_t nargs)
                      Int_t ExecuteSetCommand(Int_t)
               virtual void FixParameter(Int_t ipar)
          virtual Double_t* GetCovarianceMatrix() const
           virtual Double_t GetCovarianceMatrixElement(Int_t i, Int_t j) const
              virtual Int_t GetErrors(Int_t ipar, Double_t& eplus, Double_t& eminus, Double_t& eparab, Double_t& globcc) const
              virtual Int_t GetNumberFreeParameters() const
              virtual Int_t GetNumberTotalParameters() const
           virtual Double_t GetParameter(Int_t ipar) const
              virtual Int_t GetParameter(Int_t ipar, char* name, Double_t& value, Double_t& verr, Double_t& vlow, Double_t& vhigh) const
           virtual Double_t GetParError(Int_t ipar) const
                  Double_t* GetPL0() const
              virtual Int_t GetStats(Double_t& amin, Double_t& edm, Double_t& errdef, Int_t& nvpar, Int_t& nparx) const
           virtual Double_t GetSumLog(Int_t)
                  Double_t* GetZ() const
                       void InvertZ(Int_t)
            virtual TClass* IsA() const
             virtual Bool_t IsFixed(Int_t ipar) const
                      Int_t Minimize()
                   TFumili& operator=(const TFumili&)
               virtual void PrintResults(Int_t k, Double_t p) const
               virtual void ReleaseParameter(Int_t ipar)
                       void SetData(Double_t*, Int_t, Int_t)
               virtual void SetFitMethod(const char* name)
              virtual Int_t SetParameter(Int_t ipar, const char* parname, Double_t value, Double_t verr, Double_t vlow, Double_t vhigh)
                       void SetParNumber(Int_t ParNum)
                      Int_t SGZ()
               virtual void ShowMembers(TMemberInspector& insp, char* parent)
               virtual void Streamer(TBuffer& b)
                       void StreamerNVirtual(TBuffer& b)

### Data Members

        private:

              Int_t fMaxParam
              Int_t fMaxParam2           fMaxParam*fMaxParam
              Int_t fNlog
              Int_t fNfcn                Number of FCN calls;
              Int_t fNED1                Number of experimental vectors X=(x1,x2,...xK)
              Int_t fNED2                K - Length of vector X plus 2 (for chi2)
              Int_t fNED12               fNED1+fNED2
              Int_t fNpar                fNpar - number of parameters
              Int_t fNstepDec            fNstepDec - maximum number of step decreasing counter
              Int_t fNlimMul             fNlimMul - after fNlimMul successful iterations permits four-fold increasing of fPL
              Int_t fNmaxIter            fNmaxIter - maximum number of iterations
              Int_t fLastFixed           Last fixed parameter number
              Int_t fENDFLG              End flag of fit
              Int_t fINDFLG[5]           internal flags;
             Bool_t fGRAD                user calculated gradients
             Bool_t fWARN                warnings
             Bool_t fDEBUG               debug info
             Bool_t fLogLike             LogLikelihood flag
             Bool_t fNumericDerivatives
          Double_t* fZ0                  [fMaxParam2] Matrix of approximate second derivatives of objective function
          Double_t* fZ                   [fMaxParam2] Invers fZ0 matrix - covariance matrix
          Double_t* fGr                  [fMaxParam] Gradients of objective function
          Double_t* fParamError          [fMaxParam] Parameter errors
          Double_t* fSumLog              [fNlog]
          Double_t* fEXDA                [fNED12] experimental data poInt_ter
          Double_t* fA                   [fMaxParam] Fit parameter array
          Double_t* fPL0                 [fMaxParam] Step initial bounds
          Double_t* fPL                  [fMaxParam] Limits for parameters step. If <0, then parameter is fixed
          Double_t* fDA                  [fMaxParam] Parameter step
          Double_t* fAMX                 [fMaxParam] Maximum param value
          Double_t* fAMN                 [fMaxParam] Minimum param value
          Double_t* fR                   [fMaxParam] Correlation factors
          Double_t* fDF                  [fMaxParam] // First derivatives of theoretical function
          Double_t* fCmPar               [fMaxParam] parameters of commands
           Double_t fS                   fS - objective function value (return)
           Double_t fEPS                 fEPS - required precision of parameters. If fEPS<0 then
           Double_t fRP                  Precision of fit ( machine zero on CDC 6000) quite old yeh?
           Double_t fAKAPPA
           Double_t fGT                  Expected function change in next iteration
           TString* fANames              [fMaxParam] Parameter names
            TString fCword               Command string

## Class Description

             FUMILI
      Based on ideas, proposed by I.N. Silin
        [See NIM A440, 2000 (p431)]
     converted from FORTRAN to C  by
         Sergey Yaschenko <s.yaschenko@fz-juelich.de>

    ______________________________________________________________________________

## FUMILI minimization package

FUMILI is used to minimize Chi-square function or to search maximum of likelihood function.

Experimentally measured values $F_i$ are fitted with theoretical functions $f_i({\vec x}_i,\vec\theta\,\,)$, where ${\vec x}_i$ are coordinates, and $\vec\theta$ -- vector of parameters.

For better convergence Chi-square function has to be the following form

$$ {\chi^2\over2}={1\over2}\sum^n_{i=1}\left(f_i(\vec x_i,\vec\theta\,\,)-F_i\over\sigma_i\right)^2 \eqno(1) $$

where $\sigma_i$ are errors of measured function.

The minimum condition is

$$ {\partial\chi^2\over\partial\theta_i}=\sum^n_{j=1}{1\over\sigma^2_j}\cdot {\partial f_j\over\partial\theta_i}\left[f_j(\vec x_j,\vec\theta\,\,)-F_j\right]=0,\qquad i=1\ldots m\eqno(2) $$

where m is the quantity of parameters.

Expanding left part of (2) over parameter increments and retaining only linear terms one gets

$$ \left(\partial\chi^2\over\theta_i\right)_{\vec\theta={\vec\theta}^0} +\sum_k\left(\partial^2\chi^2\over\partial\theta_i\partial\theta_k\right)_{ \vec\theta={\vec\theta}^0}\cdot(\theta_k-\theta_k^0) = 0\eqno(3) $$

Here ${\vec\theta}_0$ is some initial value of parameters. In general case:

$$ {\partial^2\chi^2\over\partial\theta_i\partial\theta_k}= \sum^n_{j=1}{1\over\sigma^2_j}{\partial f_j\over\theta_i} {\partial f_j\over\theta_k} + \sum^n_{j=1}{(f_j - F_j)\over\sigma^2_j}\cdot {\partial^2f_j\over\partial\theta_i\partial\theta_k}\eqno(4) $$

In FUMILI algorithm for second derivatives of Chi-square approximate expression is used when last term in (4) is discarded. It is often done, not always wittingly, and sometimes causes troubles, for example, if user wants to limit parameters with positive values by writing down $\theta_i^2$ instead of $\theta_i$. FUMILI will fail if one tries minimize $\chi^2 = g^2(\vec\theta)$ where g is arbitrary function.

Approximate value is:

$${\partial^2\chi^2\over\partial\theta_i\partial\theta_k}\approx Z_{ik}= \sum^n_{j=1}{1\over\sigma^2_j}{\partial f_j\over\theta_i} {\partial f_j\over\theta_k}\eqno(5) $$

Then the equations for parameter increments are

$$\left(\partial\chi^2\over\partial\theta_i\right)_{\vec\theta={\vec\theta}^0} +\sum_k Z_{ik}\cdot(\theta_k-\theta^0_k) = 0, \qquad i=1\ldots m\eqno(6) $$

Remarkable feature of algorithm is the technique for step restriction. For an initial value of parameter ${\vec\theta}^0$ a parallelepiped $P_0$ is built with the center at ${\vec\theta}^0$ and axes parallel to coordinate axes $\theta_i$. The lengths of parallelepiped sides along i-th axis is $2b_i$, where $b_i$ is such a value that the functions $f_j(\vec\theta)$ are quasi-linear all over the parallelepiped.

FUMILI takes into account simple linear inequalities in the form: $$ \theta_i^{\rm min}\le\theta_i\le\theta^{\rm max}_i\eqno(7) $$

They form parallelepiped $P$ ($P_0$ may be deformed by $P$). Very similar step formulae are used in FUMILI for negative logarithm of the likelihood function with the same idea - linearization of functional argument.

    ______________________________________________________________________________

TFumili(Int_t maxpar)

     maxpar is the maximum number of parameters used with TFumili object

void BuildArrays()

       Allocates memory for internal arrays. Called by TFumili::TFumili

~TFumili()

     TFumili destructor

Double_t Chisquare(Int_t npar, Double_t *params) const

     return a chisquare equivalent

void Clear(Option_t *)

     Resets all parameter names, values and errors to zero

     Argument opt is ignored

     NB: this procedure doesn't reset parameter limits

void DeleteArrays()

     Deallocates memory. Called from destructor TFumili::~TFumili

void Derivatives(Double_t *DF,Double_t *fX)

     Calculates partial derivatives of theoretical function

     Input:
        fX  - vector of data point
     Output:
        DF - array of derivatives

     ARITHM.F
     Converted from CERNLIB

Int_t Eval(Int_t& npar, Double_t *grad, Double_t &fval, Double_t *par, Int_t flag)

     Evaluate the minimisation function
      Input parameters:
        npar:    number of currently variable parameters
        par:     array of (constant and variable) parameters
        flag:    Indicates what is to be calculated
        grad:    array of gradients
      Output parameters:
        fval:    The calculated function value.
        grad:    The vector of first derivatives.

     The meaning of the parameters par is of course defined by the user,
     who uses the values of those parameters to calculate his function value.
     The starting values must be specified by the user.
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     Inside FCN user has to define Z-matrix by means TFumili::GetZ
      and TFumili::Derivatives,
     set theoretical function by means of TFumili::SetUserFunc,
     but first - pass number of parameters by TFumili::SetParNumber
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     Later values are determined by Fumili as it searches for the minimum
     or performs whatever analysis is requested by the user.

     The default function calls the function specified in SetFCN

Double_t EvalTFN(Double_t * /*df*/, Double_t *X)

     Evaluate theoretical function
     df: array of partial derivatives
     X:  vector of theoretical function argument

Int_t ExecuteCommand(const char *command, Double_t *args, Int_t nargs)

      Execute MINUIT commands. MINImize, SIMplex, MIGrad and FUMili all
      will call TFumili::Minimize method.

      For full command list see
      MINUIT. Reference Manual. CERN Program Library Long Writeup D506.

      Improvement and errors calculation are not yet implemented as well
      as Monte-Carlo seeking and minimization.
      Contour commands are also unsupported.

      command   : command string
      args      : array of arguments
      nargs     : number of arguments

Int_t ExecuteSetCommand(Int_t nargs)

     Called from TFumili::ExecuteCommand in case
     of "SET xxx" and "SHOW xxx".

void FixParameter(Int_t ipar)

     Fixes parameter number ipar

Double_t* GetCovarianceMatrix() const

     return a pointer to the covariance matrix

Double_t GetCovarianceMatrixElement(Int_t i, Int_t j) const

     return element i,j from the covariance matrix

Int_t GetNumberTotalParameters() const

     return the total number of parameters (free + fixed)

Int_t GetNumberFreeParameters() const

     return the number of free parameters

Double_t GetParError(Int_t ipar) const

     return error of parameter ipar

Double_t GetParameter(Int_t ipar) const

     return current value of parameter ipar

Int_t GetParameter(Int_t ipar,char *cname,Double_t &value,Double_t &verr,Double_t &vlow, Double_t &vhigh) const

     Get various ipar parameter attributs:

     cname:    parameter name
     value:    parameter value
     verr:     parameter error
     vlow:     lower limit
     vhigh:    upper limit
     WARNING! parname must be suitably dimensionned in the calling function.

Int_t GetErrors(Int_t ipar,Double_t &eplus, Double_t &eminus, Double_t &eparab, Double_t &globcc) const

     Return errors after MINOs
     not implemented

Int_t GetStats(Double_t &amin, Double_t &edm, Double_t &errdef, Int_t &nvpar, Int_t &nparx) const

     return global fit parameters
       amin     : chisquare
       edm      : estimated distance to minimum
       errdef
       nvpar    : number of variable parameters
       nparx    : total number of parameters

Double_t GetSumLog(Int_t n)

     return Sum(log(i) i=0,n
     used by log likelihood fits

void InvertZ(Int_t n)

     Inverts packed diagonal matrix Z by square-root method.
      Matrix elements corresponding to
     fix parameters are removed.

     n: number of variable parameters

Bool_t IsFixed(Int_t ipar) const

    return kTRUE if parameter ipar is fixed, kFALSE othersise)

Int_t Minimize()

    *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
             FUMILI
      Based on ideas, proposed by I.N. Silin
        [See NIM A440, 2000 (p431)]
     converted from FORTRAN to C  by
         Sergey Yaschenko <s.yaschenko@fz-juelich.de>

    *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

     This function is called after setting theoretical function
     by means of TFumili::SetUserFunc and initializing parameters.
     Optionally one can set FCN function (see TFumili::SetFCN and TFumili::Eval)
     If FCN is undefined then user has to provide data arrays by calling
      TFumili::SetData procedure.

     TFumili::Minimize return following values:
        0  - fit is converged
       -2  - functional is not decreasing (or bad derivatives)
       -3  - error estimations are infinite
       -4  - maximum number of iterations is exceeded

void PrintResults(Int_t ikode,Double_t p) const

     Prints fit results.

     ikode is the type of printing parameters
     p is functional value

      ikode = 1   - print values, errors and limits
      ikode = 2   - print values, errors and steps
      ikode = 3   - print values, errors, steps and derivatives
      ikode = 4   - print only values and errors

void ReleaseParameter(Int_t ipar)

     Releases parameter number ipar

void SetData(Double_t *exdata,Int_t numpoints,Int_t vecsize)

     Sets pointer to data array provided by user.
     Necessary if SetFCN is not called.

     numpoints:    number of experimental points
     vecsize:      size of data point vector + 2
                   (for N-dimensional fit vecsize=N+2)
     exdata:       data array with following format

       exdata[0] = ExpValue_0     - experimental data value number 0
       exdata[1] = ExpSigma_0     - error of value number 0
       exdata[2] = X_0[0]
       exdata[3] = X_0[1]
           .........
       exdata[vecsize-1] = X_0[vecsize-3]
       exdata[vecsize]   = ExpValue_1
       exdata[vecsize+1] = ExpSigma_1
       exdata[vecsize+2] = X_1[0]
           .........
       exdata[vecsize*(numpoints-1)] = ExpValue_(numpoints-1)
           .........
       exdata[vecsize*numpoints-1] = X_(numpoints-1)[vecsize-3]

void SetFitMethod(const char *name)

     ret fit method (chisquare or loglikelihood)

Int_t SetParameter(Int_t ipar,const char *parname,Double_t value,Double_t verr,Double_t vlow, Double_t vhigh)

     Sets for prameter number ipar initial parameter value,
     name parname, initial error verr and limits vlow and vhigh
     If vlow = vhigh but not equil to zero, parameter will be fixed.
     If vlow = vhigh = 0, parameter is released and its limits are discarded

Int_t SGZ()

      Evaluates objective function ( chi-square ), gradients and
      Z-matrix using data provided by user via TFumili::SetData

### Inline Functions

              Double_t* GetPL0() const
              Double_t* GetZ() const
                   void SetParNumber(Int_t ParNum)
                TClass* Class()
                TClass* IsA() const
                   void ShowMembers(TMemberInspector& insp, char* parent)
                   void Streamer(TBuffer& b)
                   void StreamerNVirtual(TBuffer& b)
                TFumili TFumili(const TFumili&)
               TFumili& operator=(const TFumili&)

Author: Stanislav Nesterov 07/05/2003
Last update: root/fumili:$Name: $:$Id: TFumili.cxx,v 1.23 2005/06/01 07:43:19 brun Exp $