[math][minuit2] Add new version of Fumili algorithm using trust region method

cmake/modules/SearchInstalledSoftware.cmake

@@ -564,7 +564,7 @@ if(mathmore OR builtin_gsl)
       endif()
     endif()
   else()
-    set(gsl_version 2.5)
+    set(gsl_version 2.8)
     message(STATUS "Downloading and building GSL version ${gsl_version}")
     foreach(l gsl gslcblas)
       list(APPEND GSL_LIBRARIES ${CMAKE_BINARY_DIR}/lib/${CMAKE_STATIC_LIBRARY_PREFIX}${l}${CMAKE_STATIC_LIBRARY_SUFFIX})

etc/plugins/ROOT@@Math@@Minimizer/P040_GSLNLSMinimizer.C

@@ -1,5 +1,5 @@
 void P040_GSLNLSMinimizer()
 {
    gPluginMgr->AddHandler("ROOT::Math::Minimizer", "GSLMultiFit", "ROOT::Math::GSLNLSMinimizer",
-      "MathMore", "GSLNLSMinimizer(int)");
+      "MathMore", "GSLNLSMinimizer(const char * )");
 }

hist/hist/src/HFitImpl.cxx

@@ -310,7 +310,7 @@ TFitResultPtr HFit::Fit(FitObject * h1, TF1 *f1 , Foption_t & fitOption , const
    fitConfig.SetMinimizerOptions(minOption);
 
    // specific  print level options
-   if (fitOption.Verbose) fitConfig.MinimizerOptions().SetPrintLevel(3);
+   if (fitOption.Verbose)  fitConfig.MinimizerOptions().SetPrintLevel(fitOption.Verbose + 1);
    if (fitOption.Quiet)    fitConfig.MinimizerOptions().SetPrintLevel(0);
 
    // specific minimizer options depending on minimizer
@@ -787,8 +787,12 @@ void ROOT::Fit::FitOptionsMake(EFitObjectType type, const char *option, Foption_
          fitOption.User = 0;
       }
    }
-   if (opt.Contains("Q")) fitOption.Quiet   = 1;
-   if (opt.Contains("V")) {fitOption.Verbose = 1; fitOption.Quiet   = 0;}
+
+   // in case of Q and V options V has precedence
+   if (opt.Contains("VVV") || opt.Contains("DEBUG")) { fitOption.Verbose = 3; }
+   else if (opt.Contains("VV")) {fitOption.Verbose = 2; }
+   else if (opt.Contains("V")) {fitOption.Verbose = 1; }
+   else if (opt.Contains("Q")) {fitOption.Quiet   = 1; }
 
 
    if (opt.Contains("E")) fitOption.Errors  = 1;
@@ -895,7 +899,7 @@ TFitResultPtr ROOT::Fit::UnBinFit(ROOT::Fit::UnBinData * data, TF1 * fitfunc, Fo
 
    fitConfig.SetMinimizerOptions(minOption);
 
-   if (fitOption.Verbose)   fitConfig.MinimizerOptions().SetPrintLevel(3);
+   if (fitOption.Verbose)   fitConfig.MinimizerOptions().SetPrintLevel(fitOption.Verbose+1);
    if (fitOption.Quiet)     fitConfig.MinimizerOptions().SetPrintLevel(0);
 
    // more

hist/hist/src/TF1Convolution.cxx

@@ -329,7 +329,9 @@ void TF1Convolution::MakeFFTConv()
 
 Double_t TF1Convolution::EvalFFTConv(Double_t t)
 {
-   if (!fFlagGraph)  MakeFFTConv();
+   if (!fFlagGraph) {
+      MakeFFTConv();
+   }
    // if cannot make FFT use numconv
    if (fGraphConv)
       return  fGraphConv -> Eval(t);
@@ -430,10 +432,10 @@ void TF1Convolution::SetParameters(Double_t p0, Double_t p1, Double_t p2, Double
 }
 
 ////////////////////////////////////////////////////////////////////////////////
-/// Set the fraction of extra range used when doing an FFT convolution.
+/// Set the fraction of extra range used when doing an convolution.
 /// The extra range is often needed to avoid mirroring effect of the resulting convolution
 /// function at the borders.
-/// By default an extra range of 0.1 is used.
+/// By default an extra range of 0.1 is used when doing FFT and it is not use for numerical convolution
 
 void TF1Convolution::SetExtraRange(Double_t percentage)
 {

math/fumili/src/TFumili.cxx

@@ -110,6 +110,8 @@ function argument.
 #include "TList.h"
 #include "TVirtualFitter.h"
 
+#include <iostream>
+
 
 extern void H1FitChisquareFumili(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag);
 extern void H1FitLikelihoodFumili(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag);
@@ -1067,6 +1069,27 @@ Bool_t TFumili::IsFixed(Int_t ipar) const
    else                return kFALSE;
 }
 
+void PrintVector(const char * name, int n, double * x) {
+   // print vector data
+   std::cout << name << " : {";
+   for (int i  = 0; i < n; i++)
+      std::cout << "  " << x[i];
+   std::cout << " }\n";
+}
+void PrintMatrix(const char * name, int n, double * x) {
+   // print matrix data
+   std::cout << name << " : \n";
+   int index = 0;
+   for (int i  = 0; i < n; i++) {
+      for (int j  = 0; j <= i; j++) {
+         std::cout << "  " << x[index];
+         index++;
+      }
+      std::cout << std::endl;
+   }
+   std::cout << "\n";
+}
+
 ////////////////////////////////////////////////////////////////////////////////
 /// Main minimization procedure
 ///
@@ -1092,6 +1115,8 @@ Int_t TFumili::Minimize()
    //
    Int_t parn;
 
+   // I think this funciton call is needed for getting fA and fNpar
+   // and fAMX fAMN ?
    if(fFCN) {
       Eval(parn,fGr,fS,fA,9); fNfcn++;
    }
@@ -1102,7 +1127,7 @@ Int_t TFumili::Minimize()
 
    Int_t nn2, n, fixFLG,  ifix1, fi, nn3, nn1, n0;
    Double_t t1;
-   Double_t sp, t, olds=0;
+   Double_t sp, t = 0, olds=0;
    Double_t bi, aiMAX=0, amb;
    Double_t afix, sigi, akap;
    Double_t alambd, al, bm, abi, abm;
@@ -1135,15 +1160,25 @@ Int_t TFumili::Minimize()
    nn1 = 1;
    t1 = 1.;
 
+
+
+//perform iteration
 L4: // New iteration
 
+   if (fDEBUG) {
+      std::cout << "New iteration - nfcn " << fNfcn << " nn1 = " << nn1 << " fGT " << fGT << "\n";
+      PrintVector("Parameter",n,fA);
+      PrintVector("PL ",n,fPL);
+      PrintVector("Step ",n,fDA);
+   }
+
    // fS - objective function value - zero first
    fS = 0.;
    // n0 - number of variable parameters in fit
    n0 = 0;
    for( i = 0; i < n; i++) {
       fGr[i]=0.; // zero gradients
-      if (fPL0[i] > .0) {
+      if (fPL0[i] > .0) {   // fLO is negstive foa fixed parameters
          n0=n0+1;
          // new iteration - new parallelepiped
          if (fPL[i] > .0) fPL0[i]=fPL[i];
@@ -1170,15 +1205,18 @@ Int_t TFumili::Minimize()
    if (ijkl == -1) fINDFLG[0]=1;
 
    // sp - scaled on fS machine precision
+   // compute precision epsilon (fRP is 1.E-15)
    sp=fRP*TMath::Abs(fS);
 
    // save fZ-matrix
    for( i=0; i < nn0; i++) fZ0[i] = fZ[i];
    if (nn3 > 0) {
       if (nn1 <= fNstepDec) {
+         // conditions on next step point
          t=2.*(fS-olds-fGT);
          if (fINDFLG[0] == 0) {
             if (TMath::Abs(fS-olds) <= sp && -fGT <= sp) goto L19;
+            if (fDEBUG) std::cout << " factor t is " << -fGT/t << std::endl;
             if(        0.59*t < -fGT) goto L19;
             t = -fGT/t;
             if (t < 0.25 ) t = 0.25;
@@ -1187,6 +1225,8 @@ Int_t TFumili::Minimize()
          fGT = fGT*t;
          t1 = t1*t;
          nn2=0;
+         // loop on parameter
+         // change parameter value using step  fPL and fDA
          for( i = 0; i < n; i++) {
             if (fPL[i] > 0.) {
                fA[i]=fA[i]-fDA[i];
@@ -1196,12 +1236,19 @@ Int_t TFumili::Minimize()
             }
          }
          nn1=nn1+1;
+         if (fDEBUG) {
+            std::cout << "change all PL and steps by factor t " << t << std::endl;
+         }
          goto L4;
       }
    }
 
 L19:
 
+   if (fDEBUG)
+      std::cout << "exit iteration (L19) t = " << t << " fGT = " << fGT << std::endl;
+
+   // flag here should be equal to zero
    if(fINDFLG[0] != 0) {
       fENDFLG=-4;
       printf("trying to execute an illegal jump at L85\n");
@@ -1224,6 +1271,8 @@ Int_t TFumili::Minimize()
             // or vice versa then fix parameter again and increment k1 by i1
             if ((fA[i] >= fAMX[i] && fGr[i] < 0.) ||
                    (fA[i] <= fAMN[i] && fGr[i] > 0.)) {
+               if (fDEBUG)
+                  std::cout << "Fix parameters with values outside boundary and negative derivative " << std::endl;
                fPL[i] = 0.;
                k1 = k1 + i1; // i1 stands for fZ-matrix row-number multiplier
                ///  - skip this row
@@ -1248,6 +1297,7 @@ Int_t TFumili::Minimize()
    }
 
    // INVERT fZ-matrix (mconvd() procedure)
+   // fZ stores lower diagonal part of symmetrix matrix
    i1 = 1;
    l  = 1;
    for( i = 0; i < n; i++) {// extract diagonal elements to fR-vector
@@ -1259,10 +1309,15 @@ Int_t TFumili::Minimize()
    }
 
    n0 = i1 - 1;
+   if (fDEBUG)
+      PrintMatrix("Approximate Hessian matrix  ",n0,fZ);
+
    InvertZ(n0);
 
    // fZ matrix now is inverted
    if (fINDFLG[0] != 0) { // problems
+      if (fDEBUG)
+         std::cout << "Some problems inverting the matrix - go back and try reducing again" << std::endl;
       // some PLs now have negative values, try to reduce fZ-matrix again
       fINDFLG[0] = 0;
       fINDFLG[1] = 1; // errors can be infinite
@@ -1271,6 +1326,12 @@ Int_t TFumili::Minimize()
       goto L19;
    }
 
+   if (fDEBUG) {
+      PrintMatrix("Approximate Covariance matrix (inverse of Hessian ) ",n0,fZ);
+      PrintVector("Current gradient",n,fGr);
+      PrintVector("Current step",n,fDA);
+   }
+
    // ... CALCULATE THEORETICAL STEP TO MINIMUM
    i1 = 1;
    for( i = 0; i < n; i++) {
@@ -1291,6 +1352,9 @@ Int_t TFumili::Minimize()
          i1=i1+1;
       }
    }
+   if (fDEBUG)
+      PrintVector("theoretical step (Newton)",n,fDA);
+
    //          ... CHECK FOR PARAMETERS ON BOUNDARY
 
    afix=0.;
@@ -1300,6 +1364,7 @@ Int_t TFumili::Minimize()
    for( i = 0; i < n; i++)
       if (fPL[i] > .0) {
          sigi = TMath::Sqrt(TMath::Abs(fZ[l - 1])); // calculate \sqrt{Z^{-1}_{ii}}
+         // original fR is diagonal of Hessian (not inverted matrix)
          fR[i] = fR[i]*fZ[l - 1];      // Z_ii * Z^-1_ii
          if (fEPS > .0) fParamError[i]=sigi;
          if ((fA[i] >= fAMX[i] && fDA[i] > 0.) || (fA[i] <= fAMN[i]
@@ -1313,6 +1378,9 @@ Int_t TFumili::Minimize()
                ifix=i;
                ifix1=i;
             }
+            if (fDEBUG)
+               std::cout << "Parameter " << i << " is outside bounds "
+                         << akap << "  " << afix << "  " << ifix << std::endl;
          }
          i1=i1+1;
          l=l+i1;
@@ -1323,6 +1391,10 @@ Int_t TFumili::Minimize()
       fPL[ifix] = -1.;
       fixFLG = fixFLG + 1;
       fi = 0;
+
+      if (fDEBUG) {
+         std::cout << "Parameter " << ifix << " is the worst - fix it and repeat step calculation " << std::endl;
+      }
       //.. REPEAT CALCULATION OF THEORETICAL STEP AFTER fiXING EACH PARAMETER
       goto L19;
    }
@@ -1351,11 +1423,17 @@ Int_t TFumili::Minimize()
          if ( bi > bm) {
             bi = bm;
             abi = abm;
+            if (fDEBUG)
+               std::cout << "Reduce parallelepide for "
+                        << i << " from " << fPL[i] << " to " << bi << std::endl;
          }
          // if calculated step is out of bounds
          if ( TMath::Abs(fDA[i]) > bi) {
             // decrease step splitter alambdA if needed
+            // by definition al is less than 1
             al = TMath::Abs(bi/fDA[i]);
+            if (fDEBUG)
+               std::cout << "step out of bound for  " << i <<  " reduce to " << al << std::endl;
             if (alambd > al) {
                imax=i;
                aiMAX=abi;
@@ -1367,6 +1445,9 @@ Int_t TFumili::Minimize()
          if (akap > fAKAPPA) fAKAPPA=akap;
       }
    }
+
+   if (fDEBUG)
+      std::cout << "Compute now corrected step - min found correction factor " << alambd << " max of akappa " << fAKAPPA << " nn2 " << nn2 << std::endl;
    //... CALCULATE NEW CORRECTED STEP
    fGT = 0.;
    amb = 1.e18;
@@ -1376,6 +1457,8 @@ Int_t TFumili::Minimize()
       if (fPL[i] > .0) {
          if (nn2 > fNlimMul ) {
             if (TMath::Abs(fDA[i]/fPL[i]) > amb ) {
+               if (fDEBUG)
+                  std::cout << "increase parallelipid 4-times for " << i << std::endl;
                fPL[i] = 4.*fPL[i]; // increase parallelepiped
                t1=4.; // flag - that fPL was increased
             }
@@ -1386,6 +1469,10 @@ Int_t TFumili::Minimize()
          fGT = fGT + fDA[i]*fGr[i];
       }
    }
+   if (fDEBUG) {
+      PrintVector("Reduced step",n,fDA);
+      std::cout << "after applying step reduction of " << alambd << " expected function change is " << fGT << std::endl;
+   }
 
    //.. CHECK IF MINIMUM ATTAINED AND SET EXIT MODE
    // if expected fGT smaller than precision
@@ -1406,6 +1493,8 @@ Int_t TFumili::Minimize()
                for( i = 0; i < fNpar; i++) fPL[i] = fPL0[i];
                fINDFLG[1] = 0;
                // and repeat iteration
+               if (fDEBUG)
+                  std::cout << "We have fixed some params - released and repeat " << std::endl;
                goto L19;
             } else {
                if( ifix1 >= 0) {
@@ -1466,6 +1555,8 @@ Int_t TFumili::Minimize()
       }
       return fENDFLG - 1;
    }
+   if (fDEBUG)
+      std::cout << "Continue and repeat iteration " << std::endl;
    goto L3;
 }