MnFunctionCross.cxx

Go to the documentation of this file.

// @(#)root/minuit2:$Id$
// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei   2003-2005
 
/**********************************************************************
 *                                                                    *
 * Copyright (c) 2005 LCG ROOT Math team,  CERN/PH-SFT                *
 *                                                                    *
 **********************************************************************/
 
#include "Minuit2/MnFunctionCross.h"
#include "Minuit2/FunctionMinimum.h"
#include "Minuit2/MnMigrad.h"
#include "Minuit2/FCNBase.h"
#include "Minuit2/MnParabola.h"
#include "Minuit2/MnParabolaPoint.h"
#include "Minuit2/MnParabolaFactory.h"
#include "Minuit2/MnCross.h"
#include "Minuit2/MnMachinePrecision.h"
 
//#define DEBUG
#include "Minuit2/MnPrint.h"
 
 
namespace ROOT {
 
   namespace Minuit2 {
 
 
 
MnCross MnFunctionCross::operator()(const std::vector<unsigned int>& par, const std::vector<double>& pmid,
                                    const std::vector<double>& pdir, double tlr, unsigned int maxcalls) const {
   // evaluate crossing point where function is equal to MIN + UP,
   // with direction pdir from values pmid
   // tlr indicate tolerance and maxcalls maximum number of calls
 
//   double edmmax = 0.5*0.001*toler*fFCN.Up();
 
 
   unsigned int npar = par.size();
   unsigned int nfcn = 0;
   const MnMachinePrecision& prec = fState.Precision();
   // tolerance used when calling Migrad
   double mgr_tlr = 0.5 * tlr;   // to be consistent with F77 version (for default values of tlr which is 0.1)
   // other olerance values are fixed at 0.01
   tlr = 0.01;
   // convergence when F is within tlf of aim and next prediction
   // of aopt is within tla of previous value of aopt
   double up = fFCN.Up();
   // for finding the point :
   double tlf = tlr*up;
   double tla = tlr;
   unsigned int maxitr = 15;
   unsigned int ipt = 0;
   double aminsv = fFval;
   double aim = aminsv + up;
   //std::cout<<"aim= "<<aim<<std::endl;
   double aopt = 0.;
   bool limset = false;
   std::vector<double> alsb(3, 0.), flsb(3, 0.);
 
   int printLevel = MnPrint::Level();
#ifdef DEBUG
   printLevel = 3;
#endif
 
#ifdef DEBUG
   for (unsigned int i = 0; i < par.size(); ++i)
         std::cout << "MnFunctionCross for parameter  " << par[i] << " value " << pmid[i] << " dir " << pdir[i]
         << " function min = " << aminsv << " contur value aim = (fmin + up) = " << aim << std::endl;
#endif
 
 
   // find the largest allowed aulim
 
   double aulim = 100.;
   for(unsigned int i = 0; i < par.size(); i++) {
      unsigned int kex = par[i];
      if(fState.Parameter(kex).HasLimits()) {
         double zmid = pmid[i];
         double zdir = pdir[i];
         //       double zlim = 0.;
         if(zdir > 0. && fState.Parameter(kex).HasUpperLimit()) {
            double zlim = fState.Parameter(kex).UpperLimit();
            if (fabs(zdir) < fState.Precision().Eps()) {
               // we have a limit
               if (fabs(zlim-zmid) < fState.Precision().Eps() ) limset = true;
               continue;
            }
            aulim = std::min(aulim, (zlim-zmid)/zdir);
         }
         else if(zdir < 0. && fState.Parameter(kex).HasLowerLimit()) {
            double zlim = fState.Parameter(kex).LowerLimit();
            if (fabs(zdir) < fState.Precision().Eps()) {
               // we have a limit
               if (fabs(zlim-zmid) < fState.Precision().Eps() ) limset = true;
               continue;
            }
            aulim = std::min(aulim, (zlim-zmid)/zdir);
         }
      }
   }
 
#ifdef DEBUG
   std::cout<<"Largest allowed aulim "<< aulim << std::endl;
#endif
 
 
   // case of a single parameter and we are at limit
   if (limset && npar == 1) {
      if (printLevel > 1)
         std::cout << "MnFunctionCross: parameter is at limit " << pmid[0] << " delta "
         << pdir[0] << std::endl;
      return MnCross(fState, nfcn, MnCross::CrossParLimit());
   }
 
   if (aulim < aopt + tla)
      limset = true;
 
   MnMigrad migrad(fFCN, fState, MnStrategy(std::max(0, int(fStrategy.Strategy()-1))));
 
   if (printLevel > 2) {
      std::cout << "MnFunctionCross: Run Migrad fixing parameters :";
      if (npar > 1) std::cout << std::endl;
   }
   for (unsigned int i = 0; i < npar; i++) {
 
      migrad.SetValue(par[i], pmid[i]);
 
      if (printLevel > 2) {
         std::cout << "\t" << fState.Name(par[i]) << " #" << par[i] << "  to " << pmid[i] << std::endl;
      }
   }
   // find minimum with respect all the other parameters (n- npar) (npar are the fixed ones)
 
   FunctionMinimum min0 = migrad(maxcalls, mgr_tlr);
   nfcn += min0.NFcn();
 
   if (printLevel > 2) {
      MnPrint::PrintState(std::cout, min0.State(), "MnFunctionCross: Result after Migrad");
      std::cout << min0.UserState().Parameters() << std::endl;
   }
 
   // case a new minimum is found
   if (min0.Fval() < fFval - tlf) {
      // case of new minimum is found
      if (printLevel > 1) {
         std::cout << "MnFunctionCross: A new minimum is found when scanning parameter  " << par.front() << " new value = "
         << min0.Fval() << " old value : " << fFval << std::endl;
      }
      return MnCross(min0.UserState(), nfcn, MnCross::CrossNewMin());
   }
   if(min0.HasReachedCallLimit())
      return MnCross(min0.UserState(), nfcn, MnCross::CrossFcnLimit());
   if(!min0.IsValid()) return MnCross(fState, nfcn);
   if(limset == true && min0.Fval() < aim)
      return MnCross(min0.UserState(), nfcn, MnCross::CrossParLimit());
 
   ipt++;
   alsb[0] = 0.;
   flsb[0] = min0.Fval();
   flsb[0] = std::max(flsb[0], aminsv + 0.1*up);
   aopt = sqrt(up/(flsb[0]-aminsv)) - 1.;
   if(fabs(flsb[0] - aim) < tlf) return MnCross(aopt, min0.UserState(), nfcn);
 
   if(aopt > 1.) aopt = 1.;
   if(aopt < -0.5) aopt = -0.5;
   limset = false;
   if(aopt > aulim) {
      aopt = aulim;
      limset = true;
   }
#ifdef DEBUG
   std::cout << "MnFunctionCross: flsb[0] = " << flsb[0] << " aopt =  " << aopt  << std::endl;
#endif
 
   if (printLevel > 2) {
      std::cout << "MnFunctionCross: Run Migrad again (#2) : ";
      if (npar > 1) std::cout << std::endl;
   }
   for (unsigned int i = 0; i < npar; i++) {
      migrad.SetValue(par[i], pmid[i] + (aopt)*pdir[i]);
      if (printLevel > 2) {
         std::cout << "\t - parameter " << i << " fixed to " << pmid[i] + (aopt)*pdir[i] << std::endl;
      }
   }
 
   FunctionMinimum min1 = migrad(maxcalls, mgr_tlr);
   nfcn += min1.NFcn();
 
   if (printLevel > 2) {
      MnPrint::PrintState(std::cout, min1.State(), "MnFunctionCross: Result after 2nd Migrad");
      std::cout << min1.UserState().Parameters() << std::endl;
   }
 
   if (min1.Fval() < fFval - tlf) // case of new minimum found
      return MnCross(min1.UserState(), nfcn, MnCross::CrossNewMin());
   if(min1.HasReachedCallLimit())
      return MnCross(min1.UserState(), nfcn, MnCross::CrossFcnLimit());
   if(!min1.IsValid()) return MnCross(fState, nfcn);
   if(limset == true && min1.Fval() < aim)
      return MnCross(min1.UserState(), nfcn, MnCross::CrossParLimit());
 
   ipt++;
   alsb[1] = aopt;
   flsb[1] = min1.Fval();
   double dfda = (flsb[1] - flsb[0])/(alsb[1] - alsb[0]);
 
#ifdef DEBUG
   std::cout << "aopt = " << aopt << " min1Val = " << flsb[1] << " dfda = " << dfda << std::endl;
#endif
 
 
L300:
      if(dfda < 0.) {
         // looking for slope of the right sign
#ifdef DEBUG
         std::cout << "MnFunctionCross: dfda < 0 - iterate from " << ipt << " to max of " << maxitr << std::endl;
#endif
         // iterate (max times is maxitr) incrementing aopt
 
         unsigned int maxlk = maxitr - ipt;
         for(unsigned int it = 0; it < maxlk; it++) {
            alsb[0] = alsb[1];
            flsb[0] = flsb[1];
            // LM: Add + 1, looking at Fortran code it starts from 1 ( see bug #8396)
            aopt = alsb[0] + 0.2*(it+1);
            limset = false;
            if(aopt > aulim) {
               aopt = aulim;
               limset = true;
            }
            if (printLevel > 2) {
               std::cout << "MnFunctionCross: Run Migrad again (iteration " << it << " ) :  ";
               if (npar > 1) std::cout << std::endl;
            }
            for(unsigned int i = 0; i < npar; i++) {
               migrad.SetValue(par[i], pmid[i] + (aopt)*pdir[i]);
               if (printLevel > 2) {
                  std::cout << "\t - parameter " << i << " fixed to " << pmid[i] + (aopt)*pdir[i] << std::endl;
               }
            }
            min1 = migrad(maxcalls, mgr_tlr);
            nfcn += min1.NFcn();
 
            if (printLevel > 2) {
               MnPrint::PrintState(std::cout, min1.State(), "MnFunctionCross: Result after Migrad");
               std::cout << min1.UserState().Parameters() << std::endl;
            }
 
            if (min1.Fval() < fFval - tlf) // case of new minimum found
               return MnCross(min1.UserState(), nfcn, MnCross::CrossNewMin());
            if(min1.HasReachedCallLimit())
               return MnCross(min1.UserState(), nfcn, MnCross::CrossFcnLimit());
            if(!min1.IsValid()) return MnCross(fState, nfcn);
            if(limset == true && min1.Fval() < aim)
               return MnCross(min1.UserState(), nfcn, MnCross::CrossParLimit());
            ipt++;
            alsb[1] = aopt;
            flsb[1] = min1.Fval();
            dfda = (flsb[1] - flsb[0])/(alsb[1] - alsb[0]);
            //       if(dfda > 0.) goto L460;
 
#ifdef DEBUG
   std::cout << "aopt = " << aopt << " min1Val = " << flsb[1] << " dfda = " << dfda << std::endl;
#endif
 
            if(dfda > 0.) break;
         }
         if(ipt > maxitr) return MnCross(fState, nfcn);
      } //if(dfda < 0.)
 
L460:
 
      // dfda > 0: we have two points with the right slope
 
      aopt = alsb[1] + (aim-flsb[1])/dfda;
 
#ifdef DEBUG
      std::cout << "MnFunctionCross: dfda > 0 : aopt = " << aopt << std::endl;
#endif
 
   double fdist = std::min(fabs(aim  - flsb[0]), fabs(aim  - flsb[1]));
   double adist = std::min(fabs(aopt - alsb[0]), fabs(aopt - alsb[1]));
   tla = tlr;
   if(fabs(aopt) > 1.) tla = tlr*fabs(aopt);
   if(adist < tla && fdist < tlf) return MnCross(aopt, min1.UserState(), nfcn);
   if(ipt > maxitr) return MnCross(fState, nfcn);
   double bmin = std::min(alsb[0], alsb[1]) - 1.;
   if(aopt < bmin) aopt = bmin;
   double bmax = std::max(alsb[0], alsb[1]) + 1.;
   if(aopt > bmax) aopt = bmax;
 
   limset = false;
   if(aopt > aulim) {
      aopt = aulim;
      limset = true;
   }
 
   if (printLevel > 2) {
      std::cout << "MnFunctionCross: Run Migrad again (#3) : ";
      if (npar > 1) std::cout << std::endl;
   }
   for(unsigned int i = 0; i < npar; i++) {
      migrad.SetValue(par[i], pmid[i] + (aopt)*pdir[i]);
      if (printLevel > 2) {
         std::cout << "\t : parameter " << i << " fixed to " << pmid[i] + (aopt)*pdir[i] << std::endl;
      }
   }
   FunctionMinimum min2 = migrad(maxcalls, mgr_tlr);
   nfcn += min2.NFcn();
 
   if (printLevel > 2) {
      MnPrint::PrintState(std::cout, min2.State(), "MnFunctionCross: Result after Migrad (#3): ");
      std::cout << min2.UserState().Parameters() << std::endl;
   }
 
   if (min2.Fval() < fFval - tlf) // case of new minimum found
      return MnCross(min2.UserState(), nfcn, MnCross::CrossNewMin());
   if(min2.HasReachedCallLimit())
      return MnCross(min2.UserState(), nfcn, MnCross::CrossFcnLimit());
   if(!min2.IsValid()) return MnCross(fState, nfcn);
   if(limset == true && min2.Fval() < aim)
      return MnCross(min2.UserState(), nfcn, MnCross::CrossParLimit());
 
   ipt++;
   alsb[2] = aopt;
   flsb[2] = min2.Fval();
 
   // now we have three points, ask how many < AIM
 
   double ecarmn = fabs(flsb[2] - aim);
   double ecarmx = 0.;
   unsigned int ibest = 2;
   unsigned int iworst = 0;
   unsigned int noless = 0;
 
   for(unsigned int i = 0; i < 3; i++) {
      double ecart = fabs(flsb[i] - aim);
      if(ecart > ecarmx) {
         ecarmx = ecart;
         iworst = i;
      }
      if(ecart < ecarmn) {
         ecarmn = ecart;
         ibest = i;
      }
      if(flsb[i] < aim) noless++;
   }
 
#ifdef DEBUG
   std::cout << "MnFunctionCross: have three points : nless < aim  = " << noless << " ibest = " << ibest << " iworst = " << iworst << std::endl;
#endif
 
   //std::cout<<"480"<<std::endl;
 
   // at least one on each side of AIM (contour), fit a parabola
   if(noless == 1 || noless == 2) goto L500;
   // if all three are above AIM, third point must be the closest to AIM, return it
   if(noless == 0 && ibest != 2) return MnCross(fState, nfcn);
   // if all three below and third is not best then the slope has again gone negative,
   // re-iterate and look for positive slope
   if(noless == 3 && ibest != 2) {
      alsb[1] = alsb[2];
      flsb[1] = flsb[2];
#ifdef DEBUG
   std::cout << "MnFunctionCross: all three points below - look again fir positive slope " << std::endl;
#endif
      goto L300;
   }
 
   // in other case new straight line thru first two points
 
   flsb[iworst] = flsb[2];
   alsb[iworst] = alsb[2];
   dfda = (flsb[1] - flsb[0])/(alsb[1] - alsb[0]);
#ifdef DEBUG
   std::cout << "MnFunctionCross: new straight line using point 1-2 - dfda =  " << dfda << std::endl;
#endif
   goto L460;
 
L500:
 
      do {
         // do parabola fit
         MnParabola parbol = MnParabolaFactory()(MnParabolaPoint(alsb[0], flsb[0]), MnParabolaPoint(alsb[1], flsb[1]), MnParabolaPoint(alsb[2], flsb[2]));
         //   aopt = parbol.X_pos(aim);
         //std::cout<<"alsb1,2,3= "<<alsb[0]<<", "<<alsb[1]<<", "<<alsb[2]<<std::endl;
         //std::cout<<"flsb1,2,3= "<<flsb[0]<<", "<<flsb[1]<<", "<<flsb[2]<<std::endl;
 
#ifdef DEBUG
   std::cout << "MnFunctionCross: parabola fit: iteration " << ipt  << std::endl;
#endif
 
         double coeff1 = parbol.C();
         double coeff2 = parbol.B();
         double coeff3 = parbol.A();
         double determ = coeff2*coeff2 - 4.*coeff3*(coeff1 - aim);
 
#ifdef DEBUG
         std::cout << "MnFunctionCross: parabola fit: a =  " << coeff3  << " b = "
                   << coeff2 << " c = " << coeff1 << " determ = " << determ << std::endl;
#endif
         // curvature is negative
         if(determ < prec.Eps()) return MnCross(fState, nfcn);
         double rt = sqrt(determ);
         double x1 = (-coeff2 + rt)/(2.*coeff3);
         double x2 = (-coeff2 - rt)/(2.*coeff3);
         double s1 = coeff2 + 2.*x1*coeff3;
         double s2 = coeff2 + 2.*x2*coeff3;
 
#ifdef DEBUG
         std::cout << "MnFunctionCross: parabola fit: x1 =  " << x1  << " x2 = "
                   << x2 << " s1 = " << s1 << " s2 = " << s2 << std::endl;
#endif
 
#ifdef WARNINGMSG
         if(s1*s2 > 0.)   MN_INFO_MSG("MnFunctionCross problem 1");
#endif
         // find with root is the right one
         aopt = x1;
         double slope = s1;
         if(s2 > 0.) {
            aopt = x2;
            slope = s2;
         }
#ifdef DEBUG
         std::cout << "MnFunctionCross: parabola fit: aopt =  " << aopt  << " slope = "
                   << slope << std::endl;
#endif
 
         // ask if converged
         tla = tlr;
         if(fabs(aopt) > 1.) tla = tlr*fabs(aopt);
 
#ifdef DEBUG
         std::cout << "MnFunctionCross: Delta(aopt) =  " << fabs(aopt - alsb[ibest])  << " tla = "
                   << tla << "Delta(F) = " << fabs(flsb[ibest] - aim) << " tlf = " << tlf << std::endl;
#endif
 
 
         if(fabs(aopt - alsb[ibest]) < tla && fabs(flsb[ibest] - aim) < tlf)
            return MnCross(aopt, min2.UserState(), nfcn);
 
         //     if(ipt > maxitr) return MnCross();
 
         // see if proposed point is in acceptable zone between L and R
         // first find ileft, iright, iout and ibest
 
         unsigned int ileft = 3;
         unsigned int iright = 3;
         unsigned int iout = 3;
         ibest = 0;
         ecarmx = 0.;
         ecarmn = fabs(aim-flsb[0]);
         for(unsigned int i = 0; i < 3; i++) {
            double ecart = fabs(flsb[i] - aim);
            if(ecart < ecarmn) {
               ecarmn = ecart;
               ibest = i;
            }
            if(ecart > ecarmx) ecarmx = ecart;
            if(flsb[i] > aim) {
               if(iright == 3) iright = i;
               else if(flsb[i] > flsb[iright]) iout = i;
               else {
                  iout = iright;
                  iright = i;
               }
            } else if(ileft == 3) ileft = i;
            else if(flsb[i] < flsb[ileft]) iout = i;
            else {
               iout = ileft;
               ileft = i;
            }
         }
 
#ifdef DEBUG
         std::cout << "MnFunctionCross: ileft =  " << ileft  << " iright = "
                   << iright << " iout = " << iout << " ibest = " << ibest << std::endl;
#endif
 
         // avoid keeping a bad point nest time around
 
         if(ecarmx > 10.*fabs(flsb[iout] - aim))
            aopt = 0.5*(aopt + 0.5*(alsb[iright] + alsb[ileft]));
 
         // knowing ileft and iright, get acceptable window
         double smalla = 0.1*tla;
         if(slope*smalla > tlf) smalla = tlf/slope;
         double aleft = alsb[ileft] + smalla;
         double aright = alsb[iright] - smalla;
 
         // move proposed point AOPT into window if necessary
         if(aopt < aleft) aopt = aleft;
         if(aopt > aright) aopt = aright;
         if(aleft > aright) aopt = 0.5*(aleft + aright);
 
         // see if proposed point outside limits (should be impossible)
         limset = false;
         if(aopt > aulim) {
            aopt = aulim;
            limset = true;
         }
 
         // evaluate at new point aopt
         if (printLevel > 2) {
            std::cout << "MnFunctionCross: Run Migrad again at new point (#iter = " << ipt << " ): ";
            if (npar > 1) std::cout << std::endl;
         }
         for(unsigned int i = 0; i < npar; i++) {
            migrad.SetValue(par[i], pmid[i] + (aopt)*pdir[i]);
            if (printLevel > 2) {
               std::cout << "\t : parameter " << i << " fixed to " << pmid[i] + (aopt)*pdir[i] << std::endl;
            }
         }
         min2 = migrad(maxcalls, mgr_tlr);
         nfcn += min2.NFcn();
 
         if (printLevel > 2) {
            MnPrint::PrintState(std::cout, min2.State(), "MnFunctionCross: Result after new Migrad : ");
            std::cout << min2.UserState().Parameters() << std::endl;
         }
 
         if (min2.Fval() < fFval - tlf) // case of new minimum found
            return MnCross(min2.UserState(), nfcn, MnCross::CrossNewMin());
         if(min2.HasReachedCallLimit())
            return MnCross(min2.UserState(), nfcn, MnCross::CrossFcnLimit());
         if(!min2.IsValid()) return MnCross(fState, nfcn);
         if(limset == true && min2.Fval() < aim)
            return MnCross(min2.UserState(), nfcn, MnCross::CrossParLimit());
 
         ipt++;
         // replace odd point with new one (which is the best of three)
         alsb[iout] = aopt;
         flsb[iout] = min2.Fval();
         ibest = iout;
      } while(ipt < maxitr);
 
   // goto L500;
 
   return MnCross(fState, nfcn);
}
 
   }  // namespace Minuit2
 
}  // namespace ROOT