MnFunctionCross.cxx

Go to the documentation of this file.

// @(#)root/minuit2:$Id: MnFunctionCross.cxx 30097 2009-09-09 16:47:06Z moneta $
// 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

#if defined(DEBUG) || defined(WARNINGMSG)
#include "Minuit2/MnPrint.h" 
#endif


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();
   //   double 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();
   double tlf = tlr*up;
   // tolerance used when calling Migrad
   double mgr_tlr = 0.05 * up;   // to be consistent with F77 version 
   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.);

#ifdef DEBUG
   std::cout<<"MnFunctionCross for parameter  "<<par.front()<< "fmin = " << 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];
         if(fabs(zdir) < fState.Precision().Eps()) continue;
         //       double zlim = 0.;
         if(zdir > 0. && fState.Parameter(kex).HasUpperLimit()) {
            double zlim = fState.Parameter(kex).UpperLimit();
            aulim = std::min(aulim, (zlim-zmid)/zdir);
         }
         else if(zdir < 0. && fState.Parameter(kex).HasLowerLimit()) {
            double zlim = fState.Parameter(kex).LowerLimit();
            aulim = std::min(aulim, (zlim-zmid)/zdir);
         }
      }
   }

#ifdef DEBUG
   std::cout<<"Largest allowed aulim "<< aulim << std::endl;
#endif
   
   if(aulim  < aopt+tla) limset = true;

   
   MnMigrad migrad(fFCN, fState, MnStrategy(std::max(0, int(fStrategy.Strategy()-1))));
   
   for(unsigned int i = 0; i < npar; i++) {
#ifdef DEBUG
      std::cout << "MnFunctionCross: Set value for " << par[i] <<  " to " << pmid[i] << std::endl; 
#endif
      migrad.SetValue(par[i], pmid[i]);
   }
   // find minimum with respect all the other parameters (n- npar) (npar are the fixed ones)
   
   FunctionMinimum min0 = migrad(maxcalls, mgr_tlr);
   nfcn += min0.NFcn();

#ifdef DEBUG
   std::cout << "MnFunctionCross: after Migrad on n-1  minimum is " << min0 << std::endl; 
#endif
   
   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
   
   for(unsigned int i = 0; i < npar; i++) {
#ifdef DEBUG
      std::cout << "MnFunctionCross: Set new value for " << par[i] <<  " from " << pmid[i] << " to " << pmid[i] + (aopt)*pdir[i] << " aopt = " << aopt << std::endl; 
#endif
      migrad.SetValue(par[i], pmid[i] + (aopt)*pdir[i]);
   }
   
   FunctionMinimum min1 = migrad(maxcalls, mgr_tlr);
   nfcn += min1.NFcn();

#ifdef DEBUG
   std::cout << "MnFunctionCross: after Migrad on n-1  minimum is " << min1 << std::endl; 
#endif
   
   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;
            }
            for(unsigned int i = 0; i < npar; i++) {
#ifdef DEBUG
      std::cout << "MnFunctionCross: Set new value for " << par[i] <<  " to " << pmid[i] + (aopt)*pdir[i] << " aopt = " << aopt << std::endl; 
#endif
               migrad.SetValue(par[i], pmid[i] + (aopt)*pdir[i]);
            }
            min1 = migrad(maxcalls, mgr_tlr);
            nfcn += min1.NFcn();

#ifdef DEBUG
   std::cout << "MnFunctionCross: after Migrad on n-1  minimum is " << min1 << std::endl; 
   std::cout << "nfcn = " << nfcn << std::endl;
#endif
            
            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;
   }
   
   for(unsigned int i = 0; i < npar; i++) {
#ifdef DEBUG
      std::cout << "MnFunctionCross: Set new value for " << par[i] <<  " from " << pmid[i] << " to " << pmid[i] + (aopt)*pdir[i] << " aopt = " << aopt << std::endl; 
#endif
      migrad.SetValue(par[i], pmid[i] + (aopt)*pdir[i]);
   }
   FunctionMinimum min2 = migrad(maxcalls, mgr_tlr);
   nfcn += min2.NFcn();

#ifdef DEBUG
   std::cout << "MnFunctionCross: after Migrad on n-1  minimum is " << min2 << std::endl; 
   std::cout << "nfcn = " << nfcn << std::endl;
#endif
   
   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
         for(unsigned int i = 0; i < npar; i++) {
#ifdef DEBUG
      std::cout << "MnFunctionCross: Set new value for " << par[i] <<  " from " << pmid[i] << " to " << pmid[i] + (aopt)*pdir[i] << " aopt = " << aopt << std::endl; 
#endif   
            migrad.SetValue(par[i], pmid[i] + (aopt)*pdir[i]);
         }
         min2 = migrad(maxcalls, mgr_tlr);
         nfcn += min2.NFcn();
#ifdef DEBUG
   std::cout << "MnFunctionCross: after Migrad on n-1  minimum is " << min2 << std::endl; 
   std::cout << "nfcn = " << nfcn << std::endl;
#endif
         
         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

---