CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_6/src/HiggsAnalysis/CombinedLimit/src/SequentialMinimizer.cc

Go to the documentation of this file.
00001 #include "../interface/SequentialMinimizer.h"
00002 
00003 #include <cmath>
00004 #include <stdexcept>
00005 #include <memory>
00006 #include <algorithm>
00007 #include <limits>
00008 #include "TString.h"
00009 #include "RooRealVar.h"
00010 #include "RooAbsReal.h"
00011 #include "RooLinkedListIter.h"
00012 #include <Math/MinimizerOptions.h>
00013 #include <Math/Factory.h>
00014 #include <boost/foreach.hpp>
00015 #include "../interface/ProfilingTools.h"
00016 #define foreach BOOST_FOREACH
00017 
00018 #define DEBUG_ODM_printf if (0) printf
00019 //#define DEBUG_SM_printf   if (0) printf
00020 //#define DEBUGV_SM_printf  if (0) printf
00021 // #define DEBUG_ODM_printf printf
00022 #define DEBUG_SM_printf  if (fDebug > 1) printf
00023 #define DEBUGV_SM_printf if (fDebug > 2) printf  
00024 
00025 namespace { 
00026     const double GOLD_R1 = 0.61803399 ;
00027     const double GOLD_R2 = 1-0.61803399 ;
00028     const double XTOL    = 10*std::sqrt(std::numeric_limits<double>::epsilon());
00029 }
00030 void cmsmath::OneDimMinimizer::initDefault(const MinimizerContext &ctx, unsigned int idx) {
00031     init(ctx, idx, -std::numeric_limits<double>::infinity(), std::numeric_limits<double>::infinity(), 1.0, Form("#%d", idx));
00032 }
00033 
00034 
00035 bool cmsmath::OneDimMinimizer::minimize(int steps, double ytol, double xtol) 
00036 {
00037     // get initial bracket
00038     seek();
00039     
00040     bool done = doloop(steps,ytol,xtol);
00041     parabolaStep();
00042     x() = xi_[1];
00043     xstep_ = xi_[2] - xi_[0];
00044     return done;
00045 }
00046 
00047 cmsmath::OneDimMinimizer::ImproveRet cmsmath::OneDimMinimizer::improve(int steps, double ytol, double xtol, bool force) 
00048 {
00049     double x0 = x();
00050     if (x0 < xi_[0] || x0 > xi_[2]) {
00051         // could happen if somebody outside this loop modifies some parameters
00052         DEBUG_ODM_printf("ODM: ALERT: variable %s outside bounds x = [%.4f, %.4f, %.4f], x0 = %.4f\n", name_.c_str(), xi_[0], xi_[1], xi_[2], x0);
00053         x0 = x() = xi_[1]; 
00054     } else {
00055         xi_[1] = x0;
00056     }
00057     double y0 = eval();
00058     yi_[1] = y0;
00059     yi_[0] = eval(xi_[0]);
00060     yi_[2] = eval(xi_[2]);
00061     if (xtol == 0) xtol = (fabs(xi_[1])+XTOL)*XTOL;
00062 
00063     DEBUG_ODM_printf("ODM: start of improve %s x = [%.4f, %.4f, %.4f], y = [%.4f, %.4f, %.4f]\n", name_.c_str(), xi_[0], xi_[1], xi_[2], yi_[0], yi_[1], yi_[2]);
00064 
00065     if (yi_[1] <= yi_[0] && yi_[1] <= yi_[2]) {
00066         if (ytol > 0 && (std::max(yi_[2],yi_[0]) - yi_[1]) < ytol) {
00067             DEBUG_ODM_printf("ODM: immediate ytol for %s: ymin %.8f, ymax %.8f, diff %.8f\n", name_.c_str(), yi_[1], std::max(yi_[2],yi_[0]), std::max(yi_[2],yi_[0]) - yi_[1]);
00068             if (!force || parabolaStep()) return Unchanged;
00069         }
00070         if (xtol > 0 && (xi_[2] - xi_[0]) < xtol) {
00071             DEBUG_ODM_printf("ODM: immediate xtol for %s: xmin %.8f, xmax %.8f, diff %.8f\n", name_.c_str(), xi_[0], xi_[2], xi_[2] - xi_[0]);
00072             if (!force || parabolaStep()) return Unchanged;
00073         }
00074     } else {
00075         reseek();
00076     }
00077 
00078     //post-condition: always a sorted interval
00079     assert(xi_[0] < xi_[2]);
00080     assert(xi_[0] <= xi_[1] && xi_[1] <= xi_[2]);
00081     // if midpoint is not not one of the extremes, it's not higher than that extreme
00082     assert(xi_[1] == xi_[0] || yi_[1] <= yi_[0]); 
00083     assert(xi_[1] == xi_[2] || yi_[1] <= yi_[2]);
00084 
00085     bool done = doloop(steps,ytol,xtol);
00086     parabolaStep(); 
00087 
00088     //post-condition: always a sorted interval
00089     assert(xi_[0] < xi_[2]);
00090     assert(xi_[0] <= xi_[1] && xi_[1] <= xi_[2]);
00091     // if midpoint is not not one of the extremes, it's not higher than that extreme
00092     assert(xi_[1] == xi_[0] || yi_[1] <= yi_[0]); 
00093     assert(xi_[1] == xi_[2] || yi_[1] <= yi_[2]);
00094 
00095 
00096     if (ytol > 0 && fabs(yi_[1] - y0) < ytol) {
00097         DEBUG_ODM_printf("ODM: final ytol for %s: ymin(old) %.8f, ymin(new) %.8f, diff %.8f\n", name_.c_str(), y0, yi_[1], y0 -yi_[1]);
00098         if (!force) x() = x0;
00099         return Unchanged;
00100     } 
00101 
00102     if (xtol > 0 && fabs(xi_[1] - x0) < xtol) {
00103         x() = (force ? xi_[1] : x0);
00104         return Unchanged;
00105     }
00106     DEBUG_ODM_printf("ODM: doloop for %s is %s\n", name_.c_str(), done ? "Done" : "NotDone");
00107     DEBUG_ODM_printf("ODM: end of improve %s x = [%.4f, %.4f, %.4f], y = [%.4f, %.4f, %.4f]\n", name_.c_str(), xi_[0], xi_[1], xi_[2], yi_[0], yi_[1], yi_[2]);
00108     x() = xi_[1];
00109     xstep_ = xi_[2] - xi_[0];
00110     return done ? Done : NotDone; 
00111 }
00112 
00113 bool cmsmath::OneDimMinimizer::doloop(int steps, double ytol, double xtol) 
00114 {
00115     if (steps <= 0) steps = 100;
00116     for (int i = 0; i < steps; ++i) {
00117         if (xtol > 0 && (xi_[2] - xi_[0]) < xtol) {
00118             return true;
00119         }
00120         goldenBisection();
00121         if (ytol > 0 && (std::max(yi_[2],yi_[0]) - yi_[1]) < ytol) {
00122             DEBUG_ODM_printf("ODM: intermediate ytol for %s: ymin %.8f, ymax %.8f, diff %.8f\n", name_.c_str(), yi_[1], std::max(yi_[2],yi_[0]), std::max(yi_[2],yi_[0]) - yi_[1]);
00123             return true;
00124         }
00125         DEBUG_ODM_printf("ODM: step %d/%d done for %s: ymin %.8f, ymax %.8f, diff %.8f\n", i+1, steps, name_.c_str(), yi_[1], std::max(yi_[2],yi_[0]), std::max(yi_[2],yi_[0]) - yi_[1]);
00126         DEBUG_ODM_printf("ODM: %s x = [%.4f, %.4f, %.4f], y = [%.4f, %.4f, %.4f]\n", name_.c_str(), xi_[0], xi_[1], xi_[2], yi_[0], yi_[1], yi_[2]);
00127     }
00128 
00129     return false;
00130 }
00131 
00132 void cmsmath::OneDimMinimizer::seek() 
00133 {
00134     if (std::isfinite(xmax_-xmin_)) {
00135         xstep_ = std::max(xstep_, 0.2*(xmax_-xmin_));
00136     } else {
00137         xstep_ = 1.0;
00138     }
00139     reseek();
00140 }
00141 void cmsmath::OneDimMinimizer::reseek() 
00142 {
00143     double xtol2 = 2*(fabs(xi_[1])+XTOL)*XTOL;
00144     if (xstep_ < xtol2) xstep_ = xtol2;
00145     xi_[1] = x(); 
00146     yi_[1] = eval(xi_[1]);
00147     xi_[0] = std::max(xmin_, xi_[1]-xstep_);
00148     yi_[0] = (xi_[0] == xi_[1] ? yi_[1] : eval(xi_[0]));
00149     xi_[2] = std::min(xmax_, xi_[1]+xstep_);
00150     yi_[2] = (xi_[2] == xi_[1] ? yi_[1] : eval(xi_[2]));
00151     assert(xi_[0] < xi_[2]);
00152     assert(xi_[0] <= xi_[1] && xi_[1] <= xi_[2]);
00153 
00154     for (;;) {
00155         //DEBUG_ODM_printf("ODM: bracketing %s in x = [%.4f, %.4f, %.4f], y = [%.4f, %.4f, %.4f]\n", name_.c_str(), xi_[0], xi_[1], xi_[2], yi_[0], yi_[1], yi_[2]);
00156         if (yi_[0] < yi_[1]) {
00157             assign(2,1); // 2:=1
00158             assign(1,0); // 1:=0
00159             xi_[0] = std::max(xmin_, xi_[1]-xstep_);
00160             yi_[0] = (xi_[0] == xi_[1] ? yi_[1] : eval(xi_[0]));
00161         } else if (yi_[2]  < yi_[1]) {
00162             assign(0,1); // 0:=1
00163             assign(1,2); // 1:=2
00164             xi_[2] = std::min(xmax_, xi_[1]+xstep_);
00165             yi_[2] = (xi_[2] == xi_[1] ? yi_[1] : eval(xi_[2]));
00166         } else if (yi_[2] == yi_[1] && yi_[1] == yi_[0]) {
00167             // function is identical in three points --> constant?
00168             // try a scan
00169             const int nPoints = 20;
00170             double xi[nPoints], yi[nPoints];
00171             double dx = (xmax_-xmin_)/nPoints, x = xmin_ - 0.5*dx;
00172             bool isConst = true;
00173             int iFound = 0; 
00174             for (int i = 0; i < nPoints; ++i, x += dx) {
00175                 xi[i] = x; yi[i] = eval(x);
00176                 if (yi[i] != yi_[1]) isConst = false;
00177                 if (yi[i] < yi[iFound]) { iFound = i; }
00178             }
00179             if (isConst) break;
00180             xi_[0] = (iFound == 0        ? xmin_ : xi[iFound-1]);
00181             xi_[2] = (iFound > nPoints-1 ? xmax_ : xi[iFound+1]);
00182             xi_[1] = iFound; yi_[1] = yi_[iFound];
00183             break;
00184         } else {
00185             xstep_ /= 2;
00186             break;
00187         }
00188         xstep_ *= 2;
00189     }
00190     //DEBUG_ODM_printf("ODM: bracketed minimum of %s in [%.4f, %.4f, %.4f]\n", name_.c_str(), xi_[0], xi_[1], xi_[2]);
00191     //post-condition: always a sorted interval
00192     assert(xi_[0] < xi_[2]);
00193     assert(xi_[0] <= xi_[1] && xi_[1] <= xi_[2]);
00194     // if midpoint is not not one of the extremes, it's not higher than that extreme
00195     assert(xi_[1] == xi_[0] || yi_[1] <= yi_[0]); 
00196     assert(xi_[1] == xi_[2] || yi_[1] <= yi_[2]);
00197 }
00198 
00199 void cmsmath::OneDimMinimizer::goldenBisection() 
00200 {
00201     //pre-condition: always a sorted interval
00202     assert(xi_[0] < xi_[2]);
00203     assert(xi_[0] <= xi_[1] && xi_[1] <= xi_[2]);
00204     // if midpoint is not not one of the extremes, it's not higher than that extreme
00205     assert(xi_[1] == xi_[0] || yi_[1] <= yi_[0]); 
00206     assert(xi_[1] == xi_[2] || yi_[1] <= yi_[2]);
00207     if (xi_[0] == xi_[1] || xi_[1] == xi_[2]) {
00208         int isame = (xi_[0] == xi_[1] ? 0 : 2);
00210         assert(yi_[isame] <= yi_[2-isame]);
00211         xi_[1] = 0.5*(xi_[0]+xi_[2]);
00212         yi_[1] = eval(xi_[1]);
00213         if (yi_[1] < yi_[isame]) {
00214             // maximum is in the interval-
00215             // leave as is, next bisection steps will find it
00216         } else {
00217             // maximum remains on the boundary, leave both points there
00218             assign(2-isame, 1);
00219             assign(1, isame); 
00220         }
00221     } else {
00222         int inear = 2, ifar = 0;
00223         if (xi_[2]-xi_[1] > xi_[1] - xi_[0]) {
00224             inear = 0, ifar = 2;
00225         } else {
00226             inear = 2, ifar = 0;
00227         }
00228         double xc = xi_[1]*GOLD_R1 + xi_[ifar]*GOLD_R2;
00229         double yc = eval(xc);
00230         //DEBUG_ODM_printf("ODM: goldenBisection:\n\t\tfar = (%.2f,%.8f)\n\t\tnear = (%.2f,%.8f)\n\t\tcenter  = (%.2f,%.8f)\n\t\tnew  = (%.2f,%.8f)\n",
00231         //            xi_[ifar],  yi_[ifar], xi_[inear], yi_[inear], xi_[1], yi_[1], xc, yc);
00232         if (yc < yi_[1]) {   // then use 1, c, far
00233             assign(inear, 1);
00234             xi_[1] = xc; yi_[1] = yc;
00235         } else {  // then use c, 1, near
00236             xi_[ifar] = xc; yi_[ifar] = yc;
00237         }
00238     }
00239     //post-condition: always a sorted interval
00240     assert(xi_[0] < xi_[2]);
00241     assert(xi_[0] <= xi_[1] && xi_[1] <= xi_[2]);
00242     // if midpoint is not not one of the extremes, it's not higher than that extreme
00243     assert(xi_[1] == xi_[0] || yi_[1] <= yi_[0]); 
00244     assert(xi_[1] == xi_[2] || yi_[1] <= yi_[2]);
00245 
00246 }
00247 
00248 double cmsmath::OneDimMinimizer::parabolaFit() 
00249 {
00250     if (xi_[0] == xi_[1] || xi_[1] == xi_[2]) { 
00251         return xi_[1]; 
00252     }
00253     double dx0 = xi_[1] - xi_[0], dx2 = xi_[1] - xi_[2];
00254     double dy0 = yi_[1] - yi_[0], dy2 = yi_[1] - yi_[2];
00255     double num = dx0*dx0*dy2 - dx2*dx2*dy0;
00256     double den = dx0*dy2     - dx2*dy0;
00257     if (den != 0) {
00258         double x = xi_[1] - 0.5*num/den;
00259         if (xi_[0] < x && x < xi_[2]) {
00260             return x;
00261         }
00262     } 
00263     return xi_[1];
00264 }
00265 
00266 bool cmsmath::OneDimMinimizer::parabolaStep() {
00267     double xc = parabolaFit();
00268     if (xc != xi_[1]) {
00269         double yc = eval(xc);
00270         if (yc < yi_[1]) {
00271             xi_[1] = xc; 
00272             yi_[1] = yc;
00273             //post-condition: always a sorted interval
00274             assert(xi_[0] < xi_[2]);
00275             assert(xi_[0] <= xi_[1] && xi_[1] <= xi_[2]);
00276             // if midpoint is not not one of the extremes, it's not higher than that extreme
00277             assert(xi_[1] == xi_[0] || yi_[1] <= yi_[0]); 
00278             assert(xi_[1] == xi_[2] || yi_[1] <= yi_[2]);
00279             return true;
00280         }
00281     }
00282     return false;
00283 }
00284 
00285 void cmsmath::OneDimMinimizer::moveTo(double x) {
00286     if (x == xmax_) {
00287         xi_[0] = xmax_ - (xi_[2]-xi_[0]); yi_[0] = eval(xi_[0]);
00288         xi_[1] = xmax_; yi_[1] = eval(xmax_);
00289         xi_[2] = xmax_; yi_[2] = yi_[1];
00290     } else if (x == xmin_) {
00291         xi_[2] = xmin_ + (xi_[2]-xi_[0]); yi_[2] = eval(xi_[0]);
00292         xi_[1] = xmin_; yi_[1] = eval(xmin_);
00293         xi_[0] = xmin_; yi_[0] = yi_[1];
00294     } else {
00295         double dx = xi_[2] - xi_[0];
00296         xi_[1] = x; yi_[0] = eval(x);
00297         xi_[0] = std::max(xmin_, x-dx); yi_[0] = eval(xi_[0]);
00298         xi_[2] = std::min(xmax_, x+dx); yi_[2] = eval(xi_[2]);
00299     }
00300 }
00301 
00302 void cmsmath::SequentialMinimizer::SetFunction(const ROOT::Math::IMultiGenFunction & func) {
00303     DEBUG_SM_printf("SequentialMinimizer::SetFunction: nDim = %u\n", func.NDim());
00304     func_.reset(new MinimizerContext(&func));
00305     nFree_ = nDim_ = func_->x.size();
00306     // create dummy workers
00307     workers_.clear();
00308     workers_.resize(nDim_);
00309     // reset states
00310     Clear();
00311 }
00312 
00313 void cmsmath::SequentialMinimizer::Clear() {
00314     DEBUGV_SM_printf("SequentialMinimizer::Clear()\n");
00315     minValue_ = std::numeric_limits<double>::quiet_NaN();
00316     edm_      = std::numeric_limits<double>::infinity();
00317     state_ = Cleared;
00318     foreach(Worker &w, workers_) w.state = Cleared;
00319 }
00320 
00321 bool cmsmath::SequentialMinimizer::SetVariable(unsigned int ivar, const std::string & name, double val, double step) {
00322     DEBUGV_SM_printf("SequentialMinimizer::SetVariable(idx %u, name %s, val %g, step %g)\n", ivar, name.c_str(), val, step);
00323     assert(ivar < nDim_);
00324     func_->x[ivar] = val;
00325     workers_[ivar].initUnbound(*func_, ivar, step, name);
00326     workers_[ivar].state = Cleared;
00327     return true;
00328 }
00329 
00330 bool cmsmath::SequentialMinimizer::SetLimitedVariable(unsigned int ivar, const std::string & name, double val, double  step, double  lower, double  upper) {
00331     DEBUGV_SM_printf("SequentialMinimizer::SetLimitedVariable(idx %u, name %s, var %g, step %g, min %g, max %g)\n", ivar, name.c_str(), val, step, lower, upper);
00332     assert(ivar < nDim_);
00333     func_->x[ivar] = val;
00334     workers_[ivar].init(*func_, ivar, lower, upper, step, name);
00335     workers_[ivar].state = Cleared;
00336     return true;
00337 }
00338 
00339 bool cmsmath::SequentialMinimizer::SetFixedVariable(unsigned int ivar, const std::string & name, double val) {
00340     DEBUGV_SM_printf("SequentialMinimizer::SetFixedVariable(idx %u, name %s, var %g)\n", ivar, name.c_str(), val);
00341     assert(ivar < nDim_);
00342     func_->x[ivar] = val;
00343     workers_[ivar].initUnbound(*func_, ivar, 1.0, name);
00344     workers_[ivar].state = Fixed;
00345     return true;
00346 }
00347 
00348 bool cmsmath::SequentialMinimizer::Minimize() {
00349     return minimize();
00350 }
00351 
00352 bool cmsmath::SequentialMinimizer::minimize(int smallsteps) 
00353 {
00354     for (unsigned int i = 0; i < nDim_; ++i) {
00355         Worker &w = workers_[i];
00356         if (!w.isInit() || w.state == Unknown) throw std::runtime_error(Form("SequentialMinimizer::worker[%u/%u] not initialized!\n", i, nDim_));
00357         if (w.state != Fixed) {
00358             w.minimize(1); 
00359             w.state = Ready; 
00360         }
00361     }
00362     state_ = Ready;
00363     return improve(smallsteps);
00364 }
00365 
00366 bool cmsmath::SequentialMinimizer::improve(int smallsteps)
00367 {
00368     static int nFailWakeUpAttempts = runtimedef::get("SeqMinimizer_nFailWakeUpAttempts");
00369 
00370     // catch improve before minimize case
00371     if (state_ == Cleared) return minimize(smallsteps);
00372 
00373     // setup default tolerances and steps
00374     double ytol = Tolerance()/sqrt(workers_.size());
00375     int bigsteps = MaxIterations()*20;
00376 
00377     // list of done workers (latest-done on top)
00378     std::list<Worker*> doneWorkers;
00379 
00380     // start with active workers, for all except constants
00381     foreach(Worker &w, workers_) {
00382         if (w.state != Fixed) w.state = Active;
00383     }
00384 
00385     state_ = Active;
00386     for (int i = 0; i < bigsteps; ++i) {
00387         DEBUG_SM_printf("Start of loop. Strategy %d, State is %s\n",fStrategy,(state_ == Done ? "DONE" : "ACTIVE"));
00388         State newstate = Done;
00389         int oldActiveWorkers = 0, newActiveWorkers = 0;
00390         foreach(Worker &w, workers_) {
00391             OneDimMinimizer::ImproveRet iret = OneDimMinimizer::Unchanged;
00392             if (w.state == Done || w.state == Fixed) continue;
00393             iret = w.improve(smallsteps,ytol); 
00394             oldActiveWorkers++; 
00395             if (iret == OneDimMinimizer::Unchanged) {
00396                 DEBUGV_SM_printf("\tMinimized %s:  Unchanged. NLL = %.8f\n", w.cname(), func_->eval());
00397                 w.state = Done;
00398                 w.nUnaffected = 0;
00399                 doneWorkers.push_front(&w);
00400             } else {
00401                 DEBUGV_SM_printf("\tMinimized %s:  Changed. NLL = %.8f\n", w.cname(), func_->eval());
00402                 w.state = Active;
00403                 newstate = Active;
00404                 newActiveWorkers++;
00405             }
00406         }
00407         if (fStrategy >= 2 && newActiveWorkers <= 30) { // arbitrary cut-off
00408             DEBUG_SM_printf("Middle of loop. Strategy %d, active workers %d: firing full minimizer\n", fStrategy, newActiveWorkers);
00409             if (doFullMinim()) newstate = Done;
00410         }
00411         if (newstate == Done) {
00412             DEBUG_SM_printf("Middle of loop. Strategy %d, State is %s, active workers %d --> %d \n",fStrategy,(state_ == Done ? "DONE" : "ACTIVE"), oldActiveWorkers, newActiveWorkers);
00413             oldActiveWorkers = 0; newActiveWorkers = 0;
00414             double y0 = func_->eval();
00415             std::list<Worker*>::iterator it = doneWorkers.begin();
00416             // The topmost worker was added on the list just now, so by definition it's already done.
00417             // We save a reference to it, remove it from the done list, and if the loop doesn't end there we set it to active again 
00418             Worker* firstWorker = *it; 
00419             it = doneWorkers.erase(it);
00420             // Then we check all the others
00421             while( it != doneWorkers.end()) {
00422                 Worker &w = **it;
00423                 if (nFailWakeUpAttempts && w.nUnaffected >= nFailWakeUpAttempts) { ++it; continue; }
00424                 OneDimMinimizer::ImproveRet iret = w.improve(smallsteps,ytol,0,/*force=*/true);
00425                 oldActiveWorkers++;
00426                 if (iret == OneDimMinimizer::Unchanged) {
00427                     DEBUGV_SM_printf("\tMinimized %s:  Unchanged. NLL = %.8f\n", w.cname(), func_->eval());
00428                     w.nUnaffected++;
00429                     ++it;
00430                 } else {
00431                     DEBUGV_SM_printf("\tMinimized %s:  Changed. NLL = %.8f\n", w.cname(), func_->eval());
00432                     w.state = Active;
00433                     newstate = Active;
00434                     w.nUnaffected = 0;
00435                     it = doneWorkers.erase(it);
00436                     newActiveWorkers++;
00437                     if (fStrategy == 0) break;
00438                 }
00439             }
00440             if (newstate == Active) { // wake up him too
00441                 firstWorker->state = Active; 
00442                 firstWorker->nUnaffected = 0; 
00443             }
00444             double y1 = func_->eval();
00445             edm_ = y0 - y1;
00446         }
00447         DEBUG_SM_printf("End of loop. Strategy %d, State is %s, active workers %d --> %d \n",fStrategy,(state_ == Done ? "DONE" : "ACTIVE"), oldActiveWorkers, newActiveWorkers);
00448         if (state_ == Done && newstate == Done) {
00449             DEBUG_SM_printf("Converged after %d big steps\n",i);
00450             minValue_ = func_->eval();
00451             fStatus   = 0;
00452             return true;
00453         }
00454         state_ = newstate;
00455         if (func_->nCalls > MaxFunctionCalls()) break;
00456     }
00457     DEBUG_SM_printf("Failed do converge after %d big steps\n",bigsteps);
00458     fStatus   = -1;
00459     minValue_ = func_->eval();
00460     return false;
00461 }
00462 
00463 namespace cmsmath {
00464     class SubspaceMultiGenFunction : public ROOT::Math::IMultiGenFunction {
00465         public:
00466            SubspaceMultiGenFunction(const ROOT::Math::IMultiGenFunction *f, int nDim, const int *idx, double *xi) :
00467                 f_(f), nDim_(nDim), idx_(idx), x_(xi) {}
00468            virtual IBaseFunctionMultiDim * Clone() const { return new SubspaceMultiGenFunction(*this); }
00469            virtual unsigned int NDim() const { return nDim_; }
00470         private:
00471            virtual double DoEval(const double * x) const {
00472                 for (int i = 0; i < nDim_; ++i) x_[idx_[i]] = x[i];
00473                 return (*f_)(x_);
00474            }
00475            const ROOT::Math::IMultiGenFunction *f_; 
00476            const int  nDim_;
00477            const int *idx_;
00478            double    *x_;
00479     };
00480 }
00481 
00482 bool cmsmath::SequentialMinimizer::doFullMinim()
00483 {
00484     if (fullMinimizer_.get() == 0) {
00485         fullMinimizer_.reset(ROOT::Math::Factory::CreateMinimizer("Minuit2", ""));
00486         fullMinimizer_->SetTolerance(Tolerance());
00487         fullMinimizer_->SetStrategy(Strategy()-2);
00488     }
00489     subspaceIndices_.clear();
00490     for (int i = 0, n = workers_.size(); i < n; ++i) {
00491         if (workers_[i].state == Active) subspaceIndices_.push_back(i);
00492     }
00493     fullMinimizer_->Clear();
00494     SubspaceMultiGenFunction subfunc(func_->func, subspaceIndices_.size(), &subspaceIndices_[0], &func_->x[0]);
00495     fullMinimizer_->SetFunction(subfunc);
00496     for (int i = 0, n = subspaceIndices_.size(); i < n; ++i) {
00497         int j = subspaceIndices_[i];
00498         Worker &w = workers_[j];
00499         fullMinimizer_->SetLimitedVariable(i, w.name(), func_->x[j], w.step(), w.min(), w.max());
00500     }
00501     bool ok = fullMinimizer_->Minimize();
00502     if (ok) {
00503         const double *ximin = fullMinimizer_->X();
00504         // move to the right place
00505         for (int i = 0, n = subspaceIndices_.size(); i < n; ++i) {
00506             func_->x[subspaceIndices_[i]] = ximin[i];
00507         }
00508         // update all workers
00509         for (int i = 0, n = subspaceIndices_.size(); i < n; ++i) {
00510             int j = subspaceIndices_[i];
00511             Worker &w = workers_[j];
00512             w.moveTo( ximin[i] );
00513         }
00514 
00515     }
00516     return ok;
00517 }
00518 
00519 
00520 #include <TPluginManager.h>
00521 namespace {
00522     static int load_seqmin() {
00523         gPluginMgr->AddHandler("ROOT::Math::Minimizer", "SeqMinimizer", "cmsmath::SequentialMinimizer", "HiggsAnalysisCombinedLimit", "SequentialMinimizer(const char *)");
00524         return 1;
00525     }
00526     static int loaded_seqmin = load_seqmin();
00527 }