NCFitter.cxx

Go to the documentation of this file.

// $Id: NCFitter.cxx,v 1.17 2009/06/18 18:12:22 bckhouse Exp $
//
// Collection of functions and classes for fitting and contour finding.
//
// C. Backhouse 11/2007

#ifndef STANDALONE

#include "NCUtils/Fitter/NCFitter.h"

#include "NCUtils/NCUtility.h"
using NC::Utility::SQR;
using NC::Utility::ReportProgress;
#include "MessageService/MsgService.h"

#else

#include "NCFitter.h"

#include <iostream>
using namespace std;
#define MSG(x,y) cout << "NCCont: "
#define CVSID(x) ;
template<class T> inline const T SQR(const T a) {return a*a;}

#endif // standalone

#include "TStopwatch.h"
#include "TMath.h"

#include <cassert>
#include <cmath>
#include <sstream>

CVSID("$Id: NCFitter.cxx,v 1.17 2009/06/18 18:12:22 bckhouse Exp $");



namespace NC{
namespace Fitter{

  /*
  CoordInt ConvertCoord(const CoordNDim& c, std::vector<NCParameter> ps)
  {
    assert(c.size() == ps.size());
  }
  */

  CoordNDim ConvertCoord(const CoordInt& c, std::vector<NCParameter> ps)
  {
    assert(c.size() == ps.size());

    CoordNDim ret(c.size());
    for(unsigned int n = 0; n < c.size(); ++n)
      ret[n] = ps[n].ValFromStep(c[n]);

    return ret;
  }

  CoordNDim ConvertContourCoord(const CoordNDim& c, vector<NCParameter> ps)
  {
    assert(c.size() == ps.size());

    CoordNDim ret(c.size());
    for(unsigned int n = 0; n < c.size(); ++n)
      ret[n] = ps[n].ValFromStep(c[n]);

    return ret;
  }


  double Discretizer::EvalAt(const CoordInt& r) const
  {
    // TODO - this loses the ret information
    return fFunc->EvalAtEx(ConvertCoord(r, fParams), 0);
  }


  double Penalizer::EvalAtEx(const CoordNDim& r, CoordNDim* retCoord) const
  {
    CoordNDim reval = r;

    // If the point is outside the physical region - move it to the boundary
    for(unsigned int n = 0; n < fParams.size(); ++n){
      if(r[n] < fParams[n].Min()) reval[n] = fParams[n].Min();
      if(r[n] > fParams[n].Max()) reval[n] = fParams[n].Max();
    }

    double ret = fFunc->EvalAtEx(reval, retCoord);

    // If we had to move the point then apply a quadratic penalty term.
    // The factor of 1000 is fairly arbitrary.
    const double amplitude = 1e3;
    for(unsigned int n = 0; n < fParams.size(); ++n){
      const double wsqr = SQR(fParams[n].Min()-fParams[n].Max());
      assert(wsqr);
      if(r[n] < fParams[n].Min())
        ret += amplitude*SQR(fParams[n].Min()-r[n])/wsqr;
      if(r[n] > fParams[n].Max())
        ret += amplitude*SQR(r[n]-fParams[n].Max())/wsqr;
    }

    return ret;
  }


  double PenalizerInt::EvalAt(const CoordInt& r) const
  {
    CoordInt reval = r;

    // If the point is outside the physical region - move it to the boundary
    for(unsigned int n = 0; n < fParams.size(); ++n){
      if(r[n] < 0) reval[n] = 0;
      if(r[n] > fParams[n].Steps()) reval[n] = fParams[n].Steps();
    }

    double ret = fFunc->EvalAt(reval);

    // If we had to move the point then apply a quadratic penalty term.
    // The factor of 1000 is fairly arbitrary.
    const double amplitude = 1e3;
    for(unsigned int n = 0; n < fParams.size(); ++n){
      const double wsqr = SQR(fParams[n].Steps());
      assert(wsqr);
      if(r[n] < 0)
        ret += amplitude*SQR(r[n])/wsqr;
      if(r[n] > fParams[n].Steps())
        ret += amplitude*SQR(r[n]-fParams[n].Steps())/wsqr;
    }

    return ret;
  }


  double FuncCache::EvalAtEx(const CoordNDim& r, CoordNDim* retCoord) const
  {
    const std::map<CoordNDim, double>::const_iterator it = fCache.find(r);
    if(it != fCache.end()){
      ++fNhit;
      if(retCoord) *retCoord = fCacheRetCoord[r];
      return it->second;
    }

    ++fNmiss;

    //  const unsigned int maxCacheSize = int(1e7); // on the order of 10meg
    //  if(cache.size() > maxCacheSize) fCache.clear();

    CoordNDim myRetCoord;
    const double ret = fFunc->EvalAtEx(r, &myRetCoord);

    fCache[r] = ret;
    fCacheRetCoord[r] = myRetCoord;

    if(retCoord) *retCoord = myRetCoord;
    return ret;
  }

  ostream& operator<<(ostream& os, const FuncCache& c)
  {
    os << "FuncCache for callable at " << c.fFunc
       << " with " << c.fNhit+c.fNmiss << " calls "
       << "of which " << c.fNhit << " hit and "
       << c.fNmiss << " missed." << endl;
    return os;
  }


  double FuncCacheInt::EvalAtEx(const CoordInt& r, CoordInt* ret) const
  {
    const std::map<CoordInt, double>::const_iterator it = fCache.find(r);

    if(it != fCache.end()){
      ++fNhit;
      if(ret) *ret = fRetCache[r];
      return it->second;
    }

    ++fNmiss;

    //  const unsigned int maxCacheSize = int(1e7); // on the order of 10meg
    //  if(cache.size() > maxCacheSize) cache.clear();

    CoordInt myret;
    const double val = fFunc->EvalAtEx(r, &myret);
    if(ret) *ret = myret;

    fCache[r] = val;
    fRetCache[r] = myret;

    return val;
  }

  ostream& operator<<(ostream& os, const FuncCacheInt& c)
  {
    os << "FuncCacheInt for callable at " << c.fFunc
       << " with " << c.fNhit+c.fNmiss << " calls "
       << "of which " << c.fNhit << " hit and "
       << c.fNmiss << " missed." << endl;
    return os;
  }

  /*
  vector<CoordInt> FuncCacheInt::GetHistory()
  {
    vector<CoordInt> ret;
    map<CoordInt, double>::iterator it;
    for(it = cache.begin(); it != cache.end(); ++it) ret.push_back(it->first);
    return ret;
  }
  */


  //------------------------- MinFinderSimple ------------------------------

  vector<string> MinFinderSimple::sfPartNames;
  vector<double> MinFinderSimple::sfPartPositions;

  MinFinderSimple::MinFinderSimple(const std::vector<NCParameter>& ps,
                                   const bool limits,
                                   const bool reuse,
                                   const ICallableND* f)
    :fParams(ps),
     fRespectParameterLimits(limits),
     reuseOld(reuse),
     fFunc(f)
  {
  }


  CoordNDim MinFinderSimple::BeginSearchPos() const
  {
    assert(fPrevBest.size() == fParams.size() || fPrevBest.size() == 0);

    if(!reuseOld || fPrevBest.empty()){
      CoordNDim ret;
      for(unsigned int n = 0; n < fParams.size(); ++n){
        ret.push_back((fParams[n].Min()+fParams[n].Max())/2);
      }
      return ret;
    }

    return fPrevBest;
  }


  int MinFinderSimple::BeginScale() const
  {
    int initScale = 1;
    if(!reuseOld || fPrevBest.empty()){
      bool bigEnough = false;
      while(!bigEnough){
        initScale *= 2;
        bigEnough = true;
        for(unsigned int n = 0; n < fParams.size(); ++n){
          if(initScale*fParams[n].Precision() <
             (fParams[n].Max()-fParams[n].Min())/4)
            bigEnough = false;
        }
      }
    }
    return initScale;
  }

  void MinFinderSimple::SetSpacePartitions(vector<string> names,
                                           vector<double> poss)
  {
    assert(names.size() == poss.size());
    sfPartNames = names;
    sfPartPositions = poss;
  }

  // ---------------------------------------------------------------------
  double MinFinderSimple::FindMinPartitions(CoordNDim& ret,
                                            bool progress) const
  {
    assert(sfPartNames.size() == sfPartPositions.size());
    assert(!sfPartNames.empty());

    // Start with `this` as a template to copy the minimizers from
    vector<const MinFinderSimple*> minimizers;
    minimizers.push_back(this);

    for(unsigned int n = 0; n < sfPartNames.size(); ++n){
      for(unsigned int i = 0; i < fParams.size(); ++i){
        // If we find we're fitting a parameter we ought to be partitioning
        if(fParams[i].ShortName() == sfPartNames[n]){
          // Take every minimizer that already exists and split it into
          // two, one for each side of the partition
          const unsigned int M = minimizers.size();
          for(unsigned int m = 0; m < M; ++m){
            // The parameter ranges inform the fitter where to begin
            vector<NCParameter> ps1 = fParams;
            ps1[i].SetRange(fParams[i].Min(), sfPartPositions[n]);
            vector<NCParameter> ps2 = fParams;
            ps2[i].SetRange(sfPartPositions[n], fParams[i].Max());

            // Constrain the ranges in [-infinity, partition point]
            // and [partition point, +infinity]
            ICallableND* f1 = new ConstrainRange(minimizers[m]->fFunc, i,
                                                 -1e50, sfPartPositions[n]);
            ICallableND* f2 = new ConstrainRange(minimizers[m]->fFunc, i,
                                                 sfPartPositions[n], +1e50);

            MinFinderSimple* m1 = new MinFinderSimple(ps1,
                                                      fRespectParameterLimits,
                                                      reuseOld,
                                                      f1);
            MinFinderSimple* m2 = new MinFinderSimple(ps2,
                                                      fRespectParameterLimits,
                                                      reuseOld,
                                                      f2);
            minimizers.push_back(m1);
            minimizers.push_back(m2);
          } // end for m
          // The first M of the minimizers are the ones before we split them
          // in half. We don't need them anymore.
          for(unsigned int m = 0; m < M; ++m){
            if(minimizers[m] != this) delete minimizers[m];
            minimizers.erase(minimizers.begin());
          } // end for m
        } // end if matches
      } // end for i
    } // end for n

    // Minimize using all the minimizers and pick the best result. Ensure
    // that ret is filled properly too. We pass true for ignorePartitions
    // otherwise we'd end up here again instead of actually minimizing.
    double lowest = 1e50;
    for(unsigned int m = 0; m < minimizers.size(); ++m){
      CoordNDim thisRet;
      const double height = minimizers[m]->FindMin(thisRet, progress, true);
      if(height < lowest){
        lowest = height;
        ret = thisRet;
      }
      // Don't need it anymore
      if(minimizers[m] != this) delete minimizers[m];
    }

    return lowest;
  }

  // ---------------------------------------------------------------------
  double MinFinderSimple::FindMin(CoordNDim& ret,
                                  bool progress,
                                  bool ignorePartitions) const
  {
    if(!sfPartNames.empty() && !ignorePartitions)
      return FindMinPartitions(ret, progress);

    ret.resize(fParams.size());

    CoordNDim curPos = BeginSearchPos();
    assert(curPos.size() == fParams.size());

    CoordNDim junk;
    double curHeight = fFunc->EvalAtEx(curPos, &junk);
    /*
    if(curHeight > 9e5 && !fPrevBest.empty()){   // Appear to have started inside some kind of penalty zone
      fPrevBest.clear();    // So, forget where we were before
      // and try again.
      curPos = BeginSearchPos();
      curHeight = func->EvalAt(curPos);
    }
    */

    std::vector<int> curIntPos(curPos.size(), 0);

    int initScale = BeginScale();

    std::map<std::vector<int>, bool> alreadyBeen;
    unsigned int lastDir = 0;

    CoordNDim curEx;

    TStopwatch sw;

    for(int scale = initScale; scale >= 1; scale /= 2){
      if(progress) ReportProgress(1-log((float)scale)/log((float)initScale), sw);
    scaleTop:
      MSG("NCFitter", Msg::kVerbose) << "MinSimple: at "
                                     << curPos << endl;
      // Try the previous direction first
      for(int tryLast = true; tryLast >= false; --tryLast){
        // What direction to look in
        for(unsigned int dir = 0; dir < curPos.size(); ++dir){
          // In what sense
          for(int sign = -1; sign <= +1; sign += 2){
            if(tryLast && dir != lastDir) continue;

            CoordNDim newPos = curPos;
            newPos[dir] += sign*fParams[dir].Precision()*scale;

            //      if(fRespectParameterLimits && (newPos[dir] < fParams[dir].Min() || newPos[dir] > fParams[dir].Max()))
            //continue;

            std::vector<int> newIntPos = curIntPos;
            newIntPos[dir] += sign*scale;

            if(alreadyBeen[newIntPos]) continue;

            alreadyBeen[newIntPos] = true;

            CoordNDim ex;

            double heightHere = fFunc->EvalAtEx(newPos, &ex);

            if(heightHere < curHeight){
              curHeight = heightHere;
              curPos = newPos;
              curIntPos = newIntPos;
              curEx = ex;
              lastDir = dir;

              //              cout << "Moved to ";
              //              for (unsigned int n = 0; n < fParams.size(); ++n)
              //                cout << fParams[n].ShortName() << "=" << curPos[n] << " ";
              //              cout << endl;

              goto scaleTop; // Repeat all these loops with the same value of scale.
            } // end if better
          } // end for s
        } //end for d
      } // end for tryLast
    } // end for scale

    //    if(curHeight < 9e5)
    fPrevBest = curPos;

    if(curEx.empty()) ret = curPos; else ret = curEx;

    return curHeight;
  }


  //------------------------- MinFinderSimpleInt ------------------------------
  CoordInt MinFinderSimpleInt::BeginSearchPos() const
  {
    assert(fPrevBest.size() == fParams.size() || fPrevBest.empty());

    if(fPrevBest.empty()){
      CoordInt ret;
      for(unsigned int n = 0; n < fParams.size(); ++n){
        //      ret.push_back((fParams[n].Min()+fParams[n].Max())/2);
        //      ret.push_back(fParams[n].Min()+.25*(fParams[n].Max()-fParams[n].Min()));
        ret.push_back(int((fParams[n].Max()-fParams[n].Min())/
                          fParams[n].Precision())/2);
      }
      return ret;
    }

    return fPrevBest;
  }


  CoordInt MinFinderSimpleInt::BeginScale() const
  {
    CoordInt ret(fParams.size());
    for(unsigned int n = 0; n < ret.size(); ++n) ret[n] = 1;

    if(!fPrevBest.empty()) return ret;

    for(unsigned int n = 0; n < fParams.size(); ++n){
      bool bigEnough = false;
      while(!bigEnough){
        ret[n] *= 2;
        bigEnough = ret[n] > .25*fParams[n].Steps();
      }
    }

    return ret;
  }


  // ---------------------------------------------------------------------
  double MinFinderSimpleInt::FindMin(CoordInt* ret,
                                     vector<TGraph>* min_graphs) const
  {
    // We need this because eg we move one way and
    // then immediately test where we just were.
    FuncCacheInt f(fFunc);

    if(ret) ret->resize(fParams.size());

    // One extra for chisq value
    if(min_graphs) min_graphs->resize(1+fParams.size());

    CoordInt curPos = BeginSearchPos();
    double curHeight = f.EvalAt(curPos);

    // Appear to have started inside some kind of penalty zone
    if(curHeight > 9e5 && !fPrevBest.empty()){
      fPrevBest.clear();    // So, forget where we were before
      // and try again.
      curPos = BeginSearchPos();
      curHeight = f.EvalAt(curPos);
    }

    CoordInt scale = BeginScale();

    // What way did we last move?
    unsigned int lastDir = 0;
    int lastSign = 0;

    // How many times in a row we have moved in this direction, at this scale
    vector<int> numGoodMoves(fParams.size(), 0);

    while(true){
    loopTop:
      MSG("NCFitter", Msg::kVerbose) << "MinSimpleInt: at "
                                     << curPos << endl;

      // Try the previous direction first
      for(int tryLast = true; tryLast >= false; --tryLast){
        // What direction to look in
        for(unsigned int dir = 0; dir < curPos.size(); ++dir){
          // In what sense
          for(int sign = -1; sign <= +1; sign += 2){
            if(tryLast && (dir != lastDir || sign != lastSign)) continue;

            CoordInt newPos = curPos;
            newPos[dir] += sign*scale[dir];

            // TODO - wrong
            //      if(fRespectParameterLimits && (newPos[dir] < fParams[dir].Min() || newPos[dir] > fParams[dir].Max()))
            //        continue;

            const double heightHere = f.EvalAt(newPos);

            if(heightHere < curHeight){
              if(min_graphs){
                min_graphs->at(0).SetPoint(min_graphs->at(0).GetN(),
                                           min_graphs->at(0).GetN(),
                                           heightHere);

                for(unsigned int n = 0; n < newPos.size(); ++n)
                  min_graphs->at(n+1).SetPoint(min_graphs->at(n+1).GetN(),
                                               min_graphs->at(n+1).GetN(),
                                               newPos[n]);
              }

              curHeight = heightHere;
              curPos = newPos;
              lastDir = dir;
              lastSign = sign;

              ++numGoodMoves[dir];
              // As soon as we move twice in a direction,
              // then double the scale that way.
              // The tryLast stuff will make sure we try again
              // in this direction first.
              if(numGoodMoves[dir] == 2){
                numGoodMoves[dir] = 0;
                scale[dir] *= 2;
              }
              goto loopTop; // Try again from this new position
            } // end if better
          } // end for s
        } //end for d
      } // end for tryLast

      // If we got here then no direction was an improvement.
      // If we are at the minimum scale then quit.
      bool quit = true;
      for(unsigned int n = 0; n < scale.size(); ++n) if(scale[n] != 1) quit = false;
      if(quit) break;

      // Otherwise reduce the scale in every direction we can.
      for(unsigned int n = 0; n < scale.size(); ++n){
        if(scale[n] > 1){
          scale[n] /= 2;
          numGoodMoves[n] = 0;
        }
      }

    } // end infinite loop

    if(curHeight < 9e5) fPrevBest = curPos;
    if(ret) *ret = curPos;

    return curHeight;
  }

#ifndef STANDALONE
  // ---------------------------------------------------------------------
  double MinFinderMinuit(const ICallableND* chiSqrFunc,
                         const std::vector<NCParameter> params,
                               CoordNDim& ret,
                               TMinuit** ppMinu)
  {
    ret.resize(params.size());

    TMinuit* minu = new MinuitForCallable(chiSqrFunc, params.size());
    minu->SetPrintLevel(-1); // Shut Minuit up.
    minu->Command("set strategy 2"); // Be careful, instead of fast.

    for(int n = 0; n < int(params.size()); ++n)
      // NB - no parameter limits imposed.
      minu->DefineParameter(n, params[n].ShortName().c_str(),
                            (params[n].Min()+params[n].Max())/2,
                            (params[n].Max()-params[n].Min())/4,
                            0, 0);

    // Minuit tries to return a lot of parameters we don't want to keep.
    double junkd; int junki; TString junks;

    if(minu->Migrad()){
      MSG("NCFitter", Msg::kWarning) << "Migrad failed, "
                                            << "falling back to Simplex.\n";
      // Things are going wrong, we should let the user see.
      minu->SetPrintLevel(0);
      minu->mnsimp();
      minu->Migrad();
    }

    double minChiSqr;
    minu->mnstat(minChiSqr, junkd, junkd, junki, junki, junki);

    for(unsigned int n = 0; n < params.size(); ++n)
      minu->GetParameter(n, ret[n], junkd);

    if(ppMinu)
      *ppMinu = minu;
    else
      delete minu;

    return minChiSqr;
  }
#endif


  // ---------------------------------------------------------------------
  std::vector<Coord> FindContour(const ICallable2D* func,
                                 const Coord& minCoord,
                                 const double height,
                                 const Coord& precision,
                                 std::vector<TGraph>* marg_graphs)
  {
    assert(!marg_graphs || marg_graphs->empty());

    const double dx = precision.x;
    const double dy = precision.y;
    assert(dx > 0); assert(dy > 0);

    typedef pair<int, int> Box;    // first == left, second == top

    std::map<Box, bool> results;   // Is this box already in the result set?
    std::vector<Box> todo;         // Boxes to investigate
    // Boxes whose topleft corner has already been evaluated, used as cache.
    std::map<Box, double> heights;

    // Top left corner of the box containing the best fit.
    const int bfxi = int(minCoord.x/dx);
    const int bfyi = int(minCoord.y/dy);

    MSG("NCFitter", Msg::kVerbose) << "Contour: walking left" << endl;
    // Walk to the left
    for(int x = bfxi; ; --x){
      // As soon as we exceed the contour
      MSG("NCFitter", Msg::kVerbose) << "Contour: at " << x << endl;
      CoordNDim bestfit;
      if(func->EvalAtEx(Coord(x*dx, bfyi*dy), bestfit) > height){
        // Create enough graphs to save our marginalisation progress
        //(2 extra to save the two variables actually in the contour)
        if(marg_graphs)
          for(unsigned int n = 0; n < 4*bestfit.size(); ++n)
            marg_graphs->push_back(TGraph());

        // make it the first element in the todo list
        todo.push_back(Box(x, bfyi));
        // and stop going left.
        break;
      }
    }

    // NB - for some reason while writing this I decided
    // that increasing y is downward.

    std::vector<Coord> contour;

    TStopwatch sw; // starts automatically

    int gxi = -1; // position on the x axis of the marg graphs.

    while(!todo.empty()){
      Box now = todo.back();
      todo.pop_back();

      ++gxi;

      // We end up evaluating the function multiple times at each point, so are
      // relying heavily on the cacheing done in GetChiSqrFromDerived()

      MSG("NCFitter", Msg::kVerbose) << "Contour: at "
                                     << now.first << ","
                                     << now.second << endl;

      // Set by the loop below.
      double topleft, topright, bottomleft, bottomright;

      // Evaluate the function at all the corners.
      for(int cx = 0; cx <= 1; ++cx){
        for(int cy = 0; cy <= 1; ++cy){
          Box corner(now.first+cx, now.second+cy);
          const std::map<Box, double>::const_iterator it = heights.find(corner);
          double heightHere;
          if(it == heights.end()){
            CoordNDim bestfit;
            heightHere = func->EvalAtEx(Coord((now.first+cx)*dx,
                                              (now.second+cy)*dy), bestfit);

            /* TODO - hacked out marg graphs because they are broken somehow. FIX!
            if(marg_graphs){
              const int dn = cx+2*cy;

              for(unsigned int n = 0; n < bestfit.size(); ++n){
                marg_graphs->at(4*n+dn).SetPoint(marg_graphs->at(4*n+dn).
                                                 GetN(), gxi, bestfit[n]);
              }
            }
            */

            heights[corner] = heightHere;
          }
          else{
            heightHere = heights[corner];
          }
          if(cx == 0 && cy == 0) topleft = heightHere;
          if(cx == 1 && cy == 0) topright = heightHere;
          if(cx == 0 && cy == 1) bottomleft = heightHere;
          if(cx == 1 && cy == 1) bottomright = heightHere;
        }
      }

      MSG("NCFitter", Msg::kVerbose) << "Contour: tl,tr,bl,br="
                                     << topleft << " "
                                     << topright << " "
                                     << bottomleft << " "
                                     << bottomright << endl;

      if(now.first-bfxi){
        const double est_frac = (M_PI+atan2(now.second-bfyi-.01,
                                            double(now.first-bfxi)))/(2*M_PI);
        ReportProgress(est_frac, sw);
      }


      // For each side of the box: see if one corner is high and the other low.
      // If so, and if the box adjacent on this side isn't already part of the
      // results, then append the adjacent box to the todo list and linearly
      // interpolate to find where to put the contour point on that edge.

      if((topright-height)*(bottomright-height) <= 0){
        Box right(now.first+1, now.second);
        if(!results[right]){
          MSG("NCFitter", Msg::kVerbose) << "Contour: went right" << endl;
          todo.push_back(right);
          results[right] = true;
          const double rightdiff = (bottomright-topright) ?
            bottomright-topright : 1;
          contour.push_back(Coord((now.first+1)*dx,
                                  (now.second+(height-topright)/rightdiff)*dy));
        }
      }
      if( (topleft-height)*(bottomleft-height) <= 0){
        Box left(now.first-1, now.second);
        if(!results[left]){
          MSG("NCFitter", Msg::kVerbose) << "Contour: went left" << endl;
          todo.push_back(left);
          results[left] = true;
          const double leftdiff = (bottomleft-topleft)
            ? bottomleft-topleft : 1;
          contour.push_back(Coord(now.first*dx,
                                  (now.second+(height-topleft)/leftdiff)*dy));
        }
      }
      if( (bottomleft-height)*(bottomright-height) <= 0){
        Box down(now.first, now.second+1);
        if(!results[down]){
          MSG("NCFitter", Msg::kVerbose) << "Contour: went down" << endl;
          todo.push_back(down);
          results[down] = true;
          const double bottomdiff = (bottomright-bottomleft) ?
            bottomright-bottomleft : 1;
          contour.push_back(Coord((now.first+
                                   (height-bottomleft)/(bottomdiff))*dx,
                                  (now.second+1)*dy));
        }
      }
      if( (topleft-height)*(topright-height) <= 0 ){
        Box up(now.first, now.second-1);
        if(!results[up]){
          MSG("NCFitter", Msg::kVerbose) << "Contour: went up" << endl;
          todo.push_back(up);
          results[up] = true;
          const double topdiff = (topright-topleft) ? topright-topleft : 1;
          contour.push_back(Coord((now.first+(height-topleft)/topdiff)*dx,
                                  now.second*dy));
        }
      }
    }

    ReportProgress(1, sw);

    if(contour.size()) contour.push_back(contour[0]); // Close the loop.
    return contour;
  }


  //----------------------------------------------------------------------
  TMatrix FindCovarianceMatrix(const ICallableND* func,
                               const CoordNDim& minCoord,
                               const std::vector<NCParameter> params)
  {
    assert(params.size() == minCoord.size());

    TMatrix ret(params.size(), params.size());

    CoordNDim junk;
    const double bestFitVal = func->EvalAtEx(minCoord, &junk);

    // First, do the 2nd derivatives in one variable
    for(unsigned int i = 0; i < params.size(); ++i){
      const double step = params[i].Precision();
      CoordNDim plus = minCoord;
      plus[i] += step;
      CoordNDim minus = minCoord;
      minus[i] -= step;
      const double dii = (func->EvalAtEx(plus, &junk)+
                          func->EvalAtEx(minus, &junk)-
                          2*bestFitVal)/(step*step);
      ret(i, i) = dii;
    }

    // Now do all the off-diagonal derivatives
    for(unsigned int i = 0; i < params.size()-1; ++i){
      for(unsigned int j = i+1; j < params.size(); ++j){
        const double di = params[i].Precision();
        const double dj = params[j].Precision();

        CoordNDim pp = minCoord;
        pp[i] += di; pp[j] += dj;
        CoordNDim pi = minCoord;
        pi[i] += di;
        CoordNDim pj = minCoord;
        pj[j] += dj;

        const double dij = (func->EvalAtEx(pp, &junk)+bestFitVal-
                            func->EvalAtEx(pi, &junk)-
                            func->EvalAtEx(pj, &junk))/(di*dj);
        ret(i, j) = dij;
        ret(j, i) = dij;
      }
    }

    ret.Invert();
    return ret;
  }

  bool Interpolate(double val1, double val2, double target, double& ret)
  {
    if( (val1-target)*(val2-target) > 0) return false;
    const double diff = (val2-val1) ? val2-val1 : 1;
    ret = (target-val1)/diff;
    return true;
  }

  // ---------------------------------------------------------------------
  std::vector<CoordNDim> FindContourInt(const ICallableInt* func_org,
                                        const CoordInt& minCoord,
                                        const double height,
                                        std::vector<TGraph>* marg_graphs)
  {
    assert(!marg_graphs || marg_graphs->empty());
    assert(minCoord.size() == 2);

    // Put a cache in front of ourselves
    FuncCacheInt* func = new FuncCacheInt(func_org);

    std::map<CoordInt, bool> results; // Is this box already in the result set?
    std::vector<CoordInt> todo;       // Boxes to investigate

    MSG("NCFitter", Msg::kVerbose) << "Contour: walking left" << endl;
    // Walk to the left
    for(int x = minCoord.x(); ; --x){
      // As soon as we exceed the contour
      MSG("NCFitter", Msg::kVerbose) << "Contour: at " << x << endl;
      CoordInt evalAt(x, minCoord.y());
      CoordInt bestfit;
      if(func->EvalAtEx(evalAt, &bestfit) > height){
        // Create enough graphs to save our marginalisation progress
        //(2 extra to save the two variables actually in the contour)
        if(marg_graphs)
          for(unsigned int n = 0; n < 2+bestfit.size(); ++n)
            marg_graphs->push_back(TGraph());

        // make it the first element in the todo list
        todo.push_back(evalAt);
        // and stop going left.
        break;
      }
    }

    std::vector<CoordNDim> contour;

    TStopwatch sw; // starts automatically

    while(!todo.empty()){
      CoordInt now = todo.back();
      todo.pop_back();

      MSG("NCFitter", Msg::kVerbose) << "Contour: at " << now << endl;

      const double est_frac = (M_PI+
                               atan2(-double(now.y()-minCoord.y()-.01),
                                     double(now.x()-minCoord.x())))/(2*M_PI);
      ReportProgress(est_frac, sw);


      // For each side of the box: see if one corner is high and the other low.
      // If so, and if the box adjacent on this side isn't already part of the
      // results, then append the adjacent box to the todo list and linearly
      // interpolate to find where to put the contour point on that edge.

      for(int dx = -1; dx <= 1; ++dx){
        for(int dy = -1; dy <= 1; ++dy){
          // no diagonals, don't stay stationary
          if(abs(dx+dy) != 1) continue;

          // where we're trying to go
          const CoordInt newpos(now.x()+dx, now.y()+dy);

          // offset to the first corner of the box we're considering
          const int cx = int(dx > 0);
          const int cy = int(dy > 0);

          // vector to the second corner we're considering
          const int ix = abs(dy);
          const int iy = abs(dx);

          // The values at said corners. Reevaluating them is fine,
          // because 'func' is cached.
          CoordInt bf1, bf2;
          const double h1 = func->EvalAtEx(CoordInt(now.x()+cx,
                                                    now.y()+cy), &bf1);
          const double h2 = func->EvalAtEx(CoordInt(now.x()+cx+ix,
                                                    now.y()+cy+iy), &bf2);

          double interp;
          if(!results[newpos] && Interpolate(h1, h2, height, interp)){
            todo.push_back(newpos);
            results[newpos] = true;
            contour.push_back(CoordNDim(now.x()+cx+ix*interp,
                                        now.y()+cy+iy*interp));

            if(marg_graphs){
              marg_graphs->at(0).SetPoint(marg_graphs->at(0).GetN(),
                                          marg_graphs->at(0).GetN(), now.x());
              marg_graphs->at(1).SetPoint(marg_graphs->at(1).GetN(),
                                          marg_graphs->at(1).GetN(), now.y());
              for(unsigned int n = 0; n < bf1.size(); ++n){
                const double bf = bf1[n]+interp*(bf2[n]-bf1[n]);
                marg_graphs->at(2+n).SetPoint(marg_graphs->at(2+n).GetN(),
                                              marg_graphs->at(2+n).GetN(), bf);
              }
            } // end if marg_graphs
          }
        } // end for dy
      } // end for dx
    } // end while todo

    if(!contour.empty()) contour.push_back(contour[0]); // Close the loop.

    delete func; // our cache

    return contour;
  }


  // ---------------------------------------------------------------------
  TH2D* GetLowResGrid(const ICallable2D* func,
                      const Coord min,
                      const Coord max,
                      const Coord prec)
  {
    const int Nx = int((max.x-min.x)/prec.x+.5);
    const int Ny = int((max.y-min.y)/prec.y+.5);

    TH2D* ret = new TH2D("foo", "bar", Nx, min.x, max.x, Ny, min.y, max.y);

    for(int nx = 0; nx < Nx ; ++nx){
      for(int ny = 0; ny < Ny; ++ny){
        const double x = min.x+(max.x-min.x)/double(Nx)*(nx+.5);
        const double y = min.y+(max.y-min.y)/double(Ny)*(ny+.5);
        CoordNDim bestfit;
        const double f = func->EvalAtEx(Coord(x, y), bestfit);
        MSG("NCFitter", Msg::kVerbose) << "GetLowResGrid: "
                                       << nx << "," << ny
                                       << " (" << x << "," << y << ") "
                                       << "-> " << f << endl;
        ret->Fill(x, y, f);
      }
    }
    return ret;
  }


  // ---------------------------------------------------------------------
  TH2D* GridSearchInt(const ICallableInt* func,
                      const CoordInt min,
                      const CoordInt max)
  {
    assert(min.size() == 2 && max.size() == 2);

    TStopwatch sw;

    TH2D* ret = new TH2D("", "", max.x()-min.x(), min.x(), max.x(),
                                 max.y()-min.y(), min.y(), max.y());

    for(int ny = min.y(); ny < max.y(); ++ny){
      for(int nx = min.x(); nx < max.x() ; ++nx){
        const double f = func->EvalAt(CoordInt(nx, ny));
        ret->Fill(nx, ny, f);
      }
      ReportProgress(double(ny)/(max.y()-min.y()), sw);
    }

    return ret;
  }


  // ---------------------------------------------------------------------
  TH2D* ConvertGrid(const TH2D* intgrid,
                    const NCParameter& xp,
                    const NCParameter& yp)
  {
    const int nx = intgrid->GetNbinsX();
    const int ny = intgrid->GetNbinsY();

    double* xb = new double[nx+1];
    double* yb = new double[ny+1];

    for(int x = 0; x <= nx; ++x) xb[x] = xp.ValFromStep(x);
    for(int y = 0; y <= ny; ++y) yb[y] = yp.ValFromStep(y);

    TH2D* ret = new TH2D("", TString::Format(";%s;%s;",
                                             xp.LatexName().c_str(),
                                             yp.LatexName().c_str()).Data(),
                         nx, xb, ny, yb);

    delete[] xb;
    delete[] yb;

    for(int x = 0; x <= nx; ++x)
      for(int y = 0; y <= ny; ++y)
        ret->SetBinContent(x, y, intgrid->GetBinContent(x, y));

    return ret;
  }


  // ---------------------------------------------------------------------
  TGraph* CreateOneDimProjection(const ICallable2D* func,
                                 const NCParameter param,
                                 const double minVal,
                                 double& leftErr,
                                 double& rightErr,
                                 std::vector<TGraph>& marg_graphs)
  {
    // +.5 makes us round to nearest int instead of toward zero.
    assert(param.Precision());
    const int steps = int(.5+(param.Max()-param.Min())/param.Precision());
    assert(steps);

    vector<double> x, y;

    TStopwatch sw; // starts automatically

    double minproj = 1e10;

    for(int nx = 0; nx <= steps; ++nx){
      const double xhere = param.Min()+nx*param.Precision();

      ReportProgress(double(nx)/steps, sw);

      CoordNDim bestfit;
      double upVal = func->EvalAtEx(Coord(xhere, -999), bestfit);
      x.push_back(xhere);
      y.push_back(upVal);

      std::stringstream op;
      op << "CreateOneDimProjection: "
         << param.ShortName() << " = " << xhere << " chisq = " << upVal
         << " other variables: ";
      for(unsigned int n = 0; n < bestfit.size(); ++n)
            op << " " << bestfit[n];
      op << endl;
      MSG("NCFitter", Msg::kVerbose) << op.str();


      if(upVal < minproj) minproj = upVal;

      for(unsigned int n = 0; n < bestfit.size(); ++n){
        // If it's our first time through we need to create the graphs too.
        //      if(nx == 0) marg_graphs.push_back(TGraph());
        // TODO Something is broken with the marg graphs, but I'm too lazy
        // to fix it, so for now just do this.
        if(marg_graphs.size() == n) marg_graphs.push_back(TGraph());
        marg_graphs[n].Set(marg_graphs[n].GetN()+1);
        marg_graphs[n].SetPoint(marg_graphs[n].GetN()-1, xhere, bestfit[n]);
      }
    }

    MSG("NCFitter", Msg::kInfo) << "CreateOneDimProjection: "
                                << "min chisq = " << minproj
                                << " cf minimizer's " << minVal << endl;

#if 0
    cout << "Beginning backtrack test...\n";
    for(int nx = steps; nx >= 0; --nx){
      const double xhere = param.Min()+nx*param.Precision();
      CoordNDim bestfit;
      double upVal = func->EvalAt(Coord(xhere, -999), bestfit);
      //      x.push_back(xhere);
      //      y.push_back(upVal);
      const double diff = upVal - y[nx];
      if(diff < -1e-4){
        cout << "WARNING - GOT BETTER RESULT TRACKING BACKWARD, FIXING\n";
        y[nx] = upVal;
      }

      if(upVal < minproj) minproj = upVal;
      /*
      for(unsigned int n = 0; n < bestfit.size(); ++n){
        // If it's our first time through we need to create the graphs too.
        if(nx == 0) marg_graphs.push_back(TGraph());
        marg_graphs[n].Set(marg_graphs[n].GetN()+1);
        marg_graphs[n].SetPoint(marg_graphs[n].GetN()-1, xhere, bestfit[n]);
      }
      */
    }
#endif

    bool needLeftErr = true;
    bool needRightErr = true;

    // Initialize right error for downward sloping projections.
    // .5 is to be consistent with binning elsewhere.
    rightErr = param.Max()+.5*param.Precision();

    for(int nx = 0; nx <= steps; ++nx){
      y[nx] -= minproj;

      const double upVal = y[nx];

      if(needRightErr && !needLeftErr && upVal > 1){
        needRightErr = false;
        rightErr = x[nx];
      }
      if(needLeftErr && upVal < 1){
        needLeftErr = false;
        leftErr = x[nx]-param.Precision();
      }
    }

    TGraph* ret = new TGraph(steps+1, &x[0], &y[0]);

    return ret;
  }


  // ---------------------------------------------------------------------
  TGraph CreateOneDimProjectionInt(const ICallableInt* func,
                                   int minstep,
                                   int maxstep,
                                   double& leftErr,
                                   double& rightErr)
  {
    vector<double> xs, ys;

    TStopwatch sw; // starts automatically

    double minproj = 1e10;

    for(int x = minstep; x <= maxstep; ++x){
      ReportProgress(double(x-minstep)/(maxstep-minstep), sw);

      //      CoordInt bestfit;
      CoordInt pos; pos.push_back(x);
      const double upVal = func->EvalAt(pos);//, bestfit);
      xs.push_back(x);
      ys.push_back(upVal);

      if(upVal < minproj) minproj = upVal;
      /*
      for(unsigned int n = 0; n < bestfit.size(); ++n){
        // If it's our first time through we need to create the graphs too.
        if(nx == 0) marg_graphs.push_back(TGraph());
        marg_graphs[n].Set(marg_graphs[n].GetN()+1);
        marg_graphs[n].SetPoint(marg_graphs[n].GetN()-1, xhere, bestfit[n]);
      }
      */
    }

    MSG("NCFitter", Msg::kInfo) << "CreateOneDimProjection: "
                                << "min chisq = " << minproj << endl;


    bool needLeftErr = true;
    bool needRightErr = true;

    // Find the intersections using a linear interpolation
    for(int x = 0; x <= maxstep-minstep; ++x){
      ys[x] -= minproj;

      const double upVal = ys[x];

      if(needRightErr && !needLeftErr && upVal > 1){
        needRightErr = false;
        //      rightErr = minstep+x;

        rightErr = (minstep+x-1)+(1-ys[x-1])/(ys[x]-ys[x-1]);
      }
      if(needLeftErr && upVal < 1){
        needLeftErr = false;
        //      leftErr = minstep+x;

        if(x == 0)
          leftErr = minstep;
        else
          leftErr = (minstep+x-1)+(1-ys[x-1])/(ys[x]-ys[x-1]);
      }
    }

    if(needRightErr) rightErr = maxstep;

    TGraph ret(maxstep-minstep, &xs[0], &ys[0]);
    return ret;
  }

  // ---------------------------------------------------------------------
  vector<TGraph> MakeSlices(const ICallableND* func,
                            const vector<NCParameter>& ps,
                            const CoordNDim& minpos)
  {
    vector<TGraph> ret;
    for(unsigned int p = 0; p < ps.size(); ++p){
      MSG("NCFitter", Msg::kInfo) << "Taking slice in "
                                  << ps[p].ShortName()
                                  << "..." << endl;

      // +.5 makes us round to nearest int instead of toward zero.
      assert(ps[p].Precision());
      const int steps = int(.5+(ps[p].Max()-ps[p].Min())/ps[p].Precision());
      assert(steps);

      vector<double> x, y;

      TStopwatch sw; // starts automatically

      for(int nx = 0; nx <= steps; ++nx){
        const double xhere = ps[p].Min()+nx*ps[p].Precision();

        ReportProgress(double(nx)/steps, sw);

        CoordNDim evalAt = minpos;
        evalAt[p] = xhere;
        CoordNDim junk;
        double upVal = func->EvalAtEx(evalAt, &junk);
        x.push_back(xhere);
        y.push_back(upVal);
      } // end for nx

      MSG("NCFitter", Msg::kInfo) << endl;

      TGraph g = TGraph(steps+1, &x[0], &y[0]);
      g.SetName(TString::Format("slice_%s", ps[p].ShortName().c_str()));
      g.SetTitle(TString::Format("Slice;%s;#chi^{2}",
                                 ps[p].LatexName().c_str()));
      ret.push_back(g);
    }

    return ret;
  }

  // ---------------------------------------------------------------------
  vector<TH2D*> Make2DSlices(const ICallableND* func,
                             const vector<NCParameter>& ps,
                             const CoordNDim& minpos)
  {
    vector<TH2D*> ret;
    for(unsigned int p = 0; p < ps.size(); ++p){
      for(unsigned int q = 0; q < p; ++q){
        // +.5 makes us round to nearest int instead of toward zero.
        assert(ps[p].Precision());
        assert(ps[q].Precision());

        NCParameter px=ps[p];
        NCParameter py=ps[q];

        const int Nx = int((px.Max()-px.Min())/px.Precision()+.5);
        const int Ny = int((py.Max()-py.Min())/py.Precision()+.5);

        TH2D* h = new TH2D("foo", "bar", Nx, px.Min(), px.Max(),
                           Ny, py.Min(), py.Max());
        h->SetName(TString::Format("2dslice_%s_%s", px.ShortName().c_str(),
                                   py.ShortName().c_str() ));

        h->GetXaxis()->SetTitle(px.LatexName().c_str());
        h->GetYaxis()->SetTitle(py.LatexName().c_str());

        TStopwatch sw; // starts automatically

        for(int nx = 1; nx <= Nx; ++nx){
          ReportProgress(double(nx)/Nx, sw);
          const double xhere = px.Min()+nx*px.Precision();
          for(int ny = 1; ny <= Ny; ++ny){
            const double yhere = py.Min()+ny*py.Precision();
            CoordNDim evalAt = minpos;
            evalAt[p] = xhere; evalAt[q] = yhere;
            CoordNDim junk;
            h->SetBinContent(nx, ny, func->EvalAtEx(evalAt, &junk));
          } // end for ny
        } // end for nx

        ret.push_back(h);
      }
    }

    return ret;
  }


  // ---------------------------------------------------------------------
  FixVars::FixVars(const ICallableND* fun,
                   const CoordNDim& fix,
                   const std::vector<int>& vary)
    :fFunc(fun), fFixed(fix), fToVary(vary){}


  double FixVars::EvalAtEx(const CoordNDim& r, CoordNDim* ret) const
  {
    assert(fToVary.size() == r.size());
    // For every parameter in the list of what we are to vary - substitute
    // in, in order, from the parameters we were passed.
    CoordNDim fixed = fFixed;
    for(unsigned int n = 0; n < fToVary.size(); ++n)
      fixed[fToVary[n]] = r[n];
    return fFunc->EvalAtEx(fixed, ret);
  }


  // ---------------------------------------------------------------------
  FixVarsInt::FixVarsInt(const ICallableInt* fun,
                         const CoordInt& fix,
                         const std::vector<int>& vary)
    :fFunc(fun), fFixed(fix), fToVary(vary){}


  double FixVarsInt::EvalAt(const CoordInt& r) const
  {
    assert(fToVary.size() == r.size());
    // For every parameter in the list of what we are to vary - substitute
    // in, in order, from the parameters we were passed.
    CoordInt fixed = fFixed;
    for(unsigned int n = 0; n < fToVary.size(); ++n)
      fixed[fToVary[n]] = r[n];
    return fFunc->EvalAt(fixed);
  }


  // ---------------------------------------------------------------------
  double ConstrainRange::EvalAtEx(const CoordNDim& r, CoordNDim* ret) const
  {
    CoordNDim r2 = r;
    double pen = 0;
    if(r2[fParamNum] > fMax){
      // TODO - there should be some scaling based on the width,
      // but we happen to know all these variables are approx 0-1
      pen += 1e3*(r2[fParamNum]-fMax)*(r2[fParamNum]-fMax);
      r2[fParamNum] = fMax;
    }
    if(r2[fParamNum] < fMin){
      pen += 1e3*(r2[fParamNum]-fMin)*(r2[fParamNum]-fMin);
      r2[fParamNum] = fMin;
    }
    return fFunc->EvalAtEx(r2, ret)+pen;
  }


  // ---------------------------------------------------------------------
  MarginalizeSimple::MarginalizeSimple(const ICallableND* f,
                                       const std::vector<NCParameter>& p,
                                       const bool reuseOld,
                                       const int xcoord,
                                       const int ycoord)
    :fParams(p), fXpos(xcoord), fYpos(ycoord)
  {
    assert(fXpos < int(fParams.size()) && fYpos < int(fParams.size()));

    vector<NCParameter> shortParams;

    for(int n = 0; n < int(fParams.size()); ++n){
      if(n != fXpos && n != fYpos){
        fVary.push_back(n);
        shortParams.push_back(fParams[n]);
      }
    }

    CoordNDim fix;
    fix.resize(fParams.size());
    f2 = new FixVars(f, fix, fVary);

    fMinFinder = new MinFinderSimple(shortParams, true, reuseOld, f2);
  }


  // ---------------------------------------------------------------------
  double MarginalizeSimple::EvalAtEx(const Coord& r, CoordNDim& ret) const
  {
    CoordNDim fix;
    fix.resize(fParams.size());
    fix[fXpos] = r.x;
    if(fYpos > -1) fix[fYpos] = r.y;
    f2->SetFixed(fix);

    return fMinFinder->FindMin(ret);
  }


  // ---------------------------------------------------------------------
  MarginalizeSimpleInt::MarginalizeSimpleInt(const ICallableInt* f,
                                             const std::vector<NCParameter>& p,
                                             const int xcoord,
                                             const int ycoord)
    :fParams(p), fXpos(xcoord), fYpos(ycoord)
  {
    assert(fXpos < int(fParams.size()) && fYpos < int(fParams.size()));

    for(int n = 0; n < int(fParams.size()); ++n){
      if(n != fXpos && n != fYpos){
        fVary.push_back(n);
        fShortParams.push_back(fParams[n]);
      }
    }

    CoordInt fix;
    fix.resize(fParams.size());
    f2 = new FixVarsInt(f, fix, fVary);

    minFinder = new MinFinderSimpleInt(fShortParams, true, f2);
  }


  // ---------------------------------------------------------------------
  double MarginalizeSimpleInt::EvalAtEx(const CoordInt& r, CoordInt* ret) const
  {
    CoordInt fix;
    fix.resize(fParams.size());
    fix[fXpos] = r.x();
    if(fYpos > -1) fix[fYpos] = r.y();
    f2->SetFixed(fix);

    return minFinder->FindMin(ret);
  }


#ifndef STANDALONE
  Int_t MinuitForCallable::Eval(Int_t /*npar*/,
                                Double_t* /*grad*/,
                                Double_t& fval,
                                Double_t* par,
                                Int_t flag)
  {
    // If we are being asked to do something that isn't just
    // calculate the function then we don't know how...
    assert(flag == 4);
    CoordNDim r(fNumPar);
    for(int n = 0; n < fNumPar; ++n) r[n] = par[n];
    CoordNDim junk;
    fval = fFunc->EvalAtEx(r, &junk);
    return 0;
  }


  // ---------------------------------------------------------------------
  MarginalizeHybrid::MarginalizeHybrid(const ICallableND* f,
                                       const std::vector<NCParameter> p,
                                       const bool testmin,
                                       const int xcoord,
                                       const int ycoord)
    :fFunc(f), fParams(p), fTestMinimum(testmin), fXpos(xcoord), fYpos(ycoord)
  {
    assert(fXpos < int(fParams.size()) && fYpos < int(fParams.size()));
    assert(fYpos > fXpos || fYpos == -1);
    // Recycling the old Minuit instance each time is a possible cause of
    // getting stuck in the wrong minimum.
    /*
    pminu = new MinuitForCallable(fFunc, fParams.size());

    pminu->SetPrintLevel(-1);
    //  minu.Command("set strategy 2"); // Be careful, instead of fast.

    for(int n = 0; n < int(fParams.size()); ++n){
      pminu->DefineParameter(n, fParams[n].ShortName().c_str(),
                             (fParams[n].min+fParams[n].Max())/2,
                             (fParams[n].Max()-fParams[n].Min())/4, 0, 0);
    }
    */
    fLastMinPos.resize(fParams.size());
  }


  // ---------------------------------------------------------------------
  double MarginalizeHybrid::EvalAtEx(const Coord& r, CoordNDim& ret) const
  {
    //assert(pminu);

    double valAtMin;
    //    if(fTestMinimum && AlreadyAtMinimum(r, valAtMin)) return valAtMin;

    MinuitForCallable* pminu = new MinuitForCallable(fFunc, fParams.size());

    pminu->SetPrintLevel(-1);
    //    pminu->Command("set strategy 2"); // Be careful, instead of fast.

    for(int n = 0; n < int(fParams.size()); ++n){
      pminu->DefineParameter(n, fParams[n].ShortName().c_str(),
                             (fParams[n].Min()+fParams[n].Max())/2,
                             (fParams[n].Max()-fParams[n].Min())/4, 0, 0);
    }

    // Minuit tries to return parameters we don't want to keep.
    double junkd; int junki;

    // Add one to the index because of ridiculous fortran numbering scheme.
    double args1[] = {fXpos+1, r.x};
    pminu->mnexcm("set param", args1, 2, junki);
    if(fYpos > 0){
      double args2[] = {fYpos+1, r.y};
      pminu->mnexcm("set param", args2, 2, junki);
    }

    pminu->mnfixp(fXpos, junki);
    if(fYpos > 0){
      // The parameters all got renumbered when we fixed x.
      pminu->mnfixp(fYpos-1, junki);
    }

    // Simplex seems to be much faster than Migrad.
    // I think Migrad can fail when put in a bad place, and spends
    // a long time doing it.
    // However - I have had bad results using Simplex, so stick with Migrad.
    //  minu.mnsimp();
    if(pminu->Migrad()){
      //return FLT_MAX; // this messes up the contours...
      MSG("NCFitter", Msg::kVerbose) << "MinWRT: Migrad didn't converge"
                                     << endl;
    }

    pminu->mnstat(valAtMin, junkd, junkd, junki, junki, junki);
    for(unsigned int n = 0; n < fParams.size(); ++n){
      pminu->GetParameter(n, fLastMinPos[n], junkd);
      if(int(n) != fXpos && int(n) != fYpos)
            ret.push_back(fLastMinPos[n]);
    }

    TString txt = "MarginalizeHybrid: min at ";
    for(unsigned int n = 0; n < fParams.size(); ++n){
      txt += fParams[n].ShortName() + "=";
      txt += fLastMinPos[n];
      txt += +"\t";
    }
    MSG("NCFitter", Msg::kVerbose) << txt << " chisq = "
                                   << valAtMin << endl;

    // Free all the parameters, we will refix them next time round.
    pminu->mnfree(0);

    delete pminu;

    return valAtMin;
  }


  // Is the function higher one step size in each of the minimized-over
  // directions. If so, return true, and place the value of the old function
  // at the old position in these coordinates in valAtMin.
  bool MarginalizeHybrid::AlreadyAtMinimum(const Coord& r, double& valAtMin) const
  {
    CoordNDim p = fLastMinPos;   // Last place we were
    // Move to where we are now on one axis of the contour
    p[fXpos] = r.x;
    if(fYpos > 0) p[fYpos] = r.y; // If we have a second axis move there too.

    CoordNDim junk;
    valAtMin = fFunc->EvalAtEx(p, &junk);

    // What direction to take a step in.
    for(int n = 0; n < int(fParams.size()); ++n){
      // Not in the contour axis directions.
      if(n == fXpos || n == fYpos) continue;
      CoordNDim q = p;
      q[n] += fParams[n].Precision();   // Move one step positive in this axis
      if(fFunc->EvalAtEx(q, &junk) < valAtMin) return false;
      q[n] -= 2*fParams[n].Precision(); // Also try one step negative.
      if(fFunc->EvalAtEx(q, &junk) < valAtMin) return false;
    }
    MSG("NCFitter", Msg::kVerbose) << "MarginalizeHybrid: "
                                   << "Already at minimum" << endl;
    return true;
  }

#endif

  // ---------------------------------------------------------------------
  void PrintGridSearch(const ICallableND* func,
                       const std::vector<NCParameter>& params,
                             std::ostream& out)
  {
    std::vector<NCParameter> lowResParams = params;
    for(unsigned int n = 0; n < lowResParams.size(); ++n)
      lowResParams[n].SetPrecision(lowResParams[n].Precision()*10);

    std::vector<CoordNDim> evalAt;

    CoordNDim first;
    for(unsigned int n = 0; n < lowResParams.size(); ++n)
      first.push_back(lowResParams[n].Min());

    evalAt.push_back(first);

    for(unsigned int d = 0; d < lowResParams.size(); ++d){
      unsigned int maxn = evalAt.size();
      for(unsigned int n = 0; n < maxn; ++n){
        for(unsigned int nx = 1; ; ++nx){
          CoordNDim base = evalAt[n];
          base[d] += lowResParams[d].Precision()*nx;
          if(base[d] > lowResParams[d].Max()) break;
          evalAt.push_back(base);
        }
      }
    }

    for(unsigned int d = 0; d < lowResParams.size(); ++d)
      out << lowResParams[d].ShortName() << "\t";
    out << "chisq" << endl;

    for(unsigned int n = 0; n < evalAt.size(); ++n){
      for(unsigned int d = 0; d < evalAt[n].size(); ++d){
        out << evalAt[n][d]<<"\t";
      }
      CoordNDim junk;
      out << func->EvalAtEx(evalAt[n], &junk) << endl;
    }
  }


}} // end namespaces

---