NCContourFinder Namespace Reference

Classes
class	NCContourFinder::Coord
class	NCContourFinder::ICallableND
class	NCContourFinder::ICallable2D
class	NCContourFinder::FindMinSimple
class	NCContourFinder::FixVars
class	NCContourFinder::MarginalizeSimple
class	NCContourFinder::MinuitForCallable
class	NCContourFinder::MarginalizeHybrid
Typedefs
typedef std::vector< double >	CoordNDim
Functions
void	ReportProgress (double est_frac, TStopwatch &sw)
double	FindMinMinuit (const ICallableND *chiSqrFunc, const std::vector< NCParameter > params, CoordNDim &ret, TMinuit **ppMinu)
std::vector< Coord >	FindContour (const ICallable2D *func, const Coord &minCoord, const double &height, const Coord &precision, std::vector< TGraph > &marg_graphs)
TMatrix	FindCovarianceMatrix (const ICallableND *func, const CoordNDim &minCoord, const std::vector< NCParameter > params)
TH2D *	GetLowResGrid (const ICallable2D *func, const Coord min, const Coord max, const Coord prec)
TGraph *	CreateOneDimProjection (const ICallable2D *func, const NCParameter param, const double minVal, double &leftErr, double &rightErr, std::vector< TGraph > &marg_graphs)
void	PrintGridSearch (const ICallableND *func, const std::vector< NCParameter > &params, std::ostream &out)

---

Typedef Documentation

typedef std::vector<double> NCContourFinder::CoordNDim

Definition at line 26 of file NCContourFinder.h. Referenced by NCContourFinder::MarginalizeHybrid::AlreadyAtMinimum(), NCContourFinder::FindMinSimple::BeginSearchPos(), CreateOneDimProjection(), NCContourFinder::MinuitForCallable::Eval(), NCExtrapolation::FillResultHistograms(), FindContour(), FindCovarianceMatrix(), NCExtrapolation::GetBestFitPoint(), NCExtrapolationJK_fit::GetChiSqrForMap(), GetLowResGrid(), NCContourFinder::MarginalizeSimple::MarginalizeSimple(), NCExtrapolation::GetChiSqrFromDerived::operator()(), NCContourFinder::MarginalizeSimple::operator()(), NCContourFinder::FindMinSimple::operator()(), NCExtrapolationJK_fit::PredictSpectrum(), and PrintGridSearch().

---

Function Documentation

TGraph * NCContourFinder::CreateOneDimProjection (  const ICallable2D *  func, const NCParameter  param, const double  minVal, double &  leftErr, double &  rightErr, std::vector< TGraph > &  marg_graphs )

Definition at line 436 of file NCContourFinder.cxx. References CoordNDim, NCParameter::max, NCParameter::min, MSG, NCParameter::prec, ReportProgress(), and NCParameter::shortName. Referenced by NCExtrapolation::Get1DProjection(). { // +.5 makes us round to nearest int instead of toward zero. assert(param.prec); const int steps = int(.5+(param.max-param.min)/param.prec); 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.prec; ReportProgress(double(nx)/steps, sw); CoordNDim bestfit; double upVal = (*func)(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(uint n = 0; n < bestfit.size(); ++n) op << " " << bestfit[n]; op << endl; MSG("NCContourFinder", Msg::kVerbose) << op.str(); if(upVal < minproj) minproj = upVal; for(uint 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("NCContourFinder", Msg::kInfo) << "CreateOneDimProjection: min chisq = " << minproj << " cf minimizer's " << minVal << endl; 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.prec; 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.prec; } } TGraph* ret = new TGraph(steps+1, &x[0], &y[0]); return ret; }

std::vector< Coord > NCContourFinder::FindContour (  const ICallable2D *  func, const Coord &  minCoord, const double &  height, const Coord &  precision, std::vector< TGraph > &  marg_graphs )

Definition at line 203 of file NCContourFinder.cxx. References atan2(), CoordNDim, MSG, ReportProgress(), NCContourFinder::Coord::x, and NCContourFinder::Coord::y. Referenced by NCExtrapolation::GetContour(). { assert(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 std::map<Box, double> heights; // Boxes whose topleft corner has already been evaluated, used as cache. // 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("NCContourFinder", Msg::kVerbose) << "Contour: walking left" << endl; // Walk to the left for(int x = bfxi; ; --x){ // As soon as we exceed the contour MSG("NCContourFinder", Msg::kVerbose) << "Contour: at " << x << endl; CoordNDim bestfit; if((*func)(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) for(uint n = 0; n < 4*(2+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("NCContourFinder", 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)(Coord((now.first+cx)*dx, (now.second+cy)*dy), bestfit); const int dn = cx+2*cy; marg_graphs[dn].Set(marg_graphs[dn].GetN()+1); marg_graphs[dn].SetPoint(marg_graphs[dn].GetN()-1, gxi, (now.first+cx)*dx); marg_graphs[4+dn].Set(marg_graphs[4+dn].GetN()+1); marg_graphs[4+dn].SetPoint(marg_graphs[4+dn].GetN()-1, gxi, (now.second+cy)*dy); for(uint n = 0; n < bestfit.size(); ++n){ marg_graphs[8+4*n+dn].Set(marg_graphs[8+4*n+dn].GetN()+1); marg_graphs[8+4*n+dn].SetPoint(marg_graphs[8+4*n+dn].GetN()-1, 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("NCContourFinder", 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, 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("NCContourFinder", 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("NCContourFinder", 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("NCContourFinder", 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("NCContourFinder", 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)); } } } if(contour.size()) contour.push_back(contour[0]); // Close the loop. return contour; }

TMatrix NCContourFinder::FindCovarianceMatrix (  const ICallableND *  func, const CoordNDim &  minCoord, const std::vector< NCParameter >  params )

Definition at line 361 of file NCContourFinder.cxx. References CoordNDim. Referenced by NCExtrapolation::MakeDeltaChiSqrMaps(). { assert(params.size() == minCoord.size()); TMatrix ret(params.size(), params.size()); const double bestFitVal = (*func)(minCoord); // First, do the 2nd derivatives in one variable for(uint i = 0; i < params.size(); ++i){ const double step = params[i].prec; CoordNDim plus = minCoord; plus[i] += step; CoordNDim minus = minCoord; minus[i] -= step; const double dii = ((*func)(plus)+(*func)(minus)-2*bestFitVal)/(step*step); ret(i, i) = dii; } // Now do all the off-diagonal derivatives for(uint i = 0; i < params.size()-1; ++i){ for(uint j = i+1; j < params.size(); ++j){ const double di = params[i].prec; const double dj = params[j].prec; 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)(pp)+bestFitVal-(*func)(pi)-(*func)(pj))/(di*dj); ret(i, j) = dij; ret(j, i) = dij; } } ret.Invert(); return ret; }

double NCContourFinder::FindMinMinuit (  const ICallableND *  chiSqrFunc, const std::vector< NCParameter >  params, CoordNDim &  ret, TMinuit **  minu = 0 )

Definition at line 155 of file NCContourFinder.cxx. References MSG. Referenced by NCExtrapolation::MakeDeltaChiSqrMaps(). { 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("NCContourFinder", Msg::kWarning) << "Migrad failed, falling back to Simplex.\n"; minu->SetPrintLevel(0); // Things are going wrong, we should let the user see. 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; }

TH2D * NCContourFinder::GetLowResGrid (  const ICallable2D *  func, const Coord  min, const Coord  max, const Coord  prec )

Definition at line 407 of file NCContourFinder.cxx. References CoordNDim, max, min, MSG, NCContourFinder::Coord::x, and NCContourFinder::Coord::y. Referenced by NCExtrapolation::Get2DGrid(). { 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)(Coord(x, y), bestfit); MSG("NCContourFinder", Msg::kVerbose) << "GetLowResGrid: " << nx << "," << ny << " (" << x << "," << y << ") " << "-> " << f << endl; ret->Fill(x, y, f); } } return ret; }

void NCContourFinder::PrintGridSearch (  const ICallableND *  func, const std::vector< NCParameter > &  params, std::ostream &  out )

Definition at line 714 of file NCContourFinder.cxx. References base, CoordNDim, and min. Referenced by NCExtrapolation::MakeDeltaChiSqrMaps(). { std::vector<NCParameter> lowResParams = params; for(unsigned int n = 0; n < lowResParams.size(); ++n) lowResParams[n].prec *= 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].prec*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"; } out << (*func)(evalAt[n])<<endl; } }

void ReportProgress (  double  est_frac, TStopwatch &  sw )  [static]

Private utility function. Prints estimated time remaining Definition at line 26 of file NCContourFinder.cxx. References MSG. Referenced by CreateOneDimProjection(), and FindContour(). { const double elapsed = sw.RealTime(); sw.Start(kFALSE); const int remain = est_frac ? int(elapsed*(1/est_frac-1)) : 0; MSG("NCContourFinder", Msg::kInfo) << "Estimated " << int(100*est_frac) << "% complete - " << remain/3600 << "h" << (remain/60)%60 << "m" << remain%60 << "s remain.\n"; }

---