Step Namespace Reference

Namespaces
	Basis

Classes
struct	Correlator
struct	FunctionalForm
struct	Polynomial
struct	Result

Enumerations
enum	EarlyStopping { None, Chi2ndf, fTest, xVal }

Functions
std::ostream &	operator<< (std::ostream &stream, const Result &result)
double	Ftest (const Result &former, const Result &current)
double	Xval (const Result &current, const Result &former)
void	PrintAllTests (std::function< double(const Result &, const Result &)> test, std::ostream &stream=std::cout, size_t best=0)
double	identity (double x)
template<double(*)(double, int) T = Basis::Chebyshev, double(*)(double) Log = identity, double(*)(double) Exp = identity>
TF1 *	GetSmoothFit (TH1 *h, TH2 *cov, int im, int iM, unsigned int maxdegree, EarlyStopping criterion=None, double stopParam=1, UInt_t nrep=0, std::ostream &stream=std::cout)
template<double(*)(double, int) T = Basis::Chebyshev, double(*)(double) Log = identity, double(*)(double) Exp = identity>
TF1 *	GetSmoothFit (TH1 *h, TH2 *cov, unsigned int maxdegree, EarlyStopping criterion=EarlyStopping::None, double stopParam=1, UInt_t nrep=0, std::ostream &stream=std::cout)
template<double(*)(double, int) T = Basis::Chebyshev, double(*)(double) Log = identity, double(*)(double) Exp = identity>
TF1 *	GetSmoothFit (TH1 *h, int im, int iM, unsigned int maxdegree, EarlyStopping criterion=EarlyStopping::None, double stopParam=1, UInt_t nrep=0, std::ostream &stream=std::cout)
template<double(*)(double, int) T = Basis::Chebyshev, double(*)(double) Log = identity, double(*)(double) Exp = identity>
TF1 *	GetSmoothFit (TH1 *h, unsigned int maxdegree, EarlyStopping criterion=EarlyStopping::None, double stopParam=1, UInt_t nrep=0, std::ostream &stream=std::cout)

Variables
static const char *	green = "\x1B[32m"
static const char *	red = "\x1B[31m"
static const char *	bold = "\e[1m"
static const char *	def = "\e[0m"
static const auto	eps = std::numeric_limits<double>::epsilon()
static const auto	inf = std::numeric_limits<double>::infinity()
static bool	verbose = false
static std::deque< Result >	chi2s
static auto	bestResult = chi2s.begin()

Detailed Description

Smoothness Test with Expansion of Polynmials.

Enumeration Type Documentation

EarlyStopping

enum EarlyStopping

Enumerator
None	iterate until max value is reached
Chi2ndf	stop as soon as Chi2ndf has reached a plateau, or as soon as it is compatible with 1
fTest	include additional parameter as long as F-test provides a p-value beyond a certain threshold
xVal	include additional parameter as long as cross validation provides a p-value beyond a certain threshold

                   {
    None, 
    Chi2ndf, 
    fTest, 
    xVal 
};

Function Documentation

Ftest()

double Step::Ftest	(	const Result &	former,
		const Result &	current
	)

F-test

Fisher, R. A. (1922). "On the interpretation of chi2 from contingency tables, and the calculation of P". Journal of the Royal Statistical Society. 85 (1): 87–94. doi:10.2307/2340521. JSTOR 2340521.

Wikipedia

{
    using namespace std;
    if (!former.valid || !current.valid) throw -1.;
    double RSS1 = former.chi2,
           RSS2 = current.chi2;
    size_t p1 = former.X.size(),
           p2 = current.X.size();
    assert(p2 > p1);
    int ndf = current.ndf; // = nbins - p2
    double x = ((RSS1-RSS2)/(p2-p1))/(RSS2/ndf);
    double pvalue = TMath::FDistI(x, p2-p1, ndf);
    return pvalue;
}

GetSmoothFit() [1/4]

TF1* Step::GetSmoothFit	(	TH1 *	h,
		int	im,
		int	iM,
		unsigned int	maxdegree,
		EarlyStopping	criterion = EarlyStopping::None,
		double	stopParam = 1,
		UInt_t	nrep = 0,
		std::ostream &	stream = std::cout
	)

Function overloading, in case of no bin-to-bin correlations.

{
    return GetSmoothFit<T,Log,Exp>(h, nullptr, im, iM,
            maxdegree, criterion, stopParam, nrep, stream);
}

GetSmoothFit() [2/4]

TF1* Step::GetSmoothFit	(	TH1 *	h,
		TH2 *	cov,
		int	im,
		int	iM,
		unsigned int	maxdegree,
		EarlyStopping	criterion = None,
		double	stopParam = 1,
		UInt_t	nrep = 0,
		std::ostream &	stream = std::cout
	)

Step::GetSmoothFit() is the only function to call to get the smooth fit. It assumes that the histogram and the covariance matrix are already normalised to the bin width (i.e. cross section, not count).

Todo:

and in case of F-test failure?

Parameters

h	input histogram
cov	covariance matrix in histogram format
im	first bin to consider in the fit
iM	last fbin to consider in the fit
maxdegree	highest possible degree in fit
criterion	early stopping criterion
stopParam	early stopping parameter: if nrep > 0, threshold for x-validation; else #sigmas to satisfy
nrep	number of replicas
stream	stream (e.g. `cout` or file)

{
    using namespace std;
 
    assert(h != nullptr);
 
    if (verbose)
        stream << "im = " << im << ", iM = " << iM << '\n';
 
    unsigned int nbins = iM-im+1;
    if (verbose)
        stream << "nbins = " << nbins << '\n';
    assert(nbins > 2);
 
    double m = h->GetBinLowEdge(im),
           M = h->GetBinLowEdge(iM+1);
    if (verbose)
        stream << "m = " << m << ", M = " << M << '\n';
 
    Polynomial<T,Log,Exp> pol(maxdegree, m, M, stream);
    Correlator cor(h, cov, im, iM, pol, stream);
    ROOT::Math::Functor functor(cor, pol.npars);
 
    auto minimiser = ROOT::Math::Factory::CreateMinimizer("Minuit2", "Migrad");
    const int Ncalls = 1e6;
    minimiser->SetMaxFunctionCalls(Ncalls);
    minimiser->SetFunction(functor);
    minimiser->SetTolerance(0.001);
    minimiser->SetPrintLevel(verbose);
 
    double fm = Log(h->GetBinContent(im)),
           fM = Log(h->GetBinContent(iM));
    auto step = 0.001;
    double p0 = (fM+fm)/2, p1 = (fM-fm)/2;
    minimiser->SetVariable(0,"p0",p0,step);
    minimiser->SetVariable(1,"p1",p1,step);
    for (unsigned int d = 2; d <= maxdegree; ++d) {
        minimiser->SetVariable(d,Form("p%d", d),0.,step);
        minimiser->FixVariable(d);
    }
 
    // running
    chi2s.clear();
    chi2s.push_back({false, inf, nbins-2, {p0,p1}});
    bestResult = chi2s.begin();
    bool found = false;
    while (chi2s.size() <= maxdegree) {
        unsigned int degree = chi2s.size();
        minimiser->ReleaseVariable(degree);
        bool valid = minimiser->Minimize();
 
        const double * X = minimiser->X();
        assert(X != nullptr);
        double chi2 = cor(X);
 
        assert(isnormal(chi2));
        unsigned int npars = minimiser->NFree(),
                     ndf = nbins - npars;
        assert(npars == degree+1);
 
        const auto& former = chi2s.back();
        chi2s.push_back({valid, chi2, ndf, vector<double>(X,X+npars)});
        auto& current = chi2s.back();
 
        switch (minimiser->Status()) {
            case 0:  stream << current;
                     {
                         // check EarlyStopping::Chi2ndf criterion
                         bool plateau = current.x2n() > former.x2n(),
                              goodchi2 = abs(current.x2n()-1) <= stopParam*current.x2nErr();
                         if (plateau)
                             stream << red << "Chi2/ndf has increased.\n" << def;
                         else if (goodchi2)
                             stream << green << "Chi2 ~ ndf.\n" << def;
 
                         // check EarlyStopping::fTest criterion
                         double ftestRes = 0;
                         if (bestResult != chi2s.begin())
                         try {
                             ftestRes = Ftest(former, current);
                             stream << "F-test result: " << ftestRes << '\n';
                         }
                         catch (double) { stream << red << "No F-test possible.\n" << def; }
                         bool fvalid = ftestRes > stopParam; 
 
                         // check EarlyStopping::xVal criterion
                         double xvalFrac = 0;
                         if (nrep > 0) {
                             current.chi2sT.reserve(nrep);
                             current.chi2sV.reserve(nrep);
                             for (UInt_t seed = 1; seed <= nrep; ++seed) {
                                 cor.SetReplica(seed);
                                 for (unsigned int i = 0; i <= degree; ++i)
                                     minimiser->SetVariableValue(i, current.X[i]);
                                 const double * X = minimiser->X();
                                 assert(X != nullptr);
                                 current.chi2sT.push_back(cor(X));
                                 cor.SetReplica(seed+nrep);
                                 current.chi2sV.push_back(cor(X));
                             }
                             cor.UnsetReplica();
                             for (unsigned int i = 0; i <= degree; ++i)
                                 minimiser->SetVariableValue(i, current.X[i]);
                             X = minimiser->X();
                             xvalFrac = Xval(current, former);
                         }
                         bool xvalid = 1.-xvalFrac > stopParam;
                         if (nrep > 0)
                             stream << "Fraction of validation replicas with better fit Chi2/ndf = " << (xvalid ? green : red) << xvalFrac << def << '\n';
 
                         // apply early stopping criterion
                         if (!found)
                         switch (criterion) {
                             case EarlyStopping::Chi2ndf:
                                 found = plateau || goodchi2;
                                 if (plateau) --bestResult;
                                 break;
                             case EarlyStopping::fTest:
                                 found = !fvalid;
                                 break;
                             case EarlyStopping::xVal:
                                 assert(nrep > 0);
                                 found = !xvalid;
                                 break;
                             case EarlyStopping::None:
                                 found = false;
                         }
                     }
                     break;
            case 1:  stream << red << "Warning: Fit parameter covariance was made positive defined\n" << def;
                     for (unsigned int i = 0; i < npars; ++i) {
                         for (unsigned int j = 0; j < npars; ++j)
                             stream << setw(12) << minimiser->CovMatrix(i,j);
                         stream << '\n';
                     }
                     break;
            case 2:  stream << red << "Warning: Hesse is invalid\n" << def;
                     {
                         double * hessian = nullptr;
                         minimiser->GetHessianMatrix(hessian);
                         assert(hessian != nullptr);
                         for (unsigned int i = 0; i < npars; ++i) {
                             for (unsigned int j = 0; j < npars; ++j)
                                 stream << setw(12) << hessian[i*pol.npars + j];
                             stream << '\n';
                         }
                     }
                     break;
            case 3:  stream << red << "Warning: Expected Distance reached from the Minimum (EDM = " << minimiser->Edm() << ") is above max\n" << def;
                     break;
            case 4:  stream << red << "Warning: Reached call limit (Ncalls = " << Ncalls << ")\n" << def;
                     break;
            case 5:
            default: stream << red << "Warning: Unknown failure\n" << def;
        }
 
        if (!found) ++bestResult;
 
        if (verbose) {
            stream << "Correlations of the fit parameters:\n";
            for (unsigned int i = 0; i < npars; ++i) {
                for (unsigned int j = 0; j < npars; ++j)
                    stream << setw(12) << minimiser->Correlation(i,j);
                stream << '\n';
            }
        }
 
        if (current.ndf == 2) {
            stream << red << "Only 2 d.o.f. left. Stopping.\n" << def;
            break;
        }
 
        stream << "---\n";
    } // end of `while`-loop on number of parameters
    chi2s.pop_front(); // removing first entry, which is not a fit (chi2 = inf)
 
    const char * name = Form("smooth_%s", h->GetName());
    TF1 * f = bestResult->GetTF1(name, pol, m, M);
    stream << "Best result: " << f->GetTitle() << "\n-----\n";
 
    stream << bold << "Running F-test:\n" << def;
    PrintAllTests(Ftest, stream, bestResult->degree());
 
    if (nrep > 0) {
        stream << bold << "Running Cross-validation:\n" << def;
        PrintAllTests(Xval, stream, bestResult->degree());
    }
 
    return f;
}

GetSmoothFit() [3/4]

TF1* Step::GetSmoothFit	(	TH1 *	h,
		TH2 *	cov,
		unsigned int	maxdegree,
		EarlyStopping	criterion = EarlyStopping::None,
		double	stopParam = 1,
		UInt_t	nrep = 0,
		std::ostream &	stream = std::cout
	)

Function overloading, where the full range of the integral is used in the fit.

{
    int im = 1, iM = h->GetNbinsX();
    return GetSmoothFit<T,Log,Exp>(h, cov, im, iM,
            maxdegree, criterion, stopParam, nrep, stream);
}

GetSmoothFit() [4/4]

TF1* Step::GetSmoothFit	(	TH1 *	h,
		unsigned int	maxdegree,
		EarlyStopping	criterion = EarlyStopping::None,
		double	stopParam = 1,
		UInt_t	nrep = 0,
		std::ostream &	stream = std::cout
	)

Function overloading, in case of no bin-to-bin correlations and with fit in full range.

{
    return GetSmoothFit<T,Log,Exp>(h, nullptr,
            maxdegree, criterion, stopParam, nrep, stream);
}

identity()

double Step::identity

(

double

x

)

Identical operation, used as default for template arguments of Step::GetSmoothFit().

{ return x; }

operator<<()

std::ostream& Step::operator<<	(	std::ostream &	stream,
		const Result &	result
	)

Operator overloading to print result from Step::chi2s.

{
    using namespace std;
    if (result.valid) {
        stream << bold << "Chi2/ndf = " << result.chi2 << " / " << result.ndf
                            << " = " << result.x2n() << "\u00B1" << result.x2nErr()
                            << " (prob = " << result.prob() << ")\n" << def;
        size_t nrep = result.chi2sT.size();
        assert(nrep == result.chi2sV.size());
        stream << result.X.size() << " parameters: ";
        for (auto x: result.X) stream << std::setw(12) << x;
        stream << '\n';
        if (nrep > 0)
            stream << nrep <<   " training replicas: Chi2T/" <<  result.ndf                    << " = " << result.x2nT() << '\n'
                   << nrep << " validation replicas: Chi2V/" << (result.ndf + result.X.size()) << " = " << result.x2nV() << '\n';
    }
    else stream << "Not a fit\n";
    return stream;
}

PrintAllTests()

void Step::PrintAllTests	(	std::function< double(const Result &, const Result &)>	test,
		std::ostream &	stream = std::cout,
		size_t	best = 0
	)

Print all F-tests based on the direct-fit results stored in Step::chi2s.

{
    using namespace std;
 
    int nGoodFits = 0;
    for (auto result = begin(chi2s); result != end(chi2s); ++result)
        if (result->valid) ++nGoodFits;
    if (nGoodFits < 2) {
        stream << red << "Only one valid fit. No test to run.\n" << def;
        return;
    }
    stream << "   " << bold;
    for (auto result = begin(chi2s); result != prev(end(chi2s)); ++result) {
        // loop over all results except the highest one
        size_t deg = result->degree();
        if (deg == best) stream << green;
        if (!result->valid) stream << red;
        stream << setw(10) << deg << def;
    }
    stream << '\n';
    for (auto result2 = next(begin(chi2s)); result2 != end(chi2s); ++result2) {
        // loop over all results
        size_t deg2 = result2->degree();
        if (deg2 == best) stream << green;
        if (!result2->valid) stream << red;
        stream << setw(3) << deg2 << def;
        if (deg2 == best) stream << def;
        stream << setprecision(4);
        for (auto result1 = begin(chi2s); result1 != result2; ++result1) {
            // compare to all results with less parameters
            try            { stream << setw(10) << test(*result1, *result2); }
            catch (double) { stream << red << "         -" << def; }
        }
        stream << bold << '\n';
    }
    stream << def;
}

Xval()

double Step::Xval	(	const Result &	current,
		const Result &	former
	)

Cross-validation with replicas

Wikipedia

{
    if (former.chi2sV.size() != current.chi2sV.size())
        return 0;
 
    size_t nrep = former.chi2sV.size();
    if (nrep == 0) return 0;
 
    int Nf = former .ndf + former .X.size(),
        Nc = current.ndf + current.X.size();
    double nBetter = 0;
    for (size_t i = 0; i < nrep; ++i)
        if (former.chi2sV[i]/Nf < current.chi2sV[i]/Nc)
            ++nBetter;
    nBetter /= nrep;
    return nBetter;
}

Variable Documentation

bestResult

auto bestResult = chi2s.begin()

static

best result from last call of Step::GetSmoothFit()

bold

const char * bold = "\e[1m"

static

chi2s

std::deque<Result> chi2s

static

results from last call of Step::GetSmoothFit()

def

const char * def = "\e[0m"

static

eps

const auto eps = std::numeric_limits<double>::epsilon()

static

green

const char* green = "\x1B[32m"

static

inf

const auto inf = std::numeric_limits<double>::infinity()

static

red

const char * red = "\x1B[31m"

static

verbose

bool verbose = false

static