Skip to content

Commit

Permalink
[Minuit2] Migrate InitialGradientCalculator to free function
Browse files Browse the repository at this point in the history
  • Loading branch information
@guitargeek
guitargeek committed Mar 3, 2025
1 parent bc236c0 commit 962f80a
Show file tree
Hide file tree
Showing 6 changed files with 29 additions and 66 deletions.
23 changes: 1 addition & 22 deletions math/minuit2/inc/Minuit2/InitialGradientCalculator.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,30 +16,9 @@ namespace ROOT {

namespace Minuit2 {

class MnFcn;
class MnUserTransformation;
class MnMachinePrecision;

/**
Class to calculate an initial estimate of the gradient
*/
class InitialGradientCalculator : public GradientCalculator {

public:
InitialGradientCalculator(const MnFcn &fcn, const MnUserTransformation &par) : fFcn(fcn), fTransformation(par) {}

FunctionGradient operator()(const MinimumParameters &) const override;

FunctionGradient operator()(const MinimumParameters &, const FunctionGradient &) const override;

const MnFcn &Fcn() const { return fFcn; }
const MnUserTransformation &Trafo() const { return fTransformation; }
const MnMachinePrecision &Precision() const;

private:
const MnFcn &fFcn;
const MnUserTransformation &fTransformation;
};
FunctionGradient calculateInitialGradient(const MinimumParameters &, const MnUserTransformation &, double errorDef);

} // namespace Minuit2

Expand Down
1 change: 0 additions & 1 deletion math/minuit2/inc/Minuit2/Numerical2PGradientCalculator.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,6 @@ class Numerical2PGradientCalculator : public GradientCalculator {

FunctionGradient operator()(const MinimumParameters &, const FunctionGradient &) const override;

const MnFcn &Fcn() const { return fFcn; }
const MnUserTransformation &Trafo() const { return fTransformation; }
const MnMachinePrecision &Precision() const;
const MnStrategy &Strategy() const { return fStrategy; }
Expand Down
3 changes: 1 addition & 2 deletions math/minuit2/src/HessianGradientCalculator.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,8 +28,7 @@ namespace Minuit2 {
FunctionGradient HessianGradientCalculator::operator()(const MinimumParameters &par) const
{
// use initial gradient as starting point
InitialGradientCalculator gc(fFcn, fTransformation);
FunctionGradient gra = gc(par);
FunctionGradient gra = calculateInitialGradient(par, fTransformation, fFcn.ErrorDef());

return (*this)(par, gra);
}
Expand Down
50 changes: 19 additions & 31 deletions math/minuit2/src/InitialGradientCalculator.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,13 +22,13 @@ namespace ROOT {

namespace Minuit2 {

FunctionGradient InitialGradientCalculator::operator()(const MinimumParameters &par) const
/// Initial rough estimate of the gradient using the parameter step size.
FunctionGradient
calculateInitialGradient(const MinimumParameters &par, const MnUserTransformation &trafo, double errorDef)
{
// initial rough estimate of the gradient using the parameter step size

assert(par.IsValid());

unsigned int n = Trafo().VariableParameters();
unsigned int n = trafo.VariableParameters();
assert(n == par.Vec().size());

MnPrint print("InitialGradientCalculator");
Expand All @@ -38,58 +38,46 @@ FunctionGradient InitialGradientCalculator::operator()(const MinimumParameters &
MnAlgebraicVector gr(n), gr2(n), gst(n);

for (unsigned int i = 0; i < n; i++) {
unsigned int exOfIn = Trafo().ExtOfInt(i);
unsigned int exOfIn = trafo.ExtOfInt(i);

double var = par.Vec()(i);
double werr = Trafo().Parameter(exOfIn).Error();
double save1 = Trafo().Int2ext(i, var);
double werr = trafo.Parameter(exOfIn).Error();
double save1 = trafo.Int2ext(i, var);
double save2 = save1 + werr;
if (Trafo().Parameter(exOfIn).HasLimits()) {
if (Trafo().Parameter(exOfIn).HasUpperLimit() && save2 > Trafo().Parameter(exOfIn).UpperLimit())
save2 = Trafo().Parameter(exOfIn).UpperLimit();
if (trafo.Parameter(exOfIn).HasLimits()) {
if (trafo.Parameter(exOfIn).HasUpperLimit() && save2 > trafo.Parameter(exOfIn).UpperLimit())
save2 = trafo.Parameter(exOfIn).UpperLimit();
}
double var2 = Trafo().Ext2int(exOfIn, save2);
double var2 = trafo.Ext2int(exOfIn, save2);
double vplu = var2 - var;
save2 = save1 - werr;
if (Trafo().Parameter(exOfIn).HasLimits()) {
if (Trafo().Parameter(exOfIn).HasLowerLimit() && save2 < Trafo().Parameter(exOfIn).LowerLimit())
save2 = Trafo().Parameter(exOfIn).LowerLimit();
if (trafo.Parameter(exOfIn).HasLimits()) {
if (trafo.Parameter(exOfIn).HasLowerLimit() && save2 < trafo.Parameter(exOfIn).LowerLimit())
save2 = trafo.Parameter(exOfIn).LowerLimit();
}
var2 = Trafo().Ext2int(exOfIn, save2);
var2 = trafo.Ext2int(exOfIn, save2);
double vmin = var2 - var;
double gsmin = 8. * Precision().Eps2() * (std::fabs(var) + Precision().Eps2());
double gsmin = 8. * trafo.Precision().Eps2() * (std::fabs(var) + trafo.Precision().Eps2());
// protect against very small step sizes which can cause dirin to zero and then nan values in grd
double dirin = std::max(0.5 * (std::fabs(vplu) + std::fabs(vmin)), gsmin);
double g2 = 2.0 * fFcn.ErrorDef() / (dirin * dirin);
double g2 = 2.0 * errorDef / (dirin * dirin);
double gstep = std::max(gsmin, 0.1 * dirin);
double grd = g2 * dirin;
if (Trafo().Parameter(exOfIn).HasLimits()) {
if (trafo.Parameter(exOfIn).HasLimits()) {
if (gstep > 0.5)
gstep = 0.5;
}
gr(i) = grd;
gr2(i) = g2;
gst(i) = gstep;

print.Debug("Computed initial gradient for parameter", Trafo().Name(exOfIn), "value", var, "[", vmin, ",", vplu,
print.Debug("Computed initial gradient for parameter", trafo.Name(exOfIn), "value", var, "[", vmin, ",", vplu,
"]", "dirin", dirin, "grd", grd, "g2", g2);
}

return FunctionGradient(gr, gr2, gst);
}

FunctionGradient InitialGradientCalculator::operator()(const MinimumParameters &par, const FunctionGradient &) const
{
// Base class interface
return (*this)(par);
}

const MnMachinePrecision &InitialGradientCalculator::Precision() const
{
// return precision (is set in transformation class)
return fTransformation.Precision();
}

} // namespace Minuit2

} // namespace ROOT
15 changes: 7 additions & 8 deletions math/minuit2/src/Numerical2PGradientCalculator.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,8 +35,7 @@ FunctionGradient Numerical2PGradientCalculator::operator()(const MinimumParamete
{
// calculate gradient using Initial gradient calculator and from MinimumParameters object

InitialGradientCalculator gc(fFcn, fTransformation);
FunctionGradient gra = gc(par);
FunctionGradient gra = calculateInitialGradient(par, fTransformation, fFcn.ErrorDef());

return (*this)(par, gra);
}
Expand All @@ -53,7 +52,7 @@ FunctionGradient Numerical2PGradientCalculator::operator()(std::span<const doubl
par(i) = params[i];
}

double fval = Fcn()(par);
double fval = fFcn(par);

MinimumParameters minpars = MinimumParameters(par, fval);

Expand Down Expand Up @@ -81,7 +80,7 @@ operator()(const MinimumParameters &par, const FunctionGradient &Gradient) const

//print.Trace("Assumed precision eps", eps, "eps2", eps2);

double dfmin = 8. * eps2 * (std::fabs(fcnmin) + Fcn().Up());
double dfmin = 8. * eps2 * (std::fabs(fcnmin) + fFcn.Up());
double vrysml = 8. * eps * eps;
// double vrysml = std::max(1.e-4, eps2);
// std::cout<<"dfmin= "<<dfmin<<std::endl;
Expand Down Expand Up @@ -156,13 +155,13 @@ operator()(const MinimumParameters &par, const FunctionGradient &Gradient) const
stepb4 = step;
// MnAlgebraicVector pstep(n);
// pstep(i) = step;
// double fs1 = Fcn()(pstate + pstep);
// double fs2 = Fcn()(pstate - pstep);
// double fs1 = fFcn(pstate + pstep);
// double fs2 = fFcn(pstate - pstep);

x(i) = xtf + step;
double fs1 = Fcn()(x);
double fs1 = fFcn(x);
x(i) = xtf - step;
double fs2 = Fcn()(x);
double fs2 = fFcn(x);
x(i) = xtf;

double grdb4 = grd(i);
Expand Down
3 changes: 1 addition & 2 deletions math/minuit2/src/SimplexSeedGenerator.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,8 +33,7 @@ operator()(const MnFcn &fcn, const GradientCalculator &, const MnUserParameterSt
x(i) = st.IntParameters()[i];
double fcnmin = fcn(x);
MinimumParameters pa(x, fcnmin);
InitialGradientCalculator igc(fcn, st.Trafo());
FunctionGradient dgrad = igc(pa);
FunctionGradient dgrad = calculateInitialGradient(pa, st.Trafo(), fcn.ErrorDef());
MnAlgebraicSymMatrix mat(n);
double dcovar = 1.;
for (unsigned int i = 0; i < n; i++)
Expand Down

0 comments on commit 962f80a

Please sign in to comment.