Skip to content

Several fixes and improvements for the vectorized implementations #649

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
lmoneta merged 4 commits into from
Jun 14, 2017
Merged

Several fixes and improvements for the vectorized implementations #649

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
aac3008
Improve vectorized implementation.
xvallspl Jun 9, 2017
e661230
Fix warning
xvallspl Jun 13, 2017
59abd4a
Fix evaluation of vectorized functions in LogL fitting
xvallspl Jun 14, 2017
079ba08
Improve formatting
xvallspl Jun 14, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion hist/hist/src/TF1.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2690,7 +2690,7 @@ void TF1::Print(Option_t *option) const
if (fFunctor)
printf("Compiled based function: %s based on a functor object. Ndim = %d, Npar = %d\n", GetName(), GetNpar(), GetNdim());
else {
printf("Function based on a list of points from a compiled based function: %s. Ndim = %d, Npar = %d, Npx = %d\n", GetName(), GetNpar(), GetNdim(), fSave.size());
printf("Function based on a list of points from a compiled based function: %s. Ndim = %d, Npar = %d, Npx = %zu\n", GetName(), GetNpar(), GetNdim(), fSave.size());
if (fSave.empty())
Warning("Print", "Function %s is based on a list of points but list is empty", GetName());
}
Expand Down
30 changes: 27 additions & 3 deletions math/mathcore/inc/Fit/FitConfig.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -103,10 +103,34 @@ class FitConfig {
set the parameter settings from a model function.
Create always new parameter setting list from a given model function
*/
void CreateParamsSettings(const ROOT::Math::IParamMultiFunction & func);
#ifdef R__HAS_VECCORE
void CreateParamsSettings(const ROOT::Math::IParamMultiFunctionTempl<ROOT::Double_v> & func);
template <class T>
void CreateParamsSettings(const ROOT::Math::IParamMultiFunctionTempl<T> &func) {
// initialize from model function
// set the parameters values from the function
unsigned int npar = func.NPar();
const double *begin = func.Parameters();
if (begin == 0) {
fSettings = std::vector<ParameterSettings>(npar);
return;
}

fSettings.clear();
fSettings.reserve(npar);
const double *end = begin + npar;
unsigned int i = 0;
for (const double *ipar = begin; ipar != end; ++ipar) {
double val = *ipar;
double step = 0.3 * std::fabs(val); // step size is 30% of par value
// double step = 2.0*std::fabs(val); // step size is 30% of par value
if (val == 0) step = 0.3;

fSettings.push_back(ParameterSettings(func.ParameterName(i), val, step));
#ifdef DEBUG
std::cout << "FitConfig: add parameter " << func.ParameterName(i) << " val = " << val << std::endl;
#endif
i++;
}
}

/**
set the parameter settings from number of parameters and a vector of values and optionally step values. If there are not existing or number of parameters does not match existing one, create a new parameter setting list.
Expand Down
87 changes: 53 additions & 34 deletions math/mathcore/inc/Fit/FitUtil.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -223,7 +223,7 @@ namespace FitUtil {
vecCore::Load<ROOT::Double_v>(xx, x);
const double *p0 = p;
auto res = (*f)( &xx, (const double *)p0);
return res[0];
return vecCore::Get<ROOT::Double_v>(res, 0);
}
#endif

Expand Down Expand Up @@ -274,16 +274,19 @@ namespace FitUtil {
return also nPoints as the effective number of used points in the LogL evaluation
*/
void EvaluateLogLGradient(const IModelFunction & func, const UnBinData & data, const double * x, double * grad, unsigned int & nPoints);

#ifdef R__HAS_VECCORE
void EvaluateLogLGradient(const IModelFunctionTempl<ROOT::Double_v> &, const UnBinData &, const double *, double *, unsigned int & ) ;
template <class NotCompileIfScalarBackend = std::enable_if<!(std::is_same<double, ROOT::Double_v>::value)>>
void EvaluateLogLGradient(const IModelFunctionTempl<ROOT::Double_v> &, const UnBinData &, const double *, double *, unsigned int & ) {}
#endif

/**
evaluate the Poisson LogL given a model function and the data at the point x.
return also nPoints as the effective number of used points in the LogL evaluation
By default is extended, pass extedend to false if want to be not extended (MultiNomial)
*/
double EvaluatePoissonLogL(const IModelFunction & func, const BinData & data, const double * x, int iWeight, bool extended, unsigned int & nPoints, const unsigned int &executionPolicy, unsigned nChunks = 0);
double EvaluatePoissonLogL(const IModelFunction &func, const BinData &data, const double *x, int iWeight, bool extended,
unsigned int &nPoints, const unsigned int &executionPolicy, unsigned nChunks = 0);

/**
evaluate the Poisson LogL given a model function and the data at the point x.
Expand All @@ -308,9 +311,20 @@ namespace FitUtil {
is used
*/
double EvaluatePdf(const IModelFunction & func, const UnBinData & data, const double * x, unsigned int ipoint, double * g = 0);

#ifdef R__HAS_VECCORE
double EvaluatePdf(const IModelFunctionTempl<ROOT::Double_v> & func, const UnBinData & data, const double * p, unsigned int i, double *);
template <class NotCompileIfScalarBackend = std::enable_if<!(std::is_same<double, ROOT::Double_v>::value)>>
double EvaluatePdf(const IModelFunctionTempl<ROOT::Double_v> &func, const UnBinData &data, const double *p, unsigned int i, double *) {
// evaluate the pdf contribution to the generic logl function in case of bin data
// return actually the log of the pdf and its derivatives
// func.SetParameters(p);
const auto x = vecCore::FromPtr<ROOT::Double_v>(data.GetCoordComponent(i, 0));
auto fval = func(&x, p);
auto logPdf = ROOT::Math::Util::EvalLog(fval);
return vecCore::Get<ROOT::Double_v>(logPdf, 0);
}
#endif

/**
evaluate the pdf contribution to the Poisson LogL given a model function and the BinPoint data.
If the pointer g is not null evaluate also the gradient of the Poisson pdf.
Expand Down Expand Up @@ -362,7 +376,7 @@ namespace FitUtil {
auto invErrorptr = (invError != nullptr) ? invError : &ones.front();
vecCore::Load<T>(invErrorVec, invErrorptr);

const T * x;
const T *x;
if(data.NDim() > 1) {
std::vector<T> xc;
xc.resize(data.NDim());
Expand Down Expand Up @@ -480,19 +494,20 @@ namespace FitUtil {
if(data.NDim() > 1) {
std::vector<T> x(data.NDim());
for (unsigned int j = 0; j < data.NDim(); ++j)
x[j] = *data.GetCoordComponent(i, j);
vecCore::Load<T>(x[j], data.GetCoordComponent(i, j));
#ifdef USE_PARAMCACHE
fval = func ( x.data() );
fval = func(x.data());
#else
fval = func ( x.data(), p );
fval = func(x.data(), p);
#endif
// one -dim case
// one -dim case
} else {
const auto x = data.GetCoordComponent(i, 0);
T x;
vecCore::Load<T>(x, data.GetCoordComponent(i, 0));
#ifdef USE_PARAMCACHE
fval = func ( x );
fval = func(&x);
#else
fval = func ( x, p );
fval = func(&x, p);
#endif
}

Expand Down Expand Up @@ -618,14 +633,14 @@ namespace FitUtil {

// reset the number of fitting data points
// nPoints = n;
// std::cout<<", n: "<<nPoints<<std::endl;
// std::cout<<", n: "<<nPoints<<std::endl;
nPoints = 0;
return -logl;

}

static double EvalPoissonLogL(const IModelFunctionTempl<T> &func, const BinData & data, const double * p, int iWeight,
bool extended, unsigned int &nPoints, const unsigned int &executionPolicy, unsigned nChunks = 0)
static double EvalPoissonLogL(const IModelFunctionTempl<T> &func, const BinData &data, const double *p, int iWeight, bool extended,
unsigned int &nPoints, const unsigned int &executionPolicy, unsigned nChunks = 0)
{
// evaluate the Poisson Log Likelihood
// for binned likelihood fits
Expand All @@ -651,30 +666,32 @@ namespace FitUtil {
// get fit option and check case of using integral of bins
const DataOptions &fitOpt = data.Opt();
if (fitOpt.fExpErrors || fitOpt.fIntegral)
Error("FitUtil::EvaluateChi2", "The vectorized implementation doesn't support Integrals or BinVolume\n. Aborting operation.");
Error("FitUtil::EvaluateChi2",
"The vectorized implementation doesn't support Integrals or BinVolume\n. Aborting operation.");
bool useW2 = (iWeight == 2);

auto mapFunction = [&](unsigned int i) {
T y;
vecCore::Load<T>(y, data.ValuePtr(i * vecSize));
T fval{};

if(data.NDim() > 1) {
if (data.NDim() > 1) {
std::vector<T> x(data.NDim());
for (unsigned int j = 0; j < data.NDim(); ++j)
x[j] = *data.GetCoordComponent(i, j);
vecCore::Load<T>(x[j], data.GetCoordComponent(i, j));
#ifdef USE_PARAMCACHE
fval = func ( x.data() );
fval = func(x.data());
#else
fval = func ( x.data(), p );
fval = func(x.data(), p);
#endif
// one -dim case
// one -dim case
} else {
const auto x = data.GetCoordComponent(i, 0);
T x;
vecCore::Load<T>(x, data.GetCoordComponent(i, 0));
#ifdef USE_PARAMCACHE
fval = func ( x );
fval = func(&x);
#else
fval = func ( x, p );
fval = func(&x, p);
#endif
}

Expand All @@ -683,7 +700,7 @@ namespace FitUtil {
auto m = vecCore::Mask_v<T>(fval < 0.0);
vecCore::MaskedAssign<T>(fval, m, 0.0);

T nloglike{}; // negative loglikelihood
T nloglike{}; // negative loglikelihood

if (useW2) {
// apply weight correction . Effective weight is error^2/ y
Expand Down Expand Up @@ -717,18 +734,16 @@ namespace FitUtil {
// this is needed for Poisson likelihood (which are extened and not for multinomial)
// the formula below include constant term due to likelihood of saturated model (f(x) = y)
// (same formula as in Baker-Cousins paper, page 439 except a factor of 2
if (extended) nloglike = fval - y ;
if (extended) nloglike = fval - y;

vecCore::MaskedAssign<T>(nloglike, y > 0, nloglike + y * (ROOT::Math::Util::EvalLog(y) -
ROOT::Math::Util::EvalLog(fval)));
vecCore::MaskedAssign<T>(
nloglike, y > 0, nloglike + y * (ROOT::Math::Util::EvalLog(y) - ROOT::Math::Util::EvalLog(fval)));
}

return nloglike;
};

auto redFunction = [](const std::vector<T> &objs) {
return std::accumulate(objs.begin(), objs.end(), T{});
};
auto redFunction = [](const std::vector<T> &objs) { return std::accumulate(objs.begin(), objs.end(), T{}); };

T res{};
if (executionPolicy == ROOT::Fit::kSerial) {
Expand All @@ -745,7 +760,9 @@ namespace FitUtil {
// ROOT::TProcessExecutor pool;
// res = pool.MapReduce(mapFunction, ROOT::TSeq<unsigned>(0, data.Size()/vecSize), redFunction);
} else {
Error("FitUtil::Evaluate<T>::EvalPoissonLogL", "Execution policy unknown. Avalaible choices:\n 0: Serial (default)\n 1: MultiThread (requires IMT)\n");
Error(
"FitUtil::Evaluate<T>::EvalPoissonLogL",
"Execution policy unknown. Avalaible choices:\n 0: Serial (default)\n 1: MultiThread (requires IMT)\n");
}

return vecCore::Reduce(res);
Expand Down Expand Up @@ -779,14 +796,16 @@ namespace FitUtil {
// optionally the integral of function in the bin is used
return FitUtil::EvaluateChi2(func, data, p, nPoints, executionPolicy, nChunks);
}

static double EvalLogL(const IModelFunctionTempl<double> &func, const UnBinData & data, const double * p, int iWeight,
bool extended, unsigned int &nPoints, const unsigned int &executionPolicy, unsigned nChunks = 0)
{
return FitUtil::EvaluateLogL(func, data, p, iWeight, extended, nPoints, executionPolicy, nChunks);
}

static double EvalPoissonLogL(const IModelFunctionTempl<double> & func, const BinData & data, const double * p, int iWeight, bool extended,
unsigned int & nPoints, const unsigned int &executionPolicy, unsigned nChunks = 0 ) {
static double EvalPoissonLogL(const IModelFunctionTempl<double> &func, const BinData &data, const double *p, int iWeight, bool extended, unsigned int &nPoints,
const unsigned int &executionPolicy, unsigned nChunks = 0)
{
return FitUtil::EvaluatePoissonLogL(func, data, p, iWeight, extended, nPoints, executionPolicy, nChunks);
}

Expand Down
29 changes: 24 additions & 5 deletions math/mathcore/inc/Fit/Fitter.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,7 @@ Classes used for fitting (regression analysis) and estimation of parameter value
#include "Fit/FitConfig.h"
#include "Fit/FitExecutionPolicy.h"
#include "Fit/FitResult.h"
#include "Math/IParamFunctionfwd.h"
#include "Math/IParamFunction.h"
#include <memory>

namespace ROOT {
Expand Down Expand Up @@ -164,11 +164,14 @@ class Fitter {
/**
Binned Likelihood fit. Default is extended
*/
bool LikelihoodFit(const BinData & data, bool extended = true, ROOT::Fit::ExecutionPolicy executionPolicy = ROOT::Fit::kSerial) {
bool LikelihoodFit(const BinData &data, bool extended = true,
ROOT::Fit::ExecutionPolicy executionPolicy = ROOT::Fit::kSerial) {
SetData(data);
return DoBinnedLikelihoodFit(extended, executionPolicy);
}
bool LikelihoodFit(const std::shared_ptr<BinData> & data, bool extended = true, ROOT::Fit::ExecutionPolicy executionPolicy = ROOT::Fit::kSerial) {

bool LikelihoodFit(const std::shared_ptr<BinData> &data, bool extended = true,
ROOT::Fit::ExecutionPolicy executionPolicy = ROOT::Fit::kSerial) {
SetData(data);
return DoBinnedLikelihoodFit(extended, executionPolicy);
}
Expand Down Expand Up @@ -324,7 +327,8 @@ class Fitter {
Set the fitted function (model function) from a vectorized parametric function interface
*/
#ifdef R__HAS_VECCORE
void SetFunction(const IModelFunction_v & func);
template <class NotCompileIfScalarBackend = std::enable_if<!(std::is_same<double, ROOT::Double_v>::value)>>
void SetFunction(const IModelFunction_v &func);
#endif
/**
Set the fitted function from a parametric 1D function interface
Expand Down Expand Up @@ -425,7 +429,7 @@ class Fitter {
/// least square fit
bool DoLeastSquareFit(ROOT::Fit::ExecutionPolicy executionPolicy = ROOT::Fit::kSerial);
/// binned likelihood fit
bool DoBinnedLikelihoodFit( bool extended = true, ROOT::Fit::ExecutionPolicy executionPolicy = ROOT::Fit::kSerial);
bool DoBinnedLikelihoodFit(bool extended = true, ROOT::Fit::ExecutionPolicy executionPolicy = ROOT::Fit::kSerial);
/// un-binned likelihood fit
bool DoUnbinnedLikelihoodFit( bool extended = false, ROOT::Fit::ExecutionPolicy executionPolicy = ROOT::Fit::kSerial);
/// linear least square fit
Expand Down Expand Up @@ -516,6 +520,21 @@ bool Fitter::GetDataFromFCN() {
}
}

#ifdef R__HAS_VECCORE
template <class NotCompileIfScalarBackend>
void Fitter::SetFunction(const IModelFunction_v &func)
{
// set the fit model function (clone the given one and keep a copy )
// std::cout << "set a non-grad function" << std::endl;
fUseGradient = false;
fFunc_v = std::shared_ptr<IModelFunction_v>(dynamic_cast<IModelFunction_v *>(func.Clone()));
assert(fFunc_v);

// creates the parameter settings
fConfig.CreateParamsSettings(*fFunc_v);
fFunc.reset();
}
#endif

} // end namespace Fit

Expand Down
3 changes: 1 addition & 2 deletions math/mathcore/inc/Fit/LogLikelihoodFCN.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -178,8 +178,7 @@ class LogLikelihoodFCN : public BasicFCN<DerivFunType,ModelFunType,UnBinData> {

mutable std::vector<double> fGrad; // for derivatives

ROOT::Fit::ExecutionPolicy fExecutionPolicy; //Execution policy

ROOT::Fit::ExecutionPolicy fExecutionPolicy; // Execution policy
};
// define useful typedef's
// using LogLikelihoodFunction_v = LogLikelihoodFCN<ROOT::Math::IMultiGenFunction, ROOT::Math::IParametricFunctionMultiDimTempl<T>>;
Expand Down
8 changes: 5 additions & 3 deletions math/mathcore/inc/Fit/PoissonLikelihoodFCN.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -162,9 +162,11 @@ class PoissonLikelihoodFCN : public BasicFCN<DerivFunType,ModelFunType,BinData>
virtual double DoEval (const double * x) const {
this->UpdateNCalls();
#ifdef R__HAS_VECCORE
return FitUtil::Evaluate<T>::EvalPoissonLogL(BaseFCN::ModelFunction(), BaseFCN::Data(), x, fWeight, fIsExtended, fNEffPoints, fExecutionPolicy);
return FitUtil::Evaluate<T>::EvalPoissonLogL(BaseFCN::ModelFunction(), BaseFCN::Data(), x, fWeight, fIsExtended,
fNEffPoints, fExecutionPolicy);
#else
return FitUtil::EvaluatePoissonLogL(BaseFCN::ModelFunction(), BaseFCN::Data(), x, fWeight, fIsExtended, fNEffPoints, fExecutionPolicy);
return FitUtil::EvaluatePoissonLogL(BaseFCN::ModelFunction(), BaseFCN::Data(), x, fWeight, fIsExtended,
fNEffPoints, fExecutionPolicy);
#endif
}

Expand All @@ -184,7 +186,7 @@ class PoissonLikelihoodFCN : public BasicFCN<DerivFunType,ModelFunType,BinData>

mutable std::vector<double> fGrad; // for derivatives

ROOT::Fit::ExecutionPolicy fExecutionPolicy; //Execution policy
ROOT::Fit::ExecutionPolicy fExecutionPolicy; // Execution policy
};

// define useful typedef's
Expand Down
Loading