Skip to content
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

Jump to bottom

[RF] Avoid code duplication with new private Algorithms.h file #11962

Merged
merged 1 commit into from Jan 5, 2023
Merged

[RF] Avoid code duplication with new private Algorithms.h file #11962

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
11 changes: 0 additions & 11 deletions roofit/roofit/inc/RooMomentMorphFuncND.h
View file
Open in desktop
Expand Up @@ -118,17 +118,6 @@ class RooMomentMorphFuncND : public RooAbsReal {
RooAbsReal *sumFunc(const RooArgSet *nset);
CacheElem *getCache(const RooArgSet *nset) const;

template <typename T>
struct Digits {
typename std::vector<T>::const_iterator begin;
typename std::vector<T>::const_iterator end;
typename std::vector<T>::const_iterator me;
};

template <typename T>
static void cartesian_product(std::vector<std::vector<T>> &out, std::vector<std::vector<T>> &in);
template <typename Iterator>
static bool next_combination(const Iterator first, Iterator k, const Iterator last);
void findShape(const std::vector<double> &x) const;

friend class CacheElem;
Expand Down
104 changes: 104 additions & 0 deletions roofit/roofit/src/RooFit/Detail/Algorithms.h
View file
Open in desktop
@@ -0,0 +1,104 @@
/*
* Project: RooFit
*
* Copyright (c) 2022, CERN
*
* Redistribution and use in source and binary forms,
* with or without modification, are permitted according to the terms
* listed in LICENSE (http://roofit.sourceforge.net/license.txt)
*/

#ifndef RooFit_Detail_Algorithms_h
#define RooFit_Detail_Algorithms_h

#include <vector>

namespace RooFit {
namespace Detail {

//_____________________________________________________________________________
// from http://stackoverflow.com/a/5279601
template <typename T>
void cartesianProduct(std::vector<std::vector<T>> &out, std::vector<std::vector<T>> &in)
{
struct Digits {
typename std::vector<T>::const_iterator begin;
typename std::vector<T>::const_iterator end;
typename std::vector<T>::const_iterator me;
};

std::vector<Digits> vd;
vd.reserve(in.size());

for (auto it = in.begin(); it != in.end(); ++it) {
Digits d = {(*it).begin(), (*it).end(), (*it).begin()};
vd.push_back(d);
}

while (true) {
std::vector<T> result;
for (auto it = vd.begin(); it != vd.end(); ++it) {
result.push_back(*(it->me));
}
out.push_back(result);

for (auto it = vd.begin();;) {
++(it->me);
if (it->me == it->end) {
if (it + 1 == vd.end()) {
return;
} else {
it->me = it->begin;
++it;
}
} else {
break;
}
}
}
}

//_____________________________________________________________________________
// from http://stackoverflow.com/a/5097100/8747
template <typename Iterator>
bool nextCombination(const Iterator first, Iterator k, const Iterator last)
{
if ((first == last) || (first == k) || (last == k)) {
return false;
}
Iterator itr1 = first;
Iterator itr2 = last;
++itr1;
if (last == itr1) {
return false;
}
itr1 = last;
--itr1;
itr1 = k;
--itr2;
while (first != itr1) {
if (*--itr1 < *itr2) {
Iterator j = k;
while (!(*itr1 < *j))
++j;
iter_swap(itr1, j);
++itr1;
++j;
itr2 = k;
rotate(itr1, j, last);
while (last != j) {
++j;
++itr2;
}
rotate(k, itr2, last);
return true;
}
}
rotate(first, k, last);
return false;
}

} // namespace Detail
} // namespace RooFit

#endif
121 changes: 24 additions & 97 deletions roofit/roofit/src/RooMomentMorphFuncND.cxx
View file
Open in desktop
Expand Up @@ -22,6 +22,8 @@
#include "RooNumIntConfig.h"
#include "RooHistPdf.h"

#include "RooFit/Detail/Algorithms.h"

#include "TMath.h"
#include "TVector.h"
#include "TMap.h"
Expand All @@ -32,16 +34,17 @@
using namespace std;

// ClassImp for RooMomentMorphFuncND
ClassImp(RooMomentMorphFuncND)
ClassImp(RooMomentMorphFuncND);

// ClassImp for RooMomentMorphFuncND
ClassImp(RooMomentMorphFuncND::Grid2)
// ClassImp for RooMomentMorphFuncND
ClassImp(RooMomentMorphFuncND::Grid2);

//_____________________________________________________________________________
RooMomentMorphFuncND::RooMomentMorphFuncND()
: _cacheMgr(this, 10, true, true), _curNormSet(0), _M(0), _MSqr(0), _setting(RooMomentMorphFuncND::Linear), _useHorizMorph(true)
//_____________________________________________________________________________
RooMomentMorphFuncND::RooMomentMorphFuncND()
: _cacheMgr(this, 10, true, true), _curNormSet(0), _M(0), _MSqr(0), _setting(RooMomentMorphFuncND::Linear),
_useHorizMorph(true)
{
TRACE_CREATE
TRACE_CREATE;
}

//_____________________________________________________________________________
Expand All @@ -63,7 +66,7 @@ RooMomentMorphFuncND::RooMomentMorphFuncND(const char *name, const char *title,
// general initialization
initialize();

TRACE_CREATE
TRACE_CREATE;
}

//_____________________________________________________________________________
Expand Down Expand Up @@ -100,7 +103,7 @@ RooMomentMorphFuncND::RooMomentMorphFuncND(const char *name, const char *title,
// general initialization
initialize();

TRACE_CREATE
TRACE_CREATE;
}

//_____________________________________________________________________________
Expand All @@ -115,7 +118,7 @@ RooMomentMorphFuncND::RooMomentMorphFuncND(const char *name, const char *title,
TVectorD mrefpoints(mrefList.getSize());
Int_t i = 0;
for (auto *mref : mrefList) {
if (!dynamic_cast<RooAbsReal*>(mref)) {
if (!dynamic_cast<RooAbsReal *>(mref)) {
coutE(InputArguments) << "RooMomentMorphFuncND::ctor(" << GetName() << ") ERROR: mref " << mref->GetName()
<< " is not of type RooAbsReal" << endl;
throw string("RooMomentMorphFuncND::ctor() ERROR mref is not of type RooAbsReal");
Expand All @@ -124,7 +127,7 @@ RooMomentMorphFuncND::RooMomentMorphFuncND(const char *name, const char *title,
coutW(InputArguments) << "RooMomentMorphFuncND::ctor(" << GetName() << ") WARNING mref point " << i
<< " is not a constant, taking a snapshot of its value" << endl;
}
mrefpoints[i] = static_cast<RooAbsReal*>(mref)->getVal();
mrefpoints[i] = static_cast<RooAbsReal *>(mref)->getVal();
++i;
}

Expand Down Expand Up @@ -152,7 +155,7 @@ RooMomentMorphFuncND::RooMomentMorphFuncND(const char *name, const char *title,
// general initialization
initialize();

TRACE_CREATE
TRACE_CREATE;
}

//_____________________________________________________________________________
Expand All @@ -166,7 +169,7 @@ RooMomentMorphFuncND::RooMomentMorphFuncND(const RooMomentMorphFuncND &other, co
// general initialization
initialize();

TRACE_CREATE
TRACE_CREATE;
}

//_____________________________________________________________________________
Expand All @@ -177,7 +180,7 @@ RooMomentMorphFuncND::~RooMomentMorphFuncND()
if (_MSqr)
delete _MSqr;

TRACE_DESTROY
TRACE_DESTROY;
}

//_____________________________________________________________________________
Expand Down Expand Up @@ -206,41 +209,6 @@ void RooMomentMorphFuncND::initializeObservables(const RooArgList &obsList)
}
}

//_____________________________________________________________________________
// from http://stackoverflow.com/a/5279601
template <typename T>
void RooMomentMorphFuncND::cartesian_product(vector<vector<T>> &out, vector<vector<T>> &in)
{
vector<Digits<T>> vd;

for (typename vector<vector<T>>::const_iterator it = in.begin(); it != in.end(); ++it) {
Digits<T> d = {(*it).begin(), (*it).end(), (*it).begin()};
vd.push_back(d);
}

while (1) {
vector<T> result;
for (typename vector<Digits<T>>::const_iterator it = vd.begin(); it != vd.end(); ++it) {
result.push_back(*(it->me));
}
out.push_back(result);

for (typename vector<Digits<T>>::iterator it = vd.begin();;) {
++(it->me);
if (it->me == it->end) {
if (it + 1 == vd.end()) {
return;
} else {
it->me = it->begin;
++it;
}
} else {
break;
}
}
}
}

//_____________________________________________________________________________
void RooMomentMorphFuncND::initialize()
{
Expand Down Expand Up @@ -293,7 +261,7 @@ void RooMomentMorphFuncND::initialize()
}

vector<vector<int>> output;
cartesian_product(output, powers);
RooFit::Detail::cartesianProduct(output, powers);
int nCombs = output.size();

for (int k = 0; k < nPdf; ++k) {
Expand Down Expand Up @@ -325,9 +293,7 @@ RooMomentMorphFuncND::Grid2::Grid2(const RooMomentMorphFuncND::Grid2 &other)
}

//_____________________________________________________________________________
RooMomentMorphFuncND::Grid2::~Grid2()
{
}
RooMomentMorphFuncND::Grid2::~Grid2() {}

//_____________________________________________________________________________
void RooMomentMorphFuncND::Grid2::addPdf(const RooAbsReal &pdf, int bin_x)
Expand Down Expand Up @@ -585,45 +551,6 @@ const RooRealVar *RooMomentMorphFuncND::CacheElem::frac(int i) const
return (RooRealVar *)(_frac.at(i));
}

//_____________________________________________________________________________
// from http://stackoverflow.com/a/5097100/8747
template <typename Iterator>
bool RooMomentMorphFuncND::next_combination(const Iterator first, Iterator k, const Iterator last)
{
if ((first == last) || (first == k) || (last == k)) {
return false;
}
Iterator itr1 = first;
Iterator itr2 = last;
++itr1;
if (last == itr1) {
return false;
}
itr1 = last;
--itr1;
itr1 = k;
--itr2;
while (first != itr1) {
if (*--itr1 < *itr2) {
Iterator j = k;
while (!(*itr1 < *j)) ++j;
iter_swap(itr1, j);
++itr1;
++j;
itr2 = k;
rotate(itr1, j, last);
while (last != j) {
++j;
++itr2;
}
rotate(k, itr2, last);
return true;
}
}
rotate(first, k, last);
return false;
}

//_____________________________________________________________________________
void RooMomentMorphFuncND::CacheElem::calculateFractions(const RooMomentMorphFuncND &self, bool verbose) const
{
Expand Down Expand Up @@ -651,7 +578,7 @@ void RooMomentMorphFuncND::CacheElem::calculateFractions(const RooMomentMorphFun
}

vector<vector<int>> output;
cartesian_product(output, powers);
RooFit::Detail::cartesianProduct(output, powers);
int nCombs = output.size();

vector<double> deltavec(nPdf, 1.0);
Expand Down Expand Up @@ -718,7 +645,7 @@ void RooMomentMorphFuncND::CacheElem::calculateFractions(const RooMomentMorphFun
vector<double> mtmp;

// loop over parList
for(auto * m : static_range_cast<RooRealVar*>(self._parList)) {
for (auto *m : static_range_cast<RooRealVar *>(self._parList)) {
mtmp.push_back(m->getVal());
}

Expand All @@ -742,7 +669,7 @@ void RooMomentMorphFuncND::CacheElem::calculateFractions(const RooMomentMorphFun
}
deltavec[nperm + 1] = dtmp;
nperm++;
} while (next_combination(xtmp.begin(), xtmp.begin() + iperm, xtmp.end()));
} while (RooFit::Detail::nextCombination(xtmp.begin(), xtmp.begin() + iperm, xtmp.end()));
}

double origFrac1(0.), origFrac2(0.);
Expand Down Expand Up @@ -803,7 +730,7 @@ void RooMomentMorphFuncND::findShape(const vector<double> &x) const
}

vector<vector<double>> output;
cartesian_product(output, boundaries);
RooFit::Detail::cartesianProduct(output, boundaries);
_squareVec = output;

for (int isq = 0; isq < depth; isq++) {
Expand Down Expand Up @@ -848,7 +775,7 @@ void RooMomentMorphFuncND::findShape(const vector<double> &x) const
}
M(k, nperm + 1) = dtmp;
nperm++;
} while (next_combination(xtmp.begin(), xtmp.begin() + iperm, xtmp.end()));
} while (RooFit::Detail::nextCombination(xtmp.begin(), xtmp.begin() + iperm, xtmp.end()));
}
}

Expand Down