Skip to content

Commit

Permalink
Merge pull request #123 from mach3-software/feature_SIMDatFD
Browse files Browse the repository at this point in the history 
SIMD for FD spline weight calculation
  • Loading branch information
@KSkwarczynski
KSkwarczynski authored Sep 18, 2024
2 parents 6ee07d3 + e39ca71 commit da7b1fa
Show file tree
Hide file tree
Showing 2 changed files with 63 additions and 57 deletions.
60 changes: 34 additions & 26 deletions samplePDF/samplePDFFDBase.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,11 +8,13 @@ samplePDFFDBase::samplePDFFDBase(double pot, std::string mc_version, covarianceX
{
(void) pot;
(void) mc_version;
std::cout << "-------------------------------------------------------------------" <<std::endl;
std::cout << "Creating samplePDFFDBase object.." << "\n" << std::endl;

MACH3LOG_INFO("-------------------------------------------------------------------");
MACH3LOG_INFO("Creating samplePDFFDBase object..");
//ETA - safety feature so you can't pass a NULL xsec_cov
if(xsec_cov == NULL){std::cerr << "[ERROR:] You've passed me a NULL xsec covariance matrix... I need this to setup splines!" << std::endl; throw;}
if(xsec_cov == NULL){
MACH3LOG_ERROR("[ERROR:] You've passed me a NULL xsec covariance matrix... I need this to setup splines!");
throw MaCh3Exception(__FILE__ , __LINE__ );}

samplePDFFD_array = nullptr;
samplePDFFD_data = nullptr;
Expand All @@ -34,7 +36,7 @@ samplePDFFDBase::samplePDFFDBase(double pot, std::string mc_version, covarianceX

samplePDFFDBase::~samplePDFFDBase()
{
std::cout << "I'm deleting samplePDFFDBase" << std::endl;
MACH3LOG_INFO("I'm deleting samplePDFFDBase");

for (unsigned int yBin=0;yBin<(YBinEdges.size()-1);yBin++) {
if(samplePDFFD_array != nullptr){delete[] samplePDFFD_array[yBin];}
Expand Down Expand Up @@ -119,8 +121,8 @@ void samplePDFFDBase::SetupSampleBinning(){
set1DBinning(SampleXBins);
}
else if(XVarStr.length() > 0 && YVarStr.length() > 0){
std::cout << "Setting Up 2D binning" << std::endl;
std::cout << XVarStr << " : " << YVarStr << std::endl;
MACH3LOG_INFO("Setting Up 2D binning");
MACH3LOG_INFO("{} : {}", XVarStr, YVarStr);
set2DBinning(SampleXBins, SampleYBins);
}
}
Expand Down Expand Up @@ -253,7 +255,7 @@ void samplePDFFDBase::calcOscWeights(int nutype, int oscnutype, double *en, doub

if (std::isnan(w[10]))
{
std::cerr << "WARNING: ProbGPU oscillation weight returned NaN! " << w[10] << std::endl;
MACH3LOG_ERROR("WARNING: ProbGPU oscillation weight returned NaN! {}", w[10]);
}
}
#endif
Expand Down Expand Up @@ -716,13 +718,13 @@ double samplePDFFDBase::CalcXsecWeightNorm(const int iSample, const int iEvent)
//DB Adding in Oscillator class support for smeared oscillation probabilities
void samplePDFFDBase::SetOscillator(Oscillator* Osc_) {
#if defined (USE_PROB3)
std::cerr << "Atmospheric Oscillator only defined using CUDAProb3 - USE_PROB3 is defined and indicates that Prob3++/probGPU is being used" << std::endl;
throw;
MACH3LOG_ERROR("Atmospheric Oscillator only defined using CUDAProb3 - USE_PROB3 is defined and indicates that Prob3++/probGPU is being used");
throw MaCh3Exception(__FILE__ , __LINE__ );
#endif

Osc = Osc_;
std::cout << "Set Oscillator" << std::endl;

MACH3LOG_INFO("Set Oscillator");
FindEventOscBin();
}

Expand All @@ -732,12 +734,12 @@ void samplePDFFDBase::FindEventOscBin() {
MCSamples[i].osc_w_pointer[j] = Osc->retPointer(MCSamples[i].nutype,MCSamples[i].oscnutype,*(MCSamples[i].rw_etru[j]),MCSamples[i].rw_truecz[j]);
}
}
std::cout << "Set all oscillation pointers to Oscillator" << std::endl;
MACH3LOG_INFO("Set all oscillation pointers to Oscillator");
}

void samplePDFFDBase::SetXsecCov(covarianceXsec *xsec){

std::cout << "SETTING UP XSEC COV!!" << std::endl;
MACH3LOG_INFO("SETTING UP XSEC COV!!");
XsecCov = xsec;

// Get the map between the normalisation parameters index, their name, what mode they should apply to, and what target
Expand All @@ -749,9 +751,9 @@ void samplePDFFDBase::SetXsecCov(covarianceXsec *xsec){
funcParsNames = XsecCov->GetFuncParsNamesFromDetID(SampleDetID);
funcParsIndex = XsecCov->GetFuncParsIndexFromDetID(SampleDetID);

std::cout << "Found " << xsec_norms.size() << " normalisation parameters" << std::endl;
std::cout << "Found " << funcParsNames.size() << " functional parameters" << std::endl;

MACH3LOG_INFO("Found {} normalisation parameters", xsec_norms.size());
MACH3LOG_INFO("Found {} functional parameters", funcParsNames.size());
return;
}

Expand All @@ -768,8 +770,8 @@ void samplePDFFDBase::SetOscCov(covarianceOsc* osc_cov){
void samplePDFFDBase::SetupNormParameters(){

if(!XsecCov){
std::cout << "XsecCov is not setup!" << std::endl;
throw;
MACH3LOG_ERROR("XsecCov is not setup!");
throw MaCh3Exception(__FILE__ , __LINE__ );
}

// Assign xsec norm bins in MCSamples tree
Expand Down Expand Up @@ -804,7 +806,7 @@ void samplePDFFDBase::CalcXsecNormsBins(int iSample){

fdmc_base *fdobj = &MCSamples[iSample];

std::cout << "FD Object has " << fdobj->nEvents << " events" << std::endl;
MACH3LOG_INFO("FD Object has {} events", fdobj->nEvents);

for(int iEvent=0; iEvent < fdobj->nEvents; ++iEvent){
std::list< int > XsecBins = {};
Expand Down Expand Up @@ -1342,7 +1344,9 @@ void samplePDFFDBase::addData(std::vector< std::vector <double> > &data) {
dathist = NULL;
dathist2d->Reset();

if (GetNDim()!=2) {std::cerr << "Trying to set a 2D 'data' histogram in a 1D sample - Quitting" << std::endl; throw;}
if (GetNDim()!=2) {
MACH3LOG_ERROR("Trying to set a 2D 'data' histogram in a 1D sample - Quitting");
throw MaCh3Exception(__FILE__ , __LINE__ );}

for (int i = 0; i < int(dataSample2D->size()); i++) {
dathist2d->Fill(dataSample2D->at(0)[i],dataSample2D->at(1)[i]);
Expand All @@ -1363,13 +1367,15 @@ void samplePDFFDBase::addData(std::vector< std::vector <double> > &data) {
}

void samplePDFFDBase::addData(TH1D* Data) {
std::cout << "adding 1D data histogram : " << Data->GetName() << " with " << Data->Integral() << " events" << std::endl;
MACH3LOG_INFO("Adding 1D data histogram : {} with {} events", Data->GetName(), Data->Integral());
dathist2d = NULL;
dathist = Data;
dataSample = NULL;
dataSample2D = NULL;

if (GetNDim()!=1) {std::cerr << "Trying to set a 1D 'data' histogram in a 2D sample - Quitting" << std::endl; throw;}
if (GetNDim()!=1) {
MACH3LOG_ERROR("Trying to set a 1D 'data' histogram in a 2D sample - Quitting");
throw MaCh3Exception(__FILE__ , __LINE__ );}

int nXBins = XBinEdges.size()-1;
int nYBins = YBinEdges.size()-1;
Expand All @@ -1386,13 +1392,15 @@ void samplePDFFDBase::addData(TH1D* Data) {
}

void samplePDFFDBase::addData(TH2D* Data) {
std::cout << "adding 2D data histogram : " << Data->GetName() << " with " << Data->Integral() << " events" << std::endl;
MACH3LOG_INFO("Adding 2D data histogram : {} with {} events", Data->GetName(), Data->Integral());
dathist2d = Data;
dathist = NULL;
dataSample = NULL;
dataSample2D = NULL;

if (GetNDim()!=2) {std::cerr << "Trying to set a 2D 'data' histogram in a 1D sample - Quitting" << std::endl; throw;}
if (GetNDim()!=2) {
MACH3LOG_ERROR("Trying to set a 2D 'data' histogram in a 1D sample - Quitting");
throw MaCh3Exception(__FILE__ , __LINE__ );}

int nXBins = XBinEdges.size()-1;
int nYBins = YBinEdges.size()-1;
Expand Down Expand Up @@ -1459,7 +1467,7 @@ void samplePDFFDBase::fillSplineBins() {
double samplePDFFDBase::GetLikelihood()
{
if (samplePDFFD_data == NULL) {
std::cerr << "data sample is empty!" << std::endl;
MACH3LOG_ERROR("data sample is empty!");
return -1;
}

Expand Down
60 changes: 29 additions & 31 deletions splines/splineFDBase.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -305,7 +305,7 @@ void splineFDBase::CalcSplineWeights()
//****************************************
{
#ifdef MULTITHREAD
#pragma omp parallel for// schedule(dynamic)
#pragma omp parallel for simd
#endif
for (unsigned int iCoeff = 0; iCoeff < uniquecoeffindices.size(); iCoeff++)
{
Expand Down Expand Up @@ -562,8 +562,8 @@ int splineFDBase::CountNumberOfLoadedSplines(bool NonFlat, int Verbosity)

if (Verbosity > 0)
{
std::cout << "Total number of splines loaded:" << FullCounter_All << std::endl;
std::cout << "Total number of non-flat splines loaded:" << FullCounter_NonFlat << std::endl;
MACH3LOG_INFO("Total number of splines loaded: {}", FullCounter_All);
MACH3LOG_INFO("Total number of non-flat splines loaded: {}", FullCounter_NonFlat);
}

if (NonFlat)
Expand Down Expand Up @@ -743,8 +743,8 @@ void splineFDBase::PrepForReweight()
}
}
// We need to grab the maximum number of knots

std::cout << "Number of combinations of Sample, OscChan, Syst and Mode which have entirely flat response:" << nCombinations_FlatSplines << " / " << nCombinations_All << std::endl;
MACH3LOG_INFO("Number of combinations of Sample, OscChan, Syst and Mode which have entirely flat response: {} / {}", nCombinations_FlatSplines, nCombinations_All);

}

//****************************************
Expand Down Expand Up @@ -782,11 +782,9 @@ void splineFDBase::getSplineCoeff_SepMany(int splineindex, _float_* &xArray, _fl
manyArray[i*4+3] = _float_(d);

if((xArray[i] == -999) | (manyArray[i*4] == -999) | (manyArray[i*4+1] == -999) | (manyArray[i*4+2] == -999) | (manyArray[i*4+3] == -999)){
std::cerr << "*********** Bad params in getSplineCoeff_SepMany() ************"<<std::endl;
std::cerr << "pre cast to _float_ (x, y, b, c, d) = "<<x<<", "<<y<<", "<<b<<", "<<c<<", "<<d<<std::endl;
std::cerr << "post cast to float (x, y, b, c, d) = "<<xArray[i]<<", "<<manyArray[i*4]<<", "<<manyArray[i*4+1]<<", "<<manyArray[i*4+2]<<", "<<manyArray[i*4+3]<<std::endl;
std::cerr<<__FILE__<<"::"<<__LINE__<<std::endl;
std::cerr << "***************************************************************"<<std::endl;
MACH3LOG_ERROR("*********** Bad params in getSplineCoeff_SepMany() ************");
MACH3LOG_ERROR("pre cast to _float_ (x, y, b, c, d) = {}, {}, {}, {}, {}",x, y, b, c, d);
MACH3LOG_ERROR("post cast to float (x, y, b, c, d) = {}, {}, {}, {}, {}",xArray[i], manyArray[i*4], manyArray[i*4+1], manyArray[i*4+2], manyArray[i*4+3]);
throw MaCh3Exception(__FILE__ , __LINE__ );
}
}
Expand Down Expand Up @@ -821,7 +819,7 @@ int splineFDBase::getSampleIndex(std::string SampleName){
}
if (SampleIndex == -1)
{
std::cerr << "Sample name not found : "<<SampleName << std::endl;
MACH3LOG_ERROR("Sample name not found : {}", SampleName);
throw MaCh3Exception(__FILE__ , __LINE__ );
}
return SampleIndex;
Expand All @@ -848,12 +846,12 @@ void splineFDBase::PrintArrayDetails(std::string SampleName)
{
int iSample = getSampleIndex(SampleName);
int nOscChannels = indexvec[iSample].size();
std::cout << "Sample " << iSample << " has " << nOscChannels << " oscillation channels" << std::endl;

MACH3LOG_INFO("Sample {} has {} oscillation channels", iSample, nOscChannels);
for (int iOscChan = 0; iOscChan < nOscChannels; iOscChan++)
{
int nSysts = indexvec[iSample][iOscChan].size();
std::cout << "Oscillation channel " << iOscChan << " has " << nSysts << " systematics" << std::endl;
MACH3LOG_INFO("Oscillation channel {} has {} systematics", iOscChan, nSysts);
}
}

Expand Down Expand Up @@ -920,56 +918,56 @@ bool splineFDBase::isValidSplineIndex(std::string SampleName, int iOscChan, int

if (iSample < 0 || iSample >= (int)indexvec.size())
{
std::cerr << "Sample index is invalid! 0 <= Index < " << indexvec.size() << std::endl;
MACH3LOG_ERROR("Sample index is invalid! 0 <= Index < {} ", indexvec.size());
isValid = false;
}

if (iOscChan < 0 || iOscChan >= (int)indexvec[iSample].size())
{
std::cerr << "OscChan index is invalid! 0 <= Index < " << indexvec[iSample].size() << std::endl;
MACH3LOG_ERROR("OscChan index is invalid! 0 <= Index < {} ", indexvec[iSample].size());
isValid = false;
}

if (iSyst < 0 || iSyst >= (int)indexvec[iSample][iOscChan].size())
{
std::cerr << "Syst index is invalid! 0 <= Index < " << indexvec[iSample][iOscChan].size() << std::endl;
MACH3LOG_ERROR("Syst index is invalid! 0 <= Index < {} ", indexvec[iSample][iOscChan].size());
isValid = false;
}

if (iMode < 0 || iMode >= (int)indexvec[iSample][iOscChan][iSyst].size())
{
std::cerr << "Mode index is invalid! 0 <= Index < " << indexvec[iSample][iOscChan][iSyst].size() << std::endl;
MACH3LOG_ERROR("Mode index is invalid! 0 <= Index < {} ", indexvec[iSample][iOscChan][iSyst].size());
isValid = false;
}

if (iVar1 < 0 || iVar1 >= (int)indexvec[iSample][iOscChan][iSyst][iMode].size())
{
std::cerr << "Var1 index is invalid! 0 <= Index < " << indexvec[iSample][iOscChan][iSyst][iMode].size() << std::endl;
MACH3LOG_ERROR("Var1 index is invalid! 0 <= Index < {} ", indexvec[iSample][iOscChan][iSyst][iMode].size());
isValid = false;
}

if (iVar2 < 0 || iVar2 >= (int)indexvec[iSample][iOscChan][iSyst][iMode][iVar1].size())
{
std::cerr << "Var2 index is invalid! 0 <= Index < " << indexvec[iSample][iOscChan][iSyst][iMode][iVar1].size() << std::endl;
MACH3LOG_ERROR("Var2 index is invalid! 0 <= Index < {} ", indexvec[iSample][iOscChan][iSyst][iMode][iVar1].size());
isValid = false;
}

if (iVar3 < 0 || iVar3 >= (int)indexvec[iSample][iOscChan][iSyst][iMode][iVar1][iVar2].size())
{
std::cerr << "Var3 index is invalid! 0 <= Index < " << indexvec[iSample][iOscChan][iSyst][iMode][iVar1][iVar2].size() << std::endl;
MACH3LOG_ERROR("Var3 index is invalid! 0 <= Index < {} ", indexvec[iSample][iOscChan][iSyst][iMode][iVar1][iVar2].size());
isValid = false;
}

if (!isValid)
{
std::cerr << "Given iSample:" << iSample << std::endl;
std::cerr << "Given iOscChan:" << iOscChan << std::endl;
std::cerr << "Given iSyst:" << iSyst << std::endl;
std::cerr << "Given iMode:" << iMode << std::endl;
std::cerr << "Given iVar1:" << iVar1 << std::endl;
std::cerr << "Given iVar2:" << iVar2 << std::endl;
std::cerr << "Given iVar3:" << iVar3 << std::endl;
std::cerr << "Come visit me at : "<<__FILE__<<" : "<<__LINE__<<std::endl;
{
MACH3LOG_ERROR("Given iSample: {}", iSample);
MACH3LOG_ERROR("Given iOscChan: {}", iOscChan);
MACH3LOG_ERROR("Given iSyst: {}", iSyst);
MACH3LOG_ERROR("Given iMode: {}", iMode);
MACH3LOG_ERROR("Given iVar1: {}", iVar1);
MACH3LOG_ERROR("Given iVar2: {}", iVar2);
MACH3LOG_ERROR("Given iVar3: {}", iVar3);
MACH3LOG_ERROR("Come visit me at : {} : {}", __FILE__, __LINE__);
throw MaCh3Exception(__FILE__ , __LINE__ );
}

Expand All @@ -982,7 +980,7 @@ void splineFDBase::PrintBinning(TAxis *Axis)
{
const int NBins = Axis->GetNbins();
const double *BinEdges = Axis->GetXbins()->GetArray();

std::cout << "\t";
for (int iBin = 0; iBin < (NBins + 1); iBin++)
{
Expand Down

0 comments on commit da7b1fa

Please sign in to comment.