Skip to content

Commit

Permalink
Merge pull request #334 from mach3-software/feature_namespace
Browse files Browse the repository at this point in the history
Breaking: Move common values under namespace
  • Loading branch information
KSkwarczynski authored Feb 7, 2025
2 parents 3efefb8 + bb55aff commit 96d89d2
Show file tree
Hide file tree
Showing 9 changed files with 64 additions and 68 deletions.
2 changes: 1 addition & 1 deletion covariance/covarianceBase.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -697,7 +697,7 @@ double covarianceBase::GetLikelihood() {
// Default behaviour is to reject negative values + do std llh calculation
const int NOutside = CheckBounds();

if(NOutside > 0) return NOutside*_LARGE_LOGL_;
if(NOutside > 0) return NOutside*M3::_LARGE_LOGL_;

return CalcLikelihood();
}
Expand Down
27 changes: 13 additions & 14 deletions manager/Core.h
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,19 @@ namespace M3 {
return std::fma(x, y, z);
#endif
}
constexpr static const double _BAD_DOUBLE_ = -999.99;
constexpr static const int _BAD_INT_ = -999;
constexpr static const double _DEFAULT_RETURN_VAL_ = -999999.123456;

/// Some commonly used variables to which we set pointers to
constexpr static const double Unity = 1.;
constexpr static const double Zero = 0.;
constexpr static const float Unity_F = 1.;
constexpr static const float Zero_F = 0.;
constexpr static const int Unity_Int = 1;

/// Large Likelihood is used it parameter go out of physical boundary, this indicates in MCMC that such step should be removed
constexpr static const double _LARGE_LOGL_ = 1234567890.0;
}

/// KS: noexcept can help with performance but is terrible for debugging, this is meant to help easy way of of turning it on or off. In near future move this to struct or other central class.
Expand All @@ -49,20 +62,6 @@ namespace M3 {
/// Number of overflow bins in TH2Poly,
#define _TH2PolyOverflowBins_ 9

constexpr static const double _BAD_DOUBLE_ = -999.99;
constexpr static const int _BAD_INT_ = -999;
constexpr static const double _DEFAULT_RETURN_VAL_ = -999999.123456;

/// Some commonly used variables to which we set pointers to
constexpr static const double Unity = 1.;
constexpr static const double Zero = 0.;
constexpr static const float Unity_F = 1.;
constexpr static const float Zero_F = 0.;
constexpr static const int Unity_Int = 1;

/// Large Likelihood is used it parameter go out of physical boundary, this indicates in MCMC that such step should be removed
constexpr static const double _LARGE_LOGL_ = 1234567890.0;

// C++ includes
#include <sstream>
#include <fstream>
Expand Down
18 changes: 9 additions & 9 deletions mcmc/FitterBase.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -274,7 +274,7 @@ void FitterBase::StartFromPreviousFit(const std::string& FitName) {
TFile *infile = new TFile(FitName.c_str(), "READ");
TTree *posts = infile->Get<TTree>("posteriors");
int step_val = 0;
double log_val = _LARGE_LOGL_;
double log_val = M3::_LARGE_LOGL_;
posts->SetBranchAddress("step",&step_val);
posts->SetBranchAddress("LogL",&log_val);

Expand Down Expand Up @@ -303,14 +303,14 @@ void FitterBase::StartFromPreviousFit(const std::string& FitName) {
CovarianceFolder->Close();
delete CovarianceFolder;

std::vector<double> branch_vals(systematics[s]->GetNumParams(), _BAD_DOUBLE_);
std::vector<double> branch_vals(systematics[s]->GetNumParams(), M3::_BAD_DOUBLE_);
for (int i = 0; i < systematics[s]->GetNumParams(); ++i) {
posts->SetBranchAddress(systematics[s]->GetParName(i).c_str(), &branch_vals[i]);
}
posts->GetEntry(posts->GetEntries()-1);

for (int i = 0; i < systematics[s]->GetNumParams(); ++i) {
if(branch_vals[i] == _BAD_DOUBLE_)
if(branch_vals[i] == M3::_BAD_DOUBLE_)
{
MACH3LOG_ERROR("Parameter {} is unvitalised with value {}", i, branch_vals[i]);
MACH3LOG_ERROR("Please check more precisely chain you passed {}", FitName);
Expand Down Expand Up @@ -396,23 +396,23 @@ void FitterBase::DragRace(const int NLaps) {
MACH3LOG_INFO("All tests will be performed with {} threads", M3::GetNThreads());

// Reweight the MC
for(unsigned int ivs = 0; ivs < samples.size(); ivs++ )
for(unsigned int ivs = 0; ivs < samples.size(); ++ivs)
{
TStopwatch clockRace;
clockRace.Start();
for(int Lap = 0; Lap < NLaps; Lap++) {
for(int Lap = 0; Lap < NLaps; ++Lap) {
samples[ivs]->reweight();
}
clockRace.Stop();
MACH3LOG_INFO("It took {:.4f} s to reweights {} times sample: {}", clockRace.RealTime(), NLaps, samples[ivs]->GetName());
MACH3LOG_INFO("On average {:.6f}", clockRace.RealTime()/NLaps);
}

for(unsigned int ivs = 0; ivs < samples.size(); ivs++ )
for(unsigned int ivs = 0; ivs < samples.size(); ++ivs)
{
TStopwatch clockRace;
clockRace.Start();
for(int Lap = 0; Lap < NLaps; Lap++) {
for(int Lap = 0; Lap < NLaps; ++Lap) {
samples[ivs]->GetLikelihood();
}
clockRace.Stop();
Expand All @@ -427,7 +427,7 @@ void FitterBase::DragRace(const int NLaps) {
for (size_t s = 0; s < systematics.size(); ++s) {
TStopwatch clockRace;
clockRace.Start();
for(int Lap = 0; Lap < NLaps; Lap++) {
for(int Lap = 0; Lap < NLaps; ++Lap) {
systematics[s]->proposeStep();
}
clockRace.Stop();
Expand All @@ -441,7 +441,7 @@ void FitterBase::DragRace(const int NLaps) {
for (size_t s = 0; s < systematics.size(); ++s) {
TStopwatch clockRace;
clockRace.Start();
for(int Lap = 0; Lap < NLaps; Lap++) {
for(int Lap = 0; Lap < NLaps; ++Lap) {
systematics[s]->GetLikelihood();
}
clockRace.Stop();
Expand Down
9 changes: 3 additions & 6 deletions mcmc/mcmc.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -76,8 +76,7 @@ void mcmc::CheckStep() {
// *******************
// Run the Markov chain with all the systematic objects added
void mcmc::runMCMC() {
// *******************

// *******************
// Save the settings into the output file
SaveSettings();

Expand Down Expand Up @@ -125,7 +124,6 @@ void mcmc::runMCMC() {
// Do the initial reconfigure of the MCMC
void mcmc::ProposeStep() {
// *******************

// Initial likelihood
double llh = 0.0;

Expand All @@ -149,7 +147,7 @@ void mcmc::ProposeStep() {

// Check if we've hit a boundary in the systematics
// In this case we can save time by not having to reconfigure the simulation
if (llh >= _LARGE_LOGL_) {
if (llh >= M3::_LARGE_LOGL_) {
reject = true;
#ifdef DEBUG
if (debug) debugFile << "Rejecting based on boundary" << std::endl;
Expand Down Expand Up @@ -181,7 +179,7 @@ void mcmc::ProposeStep() {
} else {
for (size_t i = 0; i < samples.size(); ++i) {
// Set the sample_llh[i] to be madly high also to signify a step out of bounds
sample_llh[i] = _LARGE_LOGL_;
sample_llh[i] = M3::_LARGE_LOGL_;
#ifdef DEBUG
if (debug) debugFile << "LLH after REJECT sample " << i << " " << llh << std::endl;
#endif
Expand All @@ -196,7 +194,6 @@ void mcmc::ProposeStep() {
// Print the fit output progress
void mcmc::PrintProgress() {
// *******************

MACH3LOG_INFO("Step:\t{}/{}, current: {:.2f}, proposed: {:.2f}", step - stepStart, chainLength, logLCurr, logLProp);
MACH3LOG_INFO("Accepted/Total steps: {}/{} = {:.2f}", accCount, step - stepStart, static_cast<double>(accCount) / static_cast<double>(step - stepStart));

Expand Down
10 changes: 5 additions & 5 deletions plotting/plottingUtils/inputManager.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -130,9 +130,9 @@ double InputManager::getPostFitError(int fileNum, const std::string &paramName,
return inputFileDef.postFitErrors.at(errorType).at(paramName);
}

MACH3LOG_WARN("Didn't fnd {} post fit error for {}. Returning {}", errorType, paramName, _BAD_DOUBLE_);
MACH3LOG_WARN("Didn't fnd {} post fit error for {}. Returning {}", errorType, paramName, M3::_BAD_DOUBLE_);

return _BAD_DOUBLE_;
return M3::_BAD_DOUBLE_;
}

double InputManager::getPostFitValue(int fileNum, const std::string &paramName,
Expand All @@ -159,9 +159,9 @@ double InputManager::getPostFitValue(int fileNum, const std::string &paramName,
return inputFileDef.postFitValues.at(errorType).at(paramName);
}

MACH3LOG_WARN("Didn't fnd {} post fit value for {}. Returning {}", errorType, paramName, _BAD_DOUBLE_);
MACH3LOG_WARN("Didn't fnd {} post fit value for {}. Returning {}", errorType, paramName, M3::_BAD_DOUBLE_);

return _BAD_DOUBLE_;
return M3::_BAD_DOUBLE_;
}


Expand Down Expand Up @@ -419,7 +419,7 @@ bool InputManager::findRawChainSteps(InputFile &inputFileDef, const std::string
if ( setInputBranch )
{
// EM: should probably use MCMCProcessor for this so we can use caching, gpu etc.
inputFileDef.MCMCstepParamsMap[parameter] = new double( _BAD_DOUBLE_ ); // <- initialise the parameter step values
inputFileDef.MCMCstepParamsMap[parameter] = new double( M3::_BAD_DOUBLE_ ); // <- initialise the parameter step values
inputFileDef.posteriorTree->SetBranchAddress( branchNames[paramIdx], inputFileDef.MCMCstepParamsMap.at(parameter) );
}
break;
Expand Down
2 changes: 1 addition & 1 deletion plotting/plottingUtils/inputManager.h
Original file line number Diff line number Diff line change
Expand Up @@ -320,7 +320,7 @@ class InputManager {
{
MACH3LOG_WARN("file at index {} does not have an MCMC entry for parameter {}", fileNum, paramName);
MACH3LOG_WARN("am returning a bad float");
return _BAD_DOUBLE_;
return M3::_BAD_DOUBLE_;
}

return *_fileVec[fileNum].MCMCstepParamsMap.at(paramName);
Expand Down
2 changes: 1 addition & 1 deletion samplePDF/Structs.h
Original file line number Diff line number Diff line change
Expand Up @@ -521,7 +521,7 @@ namespace MaCh3Utils {
/// beware that in the case of anti-neutrinos the NuOscillator
/// type simply gets multiplied by -1
inline int PDGToNuOscillatorFlavour(int NuPdg){
int NuOscillatorFlavour = _BAD_INT_;
int NuOscillatorFlavour = M3::_BAD_INT_;
switch(std::abs(NuPdg)){
case NuPDG::kNue:
NuOscillatorFlavour = NuOscillator::kElectron;
Expand Down
58 changes: 29 additions & 29 deletions samplePDF/samplePDFFDBase.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -976,11 +976,11 @@ void samplePDFFDBase::FindNominalBinAndEdges1D() {

//Set x_var and y_var values based on XVarStr and YVarStr
MCSamples[mc_i].x_var[event_i] = GetPointerToKinematicParameter(XVarStr, mc_i, event_i);
//Give y_var _BAD_DOUBLE_ value for the 1D case since this won't be used
MCSamples[mc_i].y_var[event_i] = &(_BAD_DOUBLE_);
//Give y_var M3::_BAD_DOUBLE_ value for the 1D case since this won't be used
MCSamples[mc_i].y_var[event_i] = &(M3::_BAD_DOUBLE_);
int bin = _hPDF1D->FindBin(*(MCSamples[mc_i].x_var[event_i]));

double low_lower_edge = _DEFAULT_RETURN_VAL_;
double low_lower_edge = M3::_DEFAULT_RETURN_VAL_;
if (bin==0) {
low_lower_edge = _hPDF1D->GetXaxis()->GetBinLowEdge(bin);
} else {
Expand All @@ -990,7 +990,7 @@ void samplePDFFDBase::FindNominalBinAndEdges1D() {
double low_edge = _hPDF1D->GetXaxis()->GetBinLowEdge(bin);
double upper_edge = _hPDF1D->GetXaxis()->GetBinUpEdge(bin);

double upper_upper_edge = _DEFAULT_RETURN_VAL_;
double upper_upper_edge = M3::_DEFAULT_RETURN_VAL_;
if (bin<(_hPDF1D->GetNbinsX()-2)) {
upper_upper_edge = _hPDF1D->GetXaxis()->GetBinLowEdge(bin+2);
} else {
Expand All @@ -1001,10 +1001,10 @@ void samplePDFFDBase::FindNominalBinAndEdges1D() {
MCSamples[mc_i].NomXBin[event_i] = bin-1;
} else {
MCSamples[mc_i].NomXBin[event_i] = -1;
low_edge = _DEFAULT_RETURN_VAL_;
upper_edge = _DEFAULT_RETURN_VAL_;
low_lower_edge = _DEFAULT_RETURN_VAL_;
upper_upper_edge = _DEFAULT_RETURN_VAL_;
low_edge = M3::_DEFAULT_RETURN_VAL_;
upper_edge = M3::_DEFAULT_RETURN_VAL_;
low_lower_edge = M3::_DEFAULT_RETURN_VAL_;
upper_upper_edge = M3::_DEFAULT_RETURN_VAL_;
}
MCSamples[mc_i].NomYBin[event_i] = 0;

Expand Down Expand Up @@ -1108,7 +1108,7 @@ void samplePDFFDBase::FindNominalBinAndEdges2D() {
_hPDF2D->GetBinXYZ(bin, bin_x, bin_y, bin_z);
//erec is the x-axis so get GetXaxis then find the bin edges using the x bin number

double low_lower_edge = _DEFAULT_RETURN_VAL_;
double low_lower_edge = M3::_DEFAULT_RETURN_VAL_;
if (bin==0) {
low_lower_edge = _hPDF2D->GetXaxis()->GetBinLowEdge(bin_x);
} else {
Expand All @@ -1118,7 +1118,7 @@ void samplePDFFDBase::FindNominalBinAndEdges2D() {
double low_edge = _hPDF2D->GetXaxis()->GetBinLowEdge(bin_x);
double upper_edge = _hPDF2D->GetXaxis()->GetBinUpEdge(bin_x);

double upper_upper_edge = _DEFAULT_RETURN_VAL_;
double upper_upper_edge = M3::_DEFAULT_RETURN_VAL_;
if (bin<(_hPDF2D->GetNbinsX()-2)) {
upper_upper_edge = _hPDF2D->GetXaxis()->GetBinLowEdge(bin_x+2);
} else {
Expand All @@ -1129,10 +1129,10 @@ void samplePDFFDBase::FindNominalBinAndEdges2D() {
MCSamples[mc_i].NomXBin[event_i] = bin_x-1;
} else {
MCSamples[mc_i].NomXBin[event_i] = -1;
low_edge = _DEFAULT_RETURN_VAL_;
upper_edge = _DEFAULT_RETURN_VAL_;
low_lower_edge = _DEFAULT_RETURN_VAL_;
upper_upper_edge = _DEFAULT_RETURN_VAL_;
low_edge = M3::_DEFAULT_RETURN_VAL_;
upper_edge = M3::_DEFAULT_RETURN_VAL_;
low_lower_edge = M3::_DEFAULT_RETURN_VAL_;
upper_upper_edge = M3::_DEFAULT_RETURN_VAL_;
}
MCSamples[mc_i].NomYBin[event_i] = bin_y-1;
if(MCSamples[mc_i].NomYBin[event_i] < 0){
Expand Down Expand Up @@ -1313,32 +1313,32 @@ void samplePDFFDBase::SetupNuOscillator() {
for (int iEvent=0;iEvent<MCSamples[iSample].nEvents;iEvent++) {
// KS: Sry but if we use low memory we need to point to float not double...
#ifdef _LOW_MEMORY_STRUCTS_
MCSamples[iSample].osc_w_pointer[iEvent] = &Unity_F;
MCSamples[iSample].osc_w_pointer[iEvent] = &M3::Unity_F;
#else
MCSamples[iSample].osc_w_pointer[iEvent] = &Unity;
MCSamples[iSample].osc_w_pointer[iEvent] = &M3::Unity;
#endif
if (MCSamples[iSample].isNC[iEvent]) {
if (*MCSamples[iSample].nupdg[iEvent] != *MCSamples[iSample].nupdgUnosc[iEvent]) {
#ifdef _LOW_MEMORY_STRUCTS_
MCSamples[iSample].osc_w_pointer[iEvent] = &Zero_F;
MCSamples[iSample].osc_w_pointer[iEvent] = &M3::Zero_F;
#else
MCSamples[iSample].osc_w_pointer[iEvent] = &Zero;
MCSamples[iSample].osc_w_pointer[iEvent] = &M3::Zero;
#endif
} else {
#ifdef _LOW_MEMORY_STRUCTS_
MCSamples[iSample].osc_w_pointer[iEvent] = &Unity_F;
MCSamples[iSample].osc_w_pointer[iEvent] = &M3::Unity_F;
#else
MCSamples[iSample].osc_w_pointer[iEvent] = &Unity;
MCSamples[iSample].osc_w_pointer[iEvent] = &M3::Unity;
#endif
}
} else {
int InitFlav = _BAD_INT_;
int FinalFlav = _BAD_INT_;
int InitFlav = M3::_BAD_INT_;
int FinalFlav = M3::_BAD_INT_;

InitFlav = MaCh3Utils::PDGToNuOscillatorFlavour((*MCSamples[iSample].nupdgUnosc[iEvent]));
FinalFlav = MaCh3Utils::PDGToNuOscillatorFlavour((*MCSamples[iSample].nupdg[iEvent]));

if (InitFlav == _BAD_INT_ || FinalFlav == _BAD_INT_) {
if (InitFlav == M3::_BAD_INT_ || FinalFlav == M3::_BAD_INT_) {
MACH3LOG_ERROR("Something has gone wrong in the mapping between MCSamples[iSample].nutype and the enum used within NuOscillator");
MACH3LOG_ERROR("MCSamples[iSample].nupdgUnosc: {}", (*MCSamples[iSample].nupdgUnosc[iEvent]));
MACH3LOG_ERROR("InitFlav: {}", InitFlav);
Expand Down Expand Up @@ -1461,9 +1461,9 @@ void samplePDFFDBase::InitialiseSingleFDMCObject(int iSample, int nEvents_) {
fdobj->ChannelIndex = iSample;

int nEvents = fdobj->nEvents;
fdobj->x_var.resize(nEvents, &Unity);
fdobj->y_var.resize(nEvents, &Unity);
fdobj->rw_etru.resize(nEvents, &Unity);
fdobj->x_var.resize(nEvents, &M3::Unity);
fdobj->y_var.resize(nEvents, &M3::Unity);
fdobj->rw_etru.resize(nEvents, &M3::Unity);
fdobj->XBin.resize(nEvents, -1);
fdobj->YBin.resize(nEvents, -1);
fdobj->NomXBin.resize(nEvents, -1);
Expand All @@ -1472,7 +1472,7 @@ void samplePDFFDBase::InitialiseSingleFDMCObject(int iSample, int nEvents_) {
fdobj->rw_lower_lower_xbinedge.resize(nEvents, -1);
fdobj->rw_upper_xbinedge.resize(nEvents, -1);
fdobj->rw_upper_upper_xbinedge.resize(nEvents, -1);
fdobj->mode.resize(nEvents, &Unity);
fdobj->mode.resize(nEvents, &M3::Unity);
fdobj->nxsec_norm_pointers.resize(nEvents);
fdobj->xsec_norm_pointers.resize(nEvents);
fdobj->xsec_norms_bins.resize(nEvents);
Expand All @@ -1486,9 +1486,9 @@ void samplePDFFDBase::InitialiseSingleFDMCObject(int iSample, int nEvents_) {
fdobj->total_weight_pointers.resize(nEvents);
fdobj->Target.resize(nEvents, 0);
#ifdef _LOW_MEMORY_STRUCTS_
fdobj->osc_w_pointer.resize(nEvents, &Unity_F);
fdobj->osc_w_pointer.resize(nEvents, &M3::Unity_F);
#else
fdobj->osc_w_pointer.resize(nEvents, &Unity);
fdobj->osc_w_pointer.resize(nEvents, &M3::Unity);
#endif

for(int iEvent = 0 ; iEvent < fdobj->nEvents ; ++iEvent){
Expand Down
4 changes: 2 additions & 2 deletions samplePDF/samplePDFFDBase.h
Original file line number Diff line number Diff line change
Expand Up @@ -211,7 +211,7 @@ class samplePDFFDBase : public samplePDFBase
//===============================================================================

/// @brief Keep track of the dimensions of the sample binning
int nDimensions = _BAD_INT_;
int nDimensions = M3::_BAD_INT_;
/// @brief A unique ID for each sample based on powers of two for quick binary operator comparisons
std::string SampleDetID;
/// holds "TrueNeutrinoEnergy" and the strings used for the sample binning.
Expand All @@ -230,7 +230,7 @@ class samplePDFFDBase : public samplePDFBase
//What gets used in IsEventSelected, which gets set equal to user input plus
//all the vectors in StoreSelection
/// @brief the Number of selections in the
int NSelections = _BAD_INT_;
int NSelections = M3::_BAD_INT_;

/// @brief What gets pulled from config options, these are constant after loading in
/// this is of length 3: 0th index is the value, 1st is lower bound, 2nd is upper bound
Expand Down

0 comments on commit 96d89d2

Please sign in to comment.