Skip to content

Commit

Permalink
Merge branch 'master' into feature_barlowbeeston
Browse files Browse the repository at this point in the history
  • Loading branch information
@willp240
willp240 authored Nov 22, 2023
2 parents 3fd1766 + 94a91e3 commit fb8265f
Show file tree
Hide file tree
Showing 17 changed files with 746 additions and 121 deletions.
4 changes: 2 additions & 2 deletions src/core/SparseMatrix.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,12 +11,12 @@ SparseMatrix::SparseMatrix(size_t rows_, size_t cols_){
}

void
SparseMatrix::PrintDense(const std::string& prefix_=""){
SparseMatrix::PrintDense(const std::string& prefix_="") const{
fArmaMat.print(prefix_);
}

void
SparseMatrix::Print(const std::string& prefix_=""){
SparseMatrix::Print(const std::string& prefix_="") const{
arma::mat B(fArmaMat);
B.print(prefix_);

Expand Down
4 changes: 2 additions & 2 deletions src/core/SparseMatrix.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,8 +32,8 @@ class SparseMatrix{
void SetZeros();
void SetToIdentity();

void Print(const std::string&);
void PrintDense(const std::string&);
void Print(const std::string&) const;
void PrintDense(const std::string&) const;

private:
arma::sp_mat fArmaMat;
Expand Down
133 changes: 117 additions & 16 deletions src/fitutil/BinnedEDManager.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,22 +16,45 @@ BinnedEDManager::GetNDims() const{
}

double
BinnedEDManager::Probability(const Event& data_) const{
BinnedEDManager::Probability(const Event& data_) {
ReassertNorms(true);
double sum = 0;

size_t j = 0;
for(size_t i = 0; i < fWorkingPdfs.size(); i++){
sum += fNormalisations.at(i) * fWorkingPdfs[i].Probability(data_);
if (fAllowNormsFittable.at(i) == DIRECT) {
sum += fNormalisations.at(i) * fWorkingPdfs[i].Probability(data_);
j++;
} else if (fAllowNormsFittable.at(i) == INDIRECT) {
sum += fFittableNorms.at(j) * fWorkingPdfs.at(i).Probability(data_);
j++;
} else {
// In case where we don't auto-normalise the pdf,
// the normalisation is held within the pdf itself!
sum += fWorkingPdfs.at(i).Probability(data_);
}
}

return sum;
}

double
BinnedEDManager::BinProbability(size_t bin_) const{
BinnedEDManager::BinProbability(size_t bin_) {
ReassertNorms(true);
double sum = 0;
size_t j = 0;
try{
for(size_t i = 0; i < fWorkingPdfs.size(); i++){
sum += fNormalisations.at(i) * fWorkingPdfs.at(i).GetBinContent(bin_);
if (fAllowNormsFittable.at(i) == DIRECT) {
sum += fNormalisations.at(i) * fWorkingPdfs.at(i).GetBinContent(bin_);
j++;
} else if (fAllowNormsFittable.at(i) == INDIRECT) {
sum += fFittableNorms.at(j) * fWorkingPdfs.at(i).GetBinContent(bin_);
j++;
} else {
// In case where we don't auto-normalise the pdf,
// the normalisation is held within the pdf itself!
sum += fWorkingPdfs.at(i).GetBinContent(bin_);
}
}
}
catch(const std::out_of_range&){
Expand All @@ -55,8 +78,19 @@ BinnedEDManager::ApplySystematics(const SystematicManager& sysMan_){

if(!sysMan_.GetNSystematics())
return;

sysMan_.DistortEDs(fOriginalPdfs,fWorkingPdfs);
if (fAllNormsDirFittable) {
sysMan_.DistortEDs(fOriginalPdfs,fWorkingPdfs);
} else {
std::vector<double> norms(fNormalisations.size(), 0);
sysMan_.DistortEDs(fOriginalPdfs,fWorkingPdfs, &norms);
for (size_t i = 0; i < norms.size(); i++) {
if (fAllowNormsFittable.at(i) == FALSE) {
fNormalisations[i] = norms.at(i);
} else if (fAllowNormsFittable.at(i) == INDIRECT) {
fNormalisations[i] *= norms.at(i);
}
}
}
}

const BinnedED&
Expand All @@ -65,18 +99,30 @@ BinnedEDManager::GetOriginalPdf(size_t index_) const{
}

void
BinnedEDManager::AddPdf(const BinnedED& pdf_){
BinnedEDManager::AddPdf(const BinnedED& pdf_, const NormFittingStatus norm_fitting_status){
fOriginalPdfs.push_back(pdf_);
fWorkingPdfs.push_back(pdf_);
fNPdfs++;
fNormalisations.resize(fOriginalPdfs.size(), 0);
fAllowNormsFittable.push_back(norm_fitting_status);
if (norm_fitting_status) { fFittableNorms.push_back(0); }
if (norm_fitting_status != DIRECT) { fAllNormsDirFittable = false; }
RegisterParameters();
}

void
BinnedEDManager::AddPdfs(const std::vector<BinnedED>& pdfs_){
BinnedEDManager::AddPdfs(const std::vector<BinnedED>& pdfs_,
const std::vector<NormFittingStatus>* norm_fitting_statuses){
if(norm_fitting_statuses != nullptr && pdfs_.size() != norm_fitting_statuses->size()) {
throw DimensionError("BinnedEDManager: number of norm_fittable bools doesn't the number of pdfs");
}

for(size_t i = 0; i < pdfs_.size(); i++){
AddPdf(pdfs_.at(i));
if (norm_fitting_statuses != nullptr) {
AddPdf(pdfs_.at(i), norm_fitting_statuses->at(i));
} else {
AddPdf(pdfs_.at(i));
}
}
RegisterParameters();
}
Expand All @@ -92,12 +138,25 @@ BinnedEDManager::ApplyShrink(const BinnedEDShrinker& shrinker_){
return;

// only shrink if not already shrunk! FIXME: more obvious behaviour
if (!fWorkingPdfs.size() || fWorkingPdfs.at(0).GetNBins() != fOriginalPdfs.at(0).GetNBins())
if (!fWorkingPdfs.size())
return;

for (size_t i = 0; i < fWorkingPdfs.size(); i++){
fWorkingPdfs[i] = shrinker_.ShrinkDist(fWorkingPdfs.at(i));
fWorkingPdfs[i].Normalise();
// Normalise if normalisation is a fittable param, but if indirect then track any change
if (fAllowNormsFittable.at(i) == DIRECT) {
fWorkingPdfs[i] = shrinker_.ShrinkDist(fWorkingPdfs.at(i));
fWorkingPdfs[i].Normalise();
} else if (fAllowNormsFittable.at(i) == FALSE) {
fWorkingPdfs[i] = shrinker_.ShrinkDist(fWorkingPdfs.at(i));
const double integral_after = fWorkingPdfs[i].Integral();
fNormalisations[i] = integral_after;
} else {
const double integral_before = fWorkingPdfs[i].Integral();
fWorkingPdfs[i] = shrinker_.ShrinkDist(fWorkingPdfs.at(i));
const double integral_after = fWorkingPdfs[i].Integral();
if (integral_before == 0. && integral_after == 0.) { fNormalisations[i] = 0.; }
else { fNormalisations[i] *= integral_after/integral_before; }
}
}

}
Expand Down Expand Up @@ -155,7 +214,49 @@ void
BinnedEDManager::RegisterParameters(){
fParameterManager.Clear();
std::vector<std::string> parameterNames;
for(size_t i = 0; i < fOriginalPdfs.size(); i++)
parameterNames.push_back(fOriginalPdfs.at(i).GetName());
fParameterManager.AddContainer(fNormalisations, parameterNames);
if (fAllNormsDirFittable) {
for(size_t i = 0; i < fOriginalPdfs.size(); i++) {
parameterNames.push_back(fOriginalPdfs.at(i).GetName());
}
fParameterManager.AddContainer(fNormalisations, parameterNames);
} else {
for(size_t i = 0; i < fOriginalPdfs.size(); i++) {
if (fAllowNormsFittable.at(i)) {
parameterNames.push_back(fOriginalPdfs.at(i).GetName());
}
}
fParameterManager.AddContainer(fFittableNorms, parameterNames);
}
}

void BinnedEDManager::ReassertNorms(bool calcing_binprob) {
/*
* Normalisations that are fittable might have been changed
* in fFittableNorms, but not in fNormalisations!
* This method corrects any changes that may have occured,
* so that both vectors are consistent.
*/
if (!fAllNormsDirFittable) {
size_t j = 0;
for (size_t i = 0; i < fNormalisations.size(); i++) {
if (fAllowNormsFittable.at(i) == ((calcing_binprob) ? DIRECT : INDIRECT)) {
fNormalisations[i] = fFittableNorms.at(j);
j++;
}
}
}
}

void BinnedEDManager::AssertDimensions(const std::vector<std::string>& observables) {
/*
* Check whether each PDF in the manager actually has dimensions
* matching that of the input vector.
* Whenever this is not the case, the given PDF is marginalised onto
* the input vector's axes.
*/
for (auto &working_pdf : fWorkingPdfs) {
if (working_pdf.GetObservables() != observables) {
working_pdf = working_pdf.Marginalise(observables);
}
}
}
28 changes: 18 additions & 10 deletions src/fitutil/BinnedEDManager.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,15 +12,18 @@
class Event;
class SystematicManager;
class BinnedEDShrinker;

enum NormFittingStatus{ FALSE=0, DIRECT=1, INDIRECT=2 };
class BinnedEDManager : public FitComponent{
public:
BinnedEDManager() : fNPdfs(0), fName("norms") {}
BinnedEDManager() : fAllNormsDirFittable(true), fNPdfs(0), fName("norms") {}

void AddPdf(const BinnedED&);
void AddPdfs(const std::vector<BinnedED>&);
void AddPdf(const BinnedED&, const NormFittingStatus norm_fitting_status=DIRECT);
void AddPdfs(const std::vector<BinnedED>&,
const std::vector<NormFittingStatus>* norm_fitting_statuses=nullptr);

double Probability(const Event&) const;
double BinProbability(size_t) const;
double Probability(const Event&);
double BinProbability(size_t);

const std::vector<double>& GetNormalisations() const;
void SetNormalisations(const std::vector<double>& normalisations_);
Expand All @@ -32,6 +35,8 @@ class BinnedEDManager : public FitComponent{
unsigned GetNPdfs() const;
size_t GetNDims() const;

void AssertDimensions(const std::vector<std::string>& observables);
void ReassertNorms(bool calcing_binprob=false);

// Make a fittable component - i.e. rescale the binned pdfs inside to fit
void SetParameter(const std::string& name_, double value);
Expand All @@ -49,11 +54,14 @@ class BinnedEDManager : public FitComponent{


private:
ParameterManager fParameterManager;
std::vector<BinnedED> fOriginalPdfs;
std::vector<BinnedED> fWorkingPdfs;
std::vector<double> fNormalisations;
int fNPdfs;
ParameterManager fParameterManager;
std::vector<BinnedED> fOriginalPdfs;
std::vector<BinnedED> fWorkingPdfs;
std::vector<double> fNormalisations;
std::vector<NormFittingStatus> fAllowNormsFittable;
std::vector<double> fFittableNorms;
bool fAllNormsDirFittable;
int fNPdfs;
size_t fNDims;

std::string fName; // component name
Expand Down
7 changes: 6 additions & 1 deletion src/fitutil/SystematicManager.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -141,7 +141,9 @@ SystematicManager::AddDist(const BinnedED& pdf, const std::string& syss_){
}

void
SystematicManager::DistortEDs(std::vector<BinnedED>& OrignalEDs_,std::vector<BinnedED>& WorkingEDs_) const {
SystematicManager::DistortEDs(const std::vector<BinnedED>& OrignalEDs_,
std::vector<BinnedED>& WorkingEDs_,
std::vector<double>* norms) const {
for(size_t j = 0; j < WorkingEDs_.size(); j++){
const std::string name = OrignalEDs_.at(j).GetName();

Expand Down Expand Up @@ -186,6 +188,9 @@ SystematicManager::DistortEDs(std::vector<BinnedED>& OrignalEDs_,std::vector<Bin
GetTotalResponse(groupName).operator()(OrignalEDs_.at(j).GetBinContents()));
}
}
if (norms != nullptr) {
norms->operator[](j) = WorkingEDs_.at(j).Integral();
}
}
}

Expand Down
4 changes: 3 additions & 1 deletion src/fitutil/SystematicManager.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,9 @@ class SystematicManager{

const SparseMatrix& GetTotalResponse(const std::string& groupName_ = "" ) const;

void DistortEDs(std::vector<BinnedED>& OrigEDs,std::vector<BinnedED>& WorkingEDs) const;
void DistortEDs(const std::vector<BinnedED>& OrigEDs,
std::vector<BinnedED>& WorkingEDs,
std::vector<double>* norms=nullptr) const;

void Construct();

Expand Down
Loading

0 comments on commit fb8265f

Please sign in to comment.