LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
LArDiscreteProbabilityVector.cc
Go to the documentation of this file.
1 
10 
11 #include <algorithm>
12 #include <cmath>
13 #include <numeric>
14 
15 namespace lar_content
16 {
17 
18 template <typename TX, typename TY>
19 DiscreteProbabilityVector::DiscreteProbabilityVector(const InputData<TX, TY> &inputData, const TX xUpperBound, const bool useWidths) :
20  m_xUpperBound(static_cast<float>(xUpperBound)),
21  m_useWidths(useWidths),
22  m_discreteProbabilityData(this->InitialiseDiscreteProbabilityData(inputData))
23 {
24  this->VerifyCompleteData();
25 }
26 
27 //------------------------------------------------------------------------------------------------------------------------------------------
28 
29 DiscreteProbabilityVector::DiscreteProbabilityVector(const DiscreteProbabilityVector &discreteProbabilityVector, std::mt19937 &randomNumberGenerator) :
30  m_xUpperBound(discreteProbabilityVector.m_xUpperBound),
31  m_useWidths(discreteProbabilityVector.m_useWidths),
32  m_discreteProbabilityData(this->RandomiseDiscreteProbabilityData(discreteProbabilityVector, randomNumberGenerator))
33 {
34  this->VerifyCompleteData();
35 }
36 
37 //------------------------------------------------------------------------------------------------------------------------------------------
38 
39 DiscreteProbabilityVector::DiscreteProbabilityVector(const DiscreteProbabilityVector &discreteProbabilityVector, const ResamplingPoints &resamplingPoints) :
40  m_xUpperBound(discreteProbabilityVector.m_xUpperBound),
41  m_useWidths(discreteProbabilityVector.m_useWidths),
42  m_discreteProbabilityData(this->ResampleDiscreteProbabilityData(discreteProbabilityVector, resamplingPoints))
43 {
44  this->VerifyCompleteData();
45 }
46 
47 //------------------------------------------------------------------------------------------------------------------------------------------
48 
50 {
51  if (x - m_discreteProbabilityData.back().GetX() > std::numeric_limits<float>::epsilon())
52  return 1.f;
53 
54  if (x - m_discreteProbabilityData.front().GetX() < std::numeric_limits<float>::epsilon())
55  return 0.f;
56 
57  for (unsigned int iDatum = 1; iDatum < m_discreteProbabilityData.size(); ++iDatum)
58  {
59  if (x - m_discreteProbabilityData.at(iDatum).GetX() > std::numeric_limits<float>::epsilon())
60  continue;
61 
62  const float xLow(m_discreteProbabilityData.at(iDatum - 1).GetX());
63  const float yLow(m_discreteProbabilityData.at(iDatum - 1).GetCumulativeDatum());
64  const float xHigh(m_discreteProbabilityData.at(iDatum).GetX());
65  const float yHigh(m_discreteProbabilityData.at(iDatum).GetCumulativeDatum());
66 
67  if (std::fabs(xHigh - xLow) < std::numeric_limits<float>::epsilon())
68  throw pandora::StatusCodeException(pandora::STATUS_CODE_INVALID_PARAMETER);
69 
70  const float m((yHigh - yLow) / (xHigh - xLow));
71  const float c(yLow - m * xLow);
72 
73  return m * x + c;
74  }
75 
76  throw pandora::StatusCodeException(pandora::STATUS_CODE_INVALID_PARAMETER);
77 }
78 
79 //------------------------------------------------------------------------------------------------------------------------------------------
80 
81 template <typename TX, typename TY>
83 {
84  if (2 > inputData.size())
85  throw pandora::StatusCodeException(pandora::STATUS_CODE_INVALID_PARAMETER);
86 
87  std::sort(inputData.begin(), inputData.end(), DiscreteProbabilityVector::SortInputDataByX<TX, TY>);
88 
89  if (inputData.back().first - m_xUpperBound > std::numeric_limits<float>::epsilon())
90  throw pandora::StatusCodeException(pandora::STATUS_CODE_INVALID_PARAMETER);
91 
92  const float normalisation(this->CalculateNormalisation(inputData));
93 
94  if (normalisation < std::numeric_limits<float>::epsilon())
95  throw pandora::StatusCodeException(pandora::STATUS_CODE_FAILURE);
96 
97  float accumulationDatum(0.f);
98 
100  for (unsigned int iDatum = 0; iDatum < inputData.size() - 1; ++iDatum)
101  {
102  const float x(static_cast<float>(inputData.at(iDatum).first));
103  const float deltaX(static_cast<float>(inputData.at(iDatum + 1).first) - x);
104  const float densityDatum(static_cast<float>(inputData.at(iDatum).second) / normalisation);
105  accumulationDatum += densityDatum * (m_useWidths ? deltaX : 1.f);
106  data.emplace_back(DiscreteProbabilityVector::DiscreteProbabilityDatum(x, densityDatum, accumulationDatum, deltaX));
107  }
108  const float x(static_cast<float>(inputData.back().first));
109  const float deltaX(m_xUpperBound - x);
110  const float densityDatum(static_cast<float>(inputData.back().second) / normalisation);
111  accumulationDatum += densityDatum * (m_useWidths ? deltaX : 1.f);
112  data.emplace_back(DiscreteProbabilityVector::DiscreteProbabilityDatum(x, densityDatum, accumulationDatum, deltaX));
113 
114  return data;
115 }
116 
117 //------------------------------------------------------------------------------------------------------------------------------------------
118 
120  const DiscreteProbabilityVector &discreteProbabilityVector, const ResamplingPoints &resamplingPoints) const
121 {
122  if (2 > resamplingPoints.size())
123  throw pandora::StatusCodeException(pandora::STATUS_CODE_INVALID_PARAMETER);
124 
125  DiscreteProbabilityData resampledProbabilityData;
126 
127  float prevCumulativeData(0.f);
128  for (unsigned int iSample = 0; iSample < resamplingPoints.size() - 1; ++iSample)
129  {
130  const float xResampled(resamplingPoints.at(iSample));
131  const float deltaX(resamplingPoints.at(iSample + 1) - xResampled);
132 
133  if (deltaX < std::numeric_limits<float>::epsilon())
134  throw pandora::StatusCodeException(pandora::STATUS_CODE_INVALID_PARAMETER);
135 
136  const float cumulativeDatumResampled(discreteProbabilityVector.EvaluateCumulativeProbability(xResampled));
137  const float densityDatumResampled((cumulativeDatumResampled - prevCumulativeData) / (m_useWidths ? deltaX : 1.f));
138  resampledProbabilityData.emplace_back(
139  DiscreteProbabilityVector::DiscreteProbabilityDatum(xResampled, densityDatumResampled, cumulativeDatumResampled, deltaX));
140  prevCumulativeData = cumulativeDatumResampled;
141  }
142 
143  const float xResampled(resamplingPoints.back());
144  const float deltaX(m_xUpperBound - xResampled);
145 
146  if (deltaX < std::numeric_limits<float>::epsilon())
147  throw pandora::StatusCodeException(pandora::STATUS_CODE_INVALID_PARAMETER);
148 
149  const float cumulativeDatumResampled(discreteProbabilityVector.EvaluateCumulativeProbability(xResampled));
150  const float densityDatumResampled((cumulativeDatumResampled - prevCumulativeData) / (m_useWidths ? deltaX : 1.f));
151  resampledProbabilityData.emplace_back(
152  DiscreteProbabilityVector::DiscreteProbabilityDatum(xResampled, densityDatumResampled, cumulativeDatumResampled, deltaX));
153 
154  return resampledProbabilityData;
155 }
156 
157 //------------------------------------------------------------------------------------------------------------------------------------------
158 
160  const DiscreteProbabilityVector &discreteProbabilityVector, std::mt19937 &randomNumberGenerator) const
161 {
162  DiscreteProbabilityData randomisedProbabilityData;
163 
164  std::vector<unsigned int> randomisedElements(discreteProbabilityVector.GetSize());
165  std::iota(std::begin(randomisedElements), std::end(randomisedElements), 0);
166  std::shuffle(std::begin(randomisedElements), std::end(randomisedElements), randomNumberGenerator);
167 
168  float xPos(discreteProbabilityVector.GetX(0));
169  float cumulativeProbability(0.f);
170  for (unsigned int iElement = 0; iElement < discreteProbabilityVector.GetSize(); ++iElement)
171  {
172  const unsigned int randomElementIndex(randomisedElements.at(iElement));
173  const float deltaX(discreteProbabilityVector.GetWidth(randomElementIndex));
174  const float probabilityDensity(discreteProbabilityVector.GetProbabilityDensity(randomElementIndex));
175  cumulativeProbability += probabilityDensity * (m_useWidths ? deltaX : 1.f);
176  randomisedProbabilityData.emplace_back(
177  DiscreteProbabilityVector::DiscreteProbabilityDatum(xPos, probabilityDensity, cumulativeProbability, deltaX));
178  xPos += deltaX;
179  }
180 
181  return randomisedProbabilityData;
182 }
183 
184 //------------------------------------------------------------------------------------------------------------------------------------------
185 
186 template <typename TX, typename TY>
188 {
189  const float deltaX(static_cast<float>(rhs.first) - static_cast<float>(lhs.first));
190  if (std::fabs(deltaX) < std::numeric_limits<float>::epsilon())
191  return (lhs.second < rhs.second);
192 
193  return (lhs.first < rhs.first);
194 }
195 
196 //------------------------------------------------------------------------------------------------------------------------------------------
197 
198 template <typename TX, typename TY>
200 {
201  if (2 > inputData.size())
202  throw pandora::StatusCodeException(pandora::STATUS_CODE_INVALID_PARAMETER);
203 
204  float normalisation(0.f);
205 
206  for (unsigned int iDatum = 0; iDatum < inputData.size() - 1; ++iDatum)
207  {
208  const float y(static_cast<float>(inputData.at(iDatum).second));
209  normalisation +=
210  y * (m_useWidths ? static_cast<float>(inputData.at(iDatum + 1).first) - static_cast<float>(inputData.at(iDatum).first) : 1.f);
211  }
212  const float y(static_cast<float>(inputData.back().second));
213  normalisation += y * (m_useWidths ? m_xUpperBound - static_cast<float>(inputData.back().first) : 1.f);
214 
215  return normalisation;
216 }
217 
218 //------------------------------------------------------------------------------------------------------------------------------------------
219 
220 template DiscreteProbabilityVector::DiscreteProbabilityVector(const InputData<int, float> &, int const, bool const);
221 template DiscreteProbabilityVector::DiscreteProbabilityVector(const InputData<float, int> &, float const, bool const);
222 template DiscreteProbabilityVector::DiscreteProbabilityVector(const AllFloatInputData &, float const, bool const);
223 template DiscreteProbabilityVector::DiscreteProbabilityVector(const InputData<int, int> &, int const, bool const);
224 
229 
234 
239 
240 } // namespace lar_content
Float_t x
Definition: compare.C:6
std::vector< InputDatum< TX, TY >> InputData
float GetProbabilityDensity(const unsigned int index) const
Get the probability density value of the element in the vector.
void VerifyCompleteData() const
Verify the integrity of the complete probability vector.
float GetX(const unsigned int index) const
Get the x value of the element in the vector.
float EvaluateCumulativeProbability(const float x) const
Evaluate the cumulative probability at arbitrary x.
Float_t y
Definition: compare.C:6
DiscreteProbabilityData ResampleDiscreteProbabilityData(const DiscreteProbabilityVector &discreteProbabilityVector, const ResamplingPoints &resamplingPoints) const
Get a resampled probability data vector by resampling another probability data vector.
std::vector< DiscreteProbabilityDatum > DiscreteProbabilityData
float CalculateNormalisation(const InputData< TX, TY > &inputData) const
Calculate the probability normalisation.
DiscreteProbabilityData InitialiseDiscreteProbabilityData(InputData< TX, TY > inputData) const
Get a initialised probability data vector from the input data.
float GetWidth(const unsigned int index) const
Get the width of the element in the vectorr.
TFile f
Definition: plotHisto.C:6
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
Definition: StdUtils.h:77
DiscreteProbabilityData RandomiseDiscreteProbabilityData(const DiscreteProbabilityVector &discreteProbabilityVector, std::mt19937 &randomNumberGenerator) const
Get a randomised probability data vector in which the x values are unchanged, the probability density...
static bool SortInputDataByX(const InputDatum< TX, TY > &lhs, const InputDatum< TX, TY > &rhs)
Sort the input data according to their x value.
Header file for the lar discrete probability vector class.
DiscreteProbabilityVector(const InputData< TX, TY > &inputData, const TX xUpperBound, const bool useWidths)
Constructor.
DiscreteProbabilityData m_discreteProbabilityData
the probability data
float m_xUpperBound
the upper bound of the probability vector
bool m_useWidths
controls whether bin widths are used in calculations
unsigned int GetSize() const
Get the size of the probability vector.
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
Definition: StdUtils.h:69