Mantid
Loading...
Searching...
No Matches
ComptonScatteringCountRate.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4// NScD Oak Ridge National Laboratory, European Spallation Source,
5// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6// SPDX - License - Identifier: GPL - 3.0 +
12
13#include <sstream>
14
16
17using namespace CurveFitting;
18using Kernel::Logger;
19
20namespace {
22const char *CONSTRAINT_MATRIX_NAME = "IntensityConstraints";
24const char *BKGD_ORDER_ATTR_NAME = "BackgroundOrderAttr";
26Logger g_log("ComptonScatteringCountRate");
27} // namespace
28
29DECLARE_FUNCTION(ComptonScatteringCountRate)
30
31//----------------------------------------------------------------------------------------------
35 : CompositeFunction(), m_profiles(), m_fixedParamIndices(), m_cmatrix(), m_eqMatrix(), m_bkgdOrderAttr("n"),
36 m_bkgdPolyN(0), m_dataErrorRatio() {
37 // Must be a string to be able to be passed through Fit
38 declareAttribute(CONSTRAINT_MATRIX_NAME, IFunction::Attribute("", true));
39 declareAttribute(BKGD_ORDER_ATTR_NAME, IFunction::Attribute(m_bkgdOrderAttr));
40}
41
42//----------------------------------------------------------------------------------------------
43// Private methods
44//----------------------------------------------------------------------------------------------
45
52 if (name == CONSTRAINT_MATRIX_NAME)
53 parseIntensityConstraintMatrix(value.asUnquotedString());
54 else if (name == BKGD_ORDER_ATTR_NAME)
55 m_bkgdOrderAttr = value.asString();
56}
57
67 if (value.empty())
68 throw std::invalid_argument("ComptonScatteringCountRate - Empty string not allowed.");
69 std::istringstream is(value);
71}
72
73namespace {
77struct NoBkgdNorm2 {
79 NoBkgdNorm2(const Kernel::DblMatrix &cmatrix, const std::vector<double> &data)
80 : cm(cmatrix), nrows(cmatrix.numRows()), ncols(cmatrix.numCols()), rhs(data) {}
81
82 double eval(const std::vector<double> &xpt) const {
83 double norm2(0.0);
84 for (size_t i = 0; i < nrows; ++i) {
85 const double *cmRow = cm[i];
86 double cx(0.0);
87 for (size_t j = 0; j < ncols; ++j) {
88 cx += cmRow[j] * xpt[j];
89 }
90 cx -= rhs[i];
91 norm2 += cx * cx;
92 }
93 return 0.5 * norm2;
94 }
96 size_t nrows;
97 size_t ncols;
98 const std::vector<double> &rhs;
99};
103struct BkgdNorm2 {
105 BkgdNorm2(const Kernel::DblMatrix &cmatrix, const std::vector<double> &data)
106 : cm(cmatrix), nrows(cmatrix.numRows()), ncols(cmatrix.numCols()), rhs(data) {}
107
108 double eval(const size_t n, const double *xpt) const {
109 double norm2(0.0);
110 for (size_t i = 0; i < nrows; ++i) {
111 const double *cmRow = cm[i];
112 double cx(0.0);
113 for (size_t j = 0; j < n; ++j) {
114 cx += cmRow[j] * xpt[j];
115 }
116 cx *= -1.0; // our definition of cm has been multiplied by -1
117 cx -= rhs[i];
118 norm2 += cx * cx;
119 }
120 return 0.5 * norm2;
121 }
122 const Kernel::DblMatrix &cm;
123 size_t nrows;
124 size_t ncols;
125 const std::vector<double> &rhs;
126};
127} // namespace
128
150 const size_t nparams(m_cmatrix.numCols());
151
152 // Compute minimization with of Cx where the amplitudes are set to 1.0
153 std::vector<double> x0(nparams, 1);
155 // Compute the constraint matrix
156 this->updateCMatrixValues();
157
158 if (m_bkgdPolyN > 0) {
159 BkgdNorm2 objf(m_cmatrix, m_dataErrorRatio);
160 using namespace std::placeholders;
161 AugmentedLagrangianOptimizer::ObjFunction objfunc = std::bind(&BkgdNorm2::eval, objf, _1, _2);
162 AugmentedLagrangianOptimizer lsqmin(nparams, objfunc, m_eqMatrix, m_cmatrix);
163 lsqmin.minimize(x0);
164 // Set the parameters for the 'real' function calls
166 } else {
167 NoBkgdNorm2 objfunc(m_cmatrix, m_dataErrorRatio);
168 Kernel::Math::SLSQPMinimizer lsqmin(nparams, objfunc, m_eqMatrix, m_cmatrix);
169 auto res = lsqmin.minimize(x0);
170 // Set the parameters for the 'real' function calls
172 }
173}
174
179void ComptonScatteringCountRate::setFixedParameterValues(const std::vector<double> &values) {
180 assert(values.size() == m_fixedParamIndices.size());
181
182 const size_t nparams = values.size();
183 for (size_t i = 0; i < nparams; ++i) {
184 this->setParameter(m_fixedParamIndices[i], values[i], true);
185 }
186
187 if (g_log.is(Logger::Priority::PRIO_DEBUG)) {
188 g_log.debug() << "--- New Intensity Parameters ---\n";
189 for (size_t i = 0; i < nparams; ++i) {
190 g_log.debug() << "x_" << i << "=" << values[i] << "\n";
191 }
192 }
193}
194
199 // -- Compute constraint matrix from each "member function" --
200 const size_t nprofiles = m_profiles.size();
201 for (size_t i = 0, start = 0; i < nprofiles; ++i) {
202 auto *profile = m_profiles[i];
203 const size_t numFilled = profile->fillConstraintMatrix(m_cmatrix, start, m_hist->e());
204 start += numFilled;
205 }
206 // Using different minimizer for background constraints as the original is not
207 // stable
208 // The background one requires the contraints are specified as <= 0
209 if (m_bkgdPolyN > 0) {
210 m_cmatrix *= -1.0;
211 }
212
213 if (g_log.is(Logger::Priority::PRIO_DEBUG)) {
214 g_log.debug() << "--- CM ---\n";
215 for (size_t i = 0; i < m_cmatrix.numRows(); ++i) {
216 for (size_t j = 0; j < m_cmatrix.numCols(); ++j) {
217 g_log.debug() << m_cmatrix[i][j] << " ";
218 }
219 g_log.debug() << "\n";
220 }
221 }
222}
223
231void ComptonScatteringCountRate::setMatrixWorkspace(std::shared_ptr<const API::MatrixWorkspace> matrix, size_t wsIndex,
232 double startX, double endX) {
233 CompositeFunction::setMatrixWorkspace(matrix, wsIndex, startX, endX);
234
235 this->m_hist = std::make_shared<HistogramData::Histogram>(matrix->histogram(wsIndex));
236 this->wsIndex = wsIndex;
237 const auto &values = m_hist->y();
238 const auto &errors = m_hist->e();
239 // Keep the errors for the constraint calculation
240 m_dataErrorRatio.resize(errors.size());
241 std::transform(values.begin(), values.end(), errors.begin(), m_dataErrorRatio.begin(), std::divides<double>());
242
243 if (g_log.is(Kernel::Logger::Priority::PRIO_DEBUG)) {
244 g_log.debug() << "-- data/error --\n";
245 for (size_t i = 0; i < errors.size(); ++i) {
246 g_log.debug() << m_dataErrorRatio[i] << "\n";
247 }
248 }
249
252}
253
254/*
255 * Casts the pointers to ComptonProfiles to avoid overhead during fit. Also
256 * stores the indices, within the composite, of all of the intensity parameters
257 * and
258 * background parameters
259 */
261 // Cache ptrs cast to ComptonProfile functions to that we can compute the
262 // constraint matrix
263 const size_t nfuncs = this->nFunctions();
264 m_profiles.reserve(nfuncs); // won't include background
265 m_fixedParamIndices.reserve(nfuncs + 3); // guess at max size
266 bool foundBkgd(false);
267
268 for (size_t i = 0; i < nfuncs; ++i) {
269 auto func = this->getFunction(i);
270 const size_t paramsOffset = this->paramOffset(i); // offset for ith function inside composite
271
272 if (auto profile = std::dynamic_pointer_cast<ComptonProfile>(func)) {
273 this->cacheComptonProfile(profile, paramsOffset);
274 continue;
275 }
276
277 auto function1D = std::dynamic_pointer_cast<API::IFunction1D>(func);
278 if (!foundBkgd) {
279 this->cacheBackground(function1D, paramsOffset);
280 foundBkgd = true;
281 } else {
282 std::ostringstream os;
283 os << "ComptonScatteringCountRate - Invalid function member at index '" << i << "'. "
284 << "All members must be of type ComptonProfile and at most a single "
285 "1D function";
286 throw std::runtime_error(os.str());
287 }
288 }
289}
290
296void ComptonScatteringCountRate::cacheComptonProfile(const std::shared_ptr<ComptonProfile> &profile,
297 const size_t paramsOffset) {
298 m_profiles.emplace_back(profile.get());
299 auto fixedParams = profile->intensityParameterIndices();
300 for (auto fixedParam : fixedParams) {
301 const size_t indexOfFixed = paramsOffset + fixedParam;
302 this->setParameterStatus(indexOfFixed, Tied);
303 m_fixedParamIndices.emplace_back(indexOfFixed);
304 }
305}
306
312void ComptonScatteringCountRate::cacheBackground(const API::IFunction1D_sptr &function1D, const size_t paramsOffset) {
313 // Check for the order attribute
314 if (function1D->hasAttribute(m_bkgdOrderAttr)) {
315 m_bkgdPolyN = function1D->getAttribute(m_bkgdOrderAttr).asInt();
316 const auto npars = static_cast<size_t>(m_bkgdPolyN + 1); // + constant term
317 // we assume the parameters are at index 0->N on the background so we need
318 // to reverse them
319 for (size_t i = npars; i > 0; --i) // i = from npars->1
320 {
321 const size_t indexOfFixed = paramsOffset + (i - 1);
322 this->setParameterStatus(indexOfFixed, Tied);
323 m_fixedParamIndices.emplace_back(indexOfFixed);
324 }
325 } else {
326 std::ostringstream os;
327 os << "ComptonScatteringCountRate - Background function does not have "
328 "attribute named '"
329 << m_bkgdOrderAttr << "' that specifies its order. Use the '" << BKGD_ORDER_ATTR_NAME
330 << "' attribute to specify the name of the order attribute.";
331 throw std::runtime_error(os.str());
332 }
333}
334
335/*
336 * The equality constraint matrix is padded out to allow for any
337 * additional intensity constraints required each specific Compton profile.
338 * Also creates the inequality matrix
339 * @param xValues The xdata from the workspace
340 */
342 const size_t nmasses = m_profiles.size();
343
344 // Sanity check that equality constraints matrix has the same number of
345 // columns as masses or is zero-sized
346 if (m_eqMatrix.numCols() > 0 && m_eqMatrix.numCols() != nmasses) {
347 std::ostringstream os;
348 os << "ComptonScatteringCountRate - Equality constraint matrix (Aeq) has "
349 "incorrect number of columns ("
350 << m_eqMatrix.numCols() << "). The number of columns should match the number of masses (" << nmasses << ")";
351 throw std::invalid_argument(os.str());
352 }
353
355 createEqualityCM(nmasses);
356
357 if (g_log.is(Logger::Priority::PRIO_DEBUG)) {
358 g_log.debug() << "\n--- aeq ---\n";
359 for (size_t i = 0; i < m_eqMatrix.numRows(); ++i) {
360 for (size_t j = 0; j < m_eqMatrix.numCols(); ++j) {
361 g_log.debug() << m_eqMatrix[i][j] << " ";
362 }
363 g_log.debug() << "\n";
364 }
365 }
366}
367
369 // -- Constraint matrix for J(y) > 0 --
370 // The first N columns are filled with J(y) for each mass + N_h for the first
371 // mass hermite
372 // terms included + (n+1) for each termn the background of order n
373 // The background columns are filled with x^j/error where j=(1,n+1)
374 const auto &xValues = m_hist->x();
375 const auto &errors = m_hist->e();
376
377 const size_t nrows(xValues.size());
378 size_t nColsCMatrix(m_fixedParamIndices.size());
379 m_cmatrix = Kernel::DblMatrix(nrows, nColsCMatrix);
380
381 // Fill background values as they don't change at all
382 if (m_bkgdPolyN > 0) {
383 size_t polyStartCol = nColsCMatrix - m_bkgdPolyN - 1;
384 for (size_t i = 0; i < nrows; ++i) // rows
385 {
386 double *row = m_cmatrix[i];
387 const double &xi = xValues[i];
388 const double &erri = errors[i];
389 size_t polyN = m_bkgdPolyN;
390 for (size_t j = polyStartCol; j < nColsCMatrix; ++j) // cols
391 {
392 row[j] = std::pow(xi, static_cast<double>(polyN)) / erri;
393 --polyN;
394 }
395 }
396 }
397}
398
403 // -- Equality constraints --
404 // The user-specified equality matrix needs to be padded on the left with the
405 // copies of the first column
406 // until the number of cols matches the number of fixed parameters to account
407 // for the extra coefficients
408 // for the hermite term.
409 // It then needs to be padded on the right with zeroes to account for the
410 // background terms
411 auto userMatrix = m_eqMatrix; // copy original
412 const size_t nconstr = userMatrix.numRows();
413
414 const size_t nColsCMatrix(m_cmatrix.numCols());
415 size_t nFixedProfilePars(m_fixedParamIndices.size());
416 // skip background for lhs padding
417 if (m_bkgdPolyN > 0) {
418 nFixedProfilePars -= (m_bkgdPolyN + 1);
419 }
420 const size_t nExtraColsLeft = (nFixedProfilePars - nmasses);
421
422 m_eqMatrix = Kernel::DblMatrix(nconstr, nColsCMatrix); // all zeroed by default.
423 for (size_t i = 0; i < nconstr; ++i) {
424 const double *userRow = userMatrix[i];
425 double *destRow = m_eqMatrix[i];
426 for (size_t j = 0; j < nFixedProfilePars; ++j) {
427 destRow[j] = (j < nExtraColsLeft) ? userRow[0] : userRow[j - nExtraColsLeft];
428 }
429 }
430}
431
432} // namespace Mantid::CurveFitting::Functions
const Kernel::DblMatrix & cm
const std::vector< double > & rhs
double value
The value of the point.
Definition: FitMW.cpp:51
#define DECLARE_FUNCTION(classname)
Macro for declaring a new type of function to be used with the FunctionFactory.
A composite function is a function containing other functions.
void setParameter(size_t, const double &value, bool explicitlySet=true) override
Set i-th parameter.
size_t paramOffset(size_t i) const
std::size_t nFunctions() const override
Number of functions.
void setMatrixWorkspace(std::shared_ptr< const API::MatrixWorkspace > workspace, size_t wi, double startX, double endX) override
Set matrix workspace.
IFunction_sptr getFunction(std::size_t i) const override
Returns the pointer to i-th function.
void setAttribute(const std::string &name, const API::IFunction::Attribute &value) override
Set a value of a named attribute.
void setParameterStatus(size_t i, ParameterStatus status) override
Change status of parameter.
Attribute is a non-fitting parameter.
Definition: IFunction.h:282
Implements the Augmented Lagrangian optimization method of Birgin & Martinez.
boost::function< double(const size_t, const double *)> ObjFunction
Function type.
void minimize(std::vector< double > &xv) const
Perform the minimization.
Implements a specialized function that encapsulates the combination of ComptonProfile functions that ...
Kernel::DblMatrix m_cmatrix
Positivity constraints on J(y)
void setFixedParameterValues(const std::vector< double > &values)
Set the fixed parameters to the given values.
void cacheBackground(const API::IFunction1D_sptr &function1D, const size_t paramsOffset)
Cache parameters positions for background function.
std::string m_bkgdOrderAttr
Name of order attribute on background function.
void updateCMatrixValues() const
Refresh the values of the C matrix for this evaluation.
Kernel::DblMatrix m_eqMatrix
Intensity equality constraints.
std::vector< double > m_dataErrorRatio
Ratio of data & errors.
void cacheFunctions()
Cache ptrs to the individual profiles and their parameters.
std::vector< ComptonProfile * > m_profiles
Holder for non-owning functions cast as ComptonProfiles.
void parseIntensityConstraintMatrix(const std::string &value)
Takes the string & constructs the constraint matrix.
void setMatrixWorkspace(std::shared_ptr< const API::MatrixWorkspace > matrix, size_t wsIndex, double startX, double endX) override
Cache reference to workspace for use in setupForFit.
void setAttribute(const std::string &name, const Attribute &value) override
Set an attribute value (and possibly cache its value)
void cacheComptonProfile(const std::shared_ptr< ComptonProfile > &profile, const size_t paramsOffset)
Cache ptr to the individual profile and its parameters.
std::shared_ptr< HistogramData::Histogram > m_hist
The histogram of the matrix workspace being cached for use.
std::vector< size_t > m_fixedParamIndices
Store parameter indices of intensity parameters that are fixed.
void createEqualityCM(const size_t nmasses)
Set up equality constraint matrix.
void iterationStarting() override
Called by the framework just before an iteration is starting.
void debug(const std::string &msg)
Logs at debug level.
Definition: Logger.cpp:114
bool is(int level) const
Returns true if at least the given log level is set.
Definition: Logger.cpp:146
Minimize an objective function using the SLSQP optimization subroutine originally implemented by Diet...
std::vector< double > minimize(const std::vector< double > &x0) const
Perform the minimization.
size_t numRows() const
Return the number of rows in the matrix.
Definition: Matrix.h:144
size_t numCols() const
Return the number of columns in the matrix.
Definition: Matrix.h:147
Kernel::Logger g_log("ExperimentInfo")
static logger object
std::shared_ptr< IFunction1D > IFunction1D_sptr
Definition: IFunction1D.h:80
double MANTID_CURVEFITTING_DLL norm2(const DoubleFortranVector &v)
Compute the 2-norm of a vector which is a square root of the sum of squares of its elements.
Definition: TrustRegion.cpp:57
void fillFromStream(std::istream &, Kernel::Matrix< T > &, const char)
Fill a Matrix from a stream using the given separator.
Definition: Matrix.cpp:1630
Mantid::Kernel::Matrix< double > DblMatrix
Definition: Matrix.h:206