17using namespace CurveFitting;
22const char *CONSTRAINT_MATRIX_NAME =
"IntensityConstraints";
24const char *BKGD_ORDER_ATTR_NAME =
"BackgroundOrderAttr";
26Logger
g_log(
"ComptonScatteringCountRate");
35 :
CompositeFunction(), m_profiles(), m_fixedParamIndices(), m_cmatrix(), m_eqMatrix(), m_bkgdOrderAttr("
n"),
36 m_bkgdPolyN(0), m_dataErrorRatio() {
52 if (
name == CONSTRAINT_MATRIX_NAME)
54 else if (
name == BKGD_ORDER_ATTR_NAME)
68 throw std::invalid_argument(
"ComptonScatteringCountRate - Empty string not allowed.");
69 std::istringstream is(
value);
80 :
cm(cmatrix),
nrows(cmatrix.numRows()),
ncols(cmatrix.numCols()),
rhs(data) {}
82 double eval(
const std::vector<double> &xpt)
const {
84 for (
size_t i = 0; i <
nrows; ++i) {
85 const double *cmRow =
cm[i];
87 for (
size_t j = 0; j <
ncols; ++j) {
88 cx += cmRow[j] * xpt[j];
98 const std::vector<double> &
rhs;
106 :
cm(cmatrix),
nrows(cmatrix.numRows()),
ncols(cmatrix.numCols()),
rhs(data) {}
108 double eval(
const size_t n,
const double *xpt)
const {
110 for (
size_t i = 0; i <
nrows; ++i) {
111 const double *cmRow =
cm[i];
113 for (
size_t j = 0; j <
n; ++j) {
114 cx += cmRow[j] * xpt[j];
125 const std::vector<double> &
rhs;
153 std::vector<double> x0(nparams, 1);
160 using namespace std::placeholders;
182 const size_t nparams = values.size();
183 for (
size_t i = 0; i < nparams; ++i) {
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";
201 for (
size_t i = 0, start = 0; i < nprofiles; ++i) {
203 const size_t numFilled = profile->fillConstraintMatrix(
m_cmatrix, start,
m_hist->e());
213 if (
g_log.
is(Logger::Priority::PRIO_DEBUG)) {
232 double startX,
double endX) {
235 this->
m_hist = std::make_shared<HistogramData::Histogram>(matrix->histogram(
wsIndex));
237 const auto &values =
m_hist->y();
238 const auto &errors =
m_hist->e();
241 std::transform(values.begin(), values.end(), errors.begin(),
m_dataErrorRatio.begin(), std::divides<double>());
243 if (
g_log.
is(Kernel::Logger::Priority::PRIO_DEBUG)) {
245 for (
size_t i = 0; i < errors.size(); ++i) {
266 bool foundBkgd(
false);
268 for (
size_t i = 0; i < nfuncs; ++i) {
272 if (
auto profile = std::dynamic_pointer_cast<ComptonProfile>(func)) {
277 auto function1D = std::dynamic_pointer_cast<API::IFunction1D>(func);
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 "
286 throw std::runtime_error(os.str());
297 const size_t paramsOffset) {
299 auto fixedParams = profile->intensityParameterIndices();
300 for (
auto fixedParam : fixedParams) {
301 const size_t indexOfFixed = paramsOffset + fixedParam;
316 const auto npars =
static_cast<size_t>(
m_bkgdPolyN + 1);
319 for (
size_t i = npars; i > 0; --i)
321 const size_t indexOfFixed = paramsOffset + (i - 1);
326 std::ostringstream os;
327 os <<
"ComptonScatteringCountRate - Background function does not have "
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());
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());
357 if (
g_log.
is(Logger::Priority::PRIO_DEBUG)) {
374 const auto &xValues =
m_hist->x();
375 const auto &errors =
m_hist->e();
377 const size_t nrows(xValues.size());
383 size_t polyStartCol = nColsCMatrix -
m_bkgdPolyN - 1;
384 for (
size_t i = 0; i <
nrows; ++i)
387 const double &xi = xValues[i];
388 const double &erri = errors[i];
390 for (
size_t j = polyStartCol; j < nColsCMatrix; ++j)
392 row[j] = std::pow(xi,
static_cast<double>(polyN)) / erri;
412 const size_t nconstr = userMatrix.
numRows();
420 const size_t nExtraColsLeft = (nFixedProfilePars - nmasses);
423 for (
size_t i = 0; i < nconstr; ++i) {
424 const double *userRow = userMatrix[i];
426 for (
size_t j = 0; j < nFixedProfilePars; ++j) {
427 destRow[j] = (j < nExtraColsLeft) ? userRow[0] : userRow[j - nExtraColsLeft];
const Kernel::DblMatrix & cm
const std::vector< double > & rhs
double value
The value of the point.
#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.
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.
void createConstraintMatrices()
Set up the constraint matrices.
std::string name() const override
String identifier.
int m_bkgdPolyN
The order of the background.
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 createPositivityCM()
Set up positivity constraint matrix.
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.
size_t wsIndex
The workspace index being worked on.
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.
bool is(int level) const
Returns true if at least the given log level is set.
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.
size_t numCols() const
Return the number of columns in the matrix.
Kernel::Logger g_log("ExperimentInfo")
static logger object
std::shared_ptr< IFunction1D > IFunction1D_sptr
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.
void fillFromStream(std::istream &, Kernel::Matrix< T > &, const char)
Fill a Matrix from a stream using the given separator.
Mantid::Kernel::Matrix< double > DblMatrix