Commit b396588f authored by Lamar Moore's avatar Lamar Moore
Browse files

Re #15542 merge with master

parent 8ecd3c4f
......@@ -55,7 +55,7 @@ private:
void init() override;
void exec() override;
void validateSpectraIndices(std::vector<int> &v);
void validateWorkspaceIndices(std::vector<int> &v);
void validateChannelIndices(std::vector<int> &v);
std::map<int, int> findElasticPeakPositions(const std::vector<int> &,
......
......@@ -59,13 +59,15 @@ public:
return "Transforms\\Distribution";
}
protected:
/// Validate inputs
std::map<std::string, std::string> validateInputs() override;
private:
/// Initialisation code
void init() override;
/// Execution code
void exec() override;
/// Validate inputs
std::map<std::string, std::string> validateInputs() override;
};
} // namespace Algorithms
......
#ifndef MANTID_ALGORITHMS_EVENTWORKSPACEACCESS_H_
#define MANTID_ALGORITHMS_EVENTWORKSPACEACCESS_H_
#include "MantidDataObjects/EventWorkspace.h"
namespace Mantid {
namespace Algorithms {
struct EventWorkspaceAccess {
static decltype(
std::mem_fn((DataObjects::EventList &
(DataObjects::EventWorkspace::*)(const std::size_t)) &
DataObjects::EventWorkspace::getEventList)) eventList;
};
}
}
#endif /* MANTID_ALGORITHMS_EVENTWORKSPACEACCESS_H_ */
......@@ -55,12 +55,14 @@ public:
/// Algorithm's category for identification overriding a virtual method
const std::string category() const override { return "Arithmetic\\FFT"; }
protected:
/// Perform validation of inputs
std::map<std::string, std::string> validateInputs() override;
private:
// Overridden Algorithm methods
void init() override;
void exec() override;
/// Perform validation of inputs
std::map<std::string, std::string> validateInputs() override;
/// Check whether supplied values are evenly spaced
bool areBinWidthsUneven(const MantidVec &xValues) const;
};
......
......@@ -157,6 +157,7 @@ private:
/// Peak profile type
std::string m_peakType;
/// Criterias for fitting peak
std::string m_minimizer;
double m_maxChiSq;
double m_minPeakHeight;
double m_leastMaxObsY;
......
......@@ -86,13 +86,13 @@ private:
void getGeometry(API::MatrixWorkspace_const_sptr WS, specnum_t mon0Spec,
specnum_t mon1Spec, double &monitor0Dist,
double &monitor1Dist) const;
std::vector<size_t> getMonitorSpecIndexs(API::MatrixWorkspace_const_sptr WS,
specnum_t specNum1,
specnum_t specNum2) const;
std::vector<size_t> getMonitorWsIndexs(API::MatrixWorkspace_const_sptr WS,
specnum_t specNum1,
specnum_t specNum2) const;
double timeToFly(double s, double E_KE) const;
double getPeakCentre(API::MatrixWorkspace_const_sptr WS,
const int64_t monitIn, const double peakTime);
void extractSpec(int64_t specInd, double start, double end);
void extractSpec(int64_t wsInd, double start, double end);
void getPeakEstimates(double &height, int64_t &centreInd,
double &background) const;
double findHalfLoc(MantidVec::size_type startInd, const double height,
......
#ifndef MANTID_ALGORITHMS_MATRIXWORKSPACEACCESS_H_
#define MANTID_ALGORITHMS_MATRIXWORKSPACEACCESS_H_
#include "MantidAPI/MatrixWorkspace.h"
namespace Mantid {
namespace Algorithms {
struct MatrixWorkspaceAccess {
static decltype(std::mem_fn(
(std::vector<double> & (API::MatrixWorkspace::*)(const std::size_t)) &
API::MatrixWorkspace::dataX)) x;
};
}
}
#endif /* MANTID_ALGORITHMS_MATRIXWORKSPACEACCESS_H_ */
......@@ -52,13 +52,15 @@ public:
/// Algorithm's summary
const std::string summary() const override;
protected:
/// Validate the input properties
std::map<std::string, std::string> validateInputs() override;
private:
/// Initialise the algorithm's properties
void init() override;
/// Run the algorithm
void exec() override;
/// Validate the input properties
std::map<std::string, std::string> validateInputs() override;
// Calculates chi-square by solving the matrix equation A*x = b
double calculateChi(const QuadraticCoefficients &coeffs, double a,
std::vector<double> &beta);
......
......@@ -77,6 +77,8 @@ public:
void calculateQuadraticCoefficients();
private:
// Initializes the member variables related to the image
void initImageSpace(const std::vector<double> &image, double background);
// Calculates the gradient of chi
std::vector<double> calculateChiGrad() const;
// Calculates the entropy
......
......@@ -34,10 +34,6 @@ namespace Algorithms {
*/
class MANTID_ALGORITHMS_DLL MaxentEntropyPositiveValues : public MaxentEntropy {
public:
// Constructor
MaxentEntropyPositiveValues() = default;
// Destructor
virtual ~MaxentEntropyPositiveValues() = default;
// First derivative
double getDerivative(double value) override;
// Second derivative
......
......@@ -108,9 +108,9 @@ private:
/// over-reading
double m_hiFrac;
/// The index of the first spectrum to calculate
int m_minSpec;
int m_minWsIndex;
/// The index of the last spectrum to calculate
int m_maxSpec;
int m_maxWsIndex;
/// Start point for integration
double m_rangeLower;
/// End point for integration
......
......@@ -69,7 +69,7 @@ private:
API::MatrixWorkspace_sptr
setUpOutputWorkspace(const std::vector<double> &binParams) const;
// these are the steps that are run on each individual spectrum
void calculateNormalization(const size_t wavStart, const size_t specInd,
void calculateNormalization(const size_t wavStart, const size_t wsIndex,
API::MatrixWorkspace_const_sptr pixelAdj,
API::MatrixWorkspace_const_sptr wavePixelAdj,
double const *const binNorms,
......@@ -77,16 +77,16 @@ private:
const MantidVec::iterator norm,
const MantidVec::iterator normETo2) const;
void pixelWeight(API::MatrixWorkspace_const_sptr pixelAdj,
const size_t specIndex, double &weight, double &error) const;
const size_t wsIndex, double &weight, double &error) const;
void addWaveAdj(const double *c, const double *Dc, MantidVec::iterator bInOut,
MantidVec::iterator e2InOut) const;
void addWaveAdj(const double *c, const double *Dc, MantidVec::iterator bInOut,
MantidVec::iterator e2InOut, MantidVec::const_iterator,
MantidVec::const_iterator) const;
void normToMask(const size_t offSet, const size_t specIndex,
void normToMask(const size_t offSet, const size_t wsIndex,
const MantidVec::iterator theNorms,
const MantidVec::iterator errorSquared) const;
void convertWavetoQ(const size_t specInd, const bool doGravity,
void convertWavetoQ(const size_t wsInd, const bool doGravity,
const size_t offset, MantidVec::iterator Qs,
const double extraLength) const;
void getQBinPlus1(const MantidVec &OutQs, const double QToFind,
......
......@@ -48,46 +48,13 @@ public:
size_t waveLengthCutOff(API::MatrixWorkspace_const_sptr dataWS,
const double RCut, const double WCut,
const size_t specInd) const;
const size_t wsInd) const;
void outputParts(API::Algorithm *alg, API::MatrixWorkspace_sptr sumOfCounts,
API::MatrixWorkspace_sptr sumOfNormFactors);
private:
/// the experimental workspace with counts across the detector
/* API::MatrixWorkspace_const_sptr m_dataWS;
///The radius cut off, should be value of the property RadiusCut. A value of
zero here will disable the cut off and all wavelengths are used
double m_RCut;
///The wavelength cut off divided by the radius cut is used in the
calculation of the first wavelength to include, it's value is only used if
RadiusCut > 0
double m_WCutOver;
void initizeCutOffs(const double RCut, const double WCut);
API::MatrixWorkspace_sptr setUpOutputWorkspace(const std::vector<double> &
binParams) const;
//these are the steps that are run on each individual spectrum
size_t waveLengthCutOff(const size_t specInd) const;
void calculateNormalization(const size_t wavStart, const size_t specInd,
API::MatrixWorkspace_const_sptr pixelAdj, double const * const binNorms,
double const * const binNormEs, const MantidVec::iterator norm, const
MantidVec::iterator normETo2) const;
void pixelWeight(API::MatrixWorkspace_const_sptr pixelAdj, const size_t
specIndex, double & weight, double & error) const;
void addWaveAdj(const double * c, const double * Dc, MantidVec::iterator
bInOut, MantidVec::iterator e2InOut) const;
void normToMask(const size_t offSet, const size_t specIndex, const
MantidVec::iterator theNorms, const MantidVec::iterator errorSquared) const;
void convertWavetoQ(const size_t specInd, const bool doGravity, const size_t
offset, MantidVec::iterator Qs) const;
void getQBinPlus1(const MantidVec & OutQs, const double QToFind,
MantidVec::const_iterator & loc) const;
void normalize(const MantidVec & normSum, const MantidVec & normError2,
MantidVec & YOut, MantidVec & errors) const;*/
};
} // namespace Algorithms
......
#ifndef MANTID_ALGORITHMS_SPECTRUMALGORITHM_H_
#define MANTID_ALGORITHMS_SPECTRUMALGORITHM_H_
#include <tuple>
#include "MantidKernel/IndexSet.h"
#include "MantidAPI/Algorithm.h"
#include "MantidAPI/MatrixWorkspace_fwd.h"
#include "MantidAlgorithms/DllConfig.h"
namespace Mantid {
namespace DataObjects {
class EventWorkspace;
}
namespace Algorithms {
/** SpectrumAlgorithm is a base class for algorithms that work with
MatrixWorkspace.
This child class of Algorithm provides several features that make writing more
generic and more compact code for algorithms easier. In particular it
provides:
1. The method for_each() that can be used to implement loops/transformations
of spectra or event lists in a workspace.
2. A way to define generic properties to allow user specified spectrum number
ranges and list.
@author Simon Heybrock, ESS
Copyright &copy; 2016 ISIS Rutherford Appleton Laboratory, NScD Oak Ridge
National Laboratory & European Spallation Source
This file is part of Mantid.
Mantid is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 3 of the License, or
(at your option) any later version.
Mantid is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
File change history is stored at: <https://github.com/mantidproject/mantid>
Code Documentation is available at: <http://doxygen.mantidproject.org>
*/
class MANTID_ALGORITHMS_DLL SpectrumAlgorithm : public API::Algorithm {
private:
/** Helpers for for_each(), struct seq and gens with a specialization.
*
* Used for generating a sequence 0,1,2,.... for extracting arguments from
* tuple. */
template <int...> struct seq {};
// Doxygen does not like recursive types.
template <int N, int... S>
struct gens /** @cond */
: gens<N - 1, N - 1, S...> /** @endcond */ {/** @cond */
};
template <int... S> struct gens<0, S...> {/** @endcond */
typedef seq<S...> type;
};
/** Helpers for for_each(), struct contains and 2 specializations.
*
* Used to determine at compile time if parameter pack contains a certain
* type. */
template <typename Tp, typename... List> struct contains : std::true_type {};
template <typename Tp, typename Head, typename... Rest>
struct contains<Tp, Head, Rest...>
: std::conditional<std::is_same<Tp, Head>::value, std::true_type,
contains<Tp, Rest...>>::type {};
template <typename Tp> struct contains<Tp> : std::false_type {};
/// Internal implementation of for_each().
template <class... Flags, class WS, class T, int... S, class OP>
void for_each(WS &workspace, T getters, seq<S...>, const OP &operation) {
// If we get the flag Indices::FromProperty we use a potential user-defined
// range property, otherwise default to the full range of all histograms.
Kernel::IndexSet indexSet(workspace.getNumberHistograms());
if (contains<Indices::FromProperty, Flags...>())
indexSet = getWorkspaceIndexSet(workspace);
auto size = static_cast<int64_t>(indexSet.size());
API::Progress progress(this, 0.0, 1.0, size);
PARALLEL_FOR1((&workspace))
for (int64_t i = 0; i < size; ++i) {
PARALLEL_START_INTERUPT_REGION
// Note the small but for now negligible overhead from the IndexSet access
// in the case where it is not used.
operation(std::get<S>(getters)(workspace, indexSet[i])...);
progress.report(name());
PARALLEL_END_INTERUPT_REGION
}
PARALLEL_CHECK_INTERUPT_REGION
ifEventWorkspaceClearMRU(workspace);
}
/// Templated function used as 'compile-time conditional', together with
/// specialization. See there for details.
template <class WS> void ifEventWorkspaceClearMRU(const WS &workspace) {
UNUSED_ARG(workspace);
}
std::string m_indexMinPropertyName;
std::string m_indexMaxPropertyName;
std::string m_indexRangePropertyName;
protected:
/// Dummy struct holding compile-time flags to for_each().
// A strongly typed enum used in non-type variadic template arguments would
// have been the other candidate. However, combining flags from multiple enums
// would not have been possible, so all flags would have been required to be
// in the same scoped enum, and we want to avoid noisy flag definitions like
// Flags::IndicesFromProperty or Flags::SkipMasked.
struct Indices {
/// Flag: Include only indices specified via properties in for_each.
struct FromProperty {};
};
/** Provides a mechanism for looping over spectra in a workspace.
*
* This variant works with a single workspace and can to in-place modification
* of spectra or event lists in the workspace. Threading and progress
* reporting is handled internally.
* @tparam Flags Variable number of flags, see struct Indices.
* @param workspace Workspace to work with.
* @param getters Tuple of getters, pointers to methods of workspace, defaults
* defined for convenience in namespace Getters.
* @param operation Callable that is executed for all spectra (etc.). Results
* from getters are passed as arguments. */
template <class... Flags, class WS, class... Args, class OP>
void for_each(WS &workspace, std::tuple<Args...> getters,
const OP &operation) {
// This comment explains some of the rationale and mechanism for the way
// templates are used in this and other variants of for_each().
// For several variants of for_each() we require a variable number of
// arguments for more than one entity, e.g., getters for the input workspace
// and getters for the output workspace. For both of them we would thus like
// to use a variadic template. However, we cannot make all those getters
// direct arguments to for_each(), because the compiler then cannot tell
// were the boundary between the two parameter packs is. Therefore getters
// are packed into tuples.
// At some point we need to extract the getters from the tuple. To this end,
// we use the gens template, which recursively constructs a sequence of
// indices which are gen used with std::get<>.
// The Flags template parameter is used for passing in flags known at
// compile time.
for_each<Flags...>(workspace, getters,
typename gens<sizeof...(Args)>::type(), operation);
}
void declareWorkspaceIndexSetProperties(
const std::string &indexMinPropertyName = "IndexMin",
const std::string &indexMaxPropertyName = "IndexMax",
const std::string &indexRangePropertyName = "WorkspaceIndexList");
Kernel::IndexSet
getWorkspaceIndexSet(const API::MatrixWorkspace &workspace) const;
};
template <>
void SpectrumAlgorithm::ifEventWorkspaceClearMRU(
const DataObjects::EventWorkspace &workspace);
/// Typedef for a shared pointer to a SpectrumAlgorithm
typedef boost::shared_ptr<SpectrumAlgorithm> SpectrumAlgorithm_sptr;
} // namespace Algorithms
} // namespace Mantid
#endif /* MANTID_ALGORITHMS_SPECTRUMALGORITHM_H_ */
......@@ -96,9 +96,9 @@ private:
/// The output spectrum id
specnum_t m_outSpecId;
/// The spectrum to start the integration from
int m_minSpec;
int m_minWsInd;
/// The spectrum to finish the integration at
int m_maxSpec;
int m_maxWsInd;
/// Set true to keep monitors
bool m_keepMonitors;
/// numberOfSpectra in the input
......
......@@ -130,5 +130,43 @@ void ApplyDeadTimeCorr::exec() {
}
}
/**
* Validates input properties:
* - input workspace must have all bins the same size (within reasonable error)
* @returns :: map of property names to error strings (empty if no error)
*/
std::map<std::string, std::string> ApplyDeadTimeCorr::validateInputs() {
std::map<std::string, std::string> errors;
MatrixWorkspace_const_sptr inputWS = getProperty("InputWorkspace");
if (inputWS) { // in case input was a WorkspaceGroup
const MantidVec &xValues = inputWS->readX(0);
const size_t xSize = xValues.size();
if (xSize < 2) {
errors["InputWorkspace"] = "Input workspace cannot be empty";
} else {
// average bin width
const double dx =
(xValues[xSize - 1] - xValues[0]) / static_cast<double>(xSize - 1);
// Use cumulative errors
auto difference = [&xValues, &dx](size_t i) {
return std::abs((xValues[i] - xValues[0] - (double)i * dx) / dx);
};
// Check each width against dx
constexpr double tolerance = 0.5;
for (size_t i = 1; i < xValues.size() - 2; i++) {
if (difference(i) > tolerance) {
g_log.error() << "dx=" << xValues[i + 1] - xValues[i] << ' ' << dx
<< ' ' << i << std::endl;
errors["InputWorkspace"] = "Uneven bin widths in input workspace";
break;
}
}
}
}
return errors;
}
} // namespace Mantid
} // namespace Algorithms
......@@ -98,7 +98,7 @@ void ApplyTransmissionCorrection::exec() {
try {
det = inputWS->getDetector(i);
} catch (Exception::NotFoundError &) {
g_log.warning() << "Spectrum index " << i
g_log.warning() << "Workspace index " << i
<< " has no detector assigned to it - discarding"
<< std::endl;
// Catch if no detector. Next line tests whether this happened - test
......
......@@ -70,28 +70,30 @@ std::map<std::string, std::string> CalMuonDetectorPhases::validateInputs() {
API::MatrixWorkspace_const_sptr inputWS = getProperty("InputWorkspace");
// Check units, should be microseconds
Unit_const_sptr unit = inputWS->getAxis(0)->unit();
if ((unit->label().ascii() != "Microseconds") &&
(unit->label().ascii() != "microsecond")) {
result["InputWorkspace"] = "InputWorkspace units must be microseconds";
}
if (inputWS) {
// Check units, should be microseconds
Unit_const_sptr unit = inputWS->getAxis(0)->unit();
if ((unit->label().ascii() != "Microseconds") &&
(unit->label().ascii() != "microsecond")) {
result["InputWorkspace"] = "InputWorkspace units must be microseconds";
}
// Check spectra numbers are valid, if specified
int nspec = static_cast<int>(inputWS->getNumberHistograms());
std::vector<int> forward = getProperty("ForwardSpectra");
std::vector<int> backward = getProperty("BackwardSpectra");
for (int spec : forward) {
if (spec < 1 || spec > nspec) {
result["ForwardSpectra"] = "Invalid spectrum numbers in ForwardSpectra";
// Check spectra numbers are valid, if specified
int nspec = static_cast<int>(inputWS->getNumberHistograms());
std::vector<int> forward = getProperty("ForwardSpectra");
std::vector<int> backward = getProperty("BackwardSpectra");
for (int spec : forward) {
if (spec < 1 || spec > nspec) {
result["ForwardSpectra"] = "Invalid spectrum numbers in ForwardSpectra";
}
}
}
for (int spec : backward) {
if (spec < 1 || spec > nspec) {
result["BackwardSpectra"] = "Invalid spectrum numbers in BackwardSpectra";
for (int spec : backward) {
if (spec < 1 || spec > nspec) {
result["BackwardSpectra"] =
"Invalid spectrum numbers in BackwardSpectra";
}
}
}
return result;
}
//----------------------------------------------------------------------------------------------
......@@ -170,6 +172,7 @@ void CalMuonDetectorPhases::fitWorkspace(const API::MatrixWorkspace_sptr &ws,
if (!fit->isExecuted() || status != success) {
std::ostringstream error;
error << "Fit failed for spectrum " << ispec;
error << ": " << status;
throw std::runtime_error(error.str());
}
......
......@@ -98,9 +98,9 @@ void CalculateFlatBackground::exec() {
double startX, endX;
this->checkRange(startX, endX);
std::vector<int> specInds = getProperty("WorkspaceIndexList");
std::vector<int> wsInds = getProperty("WorkspaceIndexList");
// check if the user passed an empty list, if so all of spec will be processed
this->getSpecInds(specInds, numHists);
this->getWsInds(wsInds, numHists);
// Are we removing the background?
const bool removeBackground =
......@@ -131,14 +131,14 @@ void CalculateFlatBackground::exec() {
// these are used to report information to the user, one progress update for
// each percent and a report on the size of the background found
double prg(0.2), backgroundTotal(0);
const double toFitsize(static_cast<double>(specInds.size()));
const double toFitsize(static_cast<double>(wsInds.size()));
const int progStep(static_cast<int>(ceil(toFitsize / 80.0)));
// Now loop over the required spectra
std::vector<int>::const_iterator specIt;
// local cache for global variable
bool skipMonitors(m_skipMonitors);
for (specIt = specInds.begin(); specIt != specInds.end(); ++specIt) {
for (specIt = wsInds.begin(); specIt != wsInds.end(); ++specIt) {
const int currentSpec = *specIt;
try {
if (skipMonitors) { // this will fail in Windows ReleaseWithDebug info
......@@ -212,7 +212,7 @@ void CalculateFlatBackground::exec() {
}
// make regular progress reports and check for canceling the algorithm
if (static_cast<int>(specInds.end() - specInds.begin()) % progStep == 0) {
if (static_cast<int>(wsInds.end() - wsInds.begin()) % progStep == 0) {
interruption_point();
prg += (progStep * 0.7 / toFitsize);
progress(prg);
......@@ -309,8 +309,8 @@ void CalculateFlatBackground::checkRange(double &startX, double &endX) {
* @param workspaceTotal :: required to be the total number of spectra in the
* workspace