15#include "MantidHistogramData/HistogramBuilder.h"
16#include "MantidHistogramData/HistogramIterator.h"
23#include "boost/make_shared.hpp"
39using Mantid::HistogramData::HistogramBuilder;
48enum class LineDirection { horizontal, vertical };
51namespace DirectionChoices {
52const static std::string HORIZONTAL{
"Horizontal"};
53const static std::string VERTICAL{
"Vertical"};
57namespace ModeChoices {
58const static std::string AVERAGE{
"Average"};
59const static std::string SUM{
"Sum"};
64const static std::string
CENTRE{
"Centre"};
65const static std::string DIRECTION{
"Direction"};
66const static std::string END{
"End"};
68const static std::string HALF_WIDTH{
"HalfWidth"};
69const static std::string IGNORE_INFS{
"IgnoreInfs"};
70const static std::string IGNORE_NANS{
"IgnoreNans"};
71const static std::string MODE{
"Mode"};
73const static std::string START{
"Start"};
100Workspace2D_sptr makeOutput(
const MatrixWorkspace &parent,
const LineDirection direction, std::vector<double> &&Xs,
101 std::vector<double> &&Ys, std::vector<double> &&Es) {
102 HistogramBuilder builder;
103 builder.setX(std::move(Xs));
104 builder.setY(std::move(Ys));
105 builder.setE(std::move(Es));
106 builder.setDistribution(direction == LineDirection::horizontal && parent.isDistribution());
107 return create<Workspace2D>(parent, 1, builder.build());
117void setAxesAndUnits(Workspace2D &outWS,
const MatrixWorkspace &ws,
const Box &box,
const LineDirection dir) {
119 outWS.setYUnit(ws.YUnit());
120 outWS.setYUnitLabel(ws.YUnitLabel());
122 auto axisIndex = dir == LineDirection::horizontal ? 0 : 1;
123 if (ws.getAxis(axisIndex)->isSpectra()) {
124 outWS.getAxis(axisIndex)->setUnit(
"Empty");
126 outWS.getAxis(0)->setUnit(ws.getAxis(axisIndex)->unit()->unitID());
129 std::vector<double> vertBins(2);
130 vertBins.front() = dir == LineDirection::horizontal ? box.top : box.left;
131 vertBins.back() = dir == LineDirection::horizontal ? box.bottom : box.right;
132 auto outVertAxis = std::make_unique<BinEdgeAxis>(vertBins);
133 axisIndex = dir == LineDirection::horizontal ? 1 : 0;
134 if (ws.getAxis(axisIndex)->isSpectra()) {
135 outVertAxis->setUnit(
"Empty");
137 outVertAxis->setUnit(ws.getAxis(axisIndex)->unit()->unitID());
139 outWS.replaceAxis(1, std::move(outVertAxis));
151template <
typename Container>
152std::pair<size_t, size_t> startAndEnd(
const Container &bins,
const bool isBinEdges,
const double lowerLimit,
153 const double upperLimit) {
154 auto lowerIt = std::upper_bound(bins.cbegin(), bins.cend(), lowerLimit);
155 if (lowerIt == bins.cend()) {
156 throw std::runtime_error(
"Profile completely outside input workspace.");
158 if (lowerIt != bins.cbegin()) {
161 auto upperIt = std::upper_bound(lowerIt, bins.cend(), upperLimit);
162 if (upperIt == bins.cbegin()) {
163 throw std::runtime_error(
"Profile completely outside input workspace.");
165 if (isBinEdges && upperIt == bins.cend()) {
168 const auto start = std::distance(bins.cbegin(), lowerIt);
169 const auto end = std::distance(bins.cbegin(), upperIt);
170 return std::pair<size_t, size_t>{start, end};
180std::vector<double> extractVerticalBins(
const Axis &axis,
const size_t numberHistograms) {
181 if (axis.isSpectra()) {
182 std::vector<double> spectrumNumbers(numberHistograms);
183 std::iota(spectrumNumbers.begin(), spectrumNumbers.end(), 1.0);
184 return spectrumNumbers;
186 std::vector<double> bins(axis.length());
187 for (
size_t i = 0; i < bins.size(); ++i) {
188 bins[i] = axis.getValue(i);
207template <
typename Container,
typename Function>
208void profile(std::vector<double> &Xs, std::vector<double> &Ys, std::vector<double> &Es,
const MatrixWorkspace &ws,
209 const LineDirection dir,
const IndexLimits &limits,
const Container &lineBins,
const bool isBinEdges,
210 Function modeFunction,
const bool ignoreNans,
const bool ignoreInfs) {
211 const auto lineSize = limits.lineEnd - limits.lineStart;
212 Xs.resize(lineSize + (isBinEdges ? 1 : 0));
215 for (
size_t i = limits.lineStart; i < limits.lineEnd; ++i) {
216 Xs[i - limits.lineStart] = lineBins[i];
220 for (
size_t j = limits.widthStart; j < limits.widthEnd; ++j) {
221 const size_t iHor = dir == LineDirection::horizontal ? i : j;
222 const size_t iVert = dir == LineDirection::horizontal ? j : i;
223 auto histogram = ws.histogram(iVert);
224 auto iter = histogram.begin();
225 std::advance(iter, iHor);
226 const double y = iter->counts();
227 if ((ignoreNans && std::isnan(
y)) || (ignoreInfs && std::isinf(
y))) {
230 const double e = iter->countStandardDeviation();
235 const size_t nTotal = limits.widthEnd - limits.widthStart;
236 Ys[i - limits.lineStart] =
n == 0 ? std::nan(
"") : modeFunction(ySum,
n, nTotal);
237 const double e = modeFunction(std::sqrt(eSqSum),
n, nTotal);
238 Es[i - limits.lineStart] = std::isnan(e) ? 0 : e;
241 Xs.back() = lineBins[limits.lineEnd];
251double averageMode(
const double sum,
const size_t n,
const size_t )
noexcept {
252 return sum /
static_cast<double>(
n);
263double sumMode(
const double sum,
const size_t n,
const size_t nTot)
noexcept {
264 return static_cast<double>(nTot) /
static_cast<double>(
n) * sum;
271auto createMode(
const std::string &modeName)
noexcept {
272 if (modeName == ModeChoices::AVERAGE) {
278void divideByBinHeight(MatrixWorkspace &ws) {
279 const BinEdgeAxis &axis = *
static_cast<BinEdgeAxis *
>(ws.getAxis(1));
280 const auto height = axis.getMax() - axis.getMin();
307 const auto mandatoryDouble = std::make_shared<MandatoryValidator<double>>();
308 const auto positiveDouble = std::make_shared<BoundedValidator<double>>();
309 positiveDouble->setLower(0.0);
310 positiveDouble->setLowerExclusive(
true);
311 const auto mandatoryPositiveDouble = std::make_shared<CompositeValidator>();
312 mandatoryPositiveDouble->add(mandatoryDouble);
313 mandatoryPositiveDouble->add(positiveDouble);
314 const auto inputWorkspaceValidator = std::make_shared<CompositeValidator>();
315 inputWorkspaceValidator->add(std::make_shared<CommonBinsValidator>());
316 inputWorkspaceValidator->add(std::make_shared<IncreasingAxisValidator>());
319 "An input workspace.");
322 "A single histogram workspace containing the profile.");
325 "Half of the width over which to calcualte the profile.");
326 const std::array<std::string, 2> directions = {{DirectionChoices::HORIZONTAL, DirectionChoices::VERTICAL}};
327 declareProperty(PropertyNames::DIRECTION, DirectionChoices::HORIZONTAL,
331 const std::array<std::string, 2> modes = {{ModeChoices::AVERAGE, ModeChoices::SUM}};
333 "How the profile is calculated over the line width.");
334 declareProperty(PropertyNames::IGNORE_INFS,
false,
"If true, ignore infinities when calculating the profile.");
335 declareProperty(PropertyNames::IGNORE_NANS,
true,
"If true, ignore not-a-numbers when calculating the profile.");
344 const bool ignoreNans =
getProperty(PropertyNames::IGNORE_NANS);
345 const bool ignoreInfs =
getProperty(PropertyNames::IGNORE_INFS);
346 const auto &horizontalBins = ws->x(0);
347 const auto horizontalIsBinEdges = ws->isHistogramData();
348 const auto vertAxis = ws->getAxis(1);
350 const auto verticalBins = extractVerticalBins(*vertAxis, ws->getNumberHistograms());
351 const auto verticalIsBinEdges = verticalBins.size() > ws->getNumberHistograms();
352 const std::string directionString =
getProperty(PropertyNames::DIRECTION);
353 LineDirection dir{LineDirection::horizontal};
354 if (directionString == DirectionChoices::VERTICAL) {
355 dir = LineDirection::vertical;
357 const double centre =
getProperty(PropertyNames::CENTRE);
358 const double halfWidth =
getProperty(PropertyNames::HALF_WIDTH);
361 start = std::numeric_limits<double>::lowest();
365 end = std::numeric_limits<double>::max();
370 if (dir == LineDirection::horizontal) {
371 bounds.top = centre - halfWidth;
372 bounds.bottom = centre + halfWidth;
378 bounds.left = centre - halfWidth;
379 bounds.right = centre + halfWidth;
382 const auto vertInterval = startAndEnd(verticalBins, verticalIsBinEdges, bounds.top, bounds.bottom);
383 const auto horInterval = startAndEnd(horizontalBins, horizontalIsBinEdges, bounds.left, bounds.right);
385 auto mode = createMode(
getProperty(PropertyNames::MODE));
387 std::vector<double> profileYs;
388 std::vector<double> profileEs;
389 std::vector<double> Xs;
390 if (dir == LineDirection::horizontal) {
392 limits.lineStart = horInterval.first;
393 limits.lineEnd = horInterval.second;
394 limits.widthStart = vertInterval.first;
395 limits.widthEnd = vertInterval.second;
396 profile(Xs, profileYs, profileEs, *ws, dir, limits, horizontalBins, horizontalIsBinEdges, mode, ignoreNans,
400 limits.lineStart = vertInterval.first;
401 limits.lineEnd = vertInterval.second;
402 limits.widthStart = horInterval.first;
403 limits.widthEnd = horInterval.second;
404 profile(Xs, profileYs, profileEs, *ws, dir, limits, verticalBins, verticalIsBinEdges, mode, ignoreNans, ignoreInfs);
407 auto outWS = makeOutput(*ws, dir, std::move(Xs), std::move(profileYs), std::move(profileEs));
411 actualBounds.top = verticalBins[vertInterval.first];
412 actualBounds.bottom =
413 vertInterval.second < verticalBins.size() ? verticalBins[vertInterval.second] : verticalBins.back();
414 actualBounds.left = horizontalBins[horInterval.first];
416 horInterval.second < horizontalBins.size() ? horizontalBins[horInterval.second] : horizontalBins.back();
417 setAxesAndUnits(*outWS, *ws, actualBounds, dir);
418 if (dir == LineDirection::vertical && ws->isDistribution()) {
419 divideByBinHeight(*outWS);
427 std::map<std::string, std::string> issues;
428 const double start =
getProperty(PropertyNames::START);
429 const double end =
getProperty(PropertyNames::END);
431 issues[PropertyNames::START] = PropertyNames::START +
" greater than " + PropertyNames::END +
".";
434 if (ws->getAxis(1)->isText()) {
#define DECLARE_ALGORITHM(classname)
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Class to represent the axis of a workspace.
Stores numeric values that are assumed to be bin edge values.
A validator which provides a TENTATIVE check that a workspace contains common bins in each spectrum.
A validator which checks that the X axis of a workspace is increasing from left to right.
Base MatrixWorkspace Abstract Class.
A property class for workspaces.
LineProfile : Calculates a horizontal or vertical line profile over a MatrixWorkspace.
void init() override
Initialize the algorithm's properties.
int version() const override
Algorithm's version for identification.
const std::string category() const override
Algorithm's category for identification.
void exec() override
Execute the algorithm.
const std::string summary() const override
Algorithm's summary for use in the GUI and help.
std::map< std::string, std::string > validateInputs() override
Validate the algorithm's inputs.
Concrete workspace implementation.
BoundedValidator is a validator that requires the values to be between upper or lower bounds,...
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
ListValidator is a validator that requires the value of a property to be one of a defined list of pos...
Validator to check that a property is not left empty.
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
static const std::string OUTPUT_WORKSPACE
static const std::string INPUT_WORKSPACE
std::string CENTRE("Centre")
std::shared_ptr< Workspace2D > Workspace2D_sptr
shared pointer to Mantid::DataObjects::Workspace2D
std::unique_ptr< T > create(const P &parent, const IndexArg &indexArg, const HistArg &histArg)
This is the create() method that all the other create() methods call.
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
Describes the direction (within an algorithm) of a Property.
@ Input
An input workspace.
@ Output
An output workspace.