From 70bc199c71db5cb150196ca86b17ef86f0d58411 Mon Sep 17 00:00:00 2001 From: Verena Reimund <reimund@ill.eu> Date: Tue, 27 Mar 2018 13:54:15 +0200 Subject: [PATCH] Simpler interface + docs, separation histograms & point data Refs #22197 --- Framework/Algorithms/src/Stitch1D.cpp | 100 ++++++++++++----------- Framework/Algorithms/test/Stitch1DTest.h | 9 +- docs/source/algorithms/Stitch1D-v3.rst | 11 ++- 3 files changed, 63 insertions(+), 57 deletions(-) diff --git a/Framework/Algorithms/src/Stitch1D.cpp b/Framework/Algorithms/src/Stitch1D.cpp index fd605eb1513..edb105e40b4 100644 --- a/Framework/Algorithms/src/Stitch1D.cpp +++ b/Framework/Algorithms/src/Stitch1D.cpp @@ -1,14 +1,12 @@ #include "MantidAlgorithms/Stitch1D.h" #include "MantidAPI/WorkspaceProperty.h" #include "MantidAPI/MatrixWorkspace.h" -#include "MantidAPI/HistogramValidator.h" #include "MantidAPI/WorkspaceFactory.h" #include "MantidHistogramData/HistogramX.h" #include "MantidHistogramData/HistogramY.h" #include "MantidHistogramData/HistogramE.h" #include "MantidKernel/ArrayProperty.h" #include "MantidKernel/BoundedValidator.h" -#include "MantidKernel/ListValidator.h" #include "MantidKernel/MultiThreaded.h" #include "MantidKernel/PropertyWithValue.h" #include "MantidKernel/RebinParamsValidator.h" @@ -118,17 +116,13 @@ void Stitch1D::maskInPlace(int a1, int a2, MatrixWorkspace_sptr source) { /** Initialize the algorithm's properties. */ void Stitch1D::init() { - Kernel::IValidator_sptr histogramValidator = - boost::make_shared<HistogramValidator>(); - declareProperty( - make_unique<WorkspaceProperty<MatrixWorkspace>>( - "LHSWorkspace", "", Direction::Input, histogramValidator->clone()), - "LHS input workspace."); - declareProperty( - make_unique<WorkspaceProperty<MatrixWorkspace>>( - "RHSWorkspace", "", Direction::Input, histogramValidator->clone()), - "RHS input workspace."); + declareProperty(make_unique<WorkspaceProperty<MatrixWorkspace>>( + "LHSWorkspace", "", Direction::Input), + "LHS input workspace."); + declareProperty(make_unique<WorkspaceProperty<MatrixWorkspace>>( + "RHSWorkspace", "", Direction::Input), + "RHS input workspace."); declareProperty(make_unique<WorkspaceProperty<MatrixWorkspace>>( "OutputWorkspace", "", Direction::Output), "Output stitched workspace."); @@ -161,12 +155,6 @@ void Stitch1D::init() { declareProperty(make_unique<PropertyWithValue<double>>( "OutScaleFactor", Mantid::EMPTY_DBL(), Direction::Output), "The actual used value for the scaling factor."); - - const std::vector<std::string> outputXOption{"WeightedMeanOverlap", - "OriginalOverlap"}; - declareProperty("OutputXOption", "WeightedMeanOverlap", - boost::make_shared<StringListValidator>(outputXOption), - "Output option to choose how to overlap x values."); } /**Gets the start of the overlapping region @@ -395,7 +383,7 @@ MatrixWorkspace_sptr Stitch1D::singleValueWS(double val) { return outWS; } -/**finds the bins containing the ends of the overlappign region +/**finds the bins containing the ends of the overlapping region @param startOverlap :: The start of the overlapping region @param endOverlap :: The end of the overlapping region @param workspace :: The workspace to determine the overlaps inside @@ -449,6 +437,7 @@ bool Stitch1D::hasNonzeroErrors(MatrixWorkspace_sptr ws) { void Stitch1D::exec() { MatrixWorkspace_sptr rhsWS = this->getProperty("RHSWorkspace"); MatrixWorkspace_sptr lhsWS = this->getProperty("LHSWorkspace"); + const MinMaxTuple intesectionXRegion = calculateXIntersection(lhsWS, rhsWS); const double intersectionMin = intesectionXRegion.get<0>(); @@ -482,9 +471,9 @@ void Stitch1D::exec() { boost::format("Stitch1D StartOverlap is outside the available X range " "after rebinning. StartOverlap: %10.9f, X min: %10.9f") % startOverlap % xMin); - throw std::runtime_error(message); } + if (endOverlap > xMax) { std::string message = boost::str( boost::format("Stitch1D EndOverlap is outside the available X range " @@ -499,40 +488,40 @@ void Stitch1D::exec() { m_infYIndexes.resize(histogramCount); m_nanEIndexes.resize(histogramCount); m_infEIndexes.resize(histogramCount); - auto rebinnedLHS = rebin(lhsWS, params); - auto rebinnedRHS = rebin(rhsWS, params); - - boost::tuple<int, int> startEnd = - findStartEndIndexes(startOverlap, endOverlap, rebinnedLHS); const bool useManualScaleFactor = this->getProperty("UseManualScaleFactor"); double scaleFactor = 0; double errorScaleFactor = 0; + MatrixWorkspace_sptr lhs, rhs; + + // If the input workspaces are histograms ... + if (rhsWS->isHistogramData() && lhsWS->isHistogramData()) { + lhs = rebin(lhsWS, params); + rhs = rebin(rhsWS, params); + } + if (useManualScaleFactor) { double manualScaleFactor = this->getProperty("ManualScaleFactor"); MatrixWorkspace_sptr manualScaleFactorWS = singleValueWS(manualScaleFactor); - if (scaleRHS) { - rebinnedRHS = rebinnedRHS * manualScaleFactorWS; - } else { - rebinnedLHS = rebinnedLHS * manualScaleFactorWS; - } + if (scaleRHS) + rhs *= manualScaleFactorWS; + else + lhs *= manualScaleFactorWS; scaleFactor = manualScaleFactor; errorScaleFactor = manualScaleFactor; } else { - auto rhsOverlapIntegrated = - integration(rebinnedRHS, startOverlap, endOverlap); - auto lhsOverlapIntegrated = - integration(rebinnedLHS, startOverlap, endOverlap); + auto rhsOverlapIntegrated = integration(rhs, startOverlap, endOverlap); + auto lhsOverlapIntegrated = integration(lhs, startOverlap, endOverlap); MatrixWorkspace_sptr ratio; if (scaleRHS) { ratio = lhsOverlapIntegrated / rhsOverlapIntegrated; - rebinnedRHS = rebinnedRHS * ratio; + rhs = rhs * ratio; } else { ratio = rhsOverlapIntegrated / lhsOverlapIntegrated; - rebinnedLHS = rebinnedLHS * ratio; + lhs *= ratio; } scaleFactor = ratio->y(0).front(); errorScaleFactor = ratio->e(0).front(); @@ -545,29 +534,46 @@ void Stitch1D::exec() { } } + boost::tuple<int, int> startEnd = + findStartEndIndexes(startOverlap, endOverlap, lhs); + int a1 = boost::tuples::get<0>(startEnd); int a2 = boost::tuples::get<1>(startEnd); // Mask out everything BUT the overlap region as a new workspace. - MatrixWorkspace_sptr overlap1 = maskAllBut(a1, a2, rebinnedLHS); + MatrixWorkspace_sptr overlap1 = maskAllBut(a1, a2, lhs); // Mask out everything BUT the overlap region as a new workspace. - MatrixWorkspace_sptr overlap2 = maskAllBut(a1, a2, rebinnedRHS); + MatrixWorkspace_sptr overlap2 = maskAllBut(a1, a2, rhs); // Mask out everything AFTER the overlap region as a new workspace. - maskInPlace(a1 + 1, static_cast<int>(rebinnedLHS->blocksize()), rebinnedLHS); + maskInPlace(a1 + 1, static_cast<int>(lhs->blocksize()), lhs); // Mask out everything BEFORE the overlap region as a new workspace. - maskInPlace(0, a2, rebinnedRHS); + maskInPlace(0, a2, rhs); MatrixWorkspace_sptr overlapave; - if (hasNonzeroErrors(overlap1) && hasNonzeroErrors(overlap2)) { - overlapave = weightedMean(overlap1, overlap2); + + // If the input workspaces are histograms ... + if (rhsWS->isHistogramData() && lhsWS->isHistogramData()) { + if (hasNonzeroErrors(overlap1) && hasNonzeroErrors(overlap2)) { + overlapave = weightedMean(overlap1, overlap2); + } else { + g_log.information("Using un-weighted mean for Stitch1D overlap mean"); + MatrixWorkspace_sptr sum = overlap1 + overlap2; + MatrixWorkspace_sptr denominator = singleValueWS(2.0); + overlapave = sum / denominator; + } } else { - g_log.information("Using un-weighted mean for Stitch1D overlap mean"); - MatrixWorkspace_sptr sum = overlap1 + overlap2; - MatrixWorkspace_sptr denominator = singleValueWS(2.0); - overlapave = sum / denominator; + // If the input workspaces are point data ... + if (!rhsWS->isHistogramData() && !lhsWS->isHistogramData()) { + g_log.error("Missing implementation"); + throw std::invalid_argument("Missing implementation"); + } else { + g_log.error("Input workspaces must be both histograms or point data"); + throw std::invalid_argument( + "Input workspaces must be both histograms or point data"); + } } - MatrixWorkspace_sptr result = rebinnedLHS + overlapave + rebinnedRHS; + MatrixWorkspace_sptr result = lhs + overlapave + rhs; reinsertSpecialValues(result); // Provide log information about the scale factors used in the calculations. diff --git a/Framework/Algorithms/test/Stitch1DTest.h b/Framework/Algorithms/test/Stitch1DTest.h index 4c9d31d090a..3002385367b 100644 --- a/Framework/Algorithms/test/Stitch1DTest.h +++ b/Framework/Algorithms/test/Stitch1DTest.h @@ -236,7 +236,8 @@ public: void test_startoverlap_greater_than_end_overlap_throws() { TSM_ASSERT_THROWS("Should have thrown with StartOverlap < x max", do_stitch1D(this->a, this->b, this->x.back(), - this->x.front(), std::vector<double>(1., 0.2)), + this->x.front(), + std::vector<double>(1., 0.2)), std::runtime_error &); } @@ -244,7 +245,7 @@ public: auto lhs_ws = make_arbitrary_point_ws(); auto rhs_ws = make_arbitrary_histogram_ws(); TSM_ASSERT_THROWS( - "LHS WS must be histogram", + "Both ws must be either histogram or point data", do_stitch1D(lhs_ws, rhs_ws, -1., 1., std::vector<double>(1., 0.2)), std::invalid_argument &); } @@ -253,12 +254,12 @@ public: auto lhs_ws = make_arbitrary_histogram_ws(); auto rhs_ws = make_arbitrary_point_ws(); TSM_ASSERT_THROWS( - "RHS WS must be histogram", + "Both ws must be either histogram or point data", do_stitch1D(lhs_ws, rhs_ws, -1., 1., std::vector<double>(1., 0.2)), std::invalid_argument &); } - void test_stitching_uses_suppiled_params() { + void test_stitching_uses_supplied_params() { std::vector<double> params = {-0.8, 0.2, 1.0}; auto ret = do_stitch1D(this->b, this->a, -0.4, 0.4, params); diff --git a/docs/source/algorithms/Stitch1D-v3.rst b/docs/source/algorithms/Stitch1D-v3.rst index c0125ad3e89..6f451d0427d 100644 --- a/docs/source/algorithms/Stitch1D-v3.rst +++ b/docs/source/algorithms/Stitch1D-v3.rst @@ -18,9 +18,7 @@ Users can optionally provide :ref:`algm-Rebin` :literal:`Params`, otherwise they Likewise, :literal:`StartOverlap` and :literal:`EndOverlap` are optional. If not provided, then these are taken to be the region of X-axis intersection. -The option :literal:`OutputXOption` defines the calculation of x values in the overlap range of the output workspace. -The default option is :literal:`WeightedMeanOverlap`. If the input workspace contains Dx values, the output workspace will not have dx values defined. -Contrarily, the option :literal:`OriginalOverlap` will propagate the corresponding Dx values to the output workspace. +The type of the input workspaces (histogram or point data) define the calculation of x values in the overlap range of the output workspace. The algorithm workflow is as follows: @@ -46,13 +44,14 @@ The algorithm workflow is as follows: #. The weighted mean of the two workspaces in range [:literal:`StartOverlap`, :literal:`EndOverlap`] is calculated. Note that if both workspaces have zero errors, an un-weighted mean will be performed instead. -#. If :literal:`OutputXOption` is :literal:`WeightedMeanOverlap`, the output workspace will be created by summing the left-hand-side workspace (values in range +#. For histograms, the output workspace will be created by summing the left-hand-side workspace (values in range [:literal:`StartX`, :literal:`StartOverlap`], where :literal:`StartX` is the minimum X value specified via :literal:`Params` or calculated from the left-hand-side workspace) + weighted mean workspace + right-hand-side workspace (values in range [:literal:`EndOverlap`, :literal:`EndX`], where :literal:`EndX` is the maximum X value specified via :literal:`Params` or calculated from the right-hand-side workspace) multiplied by the scale factor. - If :literal:`OutputXOption` is :literal:`OriginalOverlap`, the output workspace will propagate x values and Dx values present in the overlap range to the :literal:`OutputWorkspace`. + Dx values will not be present in the output workspace. + For point data, the output workspace will propagate x values and Dx values present in the overlap range to the output workspace. #. The special values are put back in the output workspace. Note that if both the left-hand-side workspace and the right-hand-side workspace happen to have a different special value in the same bin, this bin will be set to infinite in the output workspace. @@ -83,7 +82,7 @@ Usage import numpy as np def gaussian(x, mu, sigma): - """Creates a gaussian peak centered on mu and with width sigma.""" + """Creates a Gaussian peak centered on mu and with width sigma.""" return (1/ sigma * np.sqrt(2 * np.pi)) * np.exp( - (x-mu)**2 / (2*sigma**2)) #create two histograms with a single peak in each one -- GitLab