Commit 71901b4e authored by Mathieu Tillet's avatar Mathieu Tillet
Browse files

Add an inverse power rebinning method

Allow the user to rebin with bins that becomes increasingly smaller.
parent d12a0b4f
...@@ -56,6 +56,7 @@ public: ...@@ -56,6 +56,7 @@ public:
const std::vector<std::string> seeAlso() const override { const std::vector<std::string> seeAlso() const override {
return {"RebinToWorkspace", "Rebin2D", "Rebunch", "Regroup", "RebinByPulseTimes", "RebinByTimeAtSample"}; return {"RebinToWorkspace", "Rebin2D", "Rebunch", "Regroup", "RebinByPulseTimes", "RebinByTimeAtSample"};
} }
std::map<std::string, std::string> validateInputs() override;
static std::vector<double> rebinParamsFromInput(const std::vector<double> &inParams, static std::vector<double> rebinParamsFromInput(const std::vector<double> &inParams,
const API::MatrixWorkspace &inputWS, Kernel::Logger &logger); const API::MatrixWorkspace &inputWS, Kernel::Logger &logger);
......
...@@ -78,6 +78,56 @@ std::vector<double> Rebin::rebinParamsFromInput(const std::vector<double> &inPar ...@@ -78,6 +78,56 @@ std::vector<double> Rebin::rebinParamsFromInput(const std::vector<double> &inPar
// Public methods // Public methods
//--------------------------------------------------------------------------------------------- //---------------------------------------------------------------------------------------------
/// Validate that the input properties are sane.
std::map<std::string, std::string> Rebin::validateInputs() {
std::map<std::string, std::string> helpMessages;
if (existsProperty("Power") && !isDefault("Power")) {
const double power = getProperty("Power");
if (power < 0) {
helpMessages["Power"] = "Power cannot be negative";
return helpMessages;
}
if (power > 1) {
helpMessages["Power"] = "Power cannot be more than one.";
return helpMessages;
}
// attempt to roughly guess how many bins these parameters imply
double roughEstimate = 0;
if (!isDefault("Params")) {
const std::vector<double> params = getProperty("Params");
// Five significant places of the Euler-Mascheroni constant is probably more than enough for our needs
double eulerMascheroni = 0.57721;
// Params is check by the validator first, so we can assume it is in a correct format
for (size_t i = 0; i < params.size() - 2; i += 2) {
double upperLimit = params[i + 2];
double lowerLimit = params[i];
double factor = params[i + 1];
if (factor <= 0) {
helpMessages["Params"] = "Provided width value cannot be negative for inverse power binning.";
return helpMessages;
}
if (power == 1) {
roughEstimate += std::exp((upperLimit - lowerLimit) / factor - eulerMascheroni);
} else {
roughEstimate += std::pow(((upperLimit - lowerLimit) / factor) * (1 - power) + 1, 1 / (1 - power));
}
}
}
// Prevent the user form creating too many bins
if (roughEstimate > 10000) {
helpMessages["Power"] = "This binning is expected to give more than 1000 bins.";
}
}
return helpMessages;
}
/** Initialisation method. Declares properties to be used in algorithm. /** Initialisation method. Declares properties to be used in algorithm.
* *
*/ */
...@@ -110,6 +160,9 @@ void Rebin::init() { ...@@ -110,6 +160,9 @@ void Rebin::init() {
"For logarithmic intervals, the splitting starts from the end and go back to the start, ie the bins " "For logarithmic intervals, the splitting starts from the end and go back to the start, ie the bins "
"are bigger at the start and exponentially reduces until they reach the end. For these bins, the " "are bigger at the start and exponentially reduces until they reach the end. For these bins, the "
"FullBinsOnly flag is ignored."); "FullBinsOnly flag is ignored.");
declareProperty("Power", EMPTY_DBL(),
"Splits the interval in bins whose width is width / (i ^ power). Power must be between 0 and 1.");
} }
/** Executes the rebin algorithm /** Executes the rebin algorithm
...@@ -138,6 +191,7 @@ void Rebin::exec() { ...@@ -138,6 +191,7 @@ void Rebin::exec() {
bool fullBinsOnly = getProperty("FullBinsOnly"); bool fullBinsOnly = getProperty("FullBinsOnly");
bool useReverseLog = getProperty("UseReverseLogarithmic"); bool useReverseLog = getProperty("UseReverseLogarithmic");
double power = getProperty("Power");
double xmin = 0.; double xmin = 0.;
double xmax = 0.; double xmax = 0.;
...@@ -146,7 +200,7 @@ void Rebin::exec() { ...@@ -146,7 +200,7 @@ void Rebin::exec() {
HistogramData::BinEdges XValues_new(0); HistogramData::BinEdges XValues_new(0);
// create new output X axis // create new output X axis
static_cast<void>(VectorHelper::createAxisFromRebinParams(rbParams, XValues_new.mutableRawData(), true, fullBinsOnly, static_cast<void>(VectorHelper::createAxisFromRebinParams(rbParams, XValues_new.mutableRawData(), true, fullBinsOnly,
xmin, xmax, useReverseLog)); xmin, xmax, useReverseLog, power));
// Now, determine if the input workspace is actually an EventWorkspace // Now, determine if the input workspace is actually an EventWorkspace
EventWorkspace_const_sptr eventInputWS = std::dynamic_pointer_cast<const EventWorkspace>(inputWS); EventWorkspace_const_sptr eventInputWS = std::dynamic_pointer_cast<const EventWorkspace>(inputWS);
......
...@@ -689,6 +689,93 @@ public: ...@@ -689,6 +689,93 @@ public:
AnalysisDataService::Instance().remove("test_Rebin_revLog"); AnalysisDataService::Instance().remove("test_Rebin_revLog");
} }
void test_inversePowerSquareRoot() {
// Test InversePower in a simple case of square root sum
Workspace2D_sptr test_1D = Create1DWorkspace(51);
test_1D->setDistribution(false);
AnalysisDataService::Instance().add("test_Rebin_revLog", test_1D);
Rebin rebin;
rebin.initialize();
rebin.setPropertyValue("InputWorkspace", "test_Rebin_revLog");
rebin.setPropertyValue("OutputWorkspace", "test_Rebin_revLog");
rebin.setPropertyValue("Params", "1, 1, 10");
rebin.setPropertyValue("Power", "0.5");
TS_ASSERT_THROWS_NOTHING(rebin.execute());
MatrixWorkspace_sptr out = AnalysisDataService::Instance().retrieveWS<MatrixWorkspace>("test_Rebin_revLog");
auto &outX = out->x(0);
TS_ASSERT_EQUALS(outX.size(), 28);
TS_ASSERT_DELTA(outX[0], 1, 1e-5);
TS_ASSERT_DELTA(outX[1], 2, 1e-5);
TS_ASSERT_DELTA(outX[2], 2.707106781, 1e-5);
TS_ASSERT_DELTA(outX[3], 3.28445705, 1e-5);
TS_ASSERT_DELTA(outX[27], 10, 1e-5);
AnalysisDataService::Instance().remove("test_Rebin_revLog");
}
void test_inversePowerHarmonic() {
// Test InversePower in a simple case of harmonic serie
Workspace2D_sptr test_1D = Create1DWorkspace(51);
test_1D->setDistribution(false);
AnalysisDataService::Instance().add("test_Rebin_revLog", test_1D);
Rebin rebin;
rebin.initialize();
rebin.setPropertyValue("InputWorkspace", "test_Rebin_revLog");
rebin.setPropertyValue("OutputWorkspace", "test_Rebin_revLog");
rebin.setPropertyValue("Params", "1, 1, 5");
rebin.setPropertyValue("Power", "1");
TS_ASSERT_THROWS_NOTHING(rebin.execute());
MatrixWorkspace_sptr out = AnalysisDataService::Instance().retrieveWS<MatrixWorkspace>("test_Rebin_revLog");
auto &outX = out->x(0);
TS_ASSERT_EQUALS(outX.size(), 31);
TS_ASSERT_DELTA(outX[0], 1, 1e-5);
TS_ASSERT_DELTA(outX[1], 2, 1e-5);
TS_ASSERT_DELTA(outX[2], 2.5, 1e-5);
TS_ASSERT_DELTA(outX[3], 2.8333333, 1e-5);
TS_ASSERT_DELTA(outX[30], 5, 1e-5);
AnalysisDataService::Instance().remove("test_Rebin_revLog");
}
void test_inversePowerValidateHarmonic() {
// Test that the validator which forbid breating more than 10000 bins works in a harmonic series case
Workspace2D_sptr test_1D = Create1DWorkspace(51);
test_1D->setDistribution(false);
AnalysisDataService::Instance().add("test_Rebin_revLog", test_1D);
Rebin rebin;
rebin.initialize();
rebin.setPropertyValue("InputWorkspace", "test_Rebin_revLog");
rebin.setPropertyValue("OutputWorkspace", "test_Rebin_revLog");
rebin.setPropertyValue("Params", "1, 1, 100");
rebin.setPropertyValue("Power", "1");
TS_ASSERT_THROWS(rebin.execute(), const std::runtime_error &);
AnalysisDataService::Instance().remove("test_Rebin_revLog");
}
void test_inversePowerValidateInverseSquareRoot() {
// Test that the validator which forbid breating more than 10000 bins works in an inverse square root case
// We test both because they rely on different formula to compute the expected number of bins.
Workspace2D_sptr test_1D = Create1DWorkspace(51);
test_1D->setDistribution(false);
AnalysisDataService::Instance().add("test_Rebin_revLog", test_1D);
Rebin rebin;
rebin.initialize();
rebin.setPropertyValue("InputWorkspace", "test_Rebin_revLog");
rebin.setPropertyValue("OutputWorkspace", "test_Rebin_revLog");
rebin.setPropertyValue("Params", "1, 1, 1000");
rebin.setPropertyValue("Power", "0.5");
TS_ASSERT_THROWS(rebin.execute(), const std::runtime_error &);
AnalysisDataService::Instance().remove("test_Rebin_revLog");
}
void test_parallel_cloned() { ParallelTestHelpers::runParallel(run_rebin, "Parallel::StorageMode::Cloned"); } void test_parallel_cloned() { ParallelTestHelpers::runParallel(run_rebin, "Parallel::StorageMode::Cloned"); }
void test_parallel_distributed() { void test_parallel_distributed() {
......
...@@ -28,7 +28,7 @@ int MANTID_KERNEL_DLL createAxisFromRebinParams(const std::vector<double> &param ...@@ -28,7 +28,7 @@ int MANTID_KERNEL_DLL createAxisFromRebinParams(const std::vector<double> &param
const bool resize_xnew = true, const bool full_bins_only = false, const bool resize_xnew = true, const bool full_bins_only = false,
const double xMinHint = std::nan(""), const double xMinHint = std::nan(""),
const double xMaxHint = std::nan(""), const double xMaxHint = std::nan(""),
const bool useReverseLogarithmic = false); const bool useReverseLogarithmic = false, const double power = -1);
void MANTID_KERNEL_DLL rebin(const std::vector<double> &xold, const std::vector<double> &yold, void MANTID_KERNEL_DLL rebin(const std::vector<double> &xold, const std::vector<double> &yold,
const std::vector<double> &eold, const std::vector<double> &xnew, const std::vector<double> &eold, const std::vector<double> &xnew,
......
...@@ -10,6 +10,7 @@ ...@@ -10,6 +10,7 @@
#include "MantidKernel/VectorHelper.h" #include "MantidKernel/VectorHelper.h"
#include <algorithm> #include <algorithm>
#include <boost/algorithm/string.hpp> #include <boost/algorithm/string.hpp>
#include <iostream>
#include <numeric> #include <numeric>
#include <sstream> #include <sstream>
...@@ -34,7 +35,7 @@ namespace VectorHelper { ...@@ -34,7 +35,7 @@ namespace VectorHelper {
**/ **/
int DLLExport createAxisFromRebinParams(const std::vector<double> &params, std::vector<double> &xnew, int DLLExport createAxisFromRebinParams(const std::vector<double> &params, std::vector<double> &xnew,
const bool resize_xnew, const bool full_bins_only, const double xMinHint, const bool resize_xnew, const bool full_bins_only, const double xMinHint,
const double xMaxHint, const bool useReverseLogarithmic) { const double xMaxHint, const bool useReverseLogarithmic, const double power) {
std::vector<double> tmp; std::vector<double> tmp;
const std::vector<double> &fullParams = [&params, &tmp, xMinHint, xMaxHint]() { const std::vector<double> &fullParams = [&params, &tmp, xMinHint, xMaxHint]() {
if (params.size() == 1) { if (params.size() == 1) {
...@@ -70,6 +71,10 @@ int DLLExport createAxisFromRebinParams(const std::vector<double> &params, std:: ...@@ -70,6 +71,10 @@ int DLLExport createAxisFromRebinParams(const std::vector<double> &params, std::
if (resize_xnew) if (resize_xnew)
xnew.emplace_back(xcurr); xnew.emplace_back(xcurr);
int currDiv = 1;
bool isPower = power > 0 && power <= 1;
while ((ibound <= ibounds) && (istep <= isteps)) { while ((ibound <= ibounds) && (istep <= isteps)) {
// if step is negative then it is logarithmic step // if step is negative then it is logarithmic step
bool isLogBin = (fullParams[istep] < 0.0); bool isLogBin = (fullParams[istep] < 0.0);
...@@ -80,20 +85,24 @@ int DLLExport createAxisFromRebinParams(const std::vector<double> &params, std:: ...@@ -80,20 +85,24 @@ int DLLExport createAxisFromRebinParams(const std::vector<double> &params, std::
// we are starting a new bin, but since it is a rev log, xcurr needs to be at its end // we are starting a new bin, but since it is a rev log, xcurr needs to be at its end
xcurr = fullParams[ibound]; xcurr = fullParams[ibound];
} }
if (!isPower) {
if (!isLogBin) if (!isLogBin)
xs = fullParams[istep]; xs = fullParams[istep];
else { else {
if (useReverseLogarithmic) { if (useReverseLogarithmic) {
// we go through a reverse log bin by starting from its end, and working our way back to the beginning // we go through a reverse log bin by starting from its end, and working our way back to the beginning
// this way we can define the bins in a reccuring way, and with a more obvious closeness with the usual log. // this way we can define the bins in a reccuring way, and with a more obvious closeness with the usual log.
double x0 = fullParams[ibound - 2]; double x0 = fullParams[ibound - 2];
double step = x0 + fullParams[ibound] - xcurr; double step = x0 + fullParams[ibound] - xcurr;
xs = -step * alpha; xs = -step * alpha;
} else } else
xs = xcurr * fabs(fullParams[istep]); xs = xcurr * fabs(fullParams[istep]);
}
} else {
xs = fullParams[istep] * std::pow(currDiv, -power);
++currDiv;
} }
if (fabs(xs) == 0.0) { if (fabs(xs) == 0.0) {
......
...@@ -94,6 +94,19 @@ Hence the actual *Param* string used is "0, 2, 4, 3, 10". ...@@ -94,6 +94,19 @@ Hence the actual *Param* string used is "0, 2, 4, 3, 10".
This flag is ignored if UseReverseLogarithm is checked. This flag is ignored if UseReverseLogarithm is checked.
Power option
############
If a value between 0 (excluded) and 1 (included) is provided in the Power field, the binning will follow an inverse power
pattern, each bin having a width of
.. math:: w_i = \frac{F}{i^{\mathrm{power}}
where F is the factor provided between the boundaries.
Since, even though these series diverge and will reach whatever bounds are given, they might take an exponentially slow time
to reach them, and thus an exponential number of bins, a check ensure that the total number of bins is not greater
than 10000.
.. _rebin-usage: .. _rebin-usage:
Usage Usage
...@@ -161,6 +174,26 @@ Output: ...@@ -161,6 +174,26 @@ Output:
The rebinned X values are: [ 1. 6. 8. 9.] The rebinned X values are: [ 1. 6. 8. 9.]
**Example - Inverse power rebinning:**
.. testcode:: ExInversePower
# create histogram workspace
dataX = [1,2,3,4,5,6,7,8,9,10] # or use dataX=range(1,11)
dataY = [1,2,3,4,5,6,7,8,9] # or use dataY=range(1,10)
ws = CreateWorkspace(dataX, dataY)
# rebin from min to max - 1 with square root
ws = Rebin(ws, "1, 3, 10", Power=0.5)
print("The rebinned X values are: {}".format(ws.readX(0)))
Output:
.. testoutput:: ExHistRevLog
The rebinned X values are: [ 1. 4. 6.121320344 7.853371151 9.353371151 10.]
**Example - custom two regions rebinning:** **Example - custom two regions rebinning:**
.. testcode:: ExHistCustom .. testcode:: ExHistCustom
......
...@@ -21,7 +21,7 @@ Improvements ...@@ -21,7 +21,7 @@ Improvements
- :ref:`CreateSampleWorkspace <algm-CreateSampleWorkspace>` has new property InstrumentName. - :ref:`CreateSampleWorkspace <algm-CreateSampleWorkspace>` has new property InstrumentName.
- Event nexuses produced at ILL can now be loaded using :ref:`LoadEventNexus <algm-LoadEventNexus>`. - Event nexuses produced at ILL can now be loaded using :ref:`LoadEventNexus <algm-LoadEventNexus>`.
- :ref:`Rebin <algm-Rebin>` now support reverse logarithmic binning. - :ref:`Rebin <algm-Rebin>` now support reverse logarithmic and inverse power binning.
Bugfixes Bugfixes
######## ########
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment