CalculateIqt.cpp

#include "MantidAlgorithms/CalculateIqt.h"
#include "MantidAPI/AlgorithmManager.h"
#include "MantidKernel/BoundedValidator.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidKernel/MersenneTwister.h"
#include "MantidHistogramData/HistogramY.h"
#include <math.h>

using namespace Mantid::API;
using namespace Mantid::Kernel;
using namespace Mantid::HistogramData;

namespace {
constexpr int DEFAULT_ITERATIONS = 10;
constexpr int DEFAULT_SEED = 23021997;

MatrixWorkspace_sptr rebin(MatrixWorkspace_sptr workspace,
                           const std::string &params) {
  IAlgorithm_sptr rebinAlgorithm = AlgorithmManager::Instance().create("Rebin");
  rebinAlgorithm->setChild(true);
  rebinAlgorithm->initialize();
  rebinAlgorithm->setProperty("InputWorkspace", workspace);
  rebinAlgorithm->setProperty("Params", params);
  rebinAlgorithm->execute();
  return rebinAlgorithm->getProperty("OutputWorkspace");
}

MatrixWorkspace_sptr integration(MatrixWorkspace_sptr workspace) {
  IAlgorithm_sptr integrationAlgorithm =
      AlgorithmManager::Instance().create("Integration");
  integrationAlgorithm->setChild(true);
  integrationAlgorithm->initialize();
  integrationAlgorithm->setProperty("InputWorkspace", workspace);
  integrationAlgorithm->execute();
  return integrationAlgorithm->getProperty("OutputWorkspace");
}

MatrixWorkspace_sptr convertToPointData(MatrixWorkspace_sptr workspace) {
  IAlgorithm_sptr pointDataAlgorithm =
      AlgorithmManager::Instance().create("ConvertToPointData");
  pointDataAlgorithm->setChild(true);
  pointDataAlgorithm->initialize();
  pointDataAlgorithm->setProperty("InputWorkspace", workspace);
  pointDataAlgorithm->execute();
  return pointDataAlgorithm->getProperty("OutputWorkspace");
}

MatrixWorkspace_sptr extractFFTSpectrum(MatrixWorkspace_sptr workspace) {
  IAlgorithm_sptr FFTAlgorithm =
      AlgorithmManager::Instance().create("ExtractFFTSpectrum");
  FFTAlgorithm->setChild(true);
  FFTAlgorithm->initialize();
  FFTAlgorithm->setProperty("InputWorkspace", workspace);
  FFTAlgorithm->setProperty("FFTPart", 2);
  FFTAlgorithm->execute();
  return FFTAlgorithm->getProperty("OutputWorkspace");
}

MatrixWorkspace_sptr divide(MatrixWorkspace_sptr lhsWorkspace,
                            MatrixWorkspace_sptr rhsWorkspace) {
  IAlgorithm_sptr divideAlgorithm =
      AlgorithmManager::Instance().create("Divide");
  divideAlgorithm->setChild(true);
  divideAlgorithm->initialize();
  divideAlgorithm->setProperty("LHSWorkspace", lhsWorkspace);
  divideAlgorithm->setProperty("RHSWorkspace", rhsWorkspace);
  divideAlgorithm->execute();
  return divideAlgorithm->getProperty("OutputWorkspace");
}

MatrixWorkspace_sptr cropWorkspace(MatrixWorkspace_sptr workspace,
                                   double xMax) {
  IAlgorithm_sptr cropAlgorithm =
      AlgorithmManager::Instance().create("CropWorkspace");
  cropAlgorithm->setChild(true);
  cropAlgorithm->initialize();
  cropAlgorithm->setProperty("InputWorkspace", workspace);
  cropAlgorithm->setProperty("XMax", xMax);
  cropAlgorithm->execute();
  return cropAlgorithm->getProperty("OutputWorkspace");
}

MatrixWorkspace_sptr replaceSpecialValues(MatrixWorkspace_sptr workspace) {
  IAlgorithm_sptr specialValuesAlgorithm =
      AlgorithmManager::Instance().create("ReplaceSpecialValues");
  specialValuesAlgorithm->setChild(true);
  specialValuesAlgorithm->initialize();
  specialValuesAlgorithm->setProperty("InputWorkspace", workspace);
  specialValuesAlgorithm->setProperty("InfinityValue", 0.0);
  specialValuesAlgorithm->setProperty("BigNumberThreshold", 1.0001);
  specialValuesAlgorithm->setProperty("NaNValue", 0.0);
  specialValuesAlgorithm->execute();
  return specialValuesAlgorithm->getProperty("OutputWorkspace");
}

std::string createRebinString(double minimum, double maximum, double width) {
  std::stringstream rebinStream;
  rebinStream.precision(14);
  rebinStream << minimum << ", " << width << ", " << maximum;
  return rebinStream.str();
}

void randomizeRowWithinError(HistogramY &row, const HistogramE &errors,
                             std::function<double(double)> rng) {
  for (auto i = 0u; i < row.size(); ++i) {
    auto randomValue = rng(errors[i]);
    row[i] += randomValue;
  }
}

MatrixWorkspace_sptr
randomizeWorkspaceWithinError(MatrixWorkspace_sptr workspace, const int seed) {
  MersenneTwister mTwister(seed);
  auto rng = [&mTwister](const double error) {
    return mTwister.nextValue(-error, error);
  };
  for (auto i = 0u; i < workspace->getNumberHistograms(); ++i)
    randomizeRowWithinError(workspace->mutableY(i), workspace->e(i), rng);
  return workspace;
}

double standardDeviation(const std::vector<double> &inputValues) {
  double mean = std::accumulate(inputValues.begin(), inputValues.end(), 0.0) /
                inputValues.size();
  double sumOfXMinusMeanSquared = 0;
  for (auto &x : inputValues) {
    sumOfXMinusMeanSquared += pow(x - mean, 2);
  }
  return sqrt(sumOfXMinusMeanSquared / (inputValues.size() - 1));
}

std::vector<double>
standardDeviationArray(const std::vector<std::vector<double>> &yValues) {
  std::vector<double> standardDeviations;
  auto outputSize = yValues[0].size();
  standardDeviations.reserve(outputSize);
  std::vector<double> currentRow;
  currentRow.reserve(yValues.size());

  for (auto i = 0u; i < outputSize; ++i) {
    currentRow.clear();
    for (auto &yValueArray : yValues)
      currentRow.emplace_back(yValueArray[i]);
    standardDeviations.emplace_back(standardDeviation(currentRow));
  }
  return standardDeviations;
}

MatrixWorkspace_sptr cleanOutput(MatrixWorkspace_sptr workspace) {
  auto binning = static_cast<int>(std::ceil(workspace->blocksize() / 2));
  auto binV = workspace->x(0)[binning];
  workspace = cropWorkspace(workspace, binV);
  workspace = replaceSpecialValues(workspace);
  return workspace;
}

MatrixWorkspace_sptr
normalizedFourierTransform(MatrixWorkspace_sptr workspace,
                           const std::string &rebinParams) {
  workspace = rebin(workspace, rebinParams);
  auto workspace_int = integration(workspace);
  workspace = convertToPointData(workspace);
  workspace = extractFFTSpectrum(workspace);
  workspace = divide(workspace, workspace_int);
  return workspace;
}

MatrixWorkspace_sptr calculateIqt(MatrixWorkspace_sptr workspace,
                                  MatrixWorkspace_sptr resolutionWorkspace,
                                  const std::string &rebinParams) {
  workspace = normalizedFourierTransform(workspace, rebinParams);
  return divide(workspace, resolutionWorkspace);
}

} // namespace

namespace Mantid {
namespace Algorithms {

DECLARE_ALGORITHM(CalculateIqt)

const std::string CalculateIqt::name() const { return "CalculateIqt"; }

int CalculateIqt::version() const { return 1; }

const std::vector<std::string> CalculateIqt::seeAlso() const {
  return {"TransformToIqt"};
}

const std::string CalculateIqt::category() const {
  return "Inelastic\\Indirect";
}

const std::string CalculateIqt::summary() const {
  return "Calculates I(Q,t) from S(Q,w) and computes the errors using a "
         "monte-carlo routine.";
}

void CalculateIqt::init() {

  declareProperty(
      make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input),
      "The name of the sample workspace.");
  declareProperty(make_unique<WorkspaceProperty<>>("ResolutionWorkspace", "",
                                                   Direction::Input),
                  "The name of the resolution workspace.");

  declareProperty("EnergyMin", -0.5,
                  "Minimum energy for fit. Default = -.0.5.");
  declareProperty("EnergyMax", 0.5, "Maximum energy for fit. Default = 0.5.");
  declareProperty("EnergyWidth", 0.1, "Width of energy bins for fit.");

  auto positiveInt = boost::make_shared<Kernel::BoundedValidator<int>>();
  positiveInt->setLower(1);

  declareProperty("NumberOfIterations", DEFAULT_ITERATIONS, positiveInt,
                  "Number of randomised simulations within error to run.");
  declareProperty(
      "SeedValue", DEFAULT_SEED, positiveInt,
      "Seed the random number generator for monte-carlo error calculation.");

  declareProperty(make_unique<WorkspaceProperty<>>("OutputWorkspace", "",
                                                   Direction::Output),
                  "The name to use for the output workspace.");
}

void CalculateIqt::exec() {
  const auto rebinParams = rebinParamsAsString();
  const MatrixWorkspace_sptr sampleWorkspace = getProperty("InputWorkspace");
  MatrixWorkspace_sptr resolution = getProperty("ResolutionWorkspace");
  resolution = normalizedFourierTransform(resolution, rebinParams);

  // create function which only needs sampleWorkspace as input so it can be
  // easily repeated in simulations
  auto calculateIqtFunction =
      [resolution, &rebinParams](MatrixWorkspace_sptr sample) {
        return calculateIqt(sample, resolution, rebinParams);
      };

  auto outputWorkspace = monteCarloErrorCalculation(sampleWorkspace, resolution,
                                                    calculateIqtFunction);

  outputWorkspace = cleanOutput(outputWorkspace);
  setProperty("OutputWorkspace", outputWorkspace);
}

std::string CalculateIqt::rebinParamsAsString() {
  const double e_min = getProperty("EnergyMin");
  const double e_max = getProperty("EnergyMax");
  const double e_width = getProperty("EnergyWidth");
  return createRebinString(e_min, e_max, e_width);
}

MatrixWorkspace_sptr CalculateIqt::monteCarloErrorCalculation(
    MatrixWorkspace_sptr sample, MatrixWorkspace_sptr resolution,
    const std::function<MatrixWorkspace_sptr(MatrixWorkspace_sptr)> &
        calculateIqtFunction) {
  auto outputWorkspace = calculateIqtFunction(sample);
  const unsigned int nIterations = getProperty("NumberOfIterations");
  const unsigned int seed = getProperty("SeedValue");
  std::vector<MatrixWorkspace_sptr> simulatedWorkspaces;
  simulatedWorkspaces.reserve(nIterations);
  simulatedWorkspaces.emplace_back(outputWorkspace);

  PARALLEL_FOR_IF(Kernel::threadSafe(*sample, *resolution))
  for (auto i = 0; i < nIterations - 1; ++i) {
    PARALLEL_START_INTERUPT_REGION
    auto simulatedWorkspace =
        randomizeWorkspaceWithinError(sample->clone(), seed);
    simulatedWorkspace = calculateIqtFunction(simulatedWorkspace);
    simulatedWorkspaces.emplace_back(simulatedWorkspace);
    PARALLEL_END_INTERUPT_REGION
  }
  PARALLEL_CHECK_INTERUPT_REGION
  return setErrorsToStandardDeviation(nIterations, simulatedWorkspaces,
                                      outputWorkspace);
}

MatrixWorkspace_sptr CalculateIqt::setErrorsToStandardDeviation(
    int nIterations,
    const std::vector<MatrixWorkspace_sptr> simulatedWorkspaces,
    MatrixWorkspace_sptr outputWorkspace) {
  // set errors to standard deviation of y values across simulations
  std::vector<std::vector<double>> allSimY;
  allSimY.reserve(nIterations);

  for (auto i = 0u; i < outputWorkspace->getNumberHistograms(); ++i) {
    auto &outputError = outputWorkspace->mutableE(i);
    for (auto &simWorkspace : simulatedWorkspaces)
      allSimY.emplace_back(simWorkspace->y(i).rawData());
    outputError = standardDeviationArray(allSimY);
  }
  return outputWorkspace;
}

std::map<std::string, std::string> CalculateIqt::validateInputs() {
  std::map<std::string, std::string> issues;
  const double eMin = getProperty("EnergyMin");
  const double eMax = getProperty("EnergyMax");
  if (eMin > eMax) {
    auto energy_swapped = "EnergyMin is greater than EnergyMax";
    issues["EnergyMin"] = energy_swapped;
    issues["EnergyMax"] = energy_swapped;
  }
  return issues;
}

} // namespace Algorithms
} // namespace Mantid