PageRenderTime 64ms CodeModel.GetById 41ms app.highlight 19ms RepoModel.GetById 1ms app.codeStats 0ms

/Framework/Algorithms/src/CrossCorrelate.cpp

https://github.com/mantidproject/mantid
C++ | 266 lines | 192 code | 34 blank | 40 comment | 20 complexity | 6982abefe82f4f5f074dc6bb47078dfd MD5 | raw file
  1// Mantid Repository : https://github.com/mantidproject/mantid
  2//
  3// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
  4//   NScD Oak Ridge National Laboratory, European Spallation Source,
  5//   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
  6// SPDX - License - Identifier: GPL - 3.0 +
  7//----------------------------------------------------------------------
  8// Includes
  9//----------------------------------------------------------------------
 10#include "MantidAlgorithms/CrossCorrelate.h"
 11#include "MantidAPI/HistogramValidator.h"
 12#include "MantidAPI/NumericAxis.h"
 13#include "MantidAPI/RawCountValidator.h"
 14#include "MantidAPI/SpectraAxis.h"
 15#include "MantidAPI/TextAxis.h"
 16#include "MantidAPI/WorkspaceUnitValidator.h"
 17#include "MantidDataObjects/Workspace2D.h"
 18#include "MantidDataObjects/WorkspaceCreation.h"
 19#include "MantidHistogramData/Histogram.h"
 20#include "MantidKernel/BoundedValidator.h"
 21#include "MantidKernel/CompositeValidator.h"
 22#include "MantidKernel/UnitFactory.h"
 23#include "MantidKernel/VectorHelper.h"
 24#include <boost/iterator/counting_iterator.hpp>
 25#include <numeric>
 26#include <sstream>
 27
 28namespace {
 29struct Variances {
 30  double y;
 31  double e;
 32};
 33
 34Variances subtractMean(std::vector<double> &signal, std::vector<double> &error) {
 35  double mean = std::accumulate(signal.cbegin(), signal.cend(), 0.0);
 36  double errorMeanSquared =
 37      std::accumulate(error.cbegin(), error.cend(), 0.0, Mantid::Kernel::VectorHelper::SumSquares<double>());
 38  const auto n = signal.size();
 39  mean /= static_cast<double>(n);
 40  errorMeanSquared /= static_cast<double>(n * n);
 41  double variance = 0.0, errorVariance = 0.0;
 42  auto itY = signal.begin();
 43  auto itE = error.begin();
 44  for (; itY != signal.end(); ++itY, ++itE) {
 45    (*itY) -= mean;                              // Now the vector is (y[i]-refMean)
 46    (*itE) = (*itE) * (*itE) + errorMeanSquared; // New error squared
 47    const double t = (*itY) * (*itY);            //(y[i]-refMean)^2
 48    variance += t;                               // Sum previous term
 49    errorVariance += 4.0 * t * (*itE);           // Error squared
 50  }
 51  return {variance, errorVariance};
 52}
 53
 54} // namespace
 55
 56namespace Mantid::Algorithms {
 57
 58// Register the class into the algorithm factory
 59DECLARE_ALGORITHM(CrossCorrelate)
 60
 61using namespace Kernel;
 62using namespace API;
 63using namespace DataObjects;
 64using namespace HistogramData;
 65
 66/// Initialisation method.
 67void CrossCorrelate::init() {
 68  auto wsValidator = std::make_shared<CompositeValidator>();
 69  wsValidator->add<API::WorkspaceUnitValidator>("dSpacing");
 70  wsValidator->add<API::HistogramValidator>();
 71  wsValidator->add<API::RawCountValidator>();
 72
 73  // Input and output workspaces
 74  declareProperty(
 75      std::make_unique<WorkspaceProperty<MatrixWorkspace>>("InputWorkspace", "", Direction::Input, wsValidator),
 76      "A 2D workspace with X values of d-spacing");
 77  declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>("OutputWorkspace", "", Direction::Output),
 78                  "The name of the output workspace");
 79
 80  auto mustBePositive = std::make_shared<BoundedValidator<int>>();
 81  mustBePositive->setLower(0);
 82  // Reference spectra against which cross correlation is performed
 83  declareProperty("ReferenceSpectra", 0, mustBePositive,
 84                  "The Workspace Index of the spectra to correlate all other "
 85                  "spectra against. ");
 86  // Spectra in the range [min to max] will be cross correlated to referenceSpectra.
 87  declareProperty("WorkspaceIndexMin", 0, mustBePositive,
 88                  "The workspace index of the first member of the range of "
 89                  "spectra to cross-correlate against.");
 90  declareProperty("WorkspaceIndexMax", 0, mustBePositive,
 91                  " The workspace index of the last member of the range of "
 92                  "spectra to cross-correlate against.");
 93  // Only the data in the range X_min, X_max will be used
 94  declareProperty("XMin", 0.0, "The starting point of the region to be cross correlated.");
 95  declareProperty("XMax", 0.0, "The ending point of the region to be cross correlated.");
 96  // max is .1
 97  declareProperty("MaxDSpaceShift", EMPTY_DBL(), "Optional float for maximum shift to calculate (in d-spacing)");
 98}
 99
100/** Executes the algorithm
101 *
102 *  @throw runtime_error Thrown if algorithm cannot execute
103 */
104void CrossCorrelate::exec() {
105  MatrixWorkspace_const_sptr inputWS = getProperty("InputWorkspace");
106  double maxDSpaceShift = getProperty("MaxDSpaceShift");
107  int referenceSpectra = getProperty("ReferenceSpectra");
108  double xmin = getProperty("XMin");
109  double xmax = getProperty("XMax");
110  const int wsIndexMin = getProperty("WorkspaceIndexMin");
111  const int wsIndexMax = getProperty("WorkspaceIndexMax");
112
113  const auto index_ref = static_cast<size_t>(referenceSpectra);
114
115  if (wsIndexMin >= wsIndexMax)
116    throw std::runtime_error("Must specify WorkspaceIndexMin<WorkspaceIndexMax");
117  // Get the number of spectra in range wsIndexMin to wsIndexMax
118  int numSpectra = 1 + wsIndexMax - wsIndexMin;
119  // Indexes of all spectra in range
120  std::vector<size_t> indexes(boost::make_counting_iterator(wsIndexMin), boost::make_counting_iterator(wsIndexMax + 1));
121
122  if (numSpectra == 0) {
123    std::ostringstream message;
124    message << "No spectra in range between" << wsIndexMin << " and " << wsIndexMax;
125    throw std::runtime_error(message.str());
126  }
127  // Output messageage information
128  g_log.information() << "There are " << numSpectra << " spectra in the range\n";
129
130  // checdataIndex that the data range specified madataIndexes sense
131  if (xmin >= xmax)
132    throw std::runtime_error("Must specify xmin < xmax, " + std::to_string(xmin) + " vs " + std::to_string(xmax));
133
134  // TadataIndexe a copy of  the referenceSpectra spectrum
135  auto &referenceSpectraE = inputWS->e(index_ref);
136  auto &referenceSpectraX = inputWS->x(index_ref);
137  auto &referenceSpectraY = inputWS->y(index_ref);
138  // Now checdataIndex if the range between x_min and x_max is valid
139  using std::placeholders::_1;
140  auto rangeStart =
141      std::find_if(referenceSpectraX.cbegin(), referenceSpectraX.cend(), std::bind(std::greater<double>(), _1, xmin));
142  if (rangeStart == referenceSpectraX.cend())
143    throw std::runtime_error("No data above XMin");
144  auto rangeEnd = std::find_if(rangeStart, referenceSpectraX.cend(), std::bind(std::greater<double>(), _1, xmax));
145  if (rangeStart == rangeEnd)
146    throw std::runtime_error("Range is not valid");
147
148  MantidVec::difference_type rangeStartCorrection = std::distance(referenceSpectraX.cbegin(), rangeStart);
149  MantidVec::difference_type rangeEndCorrection = std::distance(referenceSpectraX.cbegin(), rangeEnd);
150
151  const std::vector<double> referenceXVector(rangeStart, rangeEnd);
152  std::vector<double> referenceYVector(referenceSpectraY.cbegin() + rangeStartCorrection,
153                                       referenceSpectraY.cbegin() + (rangeEndCorrection - 1));
154  std::vector<double> referenceEVector(referenceSpectraE.cbegin() + rangeStartCorrection,
155                                       referenceSpectraE.cbegin() + (rangeEndCorrection - 1));
156
157  g_log.information() << "min max " << referenceXVector.front() << " " << referenceXVector.back() << '\n';
158
159  // Now start the real stuff
160  // Create a 2DWorkspace that will hold the result
161  auto numReferenceY = static_cast<int>(referenceYVector.size());
162
163  // max the shift
164  int shiftCorrection = 0;
165  if (maxDSpaceShift != EMPTY_DBL()) {
166    if (xmax - xmin < maxDSpaceShift)
167      g_log.warning() << "maxDSpaceShift(" << std::to_string(maxDSpaceShift)
168                      << ") is larger than specified range of xmin(" << xmin << ") to xmax(" << xmax
169                      << "), please make it smaller or removed it entirely!"
170                      << "\n";
171
172    // convert dspacing to bins, where maxDSpaceShift is at least 0.1
173    const auto maxBins = std::max(0.0 + maxDSpaceShift * 2, 0.1) / inputWS->getDimension(0)->getBinWidth();
174    // calc range based on max bins
175    shiftCorrection = (int)std::max(0.0, abs((-numReferenceY + 2) - (numReferenceY - 2)) - maxBins) / 2;
176  }
177
178  const int numPoints = 2 * (numReferenceY - shiftCorrection) - 3;
179  if (numPoints < 1)
180    throw std::runtime_error("Range is not valid");
181
182  MatrixWorkspace_sptr out = create<HistoWorkspace>(*inputWS, numSpectra, Points(numPoints));
183
184  const auto referenceVariance = subtractMean(referenceYVector, referenceEVector);
185
186  const double referenceNorm = 1.0 / sqrt(referenceVariance.y);
187  double referenceNormE = 0.5 * pow(referenceNorm, 3) * sqrt(referenceVariance.e);
188
189  // Now copy the other spectra
190  bool isDistribution = inputWS->isDistribution();
191
192  auto &outX = out->mutableX(0);
193  for (int i = 0; i < static_cast<int>(outX.size()); ++i) {
194    outX[i] = static_cast<double>(i - (numReferenceY - shiftCorrection) + 2);
195  }
196
197  // Initialise the progress reporting object
198  m_progress = std::make_unique<Progress>(this, 0.0, 1.0, numSpectra);
199  PARALLEL_FOR_IF(Kernel::threadSafe(*inputWS, *out))
200  for (int currentSpecIndex = 0; currentSpecIndex < numSpectra; ++currentSpecIndex) // Now loop on all spectra
201  {
202    PARALLEL_START_INTERRUPT_REGION
203    size_t wsIndex = indexes[currentSpecIndex]; // Get the ws index from the table
204    // Copy spectra info from input Workspace
205    out->getSpectrum(currentSpecIndex).copyInfoFrom(inputWS->getSpectrum(wsIndex));
206    out->setSharedX(currentSpecIndex, out->sharedX(0));
207    // Get temp referenceSpectras
208    const auto &inputXVector = inputWS->x(wsIndex);
209    const auto &inputYVector = inputWS->y(wsIndex);
210    const auto &inputEVector = inputWS->e(wsIndex);
211    // Copy Y,E data of spec(currentSpecIndex) to temp vector
212    // Now rebin on the grid of referenceSpectra
213    std::vector<double> tempY(numReferenceY);
214    std::vector<double> tempE(numReferenceY);
215
216    VectorHelper::rebin(inputXVector.rawData(), inputYVector.rawData(), inputEVector.rawData(), referenceXVector, tempY,
217                        tempE, isDistribution);
218    const auto tempVar = subtractMean(tempY, tempE);
219
220    // Calculate the normalisation constant
221    const double tempNorm = 1.0 / sqrt(tempVar.y);
222    const double tempNormE = 0.5 * pow(tempNorm, 3) * sqrt(tempVar.e);
223    const double normalisation = referenceNorm * tempNorm;
224    const double normalisationE2 = pow((referenceNorm * tempNormE), 2) + pow((tempNorm * referenceNormE), 2);
225    // Get referenceSpectr to the ouput spectrum
226    auto &outY = out->mutableY(currentSpecIndex);
227    auto &outE = out->mutableE(currentSpecIndex);
228
229    for (int dataIndex = -numReferenceY + 2 + shiftCorrection; dataIndex <= numReferenceY - 2 - shiftCorrection;
230         ++dataIndex) {
231      const int dataIndexP = abs(dataIndex);
232      double val = 0, err2 = 0, x, y, xE, yE;
233      for (int j = numReferenceY - 1 - dataIndexP; j >= 0; --j) {
234        if (dataIndex >= 0) {
235          x = referenceYVector[j];
236          y = tempY[j + dataIndexP];
237          xE = referenceEVector[j];
238          yE = tempE[j + dataIndexP];
239        } else {
240          x = tempY[j];
241          y = referenceYVector[j + dataIndexP];
242          xE = tempE[j];
243          yE = referenceEVector[j + dataIndexP];
244        }
245        val += (x * y);
246        err2 += x * x * yE + y * y * xE;
247      }
248      outY[dataIndex + numReferenceY - shiftCorrection - 2] = (val * normalisation);
249      outE[dataIndex + numReferenceY - shiftCorrection - 2] =
250          sqrt(val * val * normalisationE2 + normalisation * normalisation * err2);
251    }
252    // Update progress information
253    m_progress->report();
254    PARALLEL_END_INTERRUPT_REGION
255  }
256  PARALLEL_CHECK_INTERRUPT_REGION
257
258  out->getAxis(0)->unit() = UnitFactory::Instance().create("Label");
259  Unit_sptr unit = out->getAxis(0)->unit();
260  std::shared_ptr<Units::Label> label = std::dynamic_pointer_cast<Units::Label>(unit);
261  label->setLabel("Bins of Shift", "\\mathbb{Z}");
262
263  setProperty("OutputWorkspace", out);
264}
265
266} // namespace Mantid::Algorithms