TabulatedFunction.cpp 10.5 KB
Newer Older
1
2
3
// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4
5
//   NScD Oak Ridge National Laboratory, European Spallation Source,
//   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6
// SPDX - License - Identifier: GPL - 3.0 +
7
8
9
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
10
#include "MantidCurveFitting/Functions/TabulatedFunction.h"
11
12
#include "MantidAPI/Algorithm.h"
#include "MantidAPI/AnalysisDataService.h"
13
14
#include "MantidAPI/FunctionFactory.h"
#include "MantidAPI/MatrixWorkspace.h"
LamarMoore's avatar
LamarMoore committed
15
#include "MantidKernel/FileValidator.h"
16

17
#include <algorithm>
18
19
20
#include <cmath>
#include <fstream>
#include <sstream>
David Fairbrother's avatar
David Fairbrother committed
21
#include <utility>
22

23
24
namespace Mantid {
namespace CurveFitting {
25
26
27
28
namespace Functions {

using namespace CurveFitting;

29
using namespace Kernel;
30

31
32
33
34
35
using namespace API;

namespace {
/// static logger
Logger g_log("TabulatedFunction");
LamarMoore's avatar
LamarMoore committed
36
} // namespace
37
38
39

DECLARE_FUNCTION(TabulatedFunction)

40
41
const int TabulatedFunction::defaultIndexValue = 0;

42
/// Constructor
43
44
TabulatedFunction::TabulatedFunction()
    : m_setupFinished(false), m_explicitXY(false) {
45
  declareParameter("Scaling", 1.0, "A scaling factor");
46
  declareParameter("Shift", 0.0, "Shift in the abscissa");
47
  declareParameter("XScaling", 1.0, "Scaling factor in X");
48
49
50
  declareAttribute("FileName", Attribute("", true));
  declareAttribute("Workspace", Attribute(""));
  declareAttribute("WorkspaceIndex", Attribute(defaultIndexValue));
51
52
  declareAttribute("X", Attribute(std::vector<double>()));
  declareAttribute("Y", Attribute(std::vector<double>()));
53
54
55
}

/// Evaluate the function for a list of arguments and given scaling factor
56
57
58
void TabulatedFunction::eval(double scaling, double xshift, double xscale,
                             double *out, const double *xValues,
                             const size_t nData) const {
59
60
  if (nData == 0)
    return;
61
62
  setupData();

63
64
  if (size() == 0)
    return;
65

66
  // shift and scale the domain over which the function is defined
67
  std::vector<double> xData(m_xData);
Hahn, Steven's avatar
Hahn, Steven committed
68
69
70
  for (double &value : xData) {
    value *= xscale;
    value += xshift;
71
72
73
  }
  const double xStart = xData.front();
  const double xEnd = xData.back();
74

75
76
  if (xStart >= xValues[nData - 1] || xEnd <= xValues[0])
    return;
77
  size_t i = 0;
78
  while (i < nData - 1 && xValues[i] < xStart) {
79
80
81
82
    out[i] = 0;
    i++;
  }
  size_t j = 0;
83
  for (; i < nData; i++) {
84
85
86
87
    double xi = xValues[i];
    while (j < size() - 1 && xi > xData[j])
      j++;
    if (j > size() - 1) {
88
      out[i] = 0;
89
90
    } else {
      if (xi == xData[j]) {
91
        out[i] = m_yData[j] * scaling;
92
      } else if (xi > xData[j]) {
93
        out[i] = 0;
94
95
      } else if (j > 0) {
        double x0 = xData[j - 1];
96
        double x1 = xData[j];
97
        double y0 = m_yData[j - 1];
98
        double y1 = m_yData[j];
99
        out[i] = y0 + (y1 - y0) * (xi - x0) / (x1 - x0);
100
        out[i] *= scaling;
101
      } else {
102
103
104
105
106
107
108
109
110
111
112
113
        out[i] = 0;
      }
    }
  }
}

/**
 * Calculate the function values.
 * @param out :: The output buffer for the calculated values.
 * @param xValues :: The array of x-values.
 * @param nData :: The size of the data.
 */
114
115
void TabulatedFunction::function1D(double *out, const double *xValues,
                                   const size_t nData) const {
116
117
  const double scaling = getParameter("Scaling");
  const double xshift = getParameter("Shift");
118
119
  const double xscale = getParameter("XScaling");
  eval(scaling, xshift, xscale, out, xValues, nData);
120
121
122
123
}

/**
 * function derivatives
124
125
 * @param out :: The output Jacobian matrix: function derivatives over its
 * parameters.
126
127
128
 * @param xValues :: The function arguments
 * @param nData :: The size of xValues.
 */
129
130
131
void TabulatedFunction::functionDeriv1D(API::Jacobian *out,
                                        const double *xValues,
                                        const size_t nData) {
132
133
  const double scaling = getParameter("Scaling");
  const double xshift = getParameter("Shift");
134
  const double xscale = getParameter("XScaling");
135
  std::vector<double> tmp(nData);
136
  // derivative with respect to Scaling parameter
137
  eval(1.0, xshift, xscale, tmp.data(), xValues, nData);
138
139
  for (size_t i = 0; i < nData; ++i) {
    out->set(i, 0, tmp[i]);
140
  }
141

142
143
144
145
  const double dx =
      (xValues[nData - 1] - xValues[0]) / static_cast<double>(nData);
  std::vector<double> tmpplus(nData);
  std::vector<double> tmpminus(nData);
146
147
148
149
150

  // There is no unique definition for the partial derivative with respect
  // to the Shift parameter. Here we take the central difference,
  eval(scaling, xshift + dx, xscale, tmpplus.data(), xValues, nData);
  eval(scaling, xshift - dx, xscale, tmpminus.data(), xValues, nData);
151
152
  for (size_t i = 0; i < nData; ++i) {
    out->set(i, 1, (tmpplus[i] - tmpminus[i]) / (2 * dx));
153
  }
154
155
156
157
158
159

  eval(scaling, xshift, xscale + dx, tmpplus.data(), xValues, nData);
  eval(scaling, xshift, xscale - dx, tmpminus.data(), xValues, nData);
  for (size_t i = 0; i < nData; ++i) {
    out->set(i, 2, (tmpplus[i] - tmpminus[i]) / (2 * dx));
  }
160
161
162
}

/// Clear all data
163
void TabulatedFunction::clear() const {
164
165
  m_xData.clear();
  m_yData.clear();
166
  m_setupFinished = false;
167
168
169
170
171
172
}

/** Set a value to attribute attName
 * @param attName :: The attribute name
 * @param value :: The new value
 */
173
174
175
void TabulatedFunction::setAttribute(const std::string &attName,
                                     const IFunction::Attribute &value) {
  if (attName == "FileName") {
176
    std::string fileName = value.asUnquotedString();
177
178
    if (fileName.empty()) {
      storeAttributeValue("FileName", Attribute("", true));
179
180
      return;
    }
181
182
    FileValidator fval;
    std::string error = fval.isValid(fileName);
183
    if (error.empty()) {
184
185
186
187
188
      storeAttributeValue(attName, Attribute(fileName, true));
      storeAttributeValue("Workspace", Attribute(""));
    } else {
      // file not found
      throw Kernel::Exception::FileError(error, fileName);
189
    }
190
    load(fileName);
191
    m_setupFinished = false;
192
    m_explicitXY = false;
193
  } else if (attName == "Workspace") {
194
    std::string wsName = value.asString();
195
196
197
198
    if (!wsName.empty()) {
      storeAttributeValue(attName, value);
      storeAttributeValue("FileName", Attribute("", true));
      loadWorkspace(wsName);
199
      m_setupFinished = false;
200
      m_explicitXY = false;
201
    }
202
  } else if (attName == "X") {
203
204
205
206
207
208
    m_xData = value.asVector();
    if (m_xData.empty()) {
      m_setupFinished = false;
      m_explicitXY = false;
      if (!m_yData.empty()) {
        m_yData.clear();
209
      }
210
211
212
213
214
215
216
217
218
      return;
    }
    if (m_xData.size() != m_yData.size()) {
      m_yData.resize(m_xData.size());
    }
    storeAttributeValue("FileName", Attribute("", true));
    storeAttributeValue("Workspace", Attribute(""));
    m_setupFinished = true;
    m_explicitXY = true;
219
  } else if (attName == "Y") {
220
221
222
223
224
225
    m_yData = value.asVector();
    if (m_yData.empty()) {
      m_setupFinished = false;
      m_explicitXY = false;
      if (!m_xData.empty()) {
        m_xData.clear();
226
      }
227
228
229
230
231
232
233
234
235
      return;
    }
    if (m_xData.size() != m_yData.size()) {
      m_xData.resize(m_yData.size());
    }
    storeAttributeValue("FileName", Attribute("", true));
    storeAttributeValue("Workspace", Attribute(""));
    m_setupFinished = true;
    m_explicitXY = true;
236
237
  } else {
    IFunction::setAttribute(attName, value);
238
    m_setupFinished = false;
239
240
241
  }
}

242
/// Returns the number of attributes associated with the function
243
size_t TabulatedFunction::nAttributes() const {
244
  // additional X and Y attributes
245
  return IFunction::nAttributes();
246
247
248
}

/// Returns a list of attribute names
249
std::vector<std::string> TabulatedFunction::getAttributeNames() const {
Stephen's avatar
Stephen committed
250
  return IFunction::getAttributeNames();
251
252
253
254
}

/// Return a value of attribute attName
/// @param attName :: The attribute name
255
256
IFunction::Attribute
TabulatedFunction::getAttribute(const std::string &attName) const {
257
  if (attName == "X") {
258
    return m_explicitXY ? Attribute(m_xData) : Attribute(std::vector<double>());
259
  } else if (attName == "Y") {
260
    return m_explicitXY ? Attribute(m_yData) : Attribute(std::vector<double>());
261
262
263
  }
  return IFunction::getAttribute(attName);
}
264
265
266
267
/**
 * Load input file as a Nexus file.
 * @param fname :: The file name
 */
268
269
270
void TabulatedFunction::load(const std::string &fname) {
  IAlgorithm_sptr loadAlg =
      Mantid::API::AlgorithmFactory::Instance().create("Load", -1);
271
272
273
  loadAlg->initialize();
  loadAlg->setChild(true);
  loadAlg->setLogging(false);
274
  try {
275
    loadAlg->setPropertyValue("Filename", fname);
276
277
    loadAlg->setPropertyValue("OutputWorkspace",
                              "_TabulatedFunction_fit_data_");
278
    loadAlg->execute();
279
280
281
  } catch (std::runtime_error &) {
    throw std::runtime_error(
        "Unable to load Nexus file for TabulatedFunction function.");
282
  }
283

284
  Workspace_sptr ws = loadAlg->getProperty("OutputWorkspace");
285
  MatrixWorkspace_sptr resData =
286
      std::dynamic_pointer_cast<Mantid::API::MatrixWorkspace>(ws);
287
  loadWorkspace(resData);
288
289
290
291
292
293
}

/**
 * Load the points from a MatrixWorkspace
 * @param wsName :: The workspace to load from
 */
294
295
296
void TabulatedFunction::loadWorkspace(const std::string &wsName) const {
  auto ws = AnalysisDataService::Instance().retrieveWS<MatrixWorkspace>(wsName);
  loadWorkspace(ws);
297
298
299
300
301
302
}

/**
 * Load the points from a MatrixWorkspace
 * @param ws :: The workspace to load from
 */
303
void TabulatedFunction::loadWorkspace(
304
    std::shared_ptr<API::MatrixWorkspace> ws) const {
David Fairbrother's avatar
David Fairbrother committed
305
  m_workspace = std::move(ws);
306
  m_setupFinished = false;
307
308
}

309
/**
LamarMoore's avatar
LamarMoore committed
310
311
 * Fill in the x and y value containers (m_xData and m_yData)
 */
312
313
void TabulatedFunction::setupData() const {
  if (m_setupFinished) {
314
    if (m_xData.size() != m_yData.size()) {
315
316
      throw std::invalid_argument(this->name() +
                                  ": X and Y vectors have different sizes.");
317
    }
318
319
320
    g_log.debug() << "Re-setting isn't required.";
    return;
  }
321

322
323
324
325
326
327
328
  if (!m_workspace) {
    std::string wsName = getAttribute("Workspace").asString();
    if (wsName.empty())
      throw std::invalid_argument("Data not set for function " + this->name());
    else
      loadWorkspace(wsName);
  }
329

330
  size_t index = static_cast<size_t>(getAttribute("WorkspaceIndex").asInt());
331

332
  g_log.debug() << "Setting up " << m_workspace->getName() << " index " << index
333
                << '\n';
334

335
336
  const auto &xData = m_workspace->points(index);
  const auto &yData = m_workspace->y(index);
337
338
  m_xData.assign(xData.begin(), xData.end());
  m_yData.assign(yData.begin(), yData.end());
339

340
341
  m_workspace.reset();
  m_setupFinished = true;
342
}
343

344
} // namespace Functions
345
346
} // namespace CurveFitting
} // namespace Mantid