Unverified Commit 0ea8f261 authored by Pete Peterson's avatar Pete Peterson Committed by GitHub
Browse files

Merge pull request #22430 from martyngigg/22331_python_fitfunction_crash

Robustness improvements to Python fit functions
parents 9ff01c48 316b031c
......@@ -80,6 +80,13 @@ public:
const size_t nData);
protected:
template <typename FunctionType>
using Function1DMethod = void (FunctionType::*)(double *, const double *,
const size_t) const;
/// Calculate a numerical derivative for the 1D data
template <typename EvaluationMethod>
void calcNumericalDerivative1D(Jacobian *jacobian, EvaluationMethod func1D,
const double *xValues, const size_t nData);
/// Calculate histogram data for the given bin boundaries.
virtual void histogram1D(double *out, double left, const double *right,
const size_t nBins) const;
......
#ifndef MANTID_API_IFUNCTION1D_TCC
#define MANTID_API_IFUNCTION1D_TCC
/*
Copyright &copy; 2016 ISIS Rutherford Appleton Laboratory, NScD Oak Ridge
National Laboratory & European Spallation Source
This file is part of Mantid.
Mantid is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 3 of the License, or
(at your option) any later version.
Mantid is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
File change history is stored at: <https://github.com/mantidproject/mantid>
Code Documentation is available at: <http://doxygen.mantidproject.org>
*/
#include "MantidAPI/IFunction1D.h"
#include <cmath>
namespace Mantid {
namespace API {
template <typename EvaluationMethod>
void IFunction1D::calcNumericalDerivative1D(Jacobian *jacobian,
EvaluationMethod eval1D,
const double *xValues,
const size_t nData) {
/*
* Uses the knowledge that this is a simple 1D domain to calculate the
* derivative. The domains are not constructed and instead function1D
* is called directly.
*
* There is a similar more generic method for all functions in IFunction
* but the method takes different parameters and uses slightly different
* function calls in places making it difficult to share code. Please also
* consider that method when updating this.
*/
using std::fabs;
constexpr double epsilon(std::numeric_limits<double>::epsilon() * 100);
constexpr double stepPercentage(0.001);
constexpr double cutoff =
100.0 * std::numeric_limits<double>::min() / stepPercentage;
applyTies(); // just in case
std::vector<double> minusStep(nData), plusStep(nData);
eval1D(minusStep.data(), xValues, nData);
double step;
const size_t nParam = nParams();
for (size_t iP = 0; iP < nParam; iP++) {
if (isActive(iP)) {
const double val = activeParameter(iP);
if (fabs(val) < cutoff) {
step = epsilon;
} else {
step = val * stepPercentage;
}
const double paramPstep = val + step;
setActiveParameter(iP, paramPstep);
applyTies();
eval1D(plusStep.data(), xValues, nData);
setActiveParameter(iP, val);
applyTies();
step = paramPstep - val;
for (size_t i = 0; i < nData; i++) {
jacobian->set(i, iP, (plusStep[i] - minusStep[i]) / step);
}
}
}
}
}
}
#endif // MANTID_API_IFUNCTION1D_TCC
......@@ -70,9 +70,9 @@ public:
/// Function evaluation method to be implemented in the inherited classes
virtual void functionLocal(double *out, const double *xValues,
const size_t nData) const = 0;
/// Derivative evaluation method to be implemented in the inherited classes
virtual void functionDerivLocal(Jacobian *out, const double *xValues,
const size_t nData) = 0;
/// Derivative evaluation method. Default is to calculate numerically
virtual void functionDerivLocal(Jacobian *jacobian, const double *xValues,
const size_t nData);
/// Get name of parameter that is associated to centre.
std::string getCentreParameterName() const;
......
......@@ -906,27 +906,31 @@ std::string IFunction::descriptionOfActive(size_t i) const {
*/
void IFunction::calNumericalDeriv(const FunctionDomain &domain,
Jacobian &jacobian) {
const double minDouble = std::numeric_limits<double>::min();
const double epsilon = std::numeric_limits<double>::epsilon() * 100;
double stepPercentage = 0.001; // step percentage
double step; // real step
double cutoff = 100.0 * minDouble / stepPercentage;
size_t nParam = nParams();
/*
* There is a similar more specialized method for 1D functions in IFunction1D
* but the method takes different parameters and uses slightly different
* function calls in places making it difficult to share code. Please also
* consider that method when updating this.
*/
constexpr double epsilon = std::numeric_limits<double>::epsilon() * 100;
constexpr double stepPercentage = 0.001;
constexpr double cutoff =
100.0 * std::numeric_limits<double>::min() / stepPercentage;
const size_t nParam = nParams();
size_t nData = getValuesSize(domain);
FunctionValues minusStep(nData);
FunctionValues plusStep(nData);
// PARALLEL_CRITICAL(numeric_deriv)
{
applyTies(); // just in case
function(domain, minusStep);
}
applyTies(); // just in case
function(domain, minusStep);
if (nData == 0) {
nData = minusStep.size();
}
double step;
for (size_t iP = 0; iP < nParam; iP++) {
if (isActive(iP)) {
const double val = activeParameter(iP);
......@@ -936,16 +940,12 @@ void IFunction::calNumericalDeriv(const FunctionDomain &domain,
step = val * stepPercentage;
}
double paramPstep = val + step;
// PARALLEL_CRITICAL(numeric_deriv)
{
setActiveParameter(iP, paramPstep);
applyTies();
function(domain, plusStep);
setActiveParameter(iP, val);
applyTies();
}
const double paramPstep = val + step;
setActiveParameter(iP, paramPstep);
applyTies();
function(domain, plusStep);
setActiveParameter(iP, val);
applyTies();
step = paramPstep - val;
for (size_t i = 0; i < nData; i++) {
......
......@@ -2,6 +2,7 @@
// Includes
//----------------------------------------------------------------------
#include "MantidAPI/IFunction1D.h"
#include "MantidAPI/IFunction1D.tcc"
#include "MantidAPI/Jacobian.h"
#include "MantidAPI/IFunctionWithLocation.h"
#include "MantidAPI/IConstraint.h"
......@@ -85,8 +86,12 @@ void IFunction1D::derivative1D(double *out, const double *xValues, size_t nData,
void IFunction1D::functionDeriv1D(Jacobian *jacobian, const double *xValues,
const size_t nData) {
FunctionDomain1DView domain(xValues, nData);
this->calNumericalDeriv(domain, *jacobian);
auto evalMethod =
[this](double *out, const double *xValues, const size_t nData) {
this->function1D(out, xValues, nData);
};
this->calcNumericalDerivative1D(jacobian, std::move(evalMethod), xValues,
nData);
}
/// Calculate histogram data for the given bin boundaries.
......
......@@ -2,6 +2,7 @@
// Includes
//----------------------------------------------------------------------
#include "MantidAPI/IPeakFunction.h"
#include "MantidAPI/IFunction1D.tcc"
#include "MantidAPI/FunctionFactory.h"
#include "MantidAPI/Jacobian.h"
#include "MantidAPI/PeakFunctionIntegrator.h"
......@@ -265,5 +266,22 @@ std::pair<double, double> IPeakFunction::getDomainInterval(double level) const {
return std::make_pair(left, right);
}
/**
* Computes the function derivative numerically
* @param jacobian An output Jacobian to receive the calculated values
* @param xValues An input array of X data
* @param nData The number of X values provided
*/
void IPeakFunction::functionDerivLocal(Jacobian *jacobian,
const double *xValues,
const size_t nData) {
auto evalMethod =
[this](double *out, const double *xValues, const size_t nData) {
this->functionLocal(out, xValues, nData);
};
this->calcNumericalDerivative1D(jacobian, std::move(evalMethod), xValues,
nData);
}
} // namespace API
} // namespace Mantid
......@@ -54,6 +54,9 @@ class IFunction1DAdapter : public virtual API::ParamFunction,
public IFunctionAdapter {
#endif
public:
// Convenience type def
using Base = API::IFunction1D;
/// A constructor that looks like a Python __init__ method
IFunction1DAdapter(PyObject *self);
......@@ -76,10 +79,6 @@ public:
void functionDeriv1D(API::Jacobian *out, const double *xValues,
const size_t nData) override;
///@}
private:
/// Flag if the functionDeriv1D method is overridden (avoids multiple checks)
bool m_derivOveridden;
};
}
}
......
......@@ -38,7 +38,8 @@ namespace PythonInterface {
class IFunctionAdapter : virtual public API::IFunction {
public:
/// A constructor that looks like a Python __init__ method
IFunctionAdapter(PyObject *self);
IFunctionAdapter(PyObject *self, std::string functionMethod,
std::string derivMethod);
/// The PyObject must be supplied to construct the object
IFunctionAdapter(const IFunctionAdapter &) = delete;
......@@ -109,17 +110,26 @@ public:
void setActiveParameter(size_t i, double value) override;
protected:
/**
* Returns the PyObject that owns this wrapper, i.e. self
* @returns A pointer to self
*/
/// @returns The PyObject that owns this wrapper, i.e. self
inline PyObject *getSelf() const { return m_self; }
/// @returns True if the instance overrides the derivative method
inline bool derivativeOverridden() const { return m_derivOveridden; }
/// Evaluate the function by calling the overridden method
void evaluateFunction(double *out, const double *xValues,
const size_t nData) const;
/// Evaluate the derivative by calling the overridden method
void evaluateDerivative(API::Jacobian *out, const double *xValues,
const size_t nData) const;
private:
/// The name of the function
std::string m_name;
/// The Python portion of the object
PyObject *m_self;
/// The name of the method to evaluate the function
std::string m_functionName;
/// The name of the method to evaluate the derivative
std::string m_derivName;
/// Flag if the derivateive method is overridden (avoids multiple checks)
bool m_derivOveridden;
};
}
}
......
......@@ -26,7 +26,7 @@
// Includes
//-----------------------------------------------------------------------------
#include "MantidAPI/IPeakFunction.h"
#include "MantidPythonInterface/api/FitFunctions/IFunction1DAdapter.h"
#include "MantidPythonInterface/api/FitFunctions/IFunctionAdapter.h"
#include <boost/python/object.hpp>
......@@ -41,8 +41,11 @@ namespace PythonInterface {
*into Python.
*/
class IPeakFunctionAdapter : public API::IPeakFunction,
public IFunction1DAdapter {
public IFunctionAdapter {
public:
// Convenience typedef
using Base = API::IPeakFunction;
/// A constructor that looks like a Python __init__ method
IPeakFunctionAdapter(PyObject *self);
......@@ -86,7 +89,7 @@ public:
/// Python-type signature for above method
boost::python::object functionLocal(const boost::python::object &xvals) const;
/// Implemented base-class method
void functionDerivLocal(API::Jacobian *out, const double *xValues,
void functionDerivLocal(API::Jacobian *jacobian, const double *xValues,
const size_t nData) override;
/// Python signature
void functionDerivLocal(const boost::python::object &xvals,
......
......@@ -26,26 +26,23 @@
#pragma GCC system_header
#endif
#include <boost/python/list.hpp>
#include "MantidKernel/WarningSuppressions.h"
#include <boost/python/detail/wrap_python.hpp>
// clang-format off
GCC_DIAG_OFF(cast-qual)
// clang-format on
// Forward declare the numpy array types
struct tagPyArrayObject;
using PyArrayObject = tagPyArrayObject;
struct _PyArray_Descr;
using PyArray_Descr = _PyArray_Descr;
// See
// http://docs.scipy.org/doc/numpy/reference/c-api.array.html#PY_ARRAY_UNIQUE_SYMBOL
#define PY_ARRAY_UNIQUE_SYMBOL KERNEL_ARRAY_API
#define NO_IMPORT_ARRAY
#include <numpy/arrayobject.h>
// clang-format off
GCC_DIAG_ON(cast-qual)
// clang-format on
/**functions containing numpy macros. We put them in a separate header file to
*suppress the warning
*ISO C++ forbids casting between pointer-to-function and pointer-to-object
*/
/*
* The numpy API is a C api where pointers to functions and objects are the
* same size.
* This is not guaranteed by the C++ standard so macros in the numpy headers,
* such as PyArray_IterNew produces warnings when compiled with a C++ compiler.
*
* The only solution appears to be wrap the calls in our own functions and
* suppress the warnings.
*/
namespace Mantid {
namespace PythonInterface {
namespace Converters {
......@@ -55,9 +52,9 @@ PyObject *func_PyArray_IterNew(PyArrayObject *arr);
/// equivalent to macro PyArray_NewFromDescr
PyArrayObject *func_PyArray_NewFromDescr(int datatype, const int ndims,
Py_intptr_t *dims);
PyArrayObject *func_PyArray_NewFromDescr(const std::string &datadescr,
const int ndims, Py_intptr_t *dims);
PyArray_Descr *func_PyArray_Descr(const std::string &datadescr);
PyArrayObject *func_PyArray_NewFromDescr(const char *datadescr, const int ndims,
Py_intptr_t *dims);
PyArray_Descr *func_PyArray_Descr(const char *datadescr);
}
}
}
......
......@@ -37,6 +37,52 @@ struct UndefinedAttributeError {
virtual ~UndefinedAttributeError() = default;
};
namespace detail {
/**
* Wrapper around boost::python::call_method. If the call raises a Python error
* then this is translated to a C++ exception object inheriting from
* std::exception or std::runtime_error depending on the type of Python error.
*
* Note that this is an implementation method that does not hold the GIL and
* is only intended to be used below.
* @param obj Pointer to Python object
* @param methodName Name of the method call
* @param args A list of arguments to forward to call_method
*/
template <typename ReturnType, typename... Args>
ReturnType callMethodImpl(PyObject *obj, const char *methodName,
const Args &... args) {
try {
return boost::python::call_method<ReturnType, Args...>(obj, methodName,
args...);
} catch (boost::python::error_already_set &) {
PyObject *exception = PyErr_Occurred();
assert(exception);
if (PyErr_GivenExceptionMatches(exception, PyExc_RuntimeError)) {
throw PythonRuntimeError();
} else {
throw PythonException();
}
}
}
}
/**
* Wrapper around boost::python::call_method to acquire GIL for duration
* of call. If the call raises a Python error then this is translated to
* a C++ exception object inheriting from std::exception or std::runtime_error
* depending on the type of Python error.
* @param obj Pointer to Python object
* @param methodName Name of the method call
* @param args A list of arguments to forward to call_method
*/
template <typename ReturnType, typename... Args>
ReturnType callMethodNoCheck(PyObject *obj, const char *methodName,
const Args &... args) {
Environment::GlobalInterpreterLock gil;
return detail::callMethodImpl<ReturnType, Args...>(obj, methodName, args...);
}
/**
* Wrapper around boost::python::call_method to acquire GIL for duration
* of call. If the attribute does not exist then an UndefinedAttributeError
......@@ -52,18 +98,8 @@ ReturnType callMethod(PyObject *obj, const char *methodName,
const Args &... args) {
GlobalInterpreterLock gil;
if (Environment::typeHasAttribute(obj, methodName)) {
try {
return boost::python::call_method<ReturnType, Args...>(obj, methodName,
args...);
} catch (boost::python::error_already_set &) {
PyObject *exception = PyErr_Occurred();
assert(exception);
if (PyErr_GivenExceptionMatches(exception, PyExc_RuntimeError)) {
throw PythonRuntimeError();
} else {
throw PythonException();
}
}
return detail::callMethodImpl<ReturnType, Args...>(obj, methodName,
args...);
} else {
throw UndefinedAttributeError();
}
......
#ifndef MANTID_PYTHONINTERFACE_KERNEL_H_
#define MANTID_PYTHONINTERFACE_KERNEL_H_
/*
* Provides a function to initialize the numpy c api
* function pointer table in the kernel module. This
* is *only* required for the C++ unit tests on Windows.
*
* Normally, importing mantid.kernel into Python causes
* the internal numpy array api to be initialized using
* the _import_array call in the module startup code. This
* assumes that the extension module itself is only loaded
* by the Python dynamic loader and not the operating system
* library loader as a dependency on another executable. In
* our Python C++ unit tests the _kernel.pyd library is
* linked to the unittest executable and is therefore loaded
* by the OS. On Windows each DLL has a private symbol table
* and importing mantid.kernel as part of the unit test
* only intializes the numpy c api for that copy of the
* shared library. The unittest executable also sees a secondary
* copy from the dynamic linking. The C API pointer for this
* copy of the library also needs initializing and this function
* provides this capability.
*/
#ifdef _WIN32
#include "MantidKernel/System.h"
DLLExport void kernel_dll_import_numpy_capi_for_unittest();
#endif
#endif // MANTID_PYTHONINTERFACE_KERNEL_H_H
\ No newline at end of file
......@@ -69,42 +69,44 @@ createCompositeFunction(FunctionFactoryImpl &self, const std::string &name) {
throw std::invalid_argument(error_message);
}
//--------------------------------------------- Function registration
//------------------------------------------------
//----- Function registration -----
/// Python algorithm registration mutex in anonymous namespace (aka static)
std::recursive_mutex FUNCTION_REGISTER_MUTEX;
/**
* A free function to register a fit function from Python
* @param obj :: A Python object that should either be a class type derived from
* IFunction
* or an instance of a class type derived from IFunction
* @param classObject A Python class derived from IFunction
*/
void subscribe(FunctionFactoryImpl &self, const boost::python::object &obj) {
void subscribe(FunctionFactoryImpl &self, PyObject *classObject) {
std::lock_guard<std::recursive_mutex> lock(FUNCTION_REGISTER_MUTEX);
static PyTypeObject *baseClass = const_cast<PyTypeObject *>(
converter::registered<IFunction>::converters.to_python_target_type());
// obj could be or instance/class, check instance first
PyObject *classObject(nullptr);
if (PyObject_IsInstance(obj.ptr(), reinterpret_cast<PyObject *>(baseClass))) {
classObject = PyObject_GetAttrString(obj.ptr(), "__class__");
} else if (PyObject_IsSubclass(obj.ptr(),
reinterpret_cast<PyObject *>(baseClass))) {
classObject = obj.ptr(); // We need to ensure the type of lifetime
// management so grab the raw pointer
} else {
throw std::invalid_argument(
"Cannot register a function that does not derive from IFunction.");
// object mantidapi(handle<>(PyImport_ImportModule("mantid.api")));
// object ifunction = mantidapi.attr("IFunction");
// obj should be a class deriving from IFunction
// PyObject_IsSubclass can set the error handler if classObject
// is not a class so this needs to be checked in two stages
const bool isSubClass =
(PyObject_IsSubclass(classObject,
reinterpret_cast<PyObject *>(baseClass)) == 1);
if (PyErr_Occurred() || !isSubClass) {
throw std::invalid_argument(std::string("subscribe(): Unexpected type. "
"Expected a class derived from "
"IFunction1D or IPeakFunction, "
"found: ") +
classObject->ob_type->tp_name);
}
// Instantiator will store a reference to the class object, so increase
// reference count with borrowed template
auto classHandle = handle<>(borrowed(classObject));
auto *creator = new PythonObjectInstantiator<IFunction>(object(classHandle));
auto *creator = new PythonObjectInstantiator<IFunction>(
object(handle<>(borrowed(classObject))));
// Find the function name
// Can the function be created and initialized? It really shouldn't go in
// to the factory if not
auto func = creator->createInstance();
func->initialize();
// Takes ownership of instantiator
self.subscribe(func->name(), creator, FunctionFactoryImpl::OverwriteCurrent);
......
#include "MantidPythonInterface/api/FitFunctions/IFunction1DAdapter.h"
#include "MantidPythonInterface/kernel/Converters/WrapWithNumpy.h"
#include "MantidPythonInterface/kernel/Environment/CallMethod.h"
#include "MantidPythonInterface/kernel/Environment/WrapperHelpers.h"
#include <boost/python/class.hpp>
#define PY_ARRAY_UNIQUE_SYMBOL KERNEL_ARRAY_API
#define PY_ARRAY_UNIQUE_SYMBOL API_ARRAY_API
#define NO_IMPORT_ARRAY
#include <numpy/arrayobject.h>
......@@ -15,6 +11,7 @@
namespace Mantid {
namespace PythonInterface {
using Environment::callMethod;
using Environment::callMethodNoCheck;
using namespace boost::python;
/**
......@@ -22,10 +19,8 @@ using namespace boost::python;
* @param self A reference to the calling Python object
*/
IFunction1DAdapter::IFunction1DAdapter(PyObject *self)
: API::ParamFunction(), API::IFunction1D(), IFunctionAdapter(self),
m_derivOveridden(false) {
m_derivOveridden = Environment::typeHasAttribute(self, "functionDeriv1D");
}
: API::ParamFunction(), API::IFunction1D(),
IFunctionAdapter(self, "function1D", "functionDeriv1D") {}
/**