IPeakFunction.cpp

//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidAPI/IPeakFunction.h"
#include "MantidAPI/FunctionFactory.h"
#include "MantidAPI/Jacobian.h"
#include "MantidAPI/PeakFunctionIntegrator.h"
#include "MantidAPI/FunctionParameterDecorator.h"
#include "MantidKernel/Exception.h"

#include <boost/lexical_cast.hpp>
#include <boost/make_shared.hpp>
#include <cmath>
#include <limits>

namespace Mantid {
namespace API {

namespace {

/** A Jacobian for individual functions
 */
class PartialJacobian1 : public Jacobian {
  Jacobian *m_J; ///< pointer to the overall Jacobian
  int m_iY0;     ///< offset in the overall Jacobian for a particular function
public:
  /** Constructor
   * @param J :: A pointer to the overall Jacobian
   * @param iY0 :: The data offset for a particular function
   */
  PartialJacobian1(Jacobian *J, int iY0) : m_J(J), m_iY0(iY0) {}
  /**
   * Overridden Jacobian::set(...).
   * @param iY :: The index of the data point
   * @param iP :: The parameter index of an individual function.
   * @param value :: The derivative value
   */
  void set(size_t iY, size_t iP, double value) override {
    m_J->set(m_iY0 + iY, iP, value);
  }
  /**
   * Overridden Jacobian::get(...).
   * @param iY :: The index of the data point
   * @param iP :: The parameter index of an individual function.
   */
  double get(size_t iY, size_t iP) override { return m_J->get(m_iY0 + iY, iP); }
  /** Zero all matrix elements.
  */
  void zero() override {
    throw Kernel::Exception::NotImplementedError(
        "zero() is not implemented for PartialJacobian1");
  }
};

class TempJacobian : public Jacobian {
public:
  TempJacobian(size_t y, size_t p) : m_y(y), m_p(p), m_J(y * p) {}
  void set(size_t iY, size_t iP, double value) override {
    m_J[iY * m_p + iP] = value;
  }
  double get(size_t iY, size_t iP) override { return m_J[iY * m_p + iP]; }
  size_t maxParam(size_t iY) {
    double max = -DBL_MAX;
    size_t maxIndex = 0;
    for (size_t i = 0; i < m_p; ++i) {
      double current = get(iY, i);
      if (current > max) {
        maxIndex = i;
        max = current;
      }
    }

    return maxIndex;
  }
  void zero() override { m_J.assign(m_J.size(), 0.0); }

protected:
  size_t m_y;
  size_t m_p;
  std::vector<double> m_J;
};

/// Tolerance for determining the smallest significant value on the peak
const double PEAK_TOLERANCE = 1e-14;

} // namespace

/// "Infinite" value for the peak radius
const int IPeakFunction::g_maxPeakRadius = std::numeric_limits<int>::max();

/**
  * Constructor.
  * @param peakRadius :: The peak radius for this instance.
  */
IPeakFunction::IPeakFunction(int peakRadius) : m_peakRadius(peakRadius) {}

/**
 * General implementation of the method for all peaks. Limits the peak
 * evaluation to
 * a certain number of FWHMs around the peak centre. The outside points are set
 * to 0.
 * Calls functionLocal() to compute the actual values
 * @param out :: Output function values
 * @param xValues :: X values for data points
 * @param nData :: Number of data points
 */
void IPeakFunction::function1D(double *out, const double *xValues,
                               const size_t nData) const {
  double c = this->centre();
  double dx = fabs(m_peakRadius * this->fwhm());
  int i0 = -1;
  int n = 0;
  for (size_t i = 0; i < nData; ++i) {
    if (fabs(xValues[i] - c) < dx) {
      if (i0 < 0)
        i0 = static_cast<int>(i);
      ++n;
    } else {
      out[i] = 0.0;
    }
  }
  if (i0 < 0 || n == 0)
    return;
  this->functionLocal(out + i0, xValues + i0, n);
}

/**
 * General implementation of the method for all peaks. Calculates derivatives
 * only
 * for a range of x values limited to a certain number of FWHM around the peak
 * centre.
 * For the points outside the range all derivatives are set to 0.
 * Calls functionDerivLocal() to compute the actual values
 * @param out :: Derivatives
 * @param xValues :: X values for data points
 * @param nData :: Number of data points
 */
void IPeakFunction::functionDeriv1D(Jacobian *out, const double *xValues,
                                    const size_t nData) {
  double c = this->centre();
  double dx = fabs(m_peakRadius * this->fwhm());
  int i0 = -1;
  int n = 0;
  for (size_t i = 0; i < nData; ++i) {
    if (fabs(xValues[i] - c) < dx) {
      if (i0 < 0)
        i0 = static_cast<int>(i);
      ++n;
    } else {
      for (size_t ip = 0; ip < this->nParams(); ++ip) {
        out->set(i, ip, 0.0);
      }
    }
  }
  if (i0 < 0 || n == 0)
    return;
  PartialJacobian1 J(out, i0);
  this->functionDerivLocal(&J, xValues + i0, n);
}

void IPeakFunction::setPeakRadius(const int &r) {
  if (r > 0) {
    m_peakRadius = r;
  } else if (r == 0) {
    m_peakRadius = g_maxPeakRadius;
  }
}

/// Returns the integral intensity of the peak function, using the peak radius
/// to determine integration borders.
double IPeakFunction::intensity() const {
  auto interval = getDomainInterval();

  PeakFunctionIntegrator integrator;
  IntegrationResult result = integrator.integrate(*this, interval.first, interval.second);

  if (!result.success) {
    return 0.0;
  }

  return result.result;
}

/// Sets the integral intensity of the peak by adjusting the height.
void IPeakFunction::setIntensity(const double newIntensity) {
  double currentHeight = height();
  double currentIntensity = intensity();

  if (currentIntensity == 0.0) {
    // Try to set a different height first.
    setHeight(2.0);

    currentHeight = height();
    currentIntensity = intensity();

    // If the current intensity is still 0, there's nothing left to do.
    if (currentIntensity == 0.0) {
      throw std::invalid_argument(
          "Cannot set new intensity, not enough information available.");
    }
  }

  setHeight(newIntensity / currentIntensity * currentHeight);
}

std::string IPeakFunction::getCentreParameterName() const {
  FunctionParameterDecorator_sptr fn =
      boost::dynamic_pointer_cast<FunctionParameterDecorator>(
          FunctionFactory::Instance().createFunction("PeakParameterFunction"));

  if (!fn) {
    throw std::runtime_error(
        "PeakParameterFunction could not be created successfully.");
  }

  fn->setDecoratedFunction(this->name());

  FunctionDomain1DVector domain(std::vector<double>(4, 0.0));
  TempJacobian jacobian(4, fn->nParams());

  fn->functionDeriv(domain, jacobian);

  return parameterName(jacobian.maxParam(0));
}

/// Get the interval on which the peak has all its values above a certain
/// level. All values outside the interval are below that level.
/// @param level :: A fraction of the peak height.
/// @return A pair of doubles giving the bounds of the interval.
std::pair<double, double> IPeakFunction::getDomainInterval(double level) const {
  if (level < PEAK_TOLERANCE) {
    level = PEAK_TOLERANCE;
  }
  double left = 0.0;
  double right = 0.0;
  auto h = height();
  auto w = fwhm();
  auto c = centre();
  if (h == 0.0 || w == 0.0 || level >= 1.0) {
    return std::make_pair(c, c);
  }

  auto findBound = [this, c, h, level](double dx) {
    for(size_t i = 0; i < 100; ++i) {
      double x = c + dx;
      double y = 0.0;
      this->functionLocal(&y, &x, 1);
      if (fabs(y / h) < level) {
        return x;
      }
      dx *= 2;
    }
    return c + dx;
  };

  left = findBound(-w);
  right = findBound(w);
  return std::make_pair(left, right);
}


} // namespace API
} // namespace Mantid