PseudoVoigt.cpp

Go to the documentation of this file.

// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
//   NScD Oak Ridge National Laboratory, European Spallation Source,
//   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
#include "MantidCurveFitting/Functions/PseudoVoigt.h"
#include "MantidAPI/FunctionFactory.h"
#include "MantidCurveFitting/Constraints/BoundaryConstraint.h"
 
#include <iostream>
 
#include <cmath>
 
namespace {
 
const double LN2_DIV_PI = 2. * std::sqrt(M_LN2 / M_PI);
constexpr double LN2_M4 = 4 * M_LN2;
 
inline double cal_ag(const double gamma) { return LN2_DIV_PI / gamma; }
 
inline constexpr double cal_bg(const double gamma) { return LN2_M4 / (gamma * gamma); }
 
inline double cal_gaussian(const double ag, const double bg, const double xdiffsq) {
  return ag * std::exp(-bg * xdiffsq);
}
 
inline constexpr double cal_lorentzian(const double gamma_div_2, const double gammasq_div_4, const double xdiffsq) {
  return gamma_div_2 / (xdiffsq + gammasq_div_4) / M_PI;
}
 
} // namespace
 
namespace Mantid::CurveFitting::Functions {
 
using namespace CurveFitting;
using namespace Constraints;
using namespace std;
using namespace API;
 
namespace {
Kernel::Logger g_log("PseudoVoigt");
} // namespace
 
DECLARE_FUNCTION(PseudoVoigt)
 
 
void PseudoVoigt::init() {
  declareParameter("Mixing", 0.5);
  declareParameter("Intensity", 1.);
  declareParameter("PeakCentre", 0.);
  declareParameter("FWHM", 1.);
 
  auto mixingConstraint = std::make_unique<BoundaryConstraint>(this, "Mixing", 0.0, 1.0, true);
  mixingConstraint->setPenaltyFactor(1e9);
  addConstraint(std::move(mixingConstraint));
 
  auto fwhm_constraint = std::make_unique<BoundaryConstraint>(this, "FWHM", 1.E-20, true);
  fwhm_constraint->setPenaltyFactor(1e9);
  addConstraint(std::move(fwhm_constraint));
 
  // init the peak height setup parameters
  m_height = 1.; // peak height set by user
 
  // parameter set history: all start from an arbitary out of boundary value
  // (100 as easy). if set, order set to 0
  m_set_history_distances.resize(4, 128);
  // 0: mixing, 1: intensity, 2: fwhm, 3: height
}
 
void PseudoVoigt::functionLocal(double *out, const double *xValues, const size_t nData) const {
  // read values
  const double peak_intensity = getParameter("Intensity");
  const double x0 = getParameter("PeakCentre");
  const double gamma = fabs(getParameter("FWHM"));
  const double gFraction = getParameter("Mixing");
 
  if (gamma < 1.E-20)
    throw std::runtime_error("Pseudo-voigt has FWHM as 0. It will generate "
                             "infinity at center in the Lorentzian part.");
 
  const double lFraction = 1.0 - gFraction;
 
  // calculate constants
  const double a_g = cal_ag(gamma);
  const double b_g = cal_bg(gamma);
  const double gamma_div_2 = 0.5 * gamma;
  const double gammasq_div_4 = gamma_div_2 * gamma_div_2;
 
  for (size_t i = 0; i < nData; ++i) {
    double xDiffSquared = (xValues[i] - x0) * (xValues[i] - x0);
    out[i] = peak_intensity * (gFraction * cal_gaussian(a_g, b_g, xDiffSquared) +
                               lFraction * cal_lorentzian(gamma_div_2, gammasq_div_4, xDiffSquared));
  }
}
 
void PseudoVoigt::functionDerivLocal(Jacobian *out, const double *xValues, const size_t nData) {
  // get parameters
  double peak_intensity = getParameter("Intensity");
  double x0 = getParameter("PeakCentre");
  double gamma = getParameter("FWHM");
  double gFraction = getParameter("Mixing");
  double lFraction = 1.0 - gFraction;
 
  // calcualte constants
  const double a_g = cal_ag(gamma);
  const double b_g = cal_bg(gamma);
  const double gamma_div_2 = 0.5 * gamma;
  const double gammasq_div_4 = gamma_div_2 * gamma_div_2;
 
  // derivatives
  for (size_t i = 0; i < nData; ++i) {
    double xDiff = (xValues[i] - x0);
    double xDiffSquared = xDiff * xDiff;
 
    // calculate normalized Gaussian and Lorentzian
    double gaussian_term = cal_gaussian(a_g, b_g, xDiffSquared);
    double lorentzian_term = cal_lorentzian(gamma_div_2, gammasq_div_4, xDiffSquared);
 
    // derivative to mixing/eta (0-th)
    out->set(i, 0, peak_intensity * (gaussian_term - lorentzian_term));
 
    // derivative on intensity (1-st)
    out->set(i, 1, gFraction * gaussian_term + lFraction * lorentzian_term);
 
    // peak center: x0
    const double derive_g_x0 = 2. * b_g * xDiff * gaussian_term;
    const double derive_l_x0 = 4. * M_PI * xDiff / gamma * lorentzian_term * lorentzian_term;
    const double deriv_x0 = peak_intensity * (gFraction * derive_g_x0 + lFraction * derive_l_x0);
    out->set(i, 2, deriv_x0);
 
    // peak width: gamma or H
    const double t1 = -1. / gamma * gaussian_term;
    const double t2 = 2. * b_g * xDiffSquared * gaussian_term / gamma;
    const double t3 = lorentzian_term / gamma;
    const double t4 = -M_PI * lorentzian_term * lorentzian_term;
    const double derive_gamma = peak_intensity * (gFraction * (t1 + t2) + lFraction * (t3 + t4));
    out->set(i, 3, derive_gamma);
  }
}
 
void PseudoVoigt::setParameter(const std::string &name, const double &value, bool explicitlySet) {
  g_log.debug() << "[PV] Set " << name << " as " << value << "\n";
  API::IPeakFunction::setParameter(name, value, explicitlySet);
}
 
void PseudoVoigt::setParameter(size_t i, const double &value, bool explicitlySet) {
  API::IPeakFunction::setParameter(i, value, explicitlySet);
 
  g_log.debug() << "[PV] Set " << i << "-th parameter with value " << value << "\n";
 
  // explicitly set means that there is a chance that some parameter shall be
  // re-calculated
  // peak center shall be excluded as it does nothing to do with height, H,
  // intensity and mixing
  if (explicitlySet && i != 2) {
    g_log.debug() << "Update set history for " << i << "-th parameter"
                  << "\n";
    update_set_history(i);
    estimate_parameter_value();
  }
 
  return;
}
 
bool PseudoVoigt::estimate_parameter_value() {
  // peak center
  size_t to_calculate_index = get_parameter_to_calculate_from_set();
  g_log.debug() << "Time to calculate " << to_calculate_index << "-th parameter"
                << "\n";
 
  bool some_value_updated{true};
  if (to_calculate_index == 0) {
    // calculate mixings
    double peak_intensity = getParameter(1);
    double gamma = getParameter(3);
    double mixing = (m_height * 0.5 * M_PI * gamma / peak_intensity - 1) / (sqrt(M_PI * M_LN2) - 1);
    if (mixing < 0) {
      g_log.debug() << "Calculated mixing (eta) = " << mixing << " < 0.  Set to 0 instead"
                    << "\n";
      mixing = 0;
    } else if (mixing > 1) {
      g_log.debug() << "Calculated mixing (eta) = " << mixing << " > 1. Set to 1 instead"
                    << "\n";
      mixing = 1;
    }
 
    setParameter(0, mixing,
                 false); // no explicitly set: won't cause further re-calcualting
  } else if (to_calculate_index == 1) {
    // calculate intensity
    double eta = getParameter(0);
    double gamma = getParameter(3);
    double intensity = m_height / 2. / (1 + (sqrt(M_PI * M_LN2) - 1) * eta) * (M_PI * gamma);
    setParameter(1, intensity, false);
    g_log.debug() << "Estimate peak intensity = " << intensity << "\n";
  } else if (to_calculate_index == 3) {
    // calculate peak width
    double eta = getParameter(0);
    double intensity = getParameter(1);
 
    g_log.debug() << "Intensity = " << intensity << ", height = " << m_height << ", mixing = " << eta << "\n";
    double gamma = intensity * 2 * (1 + (sqrt(M_PI * M_LN2) - 1) * eta) / (M_PI * m_height);
    if (gamma < 1.E-10) {
      g_log.debug() << "Peak width (H or Gamma) = " << gamma << ". Set 1.E-2 instead"
                    << "\n";
    }
    setParameter(3, gamma, false);
  } else {
    // no value updated
    some_value_updated = false;
  }
 
  return some_value_updated;
}
 
void PseudoVoigt::update_set_history(size_t set_index) {
  if (set_index > 3)
    throw std::runtime_error("Parameter set up history index out of boundary");
 
  size_t prev_distance = m_set_history_distances[set_index];
  for (size_t i = 0; i < 4; ++i) {
    // no operation on any never-set or out-dated index
    if (m_set_history_distances[i] >= 3)
      continue;
 
    // only increase the historoy distance by 1 in case
    // the other parameters with those with distance less than the previous one
    // i.e., out of date or never been set
    if (m_set_history_distances[i] < prev_distance) {
      m_set_history_distances[i] += 1;
      // if out of range, don't worry about it
    }
  }
 
  // set the current one to 0 (most recent)
  m_set_history_distances[set_index] = 0;
}
 
size_t PseudoVoigt::get_parameter_to_calculate_from_set() {
  // go through setting-history-distance
  size_t index_largest_distance{0};
  size_t largest_distance{0};
  size_t num_been_set{0};
 
  for (size_t i = 0; i < m_set_history_distances.size(); ++i) {
    g_log.debug() << "distance[" << i << "] = " << m_set_history_distances[i] << "\n";
    if (m_set_history_distances[i] >= largest_distance) {
      // for those to set for: any number counts
      index_largest_distance = i;
      largest_distance = m_set_history_distances[i];
    }
    // only those been set and not out of date counted
    if (m_set_history_distances[i] < 3) {
      ++num_been_set;
    }
  }
 
  // signal for no enough index been set up
  if (num_been_set < 3)
    index_largest_distance = 128;
 
  return index_largest_distance;
}
 
double PseudoVoigt::height() const {
  double peak_intensity = getParameter("Intensity");
  double gamma = getParameter("FWHM");
  double eta = getParameter("Mixing");
  const double height = 2. * peak_intensity * (1 + (sqrt(M_PI * M_LN2) - 1) * eta) / (M_PI * gamma);
 
  return height;
}
 
void PseudoVoigt::setHeight(const double h) {
  // set height
  double prev_height = m_height;
  m_height = h;
 
  // update set history
  update_set_history(2);
  g_log.debug() << "PV: set height = " << m_height << "\n";
 
  // update other parameters' value
  bool updated_any_value = estimate_parameter_value();
  if (!updated_any_value) {
    // force to update peak intensity
    double old_intensity = getParameter("Intensity");
    // This is what Roman suggests!
    double new_intensity{1.};
    if (prev_height > 0)
      new_intensity = old_intensity * m_height / prev_height;
    else
      new_intensity = m_height; // this is a better than nothing estimation
    setParameter("Intensity", new_intensity,
                 false); // no need to update other parameters
  }
}
 
void PseudoVoigt::setFwhm(const double w) { setParameter("FWHM", w, true); }
 
} // namespace Mantid::CurveFitting::Functions