FitPeaks.cpp

//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidAlgorithms/FitPeaks.h"
#include "MantidAPI/Axis.h"
#include "MantidAPI/CompositeFunction.h"
#include "MantidAPI/CostFunctionFactory.h"
#include "MantidAPI/FuncMinimizerFactory.h"
#include "MantidAPI/FunctionFactory.h"
#include "MantidAPI/FunctionProperty.h"
#include "MantidAPI/MultiDomainFunction.h"
#include "MantidAPI/TableRow.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidAPI/WorkspaceProperty.h"
#include "MantidAlgorithms/FindPeakBackground.h"
#include "MantidDataObjects/TableWorkspace.h"
#include "MantidDataObjects/Workspace2D.h"
#include "MantidHistogramData/EstimatePolynomial.h"
#include "MantidHistogramData/HistogramIterator.h"
#include "MantidKernel/ArrayProperty.h"
#include "MantidKernel/BoundedValidator.h"
#include "MantidKernel/IValidator.h"
#include "MantidKernel/ListValidator.h"
#include "MantidKernel/StartsWithValidator.h"

#include "boost/algorithm/string.hpp"
#include "boost/algorithm/string/trim.hpp"
#include <limits>

using namespace Mantid;
using namespace Mantid::API;
using namespace Mantid::DataObjects;
using namespace Mantid::Kernel;
using Mantid::HistogramData::Histogram;
using namespace std;

const size_t MIN_EVENTS = 100;

namespace Mantid {
namespace Algorithms {

namespace FitPeaksAlgorithm {

//----------------------------------------------------------------------------------------------
/// Holds all of the fitting information for a single spectrum
PeakFitResult::PeakFitResult(size_t num_peaks, size_t num_params)
    : m_function_parameters_number(num_params) {
  // check input
  if (num_peaks == 0 || num_params == 0)
    throw std::runtime_error("No peak or no parameter error.");

  //
  m_fitted_peak_positions.resize(num_peaks,
                                 std::numeric_limits<double>::quiet_NaN());
  m_costs.resize(num_peaks, DBL_MAX);
  m_function_parameters_vector.resize(num_peaks);
  for (size_t ipeak = 0; ipeak < num_peaks; ++ipeak) {
    m_function_parameters_vector[ipeak].resize(
        num_params, std::numeric_limits<double>::quiet_NaN());
  }

  return;
}

//----------------------------------------------------------------------------------------------
size_t PeakFitResult::getNumberParameters() const {
  return m_function_parameters_number;
}

size_t PeakFitResult::getNumberPeaks() const {
  return m_function_parameters_vector.size();
}

double PeakFitResult::getParameterValue(size_t ipeak, size_t iparam) const {
  return m_function_parameters_vector[ipeak][iparam];
}

//----------------------------------------------------------------------------------------------
double PeakFitResult::getPeakPosition(size_t ipeak) const {
  return m_fitted_peak_positions[ipeak];
}

//----------------------------------------------------------------------------------------------
double PeakFitResult::getCost(size_t ipeak) const { return m_costs[ipeak]; }

//----------------------------------------------------------------------------------------------
/// set the peak fitting record/parameter for one peak
void PeakFitResult::setRecord(size_t ipeak, const double cost,
                              const double peak_position,
                              const FitFunction fit_functions) {
  // check input
  if (ipeak >= m_costs.size())
    throw std::runtime_error("Peak index is out of range.");

  // set the values
  m_costs[ipeak] = cost;

  // set peak position
  m_fitted_peak_positions[ipeak] = peak_position;

  // transfer from peak function to vector
  size_t peak_num_params = fit_functions.peakfunction->nParams();
  for (size_t ipar = 0; ipar < peak_num_params; ++ipar) {
    // peak function
    m_function_parameters_vector[ipeak][ipar] =
        fit_functions.peakfunction->getParameter(ipar);
  }
  for (size_t ipar = 0; ipar < fit_functions.bkgdfunction->nParams(); ++ipar) {
    // background function
    m_function_parameters_vector[ipeak][ipar + peak_num_params] =
        fit_functions.bkgdfunction->getParameter(ipar);
  }
}

/// The peak postition should be negative and indicates what went wrong
void PeakFitResult::setBadRecord(size_t ipeak, const double peak_position) {
  // check input
  if (ipeak >= m_costs.size())
    throw std::runtime_error("Peak index is out of range");
  if (peak_position >= 0.)
    throw std::runtime_error("Can only set negative postion for bad record");

  // set the values
  m_costs[ipeak] = DBL_MAX;

  // set peak position
  m_fitted_peak_positions[ipeak] = peak_position;

  // transfer from peak function to vector
  for (size_t ipar = 0; ipar < m_function_parameters_number; ++ipar) {
    m_function_parameters_vector[ipeak][ipar] = 0.;
  }
}
}

//----------------------------------------------------------------------------------------------
/** Get an index of a value in a sorted vector.  The index should be the item
 * with value nearest to X
  */
size_t findXIndex(const std::vector<double> &vecx, double x) {
  size_t index;
  if (x <= vecx.front()) {
    index = 0;
  } else if (x >= vecx.back()) {
    index = vecx.size() - 1;
  } else {
    vector<double>::const_iterator fiter =
        lower_bound(vecx.begin(), vecx.end(), x);
    if (fiter == vecx.end())
      throw runtime_error("It seems impossible to have this value. ");

    index = static_cast<size_t>(fiter - vecx.begin());
    if (x - vecx[index - 1] < vecx[index] - x)
      --index;
  }

  return index;
}

enum PeakFitResult { NOSIGNAL, LOWPEAK, OUTOFBOUND, GOOD };

//----------------------------------------------------------------------------------------------
FitPeaks::FitPeaks()
    : m_fitPeaksFromRight(true), m_numPeaksToFit(0), m_minPeakHeight(20.),
      m_bkgdSimga(1.), m_peakPosTolCase234(false) {}

//----------------------------------------------------------------------------------------------
/** initialize the properties
 */
void FitPeaks::init() {
  declareProperty(Kernel::make_unique<WorkspaceProperty<MatrixWorkspace>>(
                      "InputWorkspace", "", Direction::Input),
                  "Name of the input workspace for peak fitting.");
  declareProperty(Kernel::make_unique<WorkspaceProperty<MatrixWorkspace>>(
                      "OutputWorkspace", "", Direction::Output),
                  "Name of the output workspace containing peak centers for "
                  "fitting offset."
                  "The output workspace is point data."
                  "Each workspace index corresponds to a spectrum. "
                  "Each X value ranges from 0 to N-1, where N is the number of "
                  "peaks to fit. "
                  "Each Y value is the peak position obtained by peak fitting. "
                  "Negative value is used for error signals. "
                  "-1 for data is zero;  -2 for maximum value is smaller than "
                  "specified minimum value."
                  "and -3 for non-converged fitting.");

  // properties about fitting range and criteria
  declareProperty("StartWorkspaceIndex", EMPTY_INT(),
                  "Starting workspace index for fit");
  declareProperty("StopWorkspaceIndex", EMPTY_INT(),
                  "Last workspace index to fit (which is included)");

  // properties about peak positions to fit
  declareProperty(Kernel::make_unique<ArrayProperty<double>>("PeakCenters"),
                  "List of peak centers to fit against.");
  declareProperty(
      Kernel::make_unique<WorkspaceProperty<MatrixWorkspace>>(
          "PeakCentersWorkspace", "", Direction::Input, PropertyMode::Optional),
      "MatrixWorkspace containing peak centers");

  std::string peakcentergrp("Peak Positions");
  setPropertyGroup("PeakCenters", peakcentergrp);
  setPropertyGroup("PeakCentersWorkspace", peakcentergrp);

  // properties about peak profile
  std::vector<std::string> peakNames =
      FunctionFactory::Instance().getFunctionNames<API::IPeakFunction>();
  declareProperty("PeakFunction", "Gaussian",
                  boost::make_shared<StringListValidator>(peakNames));
  vector<string> bkgdtypes{"Flat", "Linear", "Quadratic"};
  declareProperty("BackgroundType", "Linear",
                  boost::make_shared<StringListValidator>(bkgdtypes),
                  "Type of Background.");

  std::string funcgroup("Function Types");
  setPropertyGroup("PeakFunction", funcgroup);
  setPropertyGroup("BackgroundType", funcgroup);

  // properties about peak range including fitting window and peak width
  // (percentage)
  declareProperty(
      Kernel::make_unique<ArrayProperty<double>>("FitWindowBoundaryList"),
      "List of left boundaries of the peak fitting window corresponding to "
      "PeakCenters.");

  declareProperty(Kernel::make_unique<WorkspaceProperty<MatrixWorkspace>>(
                      "FitPeakWindowWorkspace", "", Direction::Input,
                      PropertyMode::Optional),
                  "MatrixWorkspace for of peak windows");

  auto min = boost::make_shared<BoundedValidator<double>>();
  min->setLower(1e-3);
  // min->setUpper(1.); TODO make this a limit
  declareProperty("PeakWidthPercent", EMPTY_DBL(), min,
                  "The estimated peak width as a "
                  "percentage of the d-spacing "
                  "of the center of the peak. Value must be less than 1.");

  std::string fitrangeegrp("Peak Range Setup");
  setPropertyGroup("PeakWidthPercent", fitrangeegrp);
  setPropertyGroup("FitWindowBoundaryList", fitrangeegrp);
  setPropertyGroup("FitPeakWindowWorkspace", fitrangeegrp);

  // properties about peak parameters' names and value
  declareProperty(
      Kernel::make_unique<ArrayProperty<std::string>>("PeakParameterNames"),
      "List of peak parameters' names");
  declareProperty(
      Kernel::make_unique<ArrayProperty<double>>("PeakParameterValues"),
      "List of peak parameters' value");
  declareProperty(Kernel::make_unique<WorkspaceProperty<TableWorkspace>>(
                      "PeakParameterValueTable", "", Direction::Input,
                      PropertyMode::Optional),
                  "Name of the an optional workspace, whose each column "
                  "corresponds to given peak parameter names"
                  ", and each row corresponds to a subset of spectra.");

  std::string startvaluegrp("Starting Parameters Setup");
  setPropertyGroup("PeakParameterNames", startvaluegrp);
  setPropertyGroup("PeakParameterValues", startvaluegrp);
  setPropertyGroup("PeakParameterValueTable", startvaluegrp);

  // optimization setup
  declareProperty("FitFromRight", true,
                  "Flag for the order to fit peaks.  If true, peaks are fitted "
                  "from rightmost;"
                  "Otherwise peaks are fitted from leftmost.");

  std::vector<std::string> minimizerOptions =
      API::FuncMinimizerFactory::Instance().getKeys();
  declareProperty("Minimizer", "Levenberg-Marquardt",
                  Kernel::IValidator_sptr(
                      new Kernel::StartsWithValidator(minimizerOptions)),
                  "Minimizer to use for fitting. Minimizers available are "
                  "\"Levenberg-Marquardt\", \"Simplex\","
                  "\"Conjugate gradient (Fletcher-Reeves imp.)\", \"Conjugate "
                  "gradient (Polak-Ribiere imp.)\", \"BFGS\", and "
                  "\"Levenberg-MarquardtMD\"");

  std::array<string, 2> costFuncOptions = {{"Least squares", "Rwp"}};
  declareProperty("CostFunction", "Least squares",
                  Kernel::IValidator_sptr(
                      new Kernel::ListValidator<std::string>(costFuncOptions)),
                  "Cost functions");

  std::string optimizergrp("Optimization Setup");
  setPropertyGroup("Minimizer", optimizergrp);
  setPropertyGroup("CostFunction", optimizergrp);

  // other helping information
  declareProperty(
      "FindBackgroundSigma", 1.0,
      "Multiplier of standard deviations of the variance for convergence of "
      "peak elimination.  Default is 1.0. ");

  declareProperty("HighBackground", true,
                  "Flag whether the data has high background comparing to "
                  "peaks' intensities. "
                  "For example, vanadium peaks usually have high background.");

  declareProperty(
      Kernel::make_unique<ArrayProperty<double>>("PositionTolerance"),
      "List of tolerance on fitted peak positions against given peak positions."
      "If there is only one value given, then ");

  declareProperty("MinimumPeakHeight", 0.,
                  "Minimum peak height such that all the fitted peaks with "
                  "height under this value will be excluded.");

  declareProperty(
      "ConstrainPeakPositions", true,
      "If true peak position will be constrained by estimated positions "
      "(highest Y value position) and "
      "the peak width either estimted by observation or calculate.");

  std::string helpgrp("Additional Information");

  // additional output for reviewing
  declareProperty(Kernel::make_unique<WorkspaceProperty<API::ITableWorkspace>>(
                      "OutputPeakParametersWorkspace", "", Direction::Output),
                  "Name of workspace containing all fitted peak parameters.  "
                  "X-values are spectra/workspace index.");
  declareProperty(
      Kernel::make_unique<WorkspaceProperty<MatrixWorkspace>>(
          "FittedPeaksWorkspace", "", Direction::Output,
          PropertyMode::Optional),
      "Name of the output matrix workspace with fitted peak. "
      "This output workspace have the same dimesion as the input workspace."
      "The Y values belonged to peaks to fit are replaced by fitted value. "
      "Values of estimated background are used if peak fails to be fit.");

  std::string addoutgrp("Analysis");
  setPropertyGroup("OutputPeakParametersWorkspace", addoutgrp);
  setPropertyGroup("FittedPeaksWorkspace", addoutgrp);

  return;
}

std::map<std::string, std::string> FitPeaks::validateInputs() {
  map<std::string, std::string> issues;

  // check that the peak parameters are in parallel properties
  bool haveCommonPeakParameters(false);
  std::vector<string> suppliedParameterNames =
      getProperty("PeakParameterNames");
  std::vector<double> peakParamValues = getProperty("PeakParameterValues");
  if ((!suppliedParameterNames.empty()) || (!peakParamValues.empty())) {
    haveCommonPeakParameters = true;
    if (suppliedParameterNames.size() != peakParamValues.size()) {
      issues["PeakParameterNames"] =
          "must have same number of values as PeakParameterValues";
      issues["PeakParameterValues"] =
          "must have same number of values as PeakParameterNames";
    }
  }

  // get the information out of the table
  std::string partablename = getPropertyValue("PeakParameterValueTable");
  if (!partablename.empty()) {
    if (haveCommonPeakParameters) {
      const std::string msg = "Parameter value table and initial parameter "
                              "name/value vectors cannot be given "
                              "simultanenously.";
      issues["PeakParameterValueTable"] = msg;
      issues["PeakParameterNames"] = msg;
      issues["PeakParameterValues"] = msg;
    } else {
      m_profileStartingValueTable = getProperty("PeakParameterValueTable");
      suppliedParameterNames = m_profileStartingValueTable->getColumnNames();
    }
  }

  // check that the suggested peak parameter names exist in the peak function
  if (!suppliedParameterNames.empty()) {
    std::string peakfunctiontype = getPropertyValue("PeakFunction");
    m_peakFunction = boost::dynamic_pointer_cast<IPeakFunction>(
        API::FunctionFactory::Instance().createFunction(peakfunctiontype));

    // put the names in a vector
    std::vector<string> functionParameterNames;
    for (size_t i = 0; i < m_peakFunction->nParams(); ++i)
      functionParameterNames.push_back(m_peakFunction->parameterName(i));
    // check that the supplied names are in the function
    // it is acceptable to be missing parameters
    bool failed = false;
    for (const auto &name : suppliedParameterNames) {
      if (std::find(functionParameterNames.begin(),
                    functionParameterNames.end(),
                    name) == functionParameterNames.end()) {
        failed = true;
        break;
      }
    }
    if (failed) {
      std::string msg = "Specified invalid parameter for peak function";
      if (haveCommonPeakParameters)
        issues["PeakParameterNames"] = msg;
      else
        issues["PeakParameterValueTable"] = msg;
    }
  }

  return issues;
}

//----------------------------------------------------------------------------------------------
void FitPeaks::exec() {
  // process inputs
  processInputs();

  // create output workspaces
  generateOutputPeakPositionWS();
  generateFittedParametersValueWorkspace();
  generateCalculatedPeaksWS();

  // fit peaks
  fitPeaks();

  // set the output workspaces to properites
  processOutputs();
}

//----------------------------------------------------------------------------------------------
void FitPeaks::processInputs() {
  // input workspaces
  m_inputMatrixWS = getProperty("InputWorkspace");

  if (m_inputMatrixWS->getAxis(0)->unit()->unitID() == "dSpacing")
    m_inputIsDSpace = true;
  else
    m_inputIsDSpace = false;

  // spectra to fit
  int start_wi = getProperty("StartWorkspaceIndex");
  if (isEmpty(start_wi))
    m_startWorkspaceIndex = 0;
  else
    m_startWorkspaceIndex = static_cast<size_t>(start_wi);

  // last spectrum's workspace index, which is included
  int stop_wi = getProperty("StopWorkspaceIndex");
  if (isEmpty(stop_wi))
    m_stopWorkspaceIndex = m_inputMatrixWS->getNumberHistograms() - 1;
  else {
    m_stopWorkspaceIndex = static_cast<size_t>(stop_wi);
    if (m_stopWorkspaceIndex > m_inputMatrixWS->getNumberHistograms() - 1)
      m_stopWorkspaceIndex = m_inputMatrixWS->getNumberHistograms() - 1;
  }

  // optimizer, cost function and fitting scheme
  m_minimizer = getPropertyValue("Minimizer");
  m_costFunction = getPropertyValue("CostFunction");
  m_fitPeaksFromRight = getProperty("FitFromRight");
  m_constrainPeaksPosition = getProperty("ConstrainPeakPositions");

  // Peak centers, tolerance and fitting range
  processInputPeakCenters();
  // check
  if (m_numPeaksToFit == 0)
    throw std::runtime_error("number of peaks to fit is zero.");
  // about how to estimate the peak width
  m_peakWidthPercentage = getProperty("PeakWidthPercent");
  if (isEmpty(m_peakWidthPercentage))
    m_peakWidthPercentage = -1;
  if (m_peakWidthPercentage >= 1.) // TODO
    throw std::runtime_error("PeakWidthPercent must be less than 1");
  g_log.debug() << "peak width/value = " << m_peakWidthPercentage << "\n";

  // set up background
  m_highBackground = getProperty("HighBackground");
  m_bkgdSimga = getProperty("FindBackgroundSigma");

  // Set up peak and background functions
  processInputFunctions();

  // about peak width and other peak parameter estimating method
  if (m_peakWidthPercentage > 0.)
    m_peakWidthEstimateApproach = EstimatePeakWidth::InstrumentResolution;
  else if (m_peakFunction->name() == "Gaussian")
    m_peakWidthEstimateApproach = EstimatePeakWidth::Observation;
  else
    m_peakWidthEstimateApproach = EstimatePeakWidth::NoEstimation;
  //  m_peakWidthEstimateApproach = EstimatePeakWidth::NoEstimation;
  g_log.debug() << "Process inputs [3] peak type: " << m_peakFunction->name()
                << ", background type: " << m_bkgdFunction->name() << "\n";

  processInputPeakTolerance();
  processInputFitRanges();

  return;
}

//----------------------------------------------------------------------------------------------
/** process inputs for peak profile and background
 */
void FitPeaks::processInputFunctions() {
  // peak functions
  std::string peakfunctiontype = getPropertyValue("PeakFunction");
  m_peakFunction = boost::dynamic_pointer_cast<IPeakFunction>(
      API::FunctionFactory::Instance().createFunction(peakfunctiontype));

  // background functions
  std::string bkgdfunctiontype = getPropertyValue("BackgroundType");
  std::string bkgdname;
  if (bkgdfunctiontype == "Linear")
    bkgdname = "LinearBackground";
  else if (bkgdfunctiontype == "Flat")
    bkgdname = "FlatBackground";
  else
    bkgdname = bkgdfunctiontype;
  m_bkgdFunction = boost::dynamic_pointer_cast<IBackgroundFunction>(
      API::FunctionFactory::Instance().createFunction(bkgdname));
  if (m_highBackground)
    m_linearBackgroundFunction =
        boost::dynamic_pointer_cast<IBackgroundFunction>(
            API::FunctionFactory::Instance().createFunction(
                "LinearBackground"));
  else
    m_linearBackgroundFunction = nullptr;

  // TODO check that both parameter names and values exist
  // input peak parameters
  std::string partablename = getPropertyValue("PeakParameterValueTable");
  m_peakParamNames = getProperty("PeakParameterNames");

  if (partablename.empty() && (!m_peakParamNames.empty())) {
    // use uniform starting value of peak parameters
    m_initParamValues = getProperty("PeakParameterValues");
    // convert the parameter name in string to parameter name in integer index
    convertParametersNameToIndex();
    m_uniformProfileStartingValue = true;
  } else if ((!partablename.empty()) && m_peakParamNames.empty()) {
    // use non-uniform starting value of peak parameters
    m_uniformProfileStartingValue = false;
    m_profileStartingValueTable = getProperty(partablename);
  } else {
    // user specifies nothing
    g_log.warning("Neither parameter value table nor initial "
                  "parameter name/value vectors is specified. Fitting might "
                  "not be reliable for peak profile other than Gaussian");
  }

  return;
}

//----------------------------------------------------------------------------------------------
/** process and check for inputs about peak fitting range (i.e., window)
 * Note: What is the output of the method?
 */
void FitPeaks::processInputFitRanges() {
  // get peak fit window
  std::vector<double> peakwindow = getProperty("FitWindowBoundaryList");
  std::string peakwindowname = getPropertyValue("FitPeakWindowWorkspace");
  API::MatrixWorkspace_const_sptr peakwindowws =
      getProperty("FitPeakWindowWorkspace");

  // in most case, calculate window by instrument resolution is False
  m_calculateWindowInstrument = false;

  if ((!peakwindow.empty()) && peakwindowname.empty()) {
    // Peak windows are uniform among spectra: use vector for peak windows
    m_uniformPeakWindows = true;

    // check peak positions
    if (!m_uniformPeakPositions)
      throw std::invalid_argument(
          "Uniform peak range/window requires uniform peak positions.");
    // check size
    if (peakwindow.size() != m_numPeaksToFit * 2)
      throw std::invalid_argument(
          "Peak window vector must be twice as large as number of peaks.");

    // set up window to m_peakWindowVector
    m_peakWindowVector.resize(m_numPeaksToFit);
    for (size_t i = 0; i < m_numPeaksToFit; ++i) {
      std::vector<double> peakranges(2);
      peakranges[0] = peakwindow[i * 2];
      peakranges[1] = peakwindow[i * 2 + 1];
      // check peak window (range) against peak centers
      if ((peakranges[0] < m_peakCenters[i]) &&
          (m_peakCenters[i] < peakranges[1])) {
        // pass check: set
        m_peakWindowVector[i] = peakranges;
      } else {
        // failed
        std::stringstream errss;
        errss << "Peak " << i
              << ": user specifies an invalid range and peak center against "
              << peakranges[0] << " < " << m_peakCenters[i] << " < "
              << peakranges[1];
        throw std::invalid_argument(errss.str());
      }
    } // END-FOR
    // END for uniform peak window
  } else if (peakwindow.empty() && peakwindowws != nullptr) {
    // use matrix workspace for non-uniform peak windows
    m_peakWindowWorkspace = getProperty("FitPeakWindowWorkspace");
    m_uniformPeakWindows = false;

    // check size
    if (m_peakWindowWorkspace->getNumberHistograms() ==
        m_inputMatrixWS->getNumberHistograms())
      m_partialWindowSpectra = false;
    else if (m_peakWindowWorkspace->getNumberHistograms() ==
             (m_stopWorkspaceIndex - m_startWorkspaceIndex + 1))
      m_partialWindowSpectra = true;
    else
      throw std::invalid_argument(
          "Peak window workspace has unmatched number of spectra");

    // check range for peak windows and peak positions
    size_t window_index_start(0);
    if (m_partialWindowSpectra)
      window_index_start = m_startWorkspaceIndex;
    size_t center_index_start(0);
    if (m_partialSpectra)
      center_index_start = m_startWorkspaceIndex;

    // check each spectrum whether the window is defined with the correct size
    for (size_t wi = 0; wi < m_peakWindowWorkspace->getNumberHistograms();
         ++wi) {
      // check size
      if (m_peakWindowWorkspace->y(wi).size() != m_numPeaksToFit * 2) {
        std::stringstream errss;
        errss << "Peak window workspace index " << wi
              << " has incompatible number of fit windows (x2) "
              << m_peakWindowWorkspace->y(wi).size()
              << "with the number of peaks " << m_numPeaksToFit << " to fit.";
        throw std::invalid_argument(errss.str());
      }

      // check window range against peak center
      size_t window_index = window_index_start + wi;
      size_t center_index = window_index - center_index_start;

      for (size_t ipeak = 0; ipeak < m_numPeaksToFit; ++ipeak) {
        double left_w_bound =
            m_peakWindowWorkspace->x(wi)[ipeak * 2]; // TODO getting on y
        double right_w_bound = m_peakWindowWorkspace->x(wi)[ipeak * 2 + 1];
        double center = m_peakCenterWorkspace->x(center_index)[ipeak];
        if (!(left_w_bound < center && center < right_w_bound)) {
          std::stringstream errss;
          errss << "Workspace index " << wi
                << " has incompatible peak window (" // <<<<<<< HERE!!!!!!!!!
                << left_w_bound << ", " << right_w_bound << ") with " << ipeak
                << "-th expected peak's center " << center;
          throw std::runtime_error(errss.str());
        }
      }
    }
  } else if (peakwindow.empty()) {
    // no peak window is defined, then the peak window will be estimated by
    // delta(D)/D
    if (m_inputIsDSpace && m_peakWidthPercentage > 0)
      m_calculateWindowInstrument = true;
    else
      throw std::invalid_argument("Without definition of peak window, the "
                                  "input workspace must be in unit of dSpacing "
                                  "and Delta(D)/D must be given!");

  } else {
    // non-supported situation
    throw std::invalid_argument("One and only one of peak window array and "
                                "peak window workspace can be specified.");
  }

  return;
}

//----------------------------------------------------------------------------------------------
/** Processing peaks centers and fitting tolerance information from input.  the
 * parameters that are
 * set including
 * 1. m_peakCenters/m_peakCenterWorkspace/m_uniformPeakPositions
 * (bool)/m_partialSpectra (bool)
 * 2. m_peakPosTolerances (vector)
 * 3. m_numPeaksToFit
 */
void FitPeaks::processInputPeakCenters() {
  // peak centers
  m_peakCenters = getProperty("PeakCenters");
  API::MatrixWorkspace_const_sptr peakcenterws =
      getProperty("PeakCentersWorkspace");
  if (!peakcenterws)
    g_log.error("There is no peak center workspace");

  std::string peakpswsname = getPropertyValue("PeakCentersWorkspace");
  if ((!m_peakCenters.empty()) && peakcenterws == nullptr) {
    // peak positions are uniform among all spectra
    m_uniformPeakPositions = true;
    // number of peaks to fit!
    m_numPeaksToFit = m_peakCenters.size();
  } else if (m_peakCenters.empty() && peakcenterws != nullptr) {
    // peak positions can be different among spectra
    m_uniformPeakPositions = false;
    m_peakCenterWorkspace = getProperty("PeakCentersWorkspace");
    // number of peaks to fit!
    m_numPeaksToFit = m_peakCenterWorkspace->x(0).size();

    // check matrix worksapce for peak positions
    const size_t numhist = m_peakCenterWorkspace->getNumberHistograms();
    if (numhist == m_inputMatrixWS->size())
      m_partialSpectra = false;
    else if (numhist == m_stopWorkspaceIndex - m_startWorkspaceIndex + 1)
      m_partialSpectra = true;
    else
      throw std::invalid_argument(
          "Input peak center workspace has wrong number of spectra.");

  } else {
    std::stringstream errss;
    errss << "One and only one in 'PeakCenters' (vector) and "
             "'PeakCentersWorkspace' shall be given. "
          << "'PeakCenters' has size " << m_peakCenters.size()
          << ", and name of peak center workspace "
          << "is " << peakpswsname;
    throw std::invalid_argument(errss.str());
  }

  return;
}

//----------------------------------------------------------------------------------------------
/** Processing peak fitting tolerance information from input.  The parameters
 * that are
 * set including
 * 2. m_peakPosTolerances (vector)
 */
void FitPeaks::processInputPeakTolerance() {
  // check code integrity
  if (m_numPeaksToFit == 0)
    throw std::runtime_error("ProcessInputPeakTolerance() must be called after "
                             "ProcessInputPeakCenters()");

  // peak tolerance
  m_peakPosTolerances = getProperty("PositionTolerance");

  if (m_peakPosTolerances.empty()) {
    // case 2, 3, 4
    m_peakPosTolerances.clear();
    m_peakPosTolCase234 = true;
  } else if (m_peakPosTolerances.size() == 1) {
    // only 1 uniform peak position tolerance is defined: expand to all peaks
    double peak_tol = m_peakPosTolerances[0];
    m_peakPosTolerances.resize(m_numPeaksToFit, peak_tol);
  } else if (m_peakPosTolerances.size() != m_numPeaksToFit) {
    // not uniform but number of peaks does not match
    g_log.error() << "number of peak position tolerance "
                  << m_peakPosTolerances.size()
                  << " is not same as number of peaks " << m_numPeaksToFit
                  << "\n";
    throw std::runtime_error("Number of peak position tolerances and number of "
                             "peaks to fit are inconsistent.");
  }

  // minimum peak height: set default to zero
  m_minPeakHeight = getProperty("MinimumPeakHeight");
  if (isEmpty(m_minPeakHeight) || m_minPeakHeight < 0.)
    m_minPeakHeight = 0.;

  return;
}

//----------------------------------------------------------------------------------------------
/** Convert the input initial parameter name/value to parameter index/value for
 * faster access
 * according to the parameter name and peak profile function
 * Output: m_initParamIndexes will be set up
 */
void FitPeaks::convertParametersNameToIndex() {
  // get a map for peak profile parameter name and parameter index
  std::map<std::string, size_t> parname_index_map;
  for (size_t iparam = 0; iparam < m_peakFunction->nParams(); ++iparam)
    parname_index_map.insert(
        std::make_pair(m_peakFunction->parameterName(iparam), iparam));

  // define peak parameter names (class variable) if using table
  if (m_profileStartingValueTable)
    m_peakParamNames = m_profileStartingValueTable->getColumnNames();

  // map the input parameter names to parameter indexes
  for (const auto &paramName : m_peakParamNames) {
    std::map<std::string, size_t>::iterator locator =
        parname_index_map.find(paramName);
    if (locator != parname_index_map.end())
      m_initParamIndexes.push_back(locator->second);
    else {
      // a parameter name that is not defined in the peak profile function.  An
      // out-of-range index is thus set to this
      g_log.warning() << "Given peak parameter " << paramName
                      << " is not an allowed parameter of peak "
                         "function " << m_peakFunction->name() << "\n";
      m_initParamIndexes.push_back(m_peakFunction->nParams() * 10);
    }
  }

  return;
}

//----------------------------------------------------------------------------------------------
/** main method to fit peaks among all
 */
void FitPeaks::fitPeaks() {
  API::Progress prog(this, 0., 1.,
                     m_stopWorkspaceIndex - m_startWorkspaceIndex);

  // cppcheck-suppress syntaxError
  PRAGMA_OMP(parallel for schedule(dynamic, 1) )
  for (int wi = static_cast<int>(m_startWorkspaceIndex);
       wi <= static_cast<int>(m_stopWorkspaceIndex); ++wi) {

    PARALLEL_START_INTERUPT_REGION

    // peaks to fit
    std::vector<double> expected_peak_centers =
        getExpectedPeakPositions(static_cast<size_t>(wi));

    // initialize output for this
    size_t numfuncparams =
        m_peakFunction->nParams() + m_bkgdFunction->nParams();
    boost::shared_ptr<FitPeaksAlgorithm::PeakFitResult> fit_result =
        boost::make_shared<FitPeaksAlgorithm::PeakFitResult>(m_numPeaksToFit,
                                                             numfuncparams);

    fitSpectrumPeaks(static_cast<size_t>(wi), expected_peak_centers,
                     fit_result);

    PARALLEL_CRITICAL(FindPeaks_WriteOutput) {
      writeFitResult(static_cast<size_t>(wi), expected_peak_centers,
                     fit_result);
    }
    prog.report();

    PARALLEL_END_INTERUPT_REGION
  }

  PARALLEL_CHECK_INTERUPT_REGION
}

namespace {
double numberCounts(const Histogram &histogram) {
  double total = 0.;
  for (const auto &value : histogram.y())
    total += std::fabs(value);
  return total;
}

double numberCounts(const Histogram &histogram, const double xmin,
                    const double xmax) {
  const auto &vector_x = histogram.points();

  // determine left boundary
  std::vector<double>::const_iterator start_iter = vector_x.begin();
  if (xmin > vector_x.front())
    start_iter = std::lower_bound(vector_x.begin(), vector_x.end(), xmin);
  if (start_iter == vector_x.end())
    return 0.; // past the end of the data means nothing to integrate
  // determine right boundary
  std::vector<double>::const_iterator stop_iter = vector_x.end();
  if (xmax < vector_x.back()) // will set at end of vector if too large
    stop_iter = std::lower_bound(start_iter, stop_iter, xmax);

  // convert to indexes to sum over y
  size_t start_index = static_cast<size_t>(start_iter - vector_x.begin());
  size_t stop_index = static_cast<size_t>(stop_iter - vector_x.begin());

  // integrate
  double total = 0.;
  for (size_t i = start_index; i < stop_index; ++i)
    total += std::fabs(histogram.y()[i]);
  return total;
}
}

//----------------------------------------------------------------------------------------------
/** Fit peaks across one single spectrum
 */
void FitPeaks::fitSpectrumPeaks(
    size_t wi, const std::vector<double> &expected_peak_centers,
    boost::shared_ptr<FitPeaksAlgorithm::PeakFitResult> fit_result) {
  if (numberCounts(m_inputMatrixWS->histogram(wi)) <= m_minPeakHeight) {
    for (size_t i = 0; i < fit_result->getNumberPeaks(); ++i)
      fit_result->setBadRecord(i, -1.);
    return; // don't do anything
  }

  // Set up sub algorithm Fit for peak and background
  IAlgorithm_sptr peak_fitter; // both peak and background (combo)
  try {
    peak_fitter = createChildAlgorithm("Fit", -1, -1, false);
  } catch (Exception::NotFoundError &) {
    std::stringstream errss;
    errss << "The FitPeak algorithm requires the CurveFitting library";
    g_log.error(errss.str());
    throw std::runtime_error(errss.str());
  }

  // Clone the function
  IPeakFunction_sptr peakfunction =
      boost::dynamic_pointer_cast<API::IPeakFunction>(m_peakFunction->clone());
  IBackgroundFunction_sptr bkgdfunction =
      boost::dynamic_pointer_cast<API::IBackgroundFunction>(
          m_bkgdFunction->clone());

  // set up properties of algorithm (reference) 'Fit'
  peak_fitter->setProperty("Minimizer", m_minimizer);
  peak_fitter->setProperty("CostFunction", m_costFunction);
  peak_fitter->setProperty("CalcErrors", true);

  const double x0 = m_inputMatrixWS->histogram(wi).x().front();
  const double xf = m_inputMatrixWS->histogram(wi).x().back();
  for (size_t fit_index = 0; fit_index < m_numPeaksToFit; ++fit_index) {
    // convert fit index to peak index (in ascending order)
    size_t peak_index(fit_index);
    if (m_fitPeaksFromRight)
      peak_index = m_numPeaksToFit - fit_index - 1;

    // get expected peak position
    double expected_peak_pos = expected_peak_centers[peak_index];
    double cost(DBL_MAX);
    if (expected_peak_pos <= x0 || expected_peak_pos >= xf) {
      // out of range and there won't be any fit
      peakfunction->setIntensity(0);
      peakfunction->setCentre(expected_peak_pos);
    } else {
      // find out the peak position to fit
      std::pair<double, double> peak_window_i =
          getPeakFitWindow(wi, peak_index);

      bool observe_peak_width =
          decideToEstimatePeakWidth(fit_index == 0, peakfunction);

      // do fitting with peak and background function (no analysis at this
      // point)
      cost =
          fitIndividualPeak(wi, peak_fitter, expected_peak_pos, peak_window_i,
                            observe_peak_width, peakfunction, bkgdfunction);
    }

    // process fitting result
    FitPeaksAlgorithm::FitFunction fit_function;
    fit_function.peakfunction = peakfunction;
    fit_function.bkgdfunction = bkgdfunction;

    processSinglePeakFitResult(wi, peak_index, cost, expected_peak_centers,
                               fit_function, fit_result); // sets the record
  }

  return;
}

//----------------------------------------------------------------------------------------------
/** Decide whether to estimate peak width.  If not, then set the width related
 * peak parameters from user specified starting value
 */
bool FitPeaks::decideToEstimatePeakWidth(
    const bool firstPeakInSpectrum, API::IPeakFunction_sptr peak_function) {
  bool observe_peak_width(false);

  if (!m_initParamIndexes.empty()) {
    // user specifies starting value of peak parameters
    if (firstPeakInSpectrum) {
      // TODO just set the parameter values in a vector and loop over it
      // first peak.  using the user-specified value
      for (size_t i = 0; i < m_initParamIndexes.size(); ++i) {
        size_t param_index = m_initParamIndexes[i];
        double param_value = m_initParamValues[i];
        peak_function->setParameter(param_index, param_value);
      }
    } else {
      // using the fitted paramters from the previous fitting result
      // do noting
    }
  } else {
    // by observation
    observe_peak_width = true;
  }

  return observe_peak_width;
}

//----------------------------------------------------------------------------------------------
/** retrieve the fitted peak information from functions and set to output
 * vectors
 */
void FitPeaks::processSinglePeakFitResult(
    size_t wsindex, size_t peakindex, const double cost,
    const std::vector<double> &expected_peak_positions,
    FitPeaksAlgorithm::FitFunction fitfunction,
    boost::shared_ptr<FitPeaksAlgorithm::PeakFitResult> fit_result) {
  // determine peak position tolerance
  double postol(DBL_MAX);
  bool case23(false);
  if (m_peakPosTolCase234) {
    // peak tolerance is not defined
    if (m_numPeaksToFit == 1) {
      // case (d) one peak only
      postol = m_inputMatrixWS->histogram(wsindex).x().back() -
               m_inputMatrixWS->histogram(wsindex).x().front();