CrystalFieldFunction.cpp

#include "MantidCurveFitting/Functions/CrystalElectricField.h"
#include "MantidCurveFitting/Functions/CrystalFieldFunction.h"
#include "MantidCurveFitting/Functions/CrystalFieldHeatCapacity.h"
#include "MantidCurveFitting/Functions/CrystalFieldMagnetisation.h"
#include "MantidCurveFitting/Functions/CrystalFieldMoment.h"
#include "MantidCurveFitting/Functions/CrystalFieldPeakUtils.h"
#include "MantidCurveFitting/Functions/CrystalFieldPeaks.h"
#include "MantidCurveFitting/Functions/CrystalFieldSusceptibility.h"

#include "MantidAPI/FunctionFactory.h"
#include "MantidAPI/IConstraint.h"
#include "MantidAPI/IFunction1D.h"
#include "MantidAPI/IPeakFunction.h"
#include "MantidAPI/MultiDomainFunction.h"
#include "MantidAPI/ParameterTie.h"

#include "MantidKernel/Exception.h"

#include <boost/make_shared.hpp>
#include <boost/optional.hpp>
#include <boost/regex.hpp>
#include <iostream>
#include <limits>

namespace Mantid {
namespace CurveFitting {
namespace Functions {

using namespace CurveFitting;
using namespace Kernel;
using namespace API;

DECLARE_FUNCTION(CrystalFieldFunction)

namespace {

const std::string ION_PREFIX("ion");
const std::string SPECTRUM_PREFIX("sp");
const std::string BACKGROUND_PREFIX("bg");
const std::string PEAK_PREFIX("pk");

// Regex for names of attributes/parameters for a particular spectrum
// Example: sp1.FWHMX
const boost::regex SPECTRUM_ATTR_REGEX(SPECTRUM_PREFIX + "([0-9]+)\\.(.+)");
// Regex for names of attributes/parameters for a background
// Example: bg.A1
const boost::regex BACKGROUND_ATTR_REGEX(BACKGROUND_PREFIX + "\\.(.+)");
// Regex for names of attributes/parameters for peaks
// Example: pk1.PeakCentre
const boost::regex PEAK_ATTR_REGEX(PEAK_PREFIX + "([0-9]+)\\.(.+)");
// Regex for names of attributes/parameters for peaks
// Example: ion1.pk0.PeakCentre
const boost::regex ION_ATTR_REGEX(ION_PREFIX + "([0-9]+)\\.(.+)");
// Regex for names of attributes/parameters for physical properties
// Example: cv.ScaleFactor
const boost::regex
    PHYS_PROP_ATTR_REGEX("((ion[0-9]+\\.)?(cv|chi|mh|mt))\\.(.+)");

/// Define the source function for CrystalFieldFunction.
/// Its function() method is not needed.
class Peaks : public CrystalFieldPeaksBase, public API::IFunctionGeneral {
public:
  Peaks() : CrystalFieldPeaksBase() {}
  std::string name() const override { return "Peaks"; }
  size_t getNumberDomainColumns() const override {
    throw Exception::NotImplementedError(
        "This method is intentionally not implemented.");
  }
  size_t getNumberValuesPerArgument() const override {
    throw Exception::NotImplementedError(
        "This method is intentionally not implemented.");
  }
  void functionGeneral(const API::FunctionDomainGeneral &,
                       API::FunctionValues &) const override {
    throw Exception::NotImplementedError(
        "This method is intentionally not implemented.");
  }
  std::vector<size_t> m_IntensityScalingIdx;
  std::vector<size_t> m_PPLambdaIdxChild;
  std::vector<size_t> m_PPLambdaIdxSelf;
  /// Declare the intensity scaling parameters: one per spectrum.
  void declareIntensityScaling(size_t nSpec) {
    m_IntensityScalingIdx.clear();
    m_PPLambdaIdxChild.resize(nSpec, -1);
    m_PPLambdaIdxSelf.resize(nSpec, -1);
    for (size_t i = 0; i < nSpec; ++i) {
      auto si = std::to_string(i);
      try { // If parameter has already been declared, don't declare it.
        declareParameter("IntensityScaling" + si, 1.0,
                         "Intensity scaling factor for spectrum " + si);
      } catch (std::invalid_argument &) {
      }
      m_IntensityScalingIdx.push_back(parameterIndex("IntensityScaling" + si));
    }
  }
};

} // namespace

/// Constructor
CrystalFieldFunction::CrystalFieldFunction()
    : IFunction(), m_nControlParams(0), m_nControlSourceParams(0),
      m_dirtyTarget(true) {}

// Evaluates the function
void CrystalFieldFunction::function(const FunctionDomain &domain,
                                    FunctionValues &values) const {
  updateTargetFunction();
  if (!m_target) {
    throw std::logic_error(
        "FunctionGenerator failed to generate target function.");
  }
  m_target->function(domain, values);
}

/// Set the source function
/// @param source :: New source function.
void CrystalFieldFunction::setSource(IFunction_sptr source) const {
  m_source = source;
}

size_t CrystalFieldFunction::getNumberDomains() const {
  if (!m_target) {
    buildTargetFunction();
  }
  if (!m_target) {
    throw std::runtime_error("Failed to build target function.");
  }
  return m_target->getNumberDomains();
}

std::vector<IFunction_sptr>
CrystalFieldFunction::createEquivalentFunctions() const {
  checkTargetFunction();
  std::vector<IFunction_sptr> funs;
  auto &composite = dynamic_cast<CompositeFunction &>(*m_target);
  for (size_t i = 0; i < composite.nFunctions(); ++i) {
    auto fun = composite.getFunction(i);
    auto cfun = dynamic_cast<CompositeFunction *>(fun.get());
    if (cfun) {
      cfun->checkFunction();
    }
    funs.push_back(fun);
  }
  return funs;
}

/// Set i-th parameter
void CrystalFieldFunction::setParameter(size_t i, const double &value,
                                        bool explicitlySet) {
  checkSourceFunction();
  if (i < m_nControlParams) {
    m_control.setParameter(i, value, explicitlySet);
    m_dirtyTarget = true;
  } else if (i < m_nControlSourceParams) {
    m_source->setParameter(i - m_nControlParams, value, explicitlySet);
    m_dirtyTarget = true;
  } else {
    checkTargetFunction();
    m_target->setParameter(i - m_nControlSourceParams, value, explicitlySet);
  }
}

/// Set i-th parameter description
void CrystalFieldFunction::setParameterDescription(
    size_t i, const std::string &description) {
  checkSourceFunction();
  if (i < m_nControlParams) {
    m_control.setParameterDescription(i, description);
  } else if (i < m_nControlSourceParams) {
    m_source->setParameterDescription(i - m_nControlParams, description);
  } else {
    checkTargetFunction();
    m_target->setParameterDescription(i - m_nControlSourceParams, description);
  }
}

/// Get i-th parameter
double CrystalFieldFunction::getParameter(size_t i) const {
  checkSourceFunction();
  checkTargetFunction();
  if (i < m_nControlParams) {
    return m_control.getParameter(i);
  } else if (i < m_nControlSourceParams) {
    return m_source->getParameter(i - m_nControlParams);
  } else {
    return m_target->getParameter(i - m_nControlSourceParams);
  }
}

/// Check if function has a parameter with this name.
bool CrystalFieldFunction::hasParameter(const std::string &name) const {
  try {
    parameterIndex(name);
    return true;
  } catch (std::invalid_argument &) {
    return false;
  }
}

/// Set parameter by name.
void CrystalFieldFunction::setParameter(const std::string &name,
                                        const double &value,
                                        bool explicitlySet) {
  try {
    auto index = parameterIndex(name);
    setParameter(index, value, explicitlySet);
  } catch (std::invalid_argument &) {
    // Allow ignoring peak parameters: the peak may not exist.
    boost::smatch match;
    if (!boost::regex_search(name, match, PEAK_ATTR_REGEX)) {
      throw;
    }
  }
}

/// Set description of parameter by name.
void CrystalFieldFunction::setParameterDescription(
    const std::string &name, const std::string &description) {
  auto index = parameterIndex(name);
  setParameterDescription(index, description);
}

/// Get parameter by name.
double CrystalFieldFunction::getParameter(const std::string &name) const {
  auto index = parameterIndex(name);
  return getParameter(index);
}

/// Total number of parameters
size_t CrystalFieldFunction::nParams() const {
  if (!m_source) {
    // This method can be called on an uninitialised function (by tests for
    // example).
    // Return 0 so no exception is thrown an it should prevent attemts to access
    // parameters.
    return 0;
  }
  checkSourceFunction();
  checkTargetFunction();
  return m_nControlSourceParams + m_target->nParams();
}

/// Returns the index of a parameter with a given name
/// @param name :: Name of a parameter.
size_t CrystalFieldFunction::parameterIndex(const std::string &name) const {
  checkSourceFunction();
  checkTargetFunction();
  if (nParams() != m_mapIndices2Names.size()) {
    makeMaps();
  }
  auto found = m_mapNames2Indices.find(name);
  if (found == m_mapNames2Indices.end()) {
    throw std::invalid_argument("CrystalFieldFunction parameter not found: " +
                                name);
  }
  return found->second;
}

/// Returns the name of parameter i
std::string CrystalFieldFunction::parameterName(size_t i) const {
  if (i >= nParams()) {
    throw std::invalid_argument("CrystalFieldFunction's parameter index " +
                                std::to_string(i) + " is out of range " +
                                std::to_string(nParams()));
  }
  checkSourceFunction();
  checkTargetFunction();
  if (nParams() != m_mapIndices2Names.size()) {
    makeMaps();
  }
  return m_mapIndices2Names[i];
}

/// Returns the description of parameter i
std::string CrystalFieldFunction::parameterDescription(size_t i) const {
  checkSourceFunction();
  checkTargetFunction();
  if (i < m_nControlParams) {
    return m_control.parameterDescription(i);
  } else if (i < m_nControlSourceParams) {
    return m_source->parameterDescription(i - m_nControlParams);
  } else {
    return m_target->parameterDescription(i - m_nControlSourceParams);
  }
}

/// Checks if a parameter has been set explicitly
bool CrystalFieldFunction::isExplicitlySet(size_t i) const {
  checkSourceFunction();
  checkTargetFunction();
  if (i < m_nControlParams) {
    return m_control.isExplicitlySet(i);
  } else if (i < m_nControlSourceParams) {
    return m_source->isExplicitlySet(i - m_nControlParams);
  } else {
    return m_target->isExplicitlySet(i - m_nControlSourceParams);
  }
}

/// Get the fitting error for a parameter
double CrystalFieldFunction::getError(size_t i) const {
  checkSourceFunction();
  checkTargetFunction();
  if (i < m_nControlParams) {
    return m_control.getError(i);
  } else if (i < m_nControlSourceParams) {
    return m_source->getError(i - m_nControlParams);
  } else {
    return m_target->getError(i - m_nControlSourceParams);
  }
}

/// Set the fitting error for a parameter
void CrystalFieldFunction::setError(size_t i, double err) {
  checkSourceFunction();
  checkTargetFunction();
  if (i < m_nControlParams) {
    m_control.setError(i, err);
  } else if (i < m_nControlSourceParams) {
    m_source->setError(i - m_nControlParams, err);
  } else {
    m_target->setError(i - m_nControlSourceParams, err);
  }
}

/// Change status of parameter
void CrystalFieldFunction::setParameterStatus(
    size_t i, IFunction::ParameterStatus status) {
  checkSourceFunction();
  checkTargetFunction();
  if (i < m_nControlParams) {
    m_control.setParameterStatus(i, status);
  } else if (i < m_nControlSourceParams) {
    m_source->setParameterStatus(i - m_nControlParams, status);
  } else {
    m_target->setParameterStatus(i - m_nControlSourceParams, status);
  }
}

/// Get status of parameter
IFunction::ParameterStatus
CrystalFieldFunction::getParameterStatus(size_t i) const {
  checkSourceFunction();
  checkTargetFunction();
  if (i < m_nControlParams) {
    return m_control.getParameterStatus(i);
  } else if (i < m_nControlSourceParams) {
    return m_source->getParameterStatus(i - m_nControlParams);
  } else {
    return m_target->getParameterStatus(i - m_nControlSourceParams);
  }
}

/// Return parameter index from a parameter reference.
size_t
CrystalFieldFunction::getParameterIndex(const ParameterReference &ref) const {
  checkSourceFunction();
  checkTargetFunction();
  if (ref.getLocalFunction() == this) {
    return ref.getLocalIndex();
  }
  auto index = m_control.getParameterIndex(ref);
  if (index < m_nControlParams) {
    return index;
  }
  index = m_source->getParameterIndex(ref);
  if (index < m_source->nParams()) {
    return index + m_nControlParams;
  }
  return m_target->getParameterIndex(ref) + m_nControlSourceParams;
}

/// Set up the function for a fit.
void CrystalFieldFunction::setUpForFit() {
  checkSourceFunction();
  updateTargetFunction();
  IFunction::setUpForFit();
}

/// Declare a new parameter
void CrystalFieldFunction::declareParameter(const std::string &, double,
                                            const std::string &) {
  throw Kernel::Exception::NotImplementedError(
      "CrystalFieldFunction cannot have its own parameters.");
}

/// Build and cache the attribute names
void CrystalFieldFunction::buildAttributeNames() const {
  checkSourceFunction();
  checkTargetFunction();
  if (!m_attributeNames.empty()) {
    return;
  }
  m_attributeNames = IFunction::getAttributeNames();
  auto controlAttributeNames = m_control.getAttributeNames();
  // Lambda function that moves a attribute name from controlAttributeNames
  // to attNames.
  auto moveAttributeName = [&](const std::string &name) {
    auto iterFound = std::find(controlAttributeNames.begin(),
                               controlAttributeNames.end(), name);
    if (iterFound != controlAttributeNames.end()) {
      controlAttributeNames.erase(iterFound);
      m_attributeNames.push_back(name);
    }
  };
  // Prepend a prefix to attribute names, ignore NumDeriv attribute.
  auto prependPrefix =
      [&](const std::string &prefix, const std::vector<std::string> &names) {
        for (auto name : names) {
          if (name == "NumDeriv")
            continue;
          name.insert(name.begin(), prefix.begin(), prefix.end());
          m_attributeNames.push_back(name);
        }
      };
  // These names must appear first and in this order in the output vector
  moveAttributeName("Ions");
  moveAttributeName("Symmetries");
  moveAttributeName("Temperatures");
  moveAttributeName("Background");
  // Copy the rest of the names
  m_attributeNames.insert(m_attributeNames.end(), controlAttributeNames.begin(),
                          controlAttributeNames.end());
  // Get
  for (size_t iSpec = 0; iSpec < m_control.nFunctions(); ++iSpec) {
    std::string prefix(SPECTRUM_PREFIX);
    prefix.append(std::to_string(iSpec)).append(".");
    auto names = m_control.getFunction(iSpec)->getAttributeNames();
    for (auto &name : names) {
      name.insert(name.begin(), prefix.begin(), prefix.end());
    }
    m_attributeNames.insert(m_attributeNames.end(), names.begin(), names.end());
  }
  // Attributes of physical properties
  for (size_t iSpec = nSpectra(); iSpec < m_target->nFunctions(); ++iSpec) {
    auto fun = m_target->getFunction(iSpec).get();
    auto compositePhysProp = dynamic_cast<CompositeFunction *>(fun);
    if (compositePhysProp) {
      // Multi-site case
      std::string physPropPrefix(compositePhysProp->getFunction(0)->name());
      physPropPrefix.append(".");
      for (size_t ion = 0; ion < compositePhysProp->nFunctions(); ++ion) {
        std::string prefix(ION_PREFIX);
        prefix.append(std::to_string(ion)).append(".").append(physPropPrefix);
        auto names = compositePhysProp->getFunction(ion)->getAttributeNames();
        prependPrefix(prefix, names);
      }
    } else {
      // Single-site
      std::string prefix(fun->name());
      prefix.append(".");
      auto names = fun->getAttributeNames();
      prependPrefix(prefix, names);
    }
  }
}

/// Returns the number of attributes associated with the function
size_t CrystalFieldFunction::nAttributes() const {
  buildAttributeNames();
  return m_attributeNames.size();
}

/// Returns a list of attribute names
std::vector<std::string> CrystalFieldFunction::getAttributeNames() const {
  buildAttributeNames();
  return m_attributeNames;
}

/// Return a value of attribute attName
/// @param attName :: Name of an attribute.
IFunction::Attribute
CrystalFieldFunction::getAttribute(const std::string &attName) const {
  auto attRef = getAttributeReference(attName);
  if (attRef.first == nullptr) {
    // This will throw an exception because attribute doesn't exist
    return IFunction::getAttribute(attName);
  }
  return attRef.first->getAttribute(attRef.second);
}

/// Perform custom actions on setting certain attributes.
void CrystalFieldFunction::setAttribute(const std::string &attName,
                                        const Attribute &attr) {
  auto attRef = getAttributeReference(attName);
  if (attRef.first == nullptr) {
    // This will throw an exception because attribute doesn't exist
    IFunction::setAttribute(attName, attr);
  } else if (attRef.first == &m_control) {
    m_source.reset();
  }
  attRef.first->setAttribute(attRef.second, attr);
}

/// Check if attribute attName exists
bool CrystalFieldFunction::hasAttribute(const std::string &attName) const {
  auto attRef = getAttributeReference(attName);
  if (attRef.first == nullptr) {
    return false;
  }
  return attRef.first->hasAttribute(attRef.second);
}

/// Get a reference to an attribute.
/// @param attName :: A name of an attribute. It can be a code rather than an
/// actual name.  This method interprets the code and finds the function and
/// attribute it refers to.
/// @returns :: A pair (IFunction, attribute_name) where attribute_name is a
/// name that the IFunction has.
std::pair<API::IFunction *, std::string>
CrystalFieldFunction::getAttributeReference(const std::string &attName) const {
  boost::smatch match;
  if (boost::regex_match(attName, match, SPECTRUM_ATTR_REGEX)) {
    auto i = std::stoul(match[1]);
    auto name = match[2].str();
    if (m_control.nFunctions() == 0) {
      m_control.buildControls();
    }
    if (name == "FWHMX" || name == "FWHMY") {
      if (i < m_control.nFunctions()) {
        return std::make_pair(m_control.getFunction(i).get(), name);
      } else {
        return std::make_pair(nullptr, "");
      }
    }
    return std::make_pair(nullptr, "");
  } else if (boost::regex_match(attName, match, PHYS_PROP_ATTR_REGEX)) {
    auto prop = match[1].str();
    auto name = match[4].str();
    auto propIt = m_mapPrefixes2PhysProps.find(prop);
    if (propIt != m_mapPrefixes2PhysProps.end()) {
      return std::make_pair(propIt->second.get(), name);
    }
    return std::make_pair(nullptr, "");
  }
  return std::make_pair(&m_control, attName);
}

/// Get number of the number of spectra (excluding phys prop data).
size_t CrystalFieldFunction::nSpectra() const {
  auto nFuns = m_control.nFunctions();
  return nFuns;
}

/// Get the tie for i-th parameter
ParameterTie *CrystalFieldFunction::getTie(size_t i) const {
  checkSourceFunction();
  checkTargetFunction();
  auto tie = IFunction::getTie(i);
  if (tie) {
    return tie;
  }
  if (i < m_nControlParams) {
    tie = m_control.getTie(i);
  } else if (i < m_nControlSourceParams) {
    tie = m_source->getTie(i - m_nControlParams);
  } else {
    tie = m_target->getTie(i - m_nControlSourceParams);
  }
  return tie;
}

/// Get the i-th constraint
IConstraint *CrystalFieldFunction::getConstraint(size_t i) const {
  checkSourceFunction();
  auto constraint = IFunction::getConstraint(i);
  if (constraint == nullptr) {
    if (i < m_nControlParams) {
      constraint = m_control.getConstraint(i);
    } else if (i < m_nControlSourceParams) {
      constraint = m_source->getConstraint(i - m_nControlParams);
    } else {
      checkTargetFunction();
      constraint = m_target->getConstraint(i - m_nControlSourceParams);
    }
  }
  return constraint;
}

/// Check if the function is set up for a multi-site calculations.
/// (Multiple ions defined)
bool CrystalFieldFunction::isMultiSite() const {
  return m_control.isMultiSite();
}

/// Check if the function is set up for a multi-spectrum calculations
/// (Multiple temperatures defined)
bool CrystalFieldFunction::isMultiSpectrum() const {
  return m_control.isMultiSpectrum();
}

/// Check if the spectra have a background.
bool CrystalFieldFunction::hasBackground() const {
  if (!hasAttribute("Background")) {
    return false;
  }
  return !getAttribute("Background").isEmpty();
}

/// Check if there are peaks (there is at least one spectrum).
bool CrystalFieldFunction::hasPeaks() const { return m_control.hasPeaks(); }

/// Check if there are any phys. properties.
bool CrystalFieldFunction::hasPhysProperties() const {
  return m_control.hasPhysProperties();
}

/// Get a reference to the source function if it's composite
API::CompositeFunction &CrystalFieldFunction::compositeSource() const {
  auto composite = dynamic_cast<CompositeFunction *>(m_source.get());
  if (composite == nullptr) {
    throw std::logic_error("Source of CrystalFieldFunction is not composite.");
  }
  return *composite;
}

/// Build source function if necessary.
void CrystalFieldFunction::checkSourceFunction() const {
  if (!m_source) {
    buildSourceFunction();
  }
}

/// Build the source function
void CrystalFieldFunction::buildSourceFunction() const {
  setSource(m_control.buildSource());
  m_nControlParams = m_control.nParams();
  m_nControlSourceParams = m_nControlParams + m_source->nParams();
}

/// Update spectrum function if necessary.
void CrystalFieldFunction::checkTargetFunction() const {
  if (m_dirtyTarget) {
    updateTargetFunction();
  }
  if (!m_target) {
    throw std::logic_error(
        "CrystalFieldFunction failed to generate target function.");
  }
}

/// Uses source to calculate peak centres and intensities
/// then populates m_spectrum with peaks of type given in PeakShape attribute.
void CrystalFieldFunction::buildTargetFunction() const {
  checkSourceFunction();
  m_dirtyTarget = false;
  if (isMultiSite()) {
    buildMultiSite();
  } else {
    buildSingleSite();
  }
  m_attributeNames.clear();
}

/// Build the target function in a single site case.
void CrystalFieldFunction::buildSingleSite() const {
  if (isMultiSpectrum()) {
    buildSingleSiteMultiSpectrum();
  } else {
    buildSingleSiteSingleSpectrum();
  }
}

/// Build the target function in a multi site case.
void CrystalFieldFunction::buildMultiSite() const {
  if (isMultiSpectrum()) {
    buildMultiSiteMultiSpectrum();
  } else {
    buildMultiSiteSingleSpectrum();
  }
}

/// Build the target function in a single site - single spectrum case.
void CrystalFieldFunction::buildSingleSiteSingleSpectrum() const {
  auto spectrum = new CompositeFunction;
  m_target.reset(spectrum);
  m_target->setAttributeValue("NumDeriv", true);
  auto bkgdShape = getAttribute("Background").asUnquotedString();
  bool fixAllPeaks = getAttribute("FixAllPeaks").asBool();

  if (!bkgdShape.empty()) {
    auto background =
        API::FunctionFactory::Instance().createInitialized(bkgdShape);
    spectrum->addFunction(background);
  }

  FunctionDomainGeneral domain;
  FunctionValues values;
  m_source->function(domain, values);

  if (values.size() == 0) {
    return;
  }

  if (values.size() % 2 != 0) {
    throw std::runtime_error(
        "CrystalFieldPeaks returned odd number of values.");
  }

  auto xVec = m_control.getAttribute("FWHMX").asVector();
  auto yVec = m_control.getAttribute("FWHMY").asVector();
  auto &FWHMs = m_control.FWHMs();
  auto defaultFWHM = FWHMs.empty() ? 0.0 : FWHMs[0];

  auto fwhmVariation = getAttribute("FWHMVariation").asDouble();
  auto peakShape = getAttribute("PeakShape").asString();
  size_t nRequiredPeaks = getAttribute("NPeaks").asInt();
  CrystalFieldUtils::buildSpectrumFunction(*spectrum, peakShape, values, xVec,
                                           yVec, fwhmVariation, defaultFWHM,
                                           nRequiredPeaks, fixAllPeaks);
}

/// Build the target function in a single site - multi spectrum case.
void CrystalFieldFunction::buildSingleSiteMultiSpectrum() const {
  auto fun = new MultiDomainFunction;
  m_target.reset(fun);

  DoubleFortranVector energies;
  ComplexFortranMatrix waveFunctions;
  ComplexFortranMatrix hamiltonian;
  ComplexFortranMatrix hamiltonianZeeman;
  int nre = 0;
  auto &peakCalculator = dynamic_cast<CrystalFieldPeaksBase &>(*m_source);
  peakCalculator.calculateEigenSystem(energies, waveFunctions, hamiltonian,
                                      hamiltonianZeeman, nre);
  hamiltonian += hamiltonianZeeman;

  const auto nSpec = nSpectra();
  auto &temperatures = m_control.temperatures();
  auto &FWHMs = m_control.FWHMs();
  const bool addBackground = true;
  for (size_t i = 0; i < nSpec; ++i) {
    auto intensityScaling =
        m_control.getFunction(i)->getParameter("IntensityScaling");
    fun->addFunction(buildSpectrum(nre, energies, waveFunctions,
                                   temperatures[i], FWHMs[i], i, addBackground,
                                   intensityScaling));
    fun->setDomainIndex(i, i);
  }
  auto &physProps = m_control.physProps();
  size_t i = nSpec;
  for (auto &prop : physProps) {
    auto physPropFun =
        buildPhysprop(nre, energies, waveFunctions, hamiltonian, prop);
    fun->addFunction(physPropFun);
    fun->setDomainIndex(i, i);
    m_mapPrefixes2PhysProps[prop] = physPropFun;
    ++i;
  }
}

/// Build the target function in a multi site - single spectrum case.
void CrystalFieldFunction::buildMultiSiteSingleSpectrum() const {

  auto spectrum = new CompositeFunction;
  m_target.reset(spectrum);
  m_target->setAttributeValue("NumDeriv", true);
  auto bkgdShape = getAttribute("Background").asUnquotedString();
  bool fixAllPeaks = getAttribute("FixAllPeaks").asBool();

  if (!bkgdShape.empty()) {
    auto background =
        API::FunctionFactory::Instance().createInitialized(bkgdShape);
    spectrum->addFunction(background);
  }

  auto &FWHMs = m_control.FWHMs();
  auto defaultFWHM = FWHMs.empty() ? 0.0 : FWHMs[0];
  auto fwhmVariation = getAttribute("FWHMVariation").asDouble();
  auto peakShape = getAttribute("PeakShape").asString();
  size_t nRequiredPeaks = getAttribute("NPeaks").asInt();
  auto xVec = m_control.getAttribute("FWHMX").asVector();
  auto yVec = m_control.getAttribute("FWHMY").asVector();

  auto &compSource = compositeSource();
  for (size_t ionIndex = 0; ionIndex < compSource.nFunctions(); ++ionIndex) {
    FunctionDomainGeneral domain;
    FunctionValues values;
    compSource.getFunction(ionIndex)->function(domain, values);

    if (values.size() == 0) {
      continue;
    }

    if (values.size() % 2 != 0) {
      throw std::runtime_error(
          "CrystalFieldPeaks returned odd number of values.");
    }

    auto ionSpectrum = boost::make_shared<CompositeFunction>();
    CrystalFieldUtils::buildSpectrumFunction(
        *ionSpectrum, peakShape, values, xVec, yVec, fwhmVariation, defaultFWHM,
        nRequiredPeaks, fixAllPeaks);
    spectrum->addFunction(ionSpectrum);
  }
}

/// Build the target function in a multi site - multi spectrum case.
void CrystalFieldFunction::buildMultiSiteMultiSpectrum() const {
  auto multiDomain = new MultiDomainFunction;
  m_target.reset(multiDomain);

  const auto nSpec = nSpectra();
  std::vector<CompositeFunction *> spectra(nSpec);
  for (size_t i = 0; i < nSpec; ++i) {
    auto spectrum = boost::make_shared<CompositeFunction>();
    spectra[i] = spectrum.get();
    multiDomain->addFunction(spectrum);
    multiDomain->setDomainIndex(i, i);
  }
  auto &physProps = m_control.physProps();
  std::vector<CompositeFunction_sptr> compositePhysProps(physProps.size());
  std::generate(compositePhysProps.begin(), compositePhysProps.end(),
                []() { return boost::make_shared<CompositeFunction>(); });

  auto &compSource = compositeSource();
  for (size_t ionIndex = 0; ionIndex < compSource.nFunctions(); ++ionIndex) {
    DoubleFortranVector energies;
    ComplexFortranMatrix waveFunctions;
    ComplexFortranMatrix hamiltonian;
    ComplexFortranMatrix hamiltonianZeeman;
    int nre = 0;
    auto &peakCalculator = dynamic_cast<CrystalFieldPeaksBase &>(
        *compSource.getFunction(ionIndex));
    peakCalculator.calculateEigenSystem(energies, waveFunctions, hamiltonian,
                                        hamiltonianZeeman, nre);
    hamiltonian += hamiltonianZeeman;

    auto &temperatures = m_control.temperatures();
    auto &FWHMs = m_control.FWHMs();
    const bool addBackground = ionIndex == 0;
    auto ionIntensityScaling =
        compSource.getFunction(ionIndex)->getParameter("IntensityScaling");
    for (size_t i = 0; i < nSpec; ++i) {
      auto spectrumIntensityScaling =
          m_control.getFunction(i)->getParameter("IntensityScaling");
      spectra[i]->addFunction(buildSpectrum(
          nre, energies, waveFunctions, temperatures[i], FWHMs[i], i,
          addBackground, ionIntensityScaling * spectrumIntensityScaling));
    }

    size_t i = 0;
    for (auto &prop : physProps) {
      auto physPropFun =
          buildPhysprop(nre, energies, waveFunctions, hamiltonian, prop);
      compositePhysProps[i]->addFunction(physPropFun);
      std::string propName = "ion";
      propName.append(std::to_string(ionIndex)).append(".").append(prop);
      m_mapPrefixes2PhysProps[propName] = physPropFun;
      ++i;
    }
  }
  m_target->checkFunction();
  size_t i = nSpec;
  for (auto &propFun : compositePhysProps) {
    multiDomain->addFunction(propFun);
    multiDomain->setDomainIndex(i, i);
    ++i;
  }
}

/// Calculate excitations at given temperature.
/// @param nre :: An id of the ion.
/// @param energies :: A vector with energies.
/// @param waveFunctions :: A matrix with wave functions.
/// @param temperature :: A temperature of the spectrum.
/// @param values :: An object to receive computed excitations.
/// @param intensityScaling :: A scaling factor for the intensities.
void CrystalFieldFunction::calcExcitations(
    int nre, const DoubleFortranVector &energies,
    const ComplexFortranMatrix &waveFunctions, double temperature,
    FunctionValues &values, double intensityScaling) const {
  IntFortranVector degeneration;
  DoubleFortranVector eEnergies;
  DoubleFortranMatrix iEnergies;
  const double toleranceEnergy = getAttribute("ToleranceEnergy").asDouble();
  const double toleranceIntensity =
      getAttribute("ToleranceIntensity").asDouble();
  DoubleFortranVector eExcitations;
  DoubleFortranVector iExcitations;
  calculateIntensities(nre, energies, waveFunctions, temperature,
                       toleranceEnergy, degeneration, eEnergies, iEnergies);
  calculateExcitations(eEnergies, iEnergies, toleranceEnergy,
                       toleranceIntensity, eExcitations, iExcitations);
  const auto nPeaks = eExcitations.size();
  values.expand(2 * nPeaks);
  for (size_t i = 0; i < nPeaks; ++i) {
    values.setCalculated(i, eExcitations.get(i));
    values.setCalculated(i + nPeaks, iExcitations.get(i) * intensityScaling);
  }
}

/// Build a function for a single spectrum.
/// @param nre :: An id of the ion.
/// @param energies :: A vector with energies.
/// @param waveFunctions :: A matrix with wave functions.
/// @param temperature :: A temperature of the spectrum.
/// @param fwhm :: A full width at half maximum to set to each peak.
/// @param iSpec :: An index of the created spectrum in m_target composite
/// function.
/// @param addBackground :: An option to add a background to the spectrum.
/// @param intensityScaling :: A scaling factor for the peak intensities.
API::IFunction_sptr CrystalFieldFunction::buildSpectrum(
    int nre, const DoubleFortranVector &energies,
    const ComplexFortranMatrix &waveFunctions, double temperature, double fwhm,
    size_t iSpec, bool addBackground, double intensityScaling) const {
  FunctionValues values;
  calcExcitations(nre, energies, waveFunctions, temperature, values,
                  intensityScaling);
  const auto fwhmVariation = getAttribute("FWHMVariation").asDouble();
  const auto peakShape = getAttribute("PeakShape").asString();
  auto bkgdShape = getAttribute("Background").asUnquotedString();
  const size_t nRequiredPeaks = getAttribute("NPeaks").asInt();
  const bool fixAllPeaks = getAttribute("FixAllPeaks").asBool();

  auto spectrum = new CompositeFunction;

  if (addBackground && !bkgdShape.empty()) {
    if (bkgdShape.find("name=") != 0 && bkgdShape.front() != '(') {
      bkgdShape = "name=" + bkgdShape;
    }
    auto background =
        API::FunctionFactory::Instance().createInitialized(bkgdShape);
    spectrum->addFunction(background);
  }

  auto xVec = m_control.getFunction(iSpec)->getAttribute("FWHMX").asVector();
  auto yVec = m_control.getFunction(iSpec)->getAttribute("FWHMX").asVector();
  CrystalFieldUtils::buildSpectrumFunction(*spectrum, peakShape, values, xVec,
                                           yVec, fwhmVariation, fwhm,
                                           nRequiredPeaks, fixAllPeaks);
  return IFunction_sptr(spectrum);
}

/// Build a physical property function.
/// @param nre :: An id of the ion.
/// @param energies :: A vector with energies.
/// @param waveFunctions :: A matrix with wave functions.
/// @param hamiltonian :: A matrix with the hamiltonian.
/// @param propName :: the name of the physical property.
API::IFunction_sptr
CrystalFieldFunction::buildPhysprop(int nre,
                                    const DoubleFortranVector &energies,
                                    const ComplexFortranMatrix &waveFunctions,
                                    const ComplexFortranMatrix &hamiltonian,
                                    const std::string &propName) const {

  if (propName == "cv") { // HeatCapacity
    auto propFun = boost::make_shared<CrystalFieldHeatCapacityCalculation>();
    propFun->setEnergy(energies);
    return propFun;
  }
  if (propName == "chi") { // Susceptibility
    auto propFun = boost::make_shared<CrystalFieldSusceptibilityCalculation>();
    propFun->setEigensystem(energies, waveFunctions, nre);
    return propFun;
  }
  if (propName == "mh") { // Magnetisation
    auto propFun = boost::make_shared<CrystalFieldMagnetisationCalculation>();
    propFun->setHamiltonian(hamiltonian, nre);
    return propFun;
  }
  if (propName == "mt") { // MagneticMoment
    auto propFun = boost::make_shared<CrystalFieldMomentCalculation>();
    propFun->setHamiltonian(hamiltonian, nre);
    return propFun;
  }

  throw std::runtime_error("Physical property type not understood: " +
                           propName);
}

/// Update a physical property function.
/// @param nre :: An id of the ion.
/// @param energies :: A vector with energies.
/// @param waveFunctions :: A matrix with wave functions.
/// @param hamiltonian :: A matrix with the hamiltonian.
/// @param function :: A function to update.
void CrystalFieldFunction::updatePhysprop(
    int nre, const DoubleFortranVector &energies,
    const ComplexFortranMatrix &waveFunctions,
    const ComplexFortranMatrix &hamiltonian, API::IFunction &function) const {

  auto propName = function.name();

  if (propName == "cv") { // HeatCapacity
    auto &propFun =
        dynamic_cast<CrystalFieldHeatCapacityCalculation &>(function);
    propFun.setEnergy(energies);
  } else if (propName == "chi") { // Susceptibility
    auto &propFun =
        dynamic_cast<CrystalFieldSusceptibilityCalculation &>(function);
    propFun.setEigensystem(energies, waveFunctions, nre);
  } else if (propName == "mh") { // Magnetisation
    auto &propFun =
        dynamic_cast<CrystalFieldMagnetisationCalculation &>(function);
    propFun.setHamiltonian(hamiltonian, nre);
  } else if (propName == "mt") { // MagneticMoment
    auto &propFun = dynamic_cast<CrystalFieldMomentCalculation &>(function);