MultiDomainFunction.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 +
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidAPI/MultiDomainFunction.h"
#include "MantidAPI/CompositeDomain.h"
#include "MantidAPI/Expression.h"
#include "MantidAPI/FunctionFactory.h"
 
#include <algorithm>
#include <boost/lexical_cast.hpp>
#include <set>
 
namespace Mantid::API {
 
DECLARE_FUNCTION(MultiDomainFunction)
 
MultiDomainFunction::MultiDomainFunction() : m_nDomains(0), m_maxIndex(0) { setAttributeValue("NumDeriv", true); }
 
void MultiDomainFunction::setDomainIndex(size_t funIndex, size_t domainIndex) {
  m_domains[funIndex] = std::vector<size_t>(1, domainIndex);
  countNumberOfDomains();
}
 
void MultiDomainFunction::setDomainIndices(size_t funIndex, const std::vector<size_t> &domainIndices) {
  m_domains[funIndex] = domainIndices;
  countNumberOfDomains();
}
 
void MultiDomainFunction::countNumberOfDomains() {
  std::set<size_t> dSet;
  for (auto &domain : m_domains) {
    if (!domain.second.empty()) {
      dSet.insert(domain.second.begin(), domain.second.end());
    }
  }
  m_nDomains = dSet.size();
  m_maxIndex = dSet.empty() ? 0 : *dSet.rbegin();
}
 
void MultiDomainFunction::countValueOffsets(const CompositeDomain &domain) const {
  m_valueOffsets.clear();
  m_valueOffsets.emplace_back(0);
  for (size_t i = 0; i < domain.getNParts(); ++i) {
    const FunctionDomain &d = domain.getDomain(i);
    m_valueOffsets.emplace_back(m_valueOffsets.back() + d.size());
  }
}
 
void MultiDomainFunction::clearDomainIndices() {
  m_domains.clear();
  countNumberOfDomains();
}
 
void MultiDomainFunction::getDomainIndices(size_t funIndex, size_t nDomains, std::vector<size_t> &domains) const {
  auto it = m_domains.find(funIndex);
  if (it == m_domains.end()) { // apply to all domains
    domains.resize(nDomains);
    for (size_t i = 0; i < domains.size(); ++i) {
      domains[i] = i;
    }
  } else { // apply to selected domains
    domains.assign(it->second.begin(), it->second.end());
  }
}
 
size_t MultiDomainFunction::getNumberDomains() const { return getMaxIndex() + 1; }
 
void MultiDomainFunction::function(const FunctionDomain &domain, FunctionValues &values) const {
  // works only on CompositeDomain
  if (!dynamic_cast<const CompositeDomain *>(&domain)) {
    // make exception if a single member function is defined
    if (m_maxIndex == 0 && nFunctions() == 1) {
      getFunction(0)->function(domain, values);
      return;
    }
    throw std::invalid_argument("Non-CompositeDomain passed to MultiDomainFunction.");
  }
  const auto &cd = dynamic_cast<const CompositeDomain &>(domain);
  // domain must not have less parts than m_maxIndex
  if (cd.getNParts() <= m_maxIndex) {
    throw std::invalid_argument("CompositeDomain has too few parts (" + std::to_string(cd.getNParts()) +
                                ") for MultiDomainFunction (max index " + std::to_string(m_maxIndex) + ").");
  }
  // domain and values must be consistent
  if (cd.size() != values.size()) {
    throw std::invalid_argument("MultiDomainFunction: domain and values have different sizes.");
  }
 
  countValueOffsets(cd);
  // evaluate member functions
  values.zeroCalculated();
  for (size_t iFun = 0; iFun < nFunctions(); ++iFun) {
    // find the domains member function must be applied to
    std::vector<size_t> domains;
    getDomainIndices(iFun, cd.getNParts(), domains);
 
    for (auto const &dom : domains) {
      const FunctionDomain &d = cd.getDomain(dom);
      FunctionValues tmp(d);
      getFunction(iFun)->function(d, tmp);
      values.addToCalculated(m_valueOffsets[dom], tmp);
    }
  }
}
 
void MultiDomainFunction::functionDeriv(const FunctionDomain &domain, Jacobian &jacobian) {
  // works only on CompositeDomain
  if (!dynamic_cast<const CompositeDomain *>(&domain)) {
    // make exception if a single member function is defined
    if (m_maxIndex == 0 && nFunctions() == 1) {
      getFunction(0)->functionDeriv(domain, jacobian);
      return;
    }
    throw std::invalid_argument("Non-CompositeDomain passed to MultiDomainFunction.");
  }
 
  if (getAttribute("NumDeriv").asBool()) {
    calNumericalDeriv(domain, jacobian);
  } else {
    const auto &cd = dynamic_cast<const CompositeDomain &>(domain);
    // domain must not have less parts than m_maxIndex
    if (cd.getNParts() < m_maxIndex) {
      throw std::invalid_argument("CompositeDomain has too few parts (" + std::to_string(cd.getNParts()) +
                                  ") for MultiDomainFunction (max index " + std::to_string(m_maxIndex) + ").");
    }
 
    jacobian.zero();
    countValueOffsets(cd);
    // evaluate member functions derivatives
    for (size_t iFun = 0; iFun < nFunctions(); ++iFun) {
      // find the domains member function must be applied to
      std::vector<size_t> domains;
      getDomainIndices(iFun, cd.getNParts(), domains);
 
      for (auto const &dom : domains) {
        const FunctionDomain &d = cd.getDomain(dom);
        PartialJacobian J(&jacobian, m_valueOffsets[dom], paramOffset(iFun));
        getFunction(iFun)->functionDeriv(d, J);
      }
    }
  }
}
 
void MultiDomainFunction::iterationStarting() {
  for (size_t iFun = 0; iFun < nFunctions(); ++iFun) {
    getFunction(iFun)->iterationStarting();
  }
}
 
void MultiDomainFunction::iterationFinished() {
  for (size_t iFun = 0; iFun < nFunctions(); ++iFun) {
    getFunction(iFun)->iterationFinished();
  }
}
 
IFunction::Attribute MultiDomainFunction::getLocalAttribute(size_t funIndex, const std::string &attName) const {
  if (attName != "domains") {
    throw std::invalid_argument("MultiDomainFunction does not have attribute " + attName);
  }
  if (funIndex >= nFunctions()) {
    throw std::out_of_range("Function index is out of range.");
  }
  try {
    auto it = m_domains.at(funIndex);
    if (it.size() == 1 && it.front() == funIndex) {
      return IFunction::Attribute("i");
    } else if (!it.empty()) {
      auto out = std::to_string(it.front());
      for (auto i = it.begin() + 1; i != it.end(); ++i) {
        out += "," + std::to_string(*i);
      }
      return IFunction::Attribute(out);
    }
  } catch (const std::out_of_range &) {
    return IFunction::Attribute("All");
  }
  return IFunction::Attribute("");
}
 
void MultiDomainFunction::setLocalAttribute(size_t funIndex, const std::string &attName,
                                            const IFunction::Attribute &att) {
  if (attName != "domains") {
    throw std::invalid_argument("MultiDomainFunction does not have attribute " + attName);
  }
  if (funIndex >= nFunctions()) {
    throw std::out_of_range("Function index is out of range.");
  }
  std::string value = att.asString();
  auto it = m_domains.find(funIndex);
 
  if (value == "All") { // fit to all domains
    if (it != m_domains.end()) {
      m_domains.erase(it);
    }
    return;
  } else if (value == "i") { // fit to domain with the same index as the function
    setDomainIndex(funIndex, funIndex);
    return;
  } else if (value.empty()) { // do not fit to any domain
    setDomainIndices(funIndex, std::vector<size_t>());
    return;
  }
 
  // fit to a selection of domains
  std::vector<size_t> indx;
  Expression list;
  list.parse(value);
  if (list.name() == "+") {
    if (list.size() != 2 || list.terms()[1].operator_name() != "-") {
      throw std::runtime_error("MultiDomainFunction: attribute \"domains\" "
                               "expects two integers separated by a \"-\"");
    }
    // value looks like "a - b". a and b must be ints and define a range of
    // domain indices
    auto start = boost::lexical_cast<size_t>(list.terms()[0].str());
    size_t end = boost::lexical_cast<size_t>(list.terms()[1].str()) + 1;
    if (start >= end) {
      throw std::runtime_error("MultiDomainFunction: attribute \"domains\": wrong range limits.");
    }
    indx.resize(end - start);
    for (size_t i = start; i < end; ++i) {
      indx[i - start] = i;
    }
  } else {
    // value must be either an int or a list of ints: "a,b,c,..."
    list.toList();
    indx.reserve(list.size());
    std::transform(list.begin(), list.end(), std::back_inserter(indx),
                   [](const Mantid::API::Expression &k) { return boost::lexical_cast<size_t>(k.name()); });
  }
  setDomainIndices(funIndex, indx);
}
 
std::vector<IFunction_sptr> MultiDomainFunction::createEquivalentFunctions() const {
  size_t nDomains = m_maxIndex + 1;
  std::vector<CompositeFunction_sptr> compositeFunctions(nDomains);
  for (size_t iFun = 0; iFun < nFunctions(); ++iFun) {
    // find the domains member function must be applied to
    std::vector<size_t> domains;
    getDomainIndices(iFun, nDomains, domains);
 
    for (auto domainIndex : domains) {
      CompositeFunction_sptr cf = compositeFunctions[domainIndex];
      if (!cf) {
        // create a composite function for each domain
        cf = CompositeFunction_sptr(new CompositeFunction());
        compositeFunctions[domainIndex] = cf;
      }
      // add copies of all functions applied to j-th domain to a single
      // compositefunction
      cf->addFunction(FunctionFactory::Instance().createInitialized(getFunction(iFun)->asString()));
    }
  }
  std::vector<IFunction_sptr> outFunctions(nDomains);
  // fill in the output vector
  // check functions containing a single member and take it out of the composite
  for (size_t i = 0; i < compositeFunctions.size(); ++i) {
    auto fun = compositeFunctions[i];
    if (!fun || fun->nFunctions() == 0) {
      throw std::runtime_error("There is no function for domain " + std::to_string(i));
    }
    if (fun->nFunctions() > 1) {
      outFunctions[i] = fun;
    } else {
      outFunctions[i] = fun->getFunction(0);
    }
  }
  return outFunctions;
}
 
} // namespace Mantid::API