EigenJacobian.h

Go to the documentation of this file.

// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2010 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 +
#pragma once
 
#include "MantidAPI/IFunction.h"
#include "MantidAPI/Jacobian.h"
#include "MantidCurveFitting/EigenMatrix.h"
 
#include <stdexcept>
#include <vector>
 
namespace Mantid {
namespace CurveFitting {
class EigenJacobian : public API::Jacobian {
  EigenMatrix m_J;
  std::vector<int> m_index;
 
public:
  EigenJacobian(const API::IFunction &fun, const size_t ny) {
    m_index.resize(fun.nParams(), -1);
    size_t np = 0; // number of active parameters
    for (size_t i = 0; i < fun.nParams(); ++i) {
      m_index[i] = static_cast<int>(np);
      if (fun.isActive(i)) {
        ++np;
      }
    }
    m_J.resize(ny, np);
  }
 
  EigenMatrix &matrix() { return m_J; }
 
  map_type &getJ() { return m_J.mutator(); }
 
  void addNumberToColumn(const double &value, const size_t &iActiveP) override {
    if (iActiveP < m_J.size2()) {
      // add penalty to first and last point and every 10th point in between
      m_J.mutator().data()[iActiveP * m_J.size1()] += value;
      m_J.mutator().data()[(iActiveP + 1) * m_J.size1() - 1] += value;
      for (size_t iY = 9; iY < m_J.size1() - 1; iY += 10)
        m_J.mutator().data()[iActiveP * m_J.size1() + iY] += value;
    } else {
      throw std::runtime_error("Try to add number to column of Jacobian matrix "
                               "which does not exist.");
    }
  }
  void set(size_t iY, size_t iP, double value) override {
    int j = m_index[iP];
    if (j >= 0)
      m_J.set(iY, j, value);
  }
  double get(size_t iY, size_t iP) override {
    int j = m_index[iP];
    if (j >= 0)
      return m_J.get(iY, j);
    return 0.0;
  }
  void zero() override { m_J.zero(); }
};
 
template <class T> class JacobianImpl1 : public API::Jacobian {
public:
  T *m_J;
  std::vector<int> m_index;
 
  void setJ(T *J) { m_J = J; }
 
  void addNumberToColumn(const double &value, const size_t &iActiveP) override {
    if (iActiveP < m_J->size2()) {
      // add penalty to first and last point and every 10th point in between
      m_J->mutator().data()[iActiveP * m_J->size1()] += value;
      m_J->mutator().data()[(iActiveP + 1) * m_J->size1() - 1] += value;
      for (size_t iY = 9; iY < m_J->size1() - 1; iY += 10)
        m_J->mutator().data()[iActiveP * m_J->size1() + iY] += value;
    } else {
      throw std::runtime_error("Try to add number to column of Jacobian matrix "
                               "which does not exist.");
    }
  }
  void set(size_t iY, size_t iP, double value) override {
    if (iP >= m_index.size()) {
      // Functions are allowed to change their number of active
      // fitting parameters during the fit (example: the crystal field
      // functions). When it happens after an iteration is finished there is a
      // special exception that signals the minimizer to re-initialize itself.
      // But it can happen during a numeric derivative calculation: a field
      // parameter is incremented, a new peak may appear and it may have a free
      // fitting parameter, so the number of parameters change. You cannot throw
      // and re-initialize the minimizer at this point because the same will
      // happen again. If you ignore these "virtual" extra parameters by the
      // time the derivative calculations finish the number of parameters will
      // return to the original value. Note that the get(...) method doesn't
      // skip any extra indices because if there are any it is an error.
      return;
    }
    int j = m_index[iP];
    if (j >= 0)
      m_J->mutator().data()[j * m_J->size1() + iY] = value;
  }
  double get(size_t iY, size_t iP) override {
    int j = m_index[iP];
    if (j >= 0)
      return m_J->inspector().data()[j * m_J->size1() + iY];
    return 0.0;
  }
  void zero() override { m_J->zero(); }
};
 
} // namespace CurveFitting
} // namespace Mantid