ParameterSet.hpp

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file ParameterSet.hpp

// This file is part of the Acts project.
//
// Copyright (C) 2016-2019 CERN for the benefit of the Acts project
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at http://mozilla.org/MPL/2.0/.

#pragma once
// STL include(s)
#include <memory>
#include <optional>
#include <type_traits>
#include <utility>

// Acts includes
#include "Acts/EventData/detail/full_parameter_set.hpp"
#include "Acts/EventData/detail/initialize_parameter_set.hpp"
#include "Acts/EventData/detail/make_projection_matrix.hpp"
#include "Acts/EventData/detail/residual_calculator.hpp"
#include "Acts/EventData/detail/value_corrector.hpp"
#include "Acts/Utilities/Definitions.hpp"
#include "Acts/Utilities/ParameterDefinitions.hpp"
#include "Acts/Utilities/detail/MPL/are_sorted.hpp"
#include "Acts/Utilities/detail/MPL/are_within.hpp"
#include "Acts/Utilities/detail/MPL/at_index.hpp"
#include "Acts/Utilities/detail/MPL/get_position.hpp"
#include "Acts/Utilities/detail/MPL/is_contained.hpp"

namespace Acts {
// forward type declaration for full parameter set
using FullParameterSet = typename detail::full_parset::type;

template <ParID_t... params>
class ParameterSet {
 private:
  // local typedefs and constants
  using ParSet_t = ParameterSet<params...>;  
  static constexpr unsigned int NPars =
      sizeof...(params);  

  // static assert to check that the template parameters are consistent
  static_assert(detail::are_sorted<true, true, ParID_t, params...>::value,
                "parameter identifiers are not sorted");
  static_assert(
      detail::are_within<unsigned int, 0, eBoundParametersSize,
                         static_cast<unsigned int>(params)...>::value,
      "parameter identifiers must be greater or "
      "equal to zero and smaller than the total number of parameters");
  static_assert(NPars > 0, "number of stored parameters can not be zero");
  static_assert(
      NPars <= eBoundParametersSize,
      "number of stored parameters can not exceed number of total parameters");

 public:
  // public typedefs
  using Projection_t = ActsMatrix<ParValue_t, NPars, eBoundParametersSize>;
  using ParVector_t = ActsVector<ParValue_t, NPars>;
  using CovMatrix_t = ActsSymMatrix<ParValue_t, NPars>;

  template <typename... Tail>
  ParameterSet(std::optional<CovMatrix_t> cov,
               std::enable_if_t<sizeof...(Tail) + 1 == NPars, ParValue_t> head,
               Tail... values)
      : m_vValues(NPars) {
    if (cov) {
      m_optCovariance = std::move(*cov);
    }
    detail::initialize_parset<ParID_t, params...>::init(*this, head, values...);
  }

  ParameterSet(std::optional<CovMatrix_t> cov, const ParVector_t& values)
      : m_vValues(NPars) {
    if (cov) {
      m_optCovariance = std::move(*cov);
    }
    detail::initialize_parset<ParID_t, params...>::init(*this, values);
  }

  ParameterSet(const ParSet_t& copy)
      : m_vValues(copy.m_vValues), m_optCovariance(copy.m_optCovariance) {}

  ParameterSet(ParSet_t&& copy) : m_vValues(std::move(copy.m_vValues)) {
    if (copy.m_optCovariance) {
      m_optCovariance = std::move(*copy.m_optCovariance);
    }
  }

  ~ParameterSet() = default;

  ParSet_t& operator=(const ParSet_t& rhs) {
    m_vValues = rhs.m_vValues;
    m_optCovariance = rhs.m_optCovariance;
    return *this;
  }

  ParSet_t& operator=(ParSet_t&& rhs) {
    m_vValues = std::move(rhs.m_vValues);
    m_optCovariance = std::move(rhs.m_optCovariance);
    return *this;
  }

  friend void swap(ParSet_t& first, ParSet_t& second) noexcept {
    using std::swap;
    swap(first.m_vValues, second.m_vValues);
    swap(first.m_optCovariance, second.m_optCovariance);
  }

  template <ParID_t parameter>
  static constexpr size_t getIndex() {
    return detail::get_position<ParID_t, parameter, params...>::value;
  }

  template <size_t index>
  static constexpr ParID_t getParID() {
    return detail::at_index<ParID_t, index, params...>::value;
  }

  template <ParID_t parameter>
  ParValue_t getParameter() const {
    return m_vValues(getIndex<parameter>());
  }

  ParVector_t getParameters() const { return m_vValues; }

  template <ParID_t parameter>
  void setParameter(ParValue_t value) {
    m_vValues(getIndex<parameter>()) =
        BoundParameterType<parameter>::getValue(value);
  }

  void setParameters(const ParVector_t& values) {
    detail::initialize_parset<ParID_t, params...>::init(*this, values);
  }

  template <ParID_t parameter>
  bool contains() const {
    return detail::is_contained<ParID_t, parameter, params...>::value;
  }

  const std::optional<CovMatrix_t>& getCovariance() const {
    return m_optCovariance;
  }

  template <ParID_t parameter>
  ParValue_t getUncertainty() const {
    if (m_optCovariance) {
      size_t index = getIndex<parameter>();
      return sqrt((*m_optCovariance)(index, index));
    } else {
      return -1;
    }
  }

  void setCovariance(const CovMatrix_t& cov) { m_optCovariance = cov; }

  bool operator==(const ParSet_t& rhs) const {
    // shortcut comparison with myself
    if (&rhs == this) {
      return true;
    }

    // parameter values
    if (m_vValues != rhs.m_vValues) {
      return false;
    }
    // both have covariance matrices set
    if ((m_optCovariance && rhs.m_optCovariance) &&
        (*m_optCovariance != *rhs.m_optCovariance)) {
      return false;
    }
    // only one has a covariance matrix set
    if ((m_optCovariance && !rhs.m_optCovariance) ||
        (!m_optCovariance && rhs.m_optCovariance)) {
      return false;
    }

    return true;
  }

  bool operator!=(const ParSet_t& rhs) const { return !(*this == rhs); }

  ParVector_t project(const FullParameterSet& fullParSet) const {
    return projector() * fullParSet.getParameters();
  }


  template <
      typename T = ParSet_t,
      std::enable_if_t<not std::is_same<T, FullParameterSet>::value, int> = 0>
  ParVector_t residual(const FullParameterSet& fullParSet) const {
    return detail::residual_calculator<params...>::result(
        m_vValues, projector() * fullParSet.getParameters());
  }

  ParVector_t residual(const ParSet_t& otherParSet) const {
    return detail::residual_calculator<params...>::result(
        m_vValues, otherParSet.m_vValues);
  }

  static const ActsMatrix<ParValue_t, NPars, eBoundParametersSize> projector() {
    return sProjector;
  }

  static constexpr unsigned int size() { return NPars; }

  static void correctValues(ParVector_t& values) {
    detail::value_corrector<params...>::result(values);
  }

 private:
  ParVector_t m_vValues{ParVector_t::Zero()};  

  std::optional<CovMatrix_t> m_optCovariance{
      std::nullopt};  

  static const Projection_t sProjector;  

};

// initialize static class members
template <ParID_t... params>
constexpr unsigned int ParameterSet<params...>::NPars;

template <ParID_t... params>
const typename ParameterSet<params...>::Projection_t
    ParameterSet<params...>::sProjector = detail::make_projection_matrix<
        eBoundParametersSize, static_cast<unsigned int>(params)...>::init();
}  // namespace Acts