NumCpp  2.10.1
A Templatized Header Only C++ Implementation of the Python NumPy Library
gaussNewtonNlls.hpp
Go to the documentation of this file.
1 #pragma once
31 
32 #include <array>
33 #include <functional>
34 #include <type_traits>
35 #include <utility>
36 
40 #include "NumCpp/Core/Shape.hpp"
41 #include "NumCpp/Core/Types.hpp"
42 #include "NumCpp/Functions/rms.hpp"
43 #include "NumCpp/Linalg/inv.hpp"
44 #include "NumCpp/NdArray.hpp"
45 
46 namespace nc::linalg
47 {
48  //============================================================================
49  // Method Description:
71  template<typename dtype,
72  typename... Params,
73  std::enable_if_t<std::is_arithmetic_v<dtype>, int> = 0,
74  std::enable_if_t<all_arithmetic_v<Params...>, int> = 0,
75  std::enable_if_t<all_same_v<dtype, Params...>, int> = 0>
76  std::pair<NdArray<double>, double>
77  gaussNewtonNlls(const uint32 numIterations,
78  const NdArray<dtype>& coordinates,
79  const NdArray<dtype>& measurements,
80  const std::function<dtype(const NdArray<dtype>&, const NdArray<dtype>&)>& function,
81  const std::array<std::function<dtype(const NdArray<dtype>&, const NdArray<dtype>&)>,
82  sizeof...(Params)>& derivatives,
83  Params... initialGuess)
84  {
86 
87  const auto coordinatesShape = coordinates.shape();
88 
89  if (coordinatesShape.rows != measurements.size())
90  {
91  THROW_INVALID_ARGUMENT_ERROR("coordinates number of rows, and measurements size must be the same.");
92  }
93 
94  NdArray<double> beta = NdArray<dtype>({ initialGuess... }).template astype<double>().transpose();
95  NdArray<double> residuals(coordinatesShape.rows, 1);
96  NdArray<double> jacobian(coordinatesShape.rows, sizeof...(Params));
97 
98  const auto colSlice = coordinates.cSlice();
99  for (uint32 iteration = 1; iteration <= numIterations; ++iteration)
100  {
101  for (uint32 measIdx = 0; measIdx < coordinatesShape.rows; ++measIdx)
102  {
103  const auto coordinate = coordinates(measIdx, colSlice);
104 
105  residuals[measIdx] =
106  static_cast<double>(measurements[measIdx]) - static_cast<double>(function(coordinate, beta));
107 
108  for (uint32 paramIdx = 0; paramIdx < sizeof...(Params); ++paramIdx)
109  {
110  const auto& derivative = derivatives[paramIdx];
111  jacobian(measIdx, paramIdx) = static_cast<double>(derivative(coordinate, beta));
112  }
113  }
114 
115  // perform the gauss-newton linear algebra
116  const auto jacobianT = jacobian.transpose();
117  const auto jacobianPsuedoInverse = linalg::inv(jacobianT.dot(jacobian));
118  const auto intermediate = jacobianPsuedoInverse.dot(jacobianT);
119  const auto deltaBeta = intermediate.dot(residuals);
120  beta += deltaBeta;
121  }
122 
123  // calculate the final rms of the residuals
124  for (uint32 measIdx = 0; measIdx < coordinatesShape.rows; ++measIdx)
125  {
126  const auto coordinate = coordinates(measIdx, colSlice);
127 
128  residuals[measIdx] =
129  static_cast<double>(measurements[measIdx]) - static_cast<double>(function(coordinate, beta));
130  }
131 
132  return std::make_pair(beta.flatten(), rms(residuals).item());
133  }
134 } // namespace nc::linalg
#define THROW_INVALID_ARGUMENT_ERROR(msg)
Definition: Error.hpp:37
#define STATIC_ASSERT_ARITHMETIC(dtype)
Definition: StaticAsserts.hpp:39
Holds 1D and 2D arrays, the main work horse of the NumCpp library.
Definition: NdArrayCore.hpp:138
size_type size() const noexcept
Definition: NdArrayCore.hpp:4415
self_type transpose() const
Definition: NdArrayCore.hpp:4775
Shape shape() const noexcept
Definition: NdArrayCore.hpp:4402
Slice cSlice(index_type inStartIdx=0, size_type inStepSize=1) const
Definition: NdArrayCore.hpp:951
Definition: cholesky.hpp:41
NdArray< double > inv(const NdArray< dtype > &inArray)
Definition: inv.hpp:54
std::pair< NdArray< double >, double > gaussNewtonNlls(const uint32 numIterations, const NdArray< dtype > &coordinates, const NdArray< dtype > &measurements, const std::function< dtype(const NdArray< dtype > &, const NdArray< dtype > &)> &function, const std::array< std::function< dtype(const NdArray< dtype > &, const NdArray< dtype > &)>, sizeof...(Params)> &derivatives, Params... initialGuess)
Definition: gaussNewtonNlls.hpp:77
dtype beta(GeneratorType &generator, dtype inAlpha, dtype inBeta)
Definition: Random/beta.hpp:61
constexpr bool all_arithmetic_v
Definition: TypeTraits.hpp:67
NdArray< double > rms(const NdArray< dtype > &inArray, Axis inAxis=Axis::NONE)
Definition: rms.hpp:51
constexpr bool all_same_v
Definition: TypeTraits.hpp:101
std::uint32_t uint32
Definition: Types.hpp:40