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