NumCpp  2.7.0
A Templatized Header Only C++ Implementation of the Python NumPy Library
Poly1d.hpp
Go to the documentation of this file.
1 #pragma once
29 
34 #include "NumCpp/Core/Types.hpp"
36 #include "NumCpp/Linalg/inv.hpp"
37 #include "NumCpp/NdArray.hpp"
39 #include "NumCpp/Utils/num2str.hpp"
40 #include "NumCpp/Utils/power.hpp"
41 
42 #include <iostream>
43 #include <string>
44 #include <type_traits>
45 #include <utility>
46 #include <vector>
47 
48 namespace nc
49 {
50  namespace polynomial
51  {
52  //================================================================================
56  template<typename dtype>
57  class Poly1d
58  {
59  private:
60  STATIC_ASSERT_ARITHMETIC(dtype);
61 
62  public:
63  //============================================================================
64  // Method Description:
68  Poly1d() = default;
69 
70  //============================================================================
71  // Method Description:
78  Poly1d(const NdArray<dtype>& inValues, bool isRoots = false)
79  {
80  if (inValues.size() > DtypeInfo<uint8>::max())
81  {
82  THROW_INVALID_ARGUMENT_ERROR("can only make a polynomial of order " + utils::num2str(DtypeInfo<uint8>::max()));
83  }
84 
85  if (isRoots)
86  {
87  coefficients_.push_back(1);
88  for (auto value : inValues)
89  {
90  NdArray<dtype> coeffs = { -(value), static_cast<dtype>(1) };
91  *this *= Poly1d<dtype>(coeffs, !isRoots);
92  }
93  }
94  else
95  {
96  for (auto value : inValues)
97  {
98  coefficients_.push_back(value);
99  }
100  }
101  }
102 
103  //============================================================================
104  // Method Description:
111  double area(double a, double b) const
112  {
113  if (a > b)
114  {
115  std::swap(a, b);
116  }
117 
118  auto polyIntegral = integ();
119  return polyIntegral(b) - polyIntegral(a);
120  }
121 
122  //============================================================================
123  // Method Description:
128  template<typename dtypeOut>
130  {
131  auto newCoefficients = NdArray<dtypeOut>(1, static_cast<uint32>(coefficients_.size()));
132 
133  const auto function = [](dtype value) -> dtypeOut
134  {
135  return static_cast<dtypeOut>(value);
136  };
137 
138  stl_algorithms::transform(coefficients_.begin(), coefficients_.end(), newCoefficients.begin(), function);
139 
140  return Poly1d<dtypeOut>(newCoefficients);
141  }
142 
143  //============================================================================
144  // Method Description:
150  {
151  auto coefficientsCopy = coefficients_;
152  return NdArray<dtype>(coefficientsCopy);
153  }
154 
155  //============================================================================
156  // Method Description:
161  {
162  const auto numCoefficients = static_cast<uint32>(coefficients_.size());
163  if (numCoefficients == 0)
164  {
165  return {};
166  }
167  if (numCoefficients == 1)
168  {
169  return Poly1d<dtype>({ 0 });
170  }
171 
172  NdArray<dtype> derivativeCofficients(1, numCoefficients - 1);
173 
174  uint32 counter = 0;
175  for (uint32 i = 1; i < numCoefficients; ++i)
176  {
177  derivativeCofficients[counter++] = coefficients_[i] * i;
178  }
179 
180  return Poly1d<dtype>(derivativeCofficients);
181  }
182 
183  //============================================================================
184  // Method Description:
191  static Poly1d<double> fit(const NdArray<dtype>& xValues, const NdArray<dtype>& yValues, uint8 polyOrder)
192  {
193  const auto numMeasurements = xValues.size();
194 
195  if (yValues.size() != numMeasurements)
196  {
197  THROW_INVALID_ARGUMENT_ERROR("Input x and y arrays must be of equal size.");
198  }
199 
200  if (!xValues.isflat())
201  {
202  THROW_INVALID_ARGUMENT_ERROR("Input x must be a flattened [1, n] or [n, 1] array.");
203  }
204 
205  if (!yValues.isflat())
206  {
207  THROW_INVALID_ARGUMENT_ERROR("Input y must be a flattened [n, 1] array.");
208  }
209 
210  NdArray<double> a(numMeasurements, polyOrder + 1);
211  for (uint32 measIdx = 0; measIdx < numMeasurements; ++measIdx)
212  {
213  const auto xDouble = static_cast<double>(xValues[measIdx]);
214  for (uint8 order = 0; order < a.numCols(); ++order)
215  {
216  a(measIdx, order) = utils::power(xDouble, order);
217  }
218  }
219 
220  NdArray<double> aInv;
221  if (a.issquare())
222  {
223  aInv = linalg::inv(a);
224  }
225  else
226  {
227  // psuedo-inverse
228  auto aT = a.transpose();
229  auto aTaInv = linalg::inv(aT.dot(a));
230  aInv = aTaInv.dot(aT);
231  }
232 
233  auto x = aInv.dot(yValues.template astype<double>());
234  return Poly1d<double>(x);
235  }
236 
237  //============================================================================
238  // Method Description:
246  static Poly1d<double> fit(const NdArray<dtype>& xValues, const NdArray<dtype>& yValues,
247  const NdArray<dtype>& weights, uint8 polyOrder)
248  {
249  const auto numMeasurements = xValues.size();
250 
251  if (yValues.size() != numMeasurements)
252  {
253  THROW_INVALID_ARGUMENT_ERROR("Input x and y arrays must be of equal size.");
254  }
255 
256  if (weights.size() != numMeasurements)
257  {
258  THROW_INVALID_ARGUMENT_ERROR("Input x and weights arrays must be of equal size.");
259  }
260 
261  if (!xValues.isflat())
262  {
263  THROW_INVALID_ARGUMENT_ERROR("Input x must be a flattened [1, n] or [n, 1] array.");
264  }
265 
266  if (!yValues.isflat())
267  {
268  THROW_INVALID_ARGUMENT_ERROR("Input y must be a flattened [n, 1] array.");
269  }
270 
271  if (!weights.isflat())
272  {
273  THROW_INVALID_ARGUMENT_ERROR("Input weights must be a flattened [1, n] or [n, 1] array.");
274  }
275 
276  NdArray<double> a(numMeasurements, polyOrder + 1);
277  for (uint32 measIdx = 0; measIdx < numMeasurements; ++measIdx)
278  {
279  const auto xDouble = static_cast<double>(xValues[measIdx]);
280  for (uint8 order = 0; order < a.numCols(); ++order)
281  {
282  a(measIdx, order) = utils::power(xDouble, order);
283  }
284  }
285 
286  NdArray<double> aWeighted(a.shape());
287  NdArray<double> yWeighted(yValues.shape());
288 
289  for (uint32 measIdx = 0; measIdx < numMeasurements; ++measIdx)
290  {
291  const auto weight = static_cast<double>(weights[measIdx]);
292 
293  yWeighted[measIdx] = yValues[measIdx] * weight;
294  for (uint8 order = 0; order < a.numCols(); ++order)
295  {
296  aWeighted(measIdx, order) = a(measIdx, order) * weight;
297  }
298  }
299 
300  NdArray<double> aInv;
301  if (aWeighted.issquare())
302  {
303  aInv = linalg::inv(aWeighted);
304  }
305  else
306  {
307  // psuedo-inverse
308  auto aT = a.transpose();
309  auto aTaInv = linalg::inv(aT.dot(aWeighted));
310  aInv = aTaInv.dot(aT);
311  }
312 
313  auto x = aInv.dot(yWeighted);
314  return Poly1d<double>(x);
315  }
316 
317  //============================================================================
318  // Method Description:
323  {
324  const auto numCoefficients = static_cast<uint32>(coefficients_.size());
325  if (numCoefficients == 0)
326  {
327  return {};
328  }
329 
330  NdArray<double> integralCofficients(1, numCoefficients + 1);
331  integralCofficients[0] = 0.0;
332 
333  for (uint32 i = 0; i < numCoefficients; ++i)
334  {
335  integralCofficients[i + 1] = static_cast<double>(coefficients_[i]) / static_cast<double>(i + 1);
336  }
337 
338  return Poly1d<double>(integralCofficients);
339  }
340 
341  //============================================================================
342  // Method Description:
347  uint32 order() const noexcept
348  {
349  return static_cast<uint32>(coefficients_.size() - 1);
350  }
351 
352  //============================================================================
353  // Method Description:
357  void print() const
358  {
359  std::cout << *this << std::endl;
360  }
361 
362  //============================================================================
363  // Method Description:
368  std::string str() const
369  {
370  const auto numCoeffients = static_cast<uint32>(coefficients_.size());
371 
372  std::string repr = "Poly1d<";
373  uint32 power = 0;
374  for (auto& coefficient : coefficients_)
375  {
377  {
378  if (coefficient == 0)
379  {
380  ++power;
381  continue;
382  }
383  }
384  else
385  {
386  if (utils::essentiallyEqual(coefficient, static_cast<dtype>(0.0)))
387  {
388  ++power;
389  continue;
390  }
391  }
392 
393  repr += utils::num2str(coefficient);
394 
395  if (power > 1)
396  {
397  repr += "x^" + utils::num2str(power);
398  }
399  else if (power == 1)
400  {
401  repr += "x";
402  }
403 
404  ++power;
405 
406  if (power < numCoeffients)
407  {
408  repr += " + ";
409  }
410  }
411 
412  return repr + ">";
413  }
414 
415  //============================================================================
416  // Method Description:
422  dtype operator()(dtype inValue) const noexcept
423  {
424  dtype polyValue = 0;
425  uint8 power = 0;
426  for (auto& coefficient : coefficients_)
427  {
428  polyValue += coefficient * utils::power(inValue, power++);
429  }
430 
431  return polyValue;
432  }
433 
434  //============================================================================
435  // Method Description:
441  Poly1d<dtype> operator+(const Poly1d<dtype>& inOtherPoly) const
442  {
443  return Poly1d<dtype>(*this) += inOtherPoly;
444  }
445 
446  //============================================================================
447  // Method Description:
453  Poly1d<dtype>& operator+=(const Poly1d<dtype>& inOtherPoly)
454  {
455  if (this->coefficients_.size() < inOtherPoly.coefficients_.size())
456  {
457  for (size_t i = 0; i < coefficients_.size(); ++i)
458  {
459  coefficients_[i] += inOtherPoly.coefficients_[i];
460  }
461  for (size_t i = coefficients_.size(); i < inOtherPoly.coefficients_.size(); ++i)
462  {
463  coefficients_.push_back(inOtherPoly.coefficients_[i]);
464  }
465  }
466  else
467  {
468  for (size_t i = 0; i < inOtherPoly.coefficients_.size(); ++i)
469  {
470  coefficients_[i] += inOtherPoly.coefficients_[i];
471  }
472  }
473 
474  return *this;
475  }
476 
477  //============================================================================
478  // Method Description:
484  Poly1d<dtype> operator-(const Poly1d<dtype>& inOtherPoly) const
485  {
486  return Poly1d<dtype>(*this) -= inOtherPoly;
487  }
488 
489  //============================================================================
490  // Method Description:
496  Poly1d<dtype>& operator-=(const Poly1d<dtype>& inOtherPoly)
497  {
498  if (this->coefficients_.size() < inOtherPoly.coefficients_.size())
499  {
500  for (size_t i = 0; i < coefficients_.size(); ++i)
501  {
502  coefficients_[i] -= inOtherPoly.coefficients_[i];
503  }
504  for (size_t i = coefficients_.size(); i < inOtherPoly.coefficients_.size(); ++i)
505  {
506  coefficients_.push_back(-inOtherPoly.coefficients_[i]);
507  }
508  }
509  else
510  {
511  for (size_t i = 0; i < inOtherPoly.coefficients_.size(); ++i)
512  {
513  coefficients_[i] -= inOtherPoly.coefficients_[i];
514  }
515  }
516 
517  return *this;
518  }
519 
520  //============================================================================
521  // Method Description:
527  Poly1d<dtype> operator*(const Poly1d<dtype>& inOtherPoly) const
528  {
529  return Poly1d<dtype>(*this) *= inOtherPoly;
530  }
531 
532  //============================================================================
533  // Method Description:
539  Poly1d<dtype>& operator*=(const Poly1d<dtype>& inOtherPoly)
540  {
541  const uint32 finalCoefficientsSize = order() + inOtherPoly.order() + 1;
542  std::vector<dtype> coeffsA(finalCoefficientsSize, 0);
543  std::vector<dtype> coeffsB(finalCoefficientsSize, 0);
544  stl_algorithms::copy(coefficients_.begin(), coefficients_.end(), coeffsA.begin());
545  stl_algorithms::copy(inOtherPoly.coefficients_.cbegin(), inOtherPoly.coefficients_.cend(), coeffsB.begin());
546 
547  // now multiply out the coefficients
548  std::vector<dtype> finalCoefficients(finalCoefficientsSize, 0);
549  for (uint32 i = 0; i < finalCoefficientsSize; ++i)
550  {
551  for (uint32 k = 0; k <= i; ++k)
552  {
553  finalCoefficients[i] += coeffsA[k] * coeffsB[i - k];
554  }
555  }
556 
557  this->coefficients_ = finalCoefficients;
558  return *this;
559  }
560 
561  //============================================================================
562  // Method Description:
569  {
570  return Poly1d(*this) ^= inPower;
571  }
572 
573  //============================================================================
574  // Method Description:
581  {
582  if (inPower == 0)
583  {
584  coefficients_.clear();
585  coefficients_.push_back(1);
586  return *this;
587  }
588  if (inPower == 1)
589  {
590  return *this;
591  }
592 
593  auto thisPoly(*this);
594  for (uint32 power = 1; power < inPower; ++power)
595  {
596  *this *= thisPoly;
597  }
598 
599  return *this;
600  }
601 
602  //============================================================================
603  // Method Description:
610  friend std::ostream& operator<<(std::ostream& inOStream, const Poly1d<dtype>& inPoly)
611  {
612  inOStream << inPoly.str() << std::endl;
613  return inOStream;
614  }
615 
616  private:
617  std::vector<dtype> coefficients_{};
618  };
619  } // namespace polynomial
620 } // namespace nc
#define THROW_INVALID_ARGUMENT_ERROR(msg)
Definition: Error.hpp:36
Holds info about the dtype.
Definition: DtypeInfo.hpp:41
Holds 1D and 2D arrays, the main work horse of the NumCpp library.
Definition: NdArrayCore.hpp:72
uint32 numCols() const noexcept
Definition: NdArrayCore.hpp:3418
size_type size() const noexcept
Definition: NdArrayCore.hpp:4296
Shape shape() const noexcept
Definition: NdArrayCore.hpp:4283
bool issquare() const noexcept
Definition: NdArrayCore.hpp:2918
bool isflat() const noexcept
Definition: NdArrayCore.hpp:2855
NdArray< dtype > transpose() const
Definition: NdArrayCore.hpp:4629
NdArray< dtype > dot(const NdArray< dtype > &inOtherArray) const
Definition: NdArrayCore.hpp:2633
Definition: Poly1d.hpp:58
Poly1d< dtype > & operator*=(const Poly1d< dtype > &inOtherPoly)
Definition: Poly1d.hpp:539
Poly1d< dtype > deriv() const
Definition: Poly1d.hpp:160
Poly1d< dtypeOut > astype() const
Definition: Poly1d.hpp:129
static Poly1d< double > fit(const NdArray< dtype > &xValues, const NdArray< dtype > &yValues, const NdArray< dtype > &weights, uint8 polyOrder)
Definition: Poly1d.hpp:246
Poly1d(const NdArray< dtype > &inValues, bool isRoots=false)
Definition: Poly1d.hpp:78
Poly1d< dtype > & operator+=(const Poly1d< dtype > &inOtherPoly)
Definition: Poly1d.hpp:453
Poly1d< double > integ() const
Definition: Poly1d.hpp:322
Poly1d< dtype > operator^(uint32 inPower) const
Definition: Poly1d.hpp:568
friend std::ostream & operator<<(std::ostream &inOStream, const Poly1d< dtype > &inPoly)
Definition: Poly1d.hpp:610
Poly1d< dtype > operator+(const Poly1d< dtype > &inOtherPoly) const
Definition: Poly1d.hpp:441
Poly1d< dtype > & operator-=(const Poly1d< dtype > &inOtherPoly)
Definition: Poly1d.hpp:496
Poly1d< dtype > & operator^=(uint32 inPower)
Definition: Poly1d.hpp:580
std::string str() const
Definition: Poly1d.hpp:368
Poly1d< dtype > operator*(const Poly1d< dtype > &inOtherPoly) const
Definition: Poly1d.hpp:527
void print() const
Definition: Poly1d.hpp:357
uint32 order() const noexcept
Definition: Poly1d.hpp:347
NdArray< dtype > coefficients() const
Definition: Poly1d.hpp:149
static Poly1d< double > fit(const NdArray< dtype > &xValues, const NdArray< dtype > &yValues, uint8 polyOrder)
Definition: Poly1d.hpp:191
dtype operator()(dtype inValue) const noexcept
Definition: Poly1d.hpp:422
double area(double a, double b) const
Definition: Poly1d.hpp:111
Poly1d< dtype > operator-(const Poly1d< dtype > &inOtherPoly) const
Definition: Poly1d.hpp:484
NdArray< double > inv(const NdArray< dtype > &inArray)
Definition: inv.hpp:52
OutputIt transform(InputIt first, InputIt last, OutputIt destination, UnaryOperation unaryFunction)
Definition: StlAlgorithms.hpp:702
OutputIt copy(InputIt first, InputIt last, OutputIt destination) noexcept
Definition: StlAlgorithms.hpp:95
std::string num2str(dtype inNumber)
Definition: num2str.hpp:46
dtype power(dtype inValue, uint8 inPower) noexcept
Definition: Utils/power.hpp:48
bool essentiallyEqual(dtype inValue1, dtype inValue2) noexcept
Definition: essentiallyEqual.hpp:52
Definition: Coordinate.hpp:45
constexpr dtype power(dtype inValue, uint8 inExponent) noexcept
Definition: Functions/power.hpp:52
void swap(NdArray< dtype > &inArray1, NdArray< dtype > &inArray2) noexcept
Definition: swap.hpp:42
std::uint8_t uint8
Definition: Types.hpp:42
std::uint32_t uint32
Definition: Types.hpp:40