NumCpp  2.5.1
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:
151  {
152  auto coefficientsCopy = coefficients_;
153  return NdArray<dtype>(coefficientsCopy);
154  }
155 
156  //============================================================================
157  // Method Description:
162  {
163  const auto numCoefficients = static_cast<uint32>(coefficients_.size());
164  if (numCoefficients == 0)
165  {
166  return {};
167  }
168  if (numCoefficients == 1)
169  {
170  return Poly1d<dtype>({ 0 });
171  }
172 
173  NdArray<dtype> derivativeCofficients(1, numCoefficients - 1);
174 
175  uint32 counter = 0;
176  for (uint32 i = 1; i < numCoefficients; ++i)
177  {
178  derivativeCofficients[counter++] = coefficients_[i] * i;
179  }
180 
181  return Poly1d<dtype>(derivativeCofficients);
182  }
183 
184  //============================================================================
185  // Method Description:
192  static Poly1d<double> fit(const NdArray<dtype>& xValues, const NdArray<dtype>& yValues, uint8 polyOrder)
193  {
194  const auto numMeasurements = xValues.size();
195 
196  if (yValues.size() != numMeasurements)
197  {
198  THROW_INVALID_ARGUMENT_ERROR("Input x and y arrays must be of equal size.");
199  }
200 
201  if (!xValues.isflat())
202  {
203  THROW_INVALID_ARGUMENT_ERROR("Input x must be a flattened [1, n] or [n, 1] array.");
204  }
205 
206  if (!yValues.isflat())
207  {
208  THROW_INVALID_ARGUMENT_ERROR("Input y must be a flattened [n, 1] array.");
209  }
210 
211  NdArray<double> a(numMeasurements, polyOrder + 1);
212  for (uint32 measIdx = 0; measIdx < numMeasurements; ++measIdx)
213  {
214  const auto xDouble = static_cast<double>(xValues[measIdx]);
215  for (uint8 order = 0; order < a.numCols(); ++order)
216  {
217  a(measIdx, order) = utils::power(xDouble, order);
218  }
219  }
220 
221  NdArray<double> aInv;
222  if (a.issquare())
223  {
224  aInv = linalg::inv(a);
225  }
226  else
227  {
228  // psuedo-inverse
229  auto aT = a.transpose();
230  auto aTaInv = linalg::inv(aT.dot(a));
231  aInv = aTaInv.dot(aT);
232  }
233 
234  auto x = aInv.dot(yValues.template astype<double>());
235  return Poly1d<double>(x);
236  }
237 
238  //============================================================================
239  // Method Description:
247  static Poly1d<double> fit(const NdArray<dtype>& xValues, const NdArray<dtype>& yValues,
248  const NdArray<dtype>& weights, uint8 polyOrder)
249  {
250  const auto numMeasurements = xValues.size();
251 
252  if (yValues.size() != numMeasurements)
253  {
254  THROW_INVALID_ARGUMENT_ERROR("Input x and y arrays must be of equal size.");
255  }
256 
257  if (weights.size() != numMeasurements)
258  {
259  THROW_INVALID_ARGUMENT_ERROR("Input x and weights arrays must be of equal size.");
260  }
261 
262  if (!xValues.isflat())
263  {
264  THROW_INVALID_ARGUMENT_ERROR("Input x must be a flattened [1, n] or [n, 1] array.");
265  }
266 
267  if (!yValues.isflat())
268  {
269  THROW_INVALID_ARGUMENT_ERROR("Input y must be a flattened [n, 1] array.");
270  }
271 
272  if (!weights.isflat())
273  {
274  THROW_INVALID_ARGUMENT_ERROR("Input weights must be a flattened [1, n] or [n, 1] array.");
275  }
276 
277  NdArray<double> a(numMeasurements, polyOrder + 1);
278  for (uint32 measIdx = 0; measIdx < numMeasurements; ++measIdx)
279  {
280  const auto xDouble = static_cast<double>(xValues[measIdx]);
281  for (uint8 order = 0; order < a.numCols(); ++order)
282  {
283  a(measIdx, order) = utils::power(xDouble, order);
284  }
285  }
286 
287  NdArray<double> aWeighted(a.shape());
288  NdArray<double> yWeighted(yValues.shape());
289 
290  for (uint32 measIdx = 0; measIdx < numMeasurements; ++measIdx)
291  {
292  const auto weight = static_cast<double>(weights[measIdx]);
293 
294  yWeighted[measIdx] = yValues[measIdx] * weight;
295  for (uint8 order = 0; order < a.numCols(); ++order)
296  {
297  aWeighted(measIdx, order) = a(measIdx, order) * weight;
298  }
299  }
300 
301  NdArray<double> aInv;
302  if (aWeighted.issquare())
303  {
304  aInv = linalg::inv(aWeighted);
305  }
306  else
307  {
308  // psuedo-inverse
309  auto aT = a.transpose();
310  auto aTaInv = linalg::inv(aT.dot(aWeighted));
311  aInv = aTaInv.dot(aT);
312  }
313 
314  auto x = aInv.dot(yWeighted);
315  return Poly1d<double>(x);
316  }
317 
318  //============================================================================
319  // Method Description:
324  {
325  const auto numCoefficients = static_cast<uint32>(coefficients_.size());
326  if (numCoefficients == 0)
327  {
328  return {};
329  }
330 
331  NdArray<double> integralCofficients(1, numCoefficients + 1);
332  integralCofficients[0] = 0.0;
333 
334  for (uint32 i = 0; i < numCoefficients; ++i)
335  {
336  integralCofficients[i + 1] = static_cast<double>(coefficients_[i]) / static_cast<double>(i + 1);
337  }
338 
339  return Poly1d<double>(integralCofficients);
340  }
341 
342  //============================================================================
343  // Method Description:
349  uint32 order() const noexcept
350  {
351  return static_cast<uint32>(coefficients_.size() - 1);
352  }
353 
354  //============================================================================
355  // Method Description:
359  void print() const
360  {
361  std::cout << *this << std::endl;
362  }
363 
364  //============================================================================
365  // Method Description:
371  std::string str() const
372  {
373  const auto numCoeffients = static_cast<uint32>(coefficients_.size());
374 
375  std::string repr = "Poly1d<";
376  uint32 power = 0;
377  for (auto& coefficient : coefficients_)
378  {
380  {
381  if (coefficient == 0)
382  {
383  ++power;
384  continue;
385  }
386  }
387  else
388  {
389  if (utils::essentiallyEqual(coefficient, static_cast<dtype>(0.0)))
390  {
391  ++power;
392  continue;
393  }
394  }
395 
396  repr += utils::num2str(coefficient);
397 
398  if (power > 1)
399  {
400  repr += "x^" + utils::num2str(power);
401  }
402  else if (power == 1)
403  {
404  repr += "x";
405  }
406 
407  ++power;
408 
409  if (power < numCoeffients)
410  {
411  repr += " + ";
412  }
413  }
414 
415  return repr + ">";
416  }
417 
418  //============================================================================
419  // Method Description:
427  dtype operator()(dtype inValue) const noexcept
428  {
429  dtype polyValue = 0;
430  uint8 power = 0;
431  for (auto& coefficient : coefficients_)
432  {
433  polyValue += coefficient * utils::power(inValue, power++);
434  }
435 
436  return polyValue;
437  }
438 
439  //============================================================================
440  // Method Description:
448  Poly1d<dtype> operator+(const Poly1d<dtype>& inOtherPoly) const
449  {
450  return Poly1d<dtype>(*this) += inOtherPoly;
451  }
452 
453  //============================================================================
454  // Method Description:
462  Poly1d<dtype>& operator+=(const Poly1d<dtype>& inOtherPoly)
463  {
464  if (this->coefficients_.size() < inOtherPoly.coefficients_.size())
465  {
466  for (size_t i = 0; i < coefficients_.size(); ++i)
467  {
468  coefficients_[i] += inOtherPoly.coefficients_[i];
469  }
470  for (size_t i = coefficients_.size(); i < inOtherPoly.coefficients_.size(); ++i)
471  {
472  coefficients_.push_back(inOtherPoly.coefficients_[i]);
473  }
474  }
475  else
476  {
477  for (size_t i = 0; i < inOtherPoly.coefficients_.size(); ++i)
478  {
479  coefficients_[i] += inOtherPoly.coefficients_[i];
480  }
481  }
482 
483  return *this;
484  }
485 
486  //============================================================================
487  // Method Description:
495  Poly1d<dtype> operator-(const Poly1d<dtype>& inOtherPoly) const
496  {
497  return Poly1d<dtype>(*this) -= inOtherPoly;
498  }
499 
500  //============================================================================
501  // Method Description:
509  Poly1d<dtype>& operator-=(const Poly1d<dtype>& inOtherPoly)
510  {
511  if (this->coefficients_.size() < inOtherPoly.coefficients_.size())
512  {
513  for (size_t i = 0; i < coefficients_.size(); ++i)
514  {
515  coefficients_[i] -= inOtherPoly.coefficients_[i];
516  }
517  for (size_t i = coefficients_.size(); i < inOtherPoly.coefficients_.size(); ++i)
518  {
519  coefficients_.push_back(-inOtherPoly.coefficients_[i]);
520  }
521  }
522  else
523  {
524  for (size_t i = 0; i < inOtherPoly.coefficients_.size(); ++i)
525  {
526  coefficients_[i] -= inOtherPoly.coefficients_[i];
527  }
528  }
529 
530  return *this;
531  }
532 
533  //============================================================================
534  // Method Description:
542  Poly1d<dtype> operator*(const Poly1d<dtype>& inOtherPoly) const
543  {
544  return Poly1d<dtype>(*this) *= inOtherPoly;
545  }
546 
547  //============================================================================
548  // Method Description:
556  Poly1d<dtype>& operator*=(const Poly1d<dtype>& inOtherPoly)
557  {
558  const uint32 finalCoefficientsSize = order() + inOtherPoly.order() + 1;
559  std::vector<dtype> coeffsA(finalCoefficientsSize, 0);
560  std::vector<dtype> coeffsB(finalCoefficientsSize, 0);
561  stl_algorithms::copy(coefficients_.begin(), coefficients_.end(), coeffsA.begin());
562  stl_algorithms::copy(inOtherPoly.coefficients_.cbegin(), inOtherPoly.coefficients_.cend(), coeffsB.begin());
563 
564  // now multiply out the coefficients
565  std::vector<dtype> finalCoefficients(finalCoefficientsSize, 0);
566  for (uint32 i = 0; i < finalCoefficientsSize; ++i)
567  {
568  for (uint32 k = 0; k <= i; ++k)
569  {
570  finalCoefficients[i] += coeffsA[k] * coeffsB[i - k];
571  }
572  }
573 
574  this->coefficients_ = finalCoefficients;
575  return *this;
576  }
577 
578  //============================================================================
579  // Method Description:
588  {
589  return Poly1d(*this) ^= inPower;
590  }
591 
592  //============================================================================
593  // Method Description:
602  {
603  if (inPower == 0)
604  {
605  coefficients_.clear();
606  coefficients_.push_back(1);
607  return *this;
608  }
609  if (inPower == 1)
610  {
611  return *this;
612  }
613 
614  auto thisPoly(*this);
615  for (uint32 power = 1; power < inPower; ++power)
616  {
617  *this *= thisPoly;
618  }
619 
620  return *this;
621  }
622 
623  //============================================================================
624  // Method Description:
632  friend std::ostream& operator<<(std::ostream& inOStream, const Poly1d<dtype>& inPoly)
633  {
634  inOStream << inPoly.str() << std::endl;
635  return inOStream;
636  }
637 
638  private:
639  std::vector<dtype> coefficients_{};
640  };
641  } // namespace polynomial
642 } // 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:3602
size_type size() const noexcept
Definition: NdArrayCore.hpp:4497
Shape shape() const noexcept
Definition: NdArrayCore.hpp:4483
bool issquare() const noexcept
Definition: NdArrayCore.hpp:3088
bool isflat() const noexcept
Definition: NdArrayCore.hpp:3025
NdArray< dtype > transpose() const
Definition: NdArrayCore.hpp:4830
NdArray< dtype > dot(const NdArray< dtype > &inOtherArray) const
Definition: NdArrayCore.hpp:2788
Definition: Poly1d.hpp:58
Poly1d< dtype > & operator*=(const Poly1d< dtype > &inOtherPoly)
Definition: Poly1d.hpp:556
Poly1d< dtype > deriv() const
Definition: Poly1d.hpp:161
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:247
Poly1d(const NdArray< dtype > &inValues, bool isRoots=false)
Definition: Poly1d.hpp:78
Poly1d< dtype > & operator+=(const Poly1d< dtype > &inOtherPoly)
Definition: Poly1d.hpp:462
Poly1d< double > integ() const
Definition: Poly1d.hpp:323
Poly1d< dtype > operator^(uint32 inPower) const
Definition: Poly1d.hpp:587
friend std::ostream & operator<<(std::ostream &inOStream, const Poly1d< dtype > &inPoly)
Definition: Poly1d.hpp:632
Poly1d< dtype > operator+(const Poly1d< dtype > &inOtherPoly) const
Definition: Poly1d.hpp:448
Poly1d< dtype > & operator-=(const Poly1d< dtype > &inOtherPoly)
Definition: Poly1d.hpp:509
Poly1d< dtype > & operator^=(uint32 inPower)
Definition: Poly1d.hpp:601
std::string str() const
Definition: Poly1d.hpp:371
Poly1d< dtype > operator*(const Poly1d< dtype > &inOtherPoly) const
Definition: Poly1d.hpp:542
void print() const
Definition: Poly1d.hpp:359
uint32 order() const noexcept
Definition: Poly1d.hpp:349
NdArray< dtype > coefficients() const
Definition: Poly1d.hpp:150
static Poly1d< double > fit(const NdArray< dtype > &xValues, const NdArray< dtype > &yValues, uint8 polyOrder)
Definition: Poly1d.hpp:192
dtype operator()(dtype inValue) const noexcept
Definition: Poly1d.hpp:427
double area(double a, double b) const
Definition: Poly1d.hpp:111
Poly1d< dtype > operator-(const Poly1d< dtype > &inOtherPoly) const
Definition: Poly1d.hpp:495
NdArray< double > inv(const NdArray< dtype > &inArray)
Definition: inv.hpp:54
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:53
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