NumCpp  2.11.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 
30 #include <iostream>
31 #include <numeric>
32 #include <string>
33 #include <type_traits>
34 #include <utility>
35 #include <vector>
36 
41 #include "NumCpp/Core/Types.hpp"
43 #include "NumCpp/Linalg/inv.hpp"
44 #include "NumCpp/NdArray.hpp"
46 #include "NumCpp/Utils/num2str.hpp"
47 #include "NumCpp/Utils/power.hpp"
48 
49 namespace nc::polynomial
50 {
51  //================================================================================
55  template<typename dtype>
56  class Poly1d
57  {
58  private:
59  STATIC_ASSERT_ARITHMETIC(dtype);
60 
61  public:
62  //============================================================================
63  // Method Description:
67  Poly1d() = default;
68 
69  //============================================================================
70  // Method Description:
77  Poly1d(const NdArray<dtype>& inValues, bool isRoots = false)
78  {
79  if (inValues.size() > DtypeInfo<uint8>::max())
80  {
81  THROW_INVALID_ARGUMENT_ERROR("can only make a polynomial of order " +
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  coefficients_.resize(inValues.size());
97  stl_algorithms::copy(inValues.begin(), inValues.end(), coefficients_.begin());
98  }
99  }
100 
101  //============================================================================
102  // Method Description:
109  [[nodiscard]] double area(double a, double b) const
110  {
111  if (a > b)
112  {
113  std::swap(a, b);
114  }
115 
116  auto polyIntegral = integ();
117  return polyIntegral(b) - polyIntegral(a);
118  }
119 
120  //============================================================================
121  // Method Description:
126  template<typename dtypeOut>
128  {
129  auto newCoefficients = NdArray<dtypeOut>(1, static_cast<uint32>(coefficients_.size()));
130 
131  const auto function = [](dtype value) -> dtypeOut { return static_cast<dtypeOut>(value); };
132 
133  stl_algorithms::transform(coefficients_.begin(), coefficients_.end(), newCoefficients.begin(), function);
134 
135  return Poly1d<dtypeOut>(newCoefficients);
136  }
137 
138  //============================================================================
139  // Method Description:
144  [[nodiscard]] NdArray<dtype> coefficients() const
145  {
146  auto coefficientsCopy = coefficients_;
147  return NdArray<dtype>(coefficientsCopy);
148  }
149 
150  //============================================================================
151  // Method Description:
155  [[nodiscard]] Poly1d<dtype> deriv() const
156  {
157  const auto numCoefficients = static_cast<uint32>(coefficients_.size());
158  if (numCoefficients == 0)
159  {
160  return {};
161  }
162  if (numCoefficients == 1)
163  {
164  return Poly1d<dtype>({ 0 });
165  }
166 
167  NdArray<dtype> derivativeCofficients(1, numCoefficients - 1);
168 
169  uint32 counter = 0;
170  for (uint32 i = 1; i < numCoefficients; ++i)
171  {
172  derivativeCofficients[counter++] = coefficients_[i] * i;
173  }
174 
175  return Poly1d<dtype>(derivativeCofficients);
176  }
177 
178  //============================================================================
179  // Method Description:
185  dtype eval(dtype xValue) const noexcept
186  {
187  return operator()(xValue);
188  }
189 
190  //============================================================================
191  // Method Description:
197  NdArray<dtype> eval(const NdArray<dtype>& xValues) const noexcept
198  {
199  return operator()(xValues);
200  }
201 
202  //============================================================================
203  // Method Description:
210  static Poly1d<double> fit(const NdArray<dtype>& xValues, const NdArray<dtype>& yValues, uint8 polyOrder)
211  {
212  const auto numMeasurements = xValues.size();
213 
214  if (yValues.size() != numMeasurements)
215  {
216  THROW_INVALID_ARGUMENT_ERROR("Input x and y arrays must be of equal size.");
217  }
218 
219  if (!xValues.isflat())
220  {
221  THROW_INVALID_ARGUMENT_ERROR("Input x must be a flattened [1, n] or [n, 1] array.");
222  }
223 
224  if (!yValues.isflat())
225  {
226  THROW_INVALID_ARGUMENT_ERROR("Input y must be a flattened [n, 1] array.");
227  }
228 
229  NdArray<double> a(numMeasurements, polyOrder + 1);
230  for (uint32 measIdx = 0; measIdx < numMeasurements; ++measIdx)
231  {
232  const auto xDouble = static_cast<double>(xValues[measIdx]);
233  for (uint8 order = 0; order < a.numCols(); ++order)
234  {
235  a(measIdx, order) = utils::power(xDouble, order);
236  }
237  }
238 
239  NdArray<double> aInv;
240  if (a.issquare())
241  {
242  aInv = linalg::inv(a);
243  }
244  else
245  {
246  // psuedo-inverse
247  auto aT = a.transpose();
248  auto aTaInv = linalg::inv(aT.dot(a));
249  aInv = aTaInv.dot(aT);
250  }
251 
252  auto x = aInv.dot(yValues.template astype<double>());
253  return Poly1d<double>(x);
254  }
255 
256  //============================================================================
257  // Method Description:
265  static Poly1d<double> fit(const NdArray<dtype>& xValues,
266  const NdArray<dtype>& yValues,
267  const NdArray<dtype>& weights,
268  uint8 polyOrder)
269  {
270  const auto numMeasurements = xValues.size();
271 
272  if (yValues.size() != numMeasurements)
273  {
274  THROW_INVALID_ARGUMENT_ERROR("Input x and y arrays must be of equal size.");
275  }
276 
277  if (weights.size() != numMeasurements)
278  {
279  THROW_INVALID_ARGUMENT_ERROR("Input x and weights arrays must be of equal size.");
280  }
281 
282  if (!xValues.isflat())
283  {
284  THROW_INVALID_ARGUMENT_ERROR("Input x must be a flattened [1, n] or [n, 1] array.");
285  }
286 
287  if (!yValues.isflat())
288  {
289  THROW_INVALID_ARGUMENT_ERROR("Input y must be a flattened [n, 1] array.");
290  }
291 
292  if (!weights.isflat())
293  {
294  THROW_INVALID_ARGUMENT_ERROR("Input weights must be a flattened [1, n] or [n, 1] array.");
295  }
296 
297  NdArray<double> a(numMeasurements, polyOrder + 1);
298  for (uint32 measIdx = 0; measIdx < numMeasurements; ++measIdx)
299  {
300  const auto xDouble = static_cast<double>(xValues[measIdx]);
301  for (uint8 order = 0; order < a.numCols(); ++order)
302  {
303  a(measIdx, order) = utils::power(xDouble, order);
304  }
305  }
306 
307  NdArray<double> aWeighted(a.shape());
308  NdArray<double> yWeighted(yValues.shape());
309 
310  for (uint32 measIdx = 0; measIdx < numMeasurements; ++measIdx)
311  {
312  const auto weight = static_cast<double>(weights[measIdx]);
313 
314  yWeighted[measIdx] = yValues[measIdx] * weight;
315  for (uint8 order = 0; order < a.numCols(); ++order)
316  {
317  aWeighted(measIdx, order) = a(measIdx, order) * weight;
318  }
319  }
320 
321  NdArray<double> aInv;
322  if (aWeighted.issquare())
323  {
324  aInv = linalg::inv(aWeighted);
325  }
326  else
327  {
328  // psuedo-inverse
329  auto aT = a.transpose();
330  auto aTaInv = linalg::inv(aT.dot(aWeighted));
331  aInv = aTaInv.dot(aT);
332  }
333 
334  auto x = aInv.dot(yWeighted);
335  return Poly1d<double>(x); // NOLINT(modernize-return-braced-init-list)
336  }
337 
338  //============================================================================
339  // Method Description:
343  [[nodiscard]] Poly1d<double> integ() const
344  {
345  const auto numCoefficients = static_cast<uint32>(coefficients_.size());
346  if (numCoefficients == 0)
347  {
348  return {};
349  }
350 
351  NdArray<double> integralCofficients(1, numCoefficients + 1);
352  integralCofficients[0] = 0.;
353 
354  for (uint32 i = 0; i < numCoefficients; ++i)
355  {
356  integralCofficients[i + 1] = static_cast<double>(coefficients_[i]) / static_cast<double>(i + 1);
357  }
358 
359  return Poly1d<double>(integralCofficients); // NOLINT(modernize-return-braced-init-list)
360  }
361 
362  //============================================================================
363  // Method Description:
368  [[nodiscard]] uint32 order() const noexcept
369  {
370  return static_cast<uint32>(coefficients_.size() - 1);
371  }
372 
373  //============================================================================
374  // Method Description:
378  void print() const
379  {
380  std::cout << *this << std::endl;
381  }
382 
383  //============================================================================
384  // Method Description:
389  [[nodiscard]] std::string str() const
390  {
391  const auto numCoeffients = static_cast<uint32>(coefficients_.size());
392 
393  std::string repr = "Poly1d<";
394  uint32 power = 0;
395  for (auto& coefficient : coefficients_)
396  {
397  if (utils::essentiallyEqual(coefficient, static_cast<dtype>(0)))
398  {
399  ++power;
400  continue;
401  }
402 
403  repr += utils::num2str(coefficient);
404 
405  if (power > 1)
406  {
407  repr += "x^" + utils::num2str(power);
408  }
409  else if (power == 1)
410  {
411  repr += "x";
412  }
413 
414  ++power;
415 
416  if (power < numCoeffients)
417  {
418  repr += " + ";
419  }
420  }
421 
422  return repr + ">";
423  }
424 
425  //============================================================================
426  // Method Description:
432  dtype operator()(dtype inValue) const noexcept
433  {
434  uint8 power = 0;
435  return std::accumulate(coefficients_.begin(),
436  coefficients_.end(),
437  dtype{ 0 },
438  [&power, inValue](dtype polyValue, const auto& coefficient) noexcept -> dtype
439  { return polyValue + coefficient * utils::power(inValue, power++); });
440  }
441 
442  //============================================================================
443  // Method Description:
449  NdArray<dtype> operator()(const NdArray<dtype>& xValues) const noexcept
450  {
451  NdArray<dtype> returnArray(xValues.shape());
452 
453  stl_algorithms::transform(xValues.begin(),
454  xValues.end(),
455  returnArray.begin(),
456  [this](const auto xValue) { return this->operator()(xValue); });
457  return returnArray;
458  }
459 
460  //============================================================================
461  // Method Description:
467  Poly1d<dtype> operator+(const Poly1d<dtype>& inOtherPoly) const
468  {
469  return Poly1d<dtype>(*this) += inOtherPoly;
470  }
471 
472  //============================================================================
473  // Method Description:
480  {
481  if (this->coefficients_.size() < inOtherPoly.coefficients_.size())
482  {
483  for (size_t i = 0; i < coefficients_.size(); ++i)
484  {
485  coefficients_[i] += inOtherPoly.coefficients_[i];
486  }
487  for (size_t i = coefficients_.size(); i < inOtherPoly.coefficients_.size(); ++i)
488  {
489  coefficients_.push_back(inOtherPoly.coefficients_[i]);
490  }
491  }
492  else
493  {
494  for (size_t i = 0; i < inOtherPoly.coefficients_.size(); ++i)
495  {
496  coefficients_[i] += inOtherPoly.coefficients_[i];
497  }
498  }
499 
500  return *this;
501  }
502 
503  //============================================================================
504  // Method Description:
510  Poly1d<dtype> operator-(const Poly1d<dtype>& inOtherPoly) const
511  {
512  return Poly1d<dtype>(*this) -= inOtherPoly;
513  }
514 
515  //============================================================================
516  // Method Description:
523  {
524  if (this->coefficients_.size() < inOtherPoly.coefficients_.size())
525  {
526  for (size_t i = 0; i < coefficients_.size(); ++i)
527  {
528  coefficients_[i] -= inOtherPoly.coefficients_[i];
529  }
530  for (size_t i = coefficients_.size(); i < inOtherPoly.coefficients_.size(); ++i)
531  {
532  coefficients_.push_back(-inOtherPoly.coefficients_[i]);
533  }
534  }
535  else
536  {
537  for (size_t i = 0; i < inOtherPoly.coefficients_.size(); ++i)
538  {
539  coefficients_[i] -= inOtherPoly.coefficients_[i];
540  }
541  }
542 
543  return *this;
544  }
545 
546  //============================================================================
547  // Method Description:
553  Poly1d<dtype> operator*(const Poly1d<dtype>& inOtherPoly) const
554  {
555  return Poly1d<dtype>(*this) *= inOtherPoly;
556  }
557 
558  //============================================================================
559  // Method Description:
566  {
567  const uint32 finalCoefficientsSize = order() + inOtherPoly.order() + 1;
568  std::vector<dtype> coeffsA(finalCoefficientsSize, 0);
569  std::vector<dtype> coeffsB(finalCoefficientsSize, 0);
570  stl_algorithms::copy(coefficients_.begin(), coefficients_.end(), coeffsA.begin());
571  stl_algorithms::copy(inOtherPoly.coefficients_.cbegin(), inOtherPoly.coefficients_.cend(), coeffsB.begin());
572 
573  // now multiply out the coefficients
574  std::vector<dtype> finalCoefficients(finalCoefficientsSize, 0);
575  for (uint32 i = 0; i < finalCoefficientsSize; ++i)
576  {
577  for (uint32 k = 0; k <= i; ++k)
578  {
579  finalCoefficients[i] += coeffsA[k] * coeffsB[i - k];
580  }
581  }
582 
583  this->coefficients_ = finalCoefficients;
584  return *this;
585  }
586 
587  //============================================================================
588  // Method Description:
595  {
596  return Poly1d(*this) ^= inPower;
597  }
598 
599  //============================================================================
600  // Method Description:
607  {
608  if (inPower == 0)
609  {
610  coefficients_.clear();
611  coefficients_.push_back(1);
612  return *this;
613  }
614  if (inPower == 1)
615  {
616  return *this;
617  }
618 
619  auto thisPoly(*this);
620  for (uint32 power = 1; power < inPower; ++power)
621  {
622  *this *= thisPoly;
623  }
624 
625  return *this;
626  }
627 
628  //============================================================================
629  // Method Description:
636  friend std::ostream& operator<<(std::ostream& inOStream, const Poly1d<dtype>& inPoly)
637  {
638  inOStream << inPoly.str() << std::endl;
639  return inOStream;
640  }
641 
642  private:
643  std::vector<dtype> coefficients_{};
644  };
645 } // namespace nc::polynomial
#define THROW_INVALID_ARGUMENT_ERROR(msg)
Definition: Error.hpp:37
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:138
size_type size() const noexcept
Definition: NdArrayCore.hpp:4477
iterator end() noexcept
Definition: NdArrayCore.hpp:1576
self_type transpose() const
Definition: NdArrayCore.hpp:4837
bool issquare() const noexcept
Definition: NdArrayCore.hpp:2962
bool isflat() const noexcept
Definition: NdArrayCore.hpp:2898
self_type dot(const self_type &inOtherArray) const
Definition: NdArrayCore.hpp:2672
size_type numCols() const noexcept
Definition: NdArrayCore.hpp:3418
iterator begin() noexcept
Definition: NdArrayCore.hpp:1268
const Shape & shape() const noexcept
Definition: NdArrayCore.hpp:4464
Definition: Poly1d.hpp:57
Poly1d< dtype > & operator*=(const Poly1d< dtype > &inOtherPoly)
Definition: Poly1d.hpp:565
Poly1d< dtype > deriv() const
Definition: Poly1d.hpp:155
Poly1d< dtypeOut > astype() const
Definition: Poly1d.hpp:127
static Poly1d< double > fit(const NdArray< dtype > &xValues, const NdArray< dtype > &yValues, const NdArray< dtype > &weights, uint8 polyOrder)
Definition: Poly1d.hpp:265
dtype eval(dtype xValue) const noexcept
Definition: Poly1d.hpp:185
Poly1d(const NdArray< dtype > &inValues, bool isRoots=false)
Definition: Poly1d.hpp:77
Poly1d< dtype > & operator+=(const Poly1d< dtype > &inOtherPoly)
Definition: Poly1d.hpp:479
Poly1d< double > integ() const
Definition: Poly1d.hpp:343
Poly1d< dtype > operator^(uint32 inPower) const
Definition: Poly1d.hpp:594
friend std::ostream & operator<<(std::ostream &inOStream, const Poly1d< dtype > &inPoly)
Definition: Poly1d.hpp:636
Poly1d< dtype > operator+(const Poly1d< dtype > &inOtherPoly) const
Definition: Poly1d.hpp:467
NdArray< dtype > operator()(const NdArray< dtype > &xValues) const noexcept
Definition: Poly1d.hpp:449
Poly1d< dtype > & operator-=(const Poly1d< dtype > &inOtherPoly)
Definition: Poly1d.hpp:522
Poly1d< dtype > & operator^=(uint32 inPower)
Definition: Poly1d.hpp:606
std::string str() const
Definition: Poly1d.hpp:389
Poly1d< dtype > operator*(const Poly1d< dtype > &inOtherPoly) const
Definition: Poly1d.hpp:553
void print() const
Definition: Poly1d.hpp:378
uint32 order() const noexcept
Definition: Poly1d.hpp:368
NdArray< dtype > coefficients() const
Definition: Poly1d.hpp:144
static Poly1d< double > fit(const NdArray< dtype > &xValues, const NdArray< dtype > &yValues, uint8 polyOrder)
Definition: Poly1d.hpp:210
NdArray< dtype > eval(const NdArray< dtype > &xValues) const noexcept
Definition: Poly1d.hpp:197
dtype operator()(dtype inValue) const noexcept
Definition: Poly1d.hpp:432
double area(double a, double b) const
Definition: Poly1d.hpp:109
Poly1d< dtype > operator-(const Poly1d< dtype > &inOtherPoly) const
Definition: Poly1d.hpp:510
NdArray< double > inv(const NdArray< dtype > &inArray)
Definition: inv.hpp:54
Definition: chebyshev_t.hpp:39
OutputIt transform(InputIt first, InputIt last, OutputIt destination, UnaryOperation unaryFunction)
Definition: StlAlgorithms.hpp:775
OutputIt copy(InputIt first, InputIt last, OutputIt destination) noexcept
Definition: StlAlgorithms.hpp:97
std::string num2str(dtype inNumber)
Definition: num2str.hpp:44
dtype power(dtype inValue, uint8 inPower) noexcept
Definition: Utils/power.hpp:46
bool essentiallyEqual(dtype inValue1, dtype inValue2) noexcept
Definition: essentiallyEqual.hpp:48
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