NumCpp  2.10.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 
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:
186  static Poly1d<double> fit(const NdArray<dtype>& xValues, const NdArray<dtype>& yValues, uint8 polyOrder)
187  {
188  const auto numMeasurements = xValues.size();
189 
190  if (yValues.size() != numMeasurements)
191  {
192  THROW_INVALID_ARGUMENT_ERROR("Input x and y arrays must be of equal size.");
193  }
194 
195  if (!xValues.isflat())
196  {
197  THROW_INVALID_ARGUMENT_ERROR("Input x must be a flattened [1, n] or [n, 1] array.");
198  }
199 
200  if (!yValues.isflat())
201  {
202  THROW_INVALID_ARGUMENT_ERROR("Input y must be a flattened [n, 1] array.");
203  }
204 
205  NdArray<double> a(numMeasurements, polyOrder + 1);
206  for (uint32 measIdx = 0; measIdx < numMeasurements; ++measIdx)
207  {
208  const auto xDouble = static_cast<double>(xValues[measIdx]);
209  for (uint8 order = 0; order < a.numCols(); ++order)
210  {
211  a(measIdx, order) = utils::power(xDouble, order);
212  }
213  }
214 
215  NdArray<double> aInv;
216  if (a.issquare())
217  {
218  aInv = linalg::inv(a);
219  }
220  else
221  {
222  // psuedo-inverse
223  auto aT = a.transpose();
224  auto aTaInv = linalg::inv(aT.dot(a));
225  aInv = aTaInv.dot(aT);
226  }
227 
228  auto x = aInv.dot(yValues.template astype<double>());
229  return Poly1d<double>(x);
230  }
231 
232  //============================================================================
233  // Method Description:
241  static Poly1d<double> fit(const NdArray<dtype>& xValues,
242  const NdArray<dtype>& yValues,
243  const NdArray<dtype>& weights,
244  uint8 polyOrder)
245  {
246  const auto numMeasurements = xValues.size();
247 
248  if (yValues.size() != numMeasurements)
249  {
250  THROW_INVALID_ARGUMENT_ERROR("Input x and y arrays must be of equal size.");
251  }
252 
253  if (weights.size() != numMeasurements)
254  {
255  THROW_INVALID_ARGUMENT_ERROR("Input x and weights arrays must be of equal size.");
256  }
257 
258  if (!xValues.isflat())
259  {
260  THROW_INVALID_ARGUMENT_ERROR("Input x must be a flattened [1, n] or [n, 1] array.");
261  }
262 
263  if (!yValues.isflat())
264  {
265  THROW_INVALID_ARGUMENT_ERROR("Input y must be a flattened [n, 1] array.");
266  }
267 
268  if (!weights.isflat())
269  {
270  THROW_INVALID_ARGUMENT_ERROR("Input weights must be a flattened [1, n] or [n, 1] array.");
271  }
272 
273  NdArray<double> a(numMeasurements, polyOrder + 1);
274  for (uint32 measIdx = 0; measIdx < numMeasurements; ++measIdx)
275  {
276  const auto xDouble = static_cast<double>(xValues[measIdx]);
277  for (uint8 order = 0; order < a.numCols(); ++order)
278  {
279  a(measIdx, order) = utils::power(xDouble, order);
280  }
281  }
282 
283  NdArray<double> aWeighted(a.shape());
284  NdArray<double> yWeighted(yValues.shape());
285 
286  for (uint32 measIdx = 0; measIdx < numMeasurements; ++measIdx)
287  {
288  const auto weight = static_cast<double>(weights[measIdx]);
289 
290  yWeighted[measIdx] = yValues[measIdx] * weight;
291  for (uint8 order = 0; order < a.numCols(); ++order)
292  {
293  aWeighted(measIdx, order) = a(measIdx, order) * weight;
294  }
295  }
296 
297  NdArray<double> aInv;
298  if (aWeighted.issquare())
299  {
300  aInv = linalg::inv(aWeighted);
301  }
302  else
303  {
304  // psuedo-inverse
305  auto aT = a.transpose();
306  auto aTaInv = linalg::inv(aT.dot(aWeighted));
307  aInv = aTaInv.dot(aT);
308  }
309 
310  auto x = aInv.dot(yWeighted);
311  return Poly1d<double>(x); // NOLINT(modernize-return-braced-init-list)
312  }
313 
314  //============================================================================
315  // Method Description:
319  [[nodiscard]] Poly1d<double> integ() const
320  {
321  const auto numCoefficients = static_cast<uint32>(coefficients_.size());
322  if (numCoefficients == 0)
323  {
324  return {};
325  }
326 
327  NdArray<double> integralCofficients(1, numCoefficients + 1);
328  integralCofficients[0] = 0.;
329 
330  for (uint32 i = 0; i < numCoefficients; ++i)
331  {
332  integralCofficients[i + 1] = static_cast<double>(coefficients_[i]) / static_cast<double>(i + 1);
333  }
334 
335  return Poly1d<double>(integralCofficients); // NOLINT(modernize-return-braced-init-list)
336  }
337 
338  //============================================================================
339  // Method Description:
344  [[nodiscard]] uint32 order() const noexcept
345  {
346  return static_cast<uint32>(coefficients_.size() - 1);
347  }
348 
349  //============================================================================
350  // Method Description:
354  void print() const
355  {
356  std::cout << *this << std::endl;
357  }
358 
359  //============================================================================
360  // Method Description:
365  [[nodiscard]] std::string str() const
366  {
367  const auto numCoeffients = static_cast<uint32>(coefficients_.size());
368 
369  std::string repr = "Poly1d<";
370  uint32 power = 0;
371  for (auto& coefficient : coefficients_)
372  {
373  if (utils::essentiallyEqual(coefficient, static_cast<dtype>(0)))
374  {
375  ++power;
376  continue;
377  }
378 
379  repr += utils::num2str(coefficient);
380 
381  if (power > 1)
382  {
383  repr += "x^" + utils::num2str(power);
384  }
385  else if (power == 1)
386  {
387  repr += "x";
388  }
389 
390  ++power;
391 
392  if (power < numCoeffients)
393  {
394  repr += " + ";
395  }
396  }
397 
398  return repr + ">";
399  }
400 
401  //============================================================================
402  // Method Description:
408  dtype operator()(dtype inValue) const noexcept
409  {
410  uint8 power = 0;
411  return std::accumulate(coefficients_.begin(),
412  coefficients_.end(),
413  dtype{ 0 },
414  [&power, inValue](dtype polyValue, const auto& coefficient) noexcept -> dtype
415  { return polyValue + coefficient * utils::power(inValue, power++); });
416  }
417 
418  //============================================================================
419  // Method Description:
425  Poly1d<dtype> operator+(const Poly1d<dtype>& inOtherPoly) const
426  {
427  return Poly1d<dtype>(*this) += inOtherPoly;
428  }
429 
430  //============================================================================
431  // Method Description:
438  {
439  if (this->coefficients_.size() < inOtherPoly.coefficients_.size())
440  {
441  for (size_t i = 0; i < coefficients_.size(); ++i)
442  {
443  coefficients_[i] += inOtherPoly.coefficients_[i];
444  }
445  for (size_t i = coefficients_.size(); i < inOtherPoly.coefficients_.size(); ++i)
446  {
447  coefficients_.push_back(inOtherPoly.coefficients_[i]);
448  }
449  }
450  else
451  {
452  for (size_t i = 0; i < inOtherPoly.coefficients_.size(); ++i)
453  {
454  coefficients_[i] += inOtherPoly.coefficients_[i];
455  }
456  }
457 
458  return *this;
459  }
460 
461  //============================================================================
462  // Method Description:
468  Poly1d<dtype> operator-(const Poly1d<dtype>& inOtherPoly) const
469  {
470  return Poly1d<dtype>(*this) -= inOtherPoly;
471  }
472 
473  //============================================================================
474  // Method Description:
481  {
482  if (this->coefficients_.size() < inOtherPoly.coefficients_.size())
483  {
484  for (size_t i = 0; i < coefficients_.size(); ++i)
485  {
486  coefficients_[i] -= inOtherPoly.coefficients_[i];
487  }
488  for (size_t i = coefficients_.size(); i < inOtherPoly.coefficients_.size(); ++i)
489  {
490  coefficients_.push_back(-inOtherPoly.coefficients_[i]);
491  }
492  }
493  else
494  {
495  for (size_t i = 0; i < inOtherPoly.coefficients_.size(); ++i)
496  {
497  coefficients_[i] -= inOtherPoly.coefficients_[i];
498  }
499  }
500 
501  return *this;
502  }
503 
504  //============================================================================
505  // Method Description:
511  Poly1d<dtype> operator*(const Poly1d<dtype>& inOtherPoly) const
512  {
513  return Poly1d<dtype>(*this) *= inOtherPoly;
514  }
515 
516  //============================================================================
517  // Method Description:
524  {
525  const uint32 finalCoefficientsSize = order() + inOtherPoly.order() + 1;
526  std::vector<dtype> coeffsA(finalCoefficientsSize, 0);
527  std::vector<dtype> coeffsB(finalCoefficientsSize, 0);
528  stl_algorithms::copy(coefficients_.begin(), coefficients_.end(), coeffsA.begin());
529  stl_algorithms::copy(inOtherPoly.coefficients_.cbegin(), inOtherPoly.coefficients_.cend(), coeffsB.begin());
530 
531  // now multiply out the coefficients
532  std::vector<dtype> finalCoefficients(finalCoefficientsSize, 0);
533  for (uint32 i = 0; i < finalCoefficientsSize; ++i)
534  {
535  for (uint32 k = 0; k <= i; ++k)
536  {
537  finalCoefficients[i] += coeffsA[k] * coeffsB[i - k];
538  }
539  }
540 
541  this->coefficients_ = finalCoefficients;
542  return *this;
543  }
544 
545  //============================================================================
546  // Method Description:
553  {
554  return Poly1d(*this) ^= inPower;
555  }
556 
557  //============================================================================
558  // Method Description:
565  {
566  if (inPower == 0)
567  {
568  coefficients_.clear();
569  coefficients_.push_back(1);
570  return *this;
571  }
572  if (inPower == 1)
573  {
574  return *this;
575  }
576 
577  auto thisPoly(*this);
578  for (uint32 power = 1; power < inPower; ++power)
579  {
580  *this *= thisPoly;
581  }
582 
583  return *this;
584  }
585 
586  //============================================================================
587  // Method Description:
594  friend std::ostream& operator<<(std::ostream& inOStream, const Poly1d<dtype>& inPoly)
595  {
596  inOStream << inPoly.str() << std::endl;
597  return inOStream;
598  }
599 
600  private:
601  std::vector<dtype> coefficients_{};
602  };
603 } // 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:4415
iterator end() noexcept
Definition: NdArrayCore.hpp:1566
self_type transpose() const
Definition: NdArrayCore.hpp:4775
Shape shape() const noexcept
Definition: NdArrayCore.hpp:4402
bool issquare() const noexcept
Definition: NdArrayCore.hpp:2932
bool isflat() const noexcept
Definition: NdArrayCore.hpp:2868
self_type dot(const self_type &inOtherArray) const
Definition: NdArrayCore.hpp:2642
size_type numCols() const noexcept
Definition: NdArrayCore.hpp:3388
iterator begin() noexcept
Definition: NdArrayCore.hpp:1258
Definition: Poly1d.hpp:57
Poly1d< dtype > & operator*=(const Poly1d< dtype > &inOtherPoly)
Definition: Poly1d.hpp:523
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:241
Poly1d(const NdArray< dtype > &inValues, bool isRoots=false)
Definition: Poly1d.hpp:77
Poly1d< dtype > & operator+=(const Poly1d< dtype > &inOtherPoly)
Definition: Poly1d.hpp:437
Poly1d< double > integ() const
Definition: Poly1d.hpp:319
Poly1d< dtype > operator^(uint32 inPower) const
Definition: Poly1d.hpp:552
friend std::ostream & operator<<(std::ostream &inOStream, const Poly1d< dtype > &inPoly)
Definition: Poly1d.hpp:594
Poly1d< dtype > operator+(const Poly1d< dtype > &inOtherPoly) const
Definition: Poly1d.hpp:425
Poly1d< dtype > & operator-=(const Poly1d< dtype > &inOtherPoly)
Definition: Poly1d.hpp:480
Poly1d< dtype > & operator^=(uint32 inPower)
Definition: Poly1d.hpp:564
std::string str() const
Definition: Poly1d.hpp:365
Poly1d< dtype > operator*(const Poly1d< dtype > &inOtherPoly) const
Definition: Poly1d.hpp:511
void print() const
Definition: Poly1d.hpp:354
uint32 order() const noexcept
Definition: Poly1d.hpp:344
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:186
dtype operator()(dtype inValue) const noexcept
Definition: Poly1d.hpp:408
double area(double a, double b) const
Definition: Poly1d.hpp:109
Poly1d< dtype > operator-(const Poly1d< dtype > &inOtherPoly) const
Definition: Poly1d.hpp:468
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