NumCpp  2.9.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 <string>
32 #include <type_traits>
33 #include <utility>
34 #include <vector>
35 
40 #include "NumCpp/Core/Types.hpp"
42 #include "NumCpp/Linalg/inv.hpp"
43 #include "NumCpp/NdArray.hpp"
45 #include "NumCpp/Utils/num2str.hpp"
46 #include "NumCpp/Utils/power.hpp"
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 " +
84  }
85 
86  if (isRoots)
87  {
88  coefficients_.push_back(1);
89  for (auto value : inValues)
90  {
91  NdArray<dtype> coeffs = { -(value), static_cast<dtype>(1) };
92  *this *= Poly1d<dtype>(coeffs, !isRoots);
93  }
94  }
95  else
96  {
97  for (auto value : inValues)
98  {
99  coefficients_.push_back(value);
100  }
101  }
102  }
103 
104  //============================================================================
105  // Method Description:
112  double area(double a, double b) const
113  {
114  if (a > b)
115  {
116  std::swap(a, b);
117  }
118 
119  auto polyIntegral = integ();
120  return polyIntegral(b) - polyIntegral(a);
121  }
122 
123  //============================================================================
124  // Method Description:
129  template<typename dtypeOut>
131  {
132  auto newCoefficients = NdArray<dtypeOut>(1, static_cast<uint32>(coefficients_.size()));
133 
134  const auto function = [](dtype value) -> dtypeOut { return static_cast<dtypeOut>(value); };
135 
136  stl_algorithms::transform(coefficients_.begin(),
137  coefficients_.end(),
138  newCoefficients.begin(),
139  function);
140 
141  return Poly1d<dtypeOut>(newCoefficients);
142  }
143 
144  //============================================================================
145  // 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,
248  const NdArray<dtype>& yValues,
249  const NdArray<dtype>& weights,
250  uint8 polyOrder)
251  {
252  const auto numMeasurements = xValues.size();
253 
254  if (yValues.size() != numMeasurements)
255  {
256  THROW_INVALID_ARGUMENT_ERROR("Input x and y arrays must be of equal size.");
257  }
258 
259  if (weights.size() != numMeasurements)
260  {
261  THROW_INVALID_ARGUMENT_ERROR("Input x and weights arrays must be of equal size.");
262  }
263 
264  if (!xValues.isflat())
265  {
266  THROW_INVALID_ARGUMENT_ERROR("Input x must be a flattened [1, n] or [n, 1] array.");
267  }
268 
269  if (!yValues.isflat())
270  {
271  THROW_INVALID_ARGUMENT_ERROR("Input y must be a flattened [n, 1] array.");
272  }
273 
274  if (!weights.isflat())
275  {
276  THROW_INVALID_ARGUMENT_ERROR("Input weights must be a flattened [1, n] or [n, 1] array.");
277  }
278 
279  NdArray<double> a(numMeasurements, polyOrder + 1);
280  for (uint32 measIdx = 0; measIdx < numMeasurements; ++measIdx)
281  {
282  const auto xDouble = static_cast<double>(xValues[measIdx]);
283  for (uint8 order = 0; order < a.numCols(); ++order)
284  {
285  a(measIdx, order) = utils::power(xDouble, order);
286  }
287  }
288 
289  NdArray<double> aWeighted(a.shape());
290  NdArray<double> yWeighted(yValues.shape());
291 
292  for (uint32 measIdx = 0; measIdx < numMeasurements; ++measIdx)
293  {
294  const auto weight = static_cast<double>(weights[measIdx]);
295 
296  yWeighted[measIdx] = yValues[measIdx] * weight;
297  for (uint8 order = 0; order < a.numCols(); ++order)
298  {
299  aWeighted(measIdx, order) = a(measIdx, order) * weight;
300  }
301  }
302 
303  NdArray<double> aInv;
304  if (aWeighted.issquare())
305  {
306  aInv = linalg::inv(aWeighted);
307  }
308  else
309  {
310  // psuedo-inverse
311  auto aT = a.transpose();
312  auto aTaInv = linalg::inv(aT.dot(aWeighted));
313  aInv = aTaInv.dot(aT);
314  }
315 
316  auto x = aInv.dot(yWeighted);
317  return Poly1d<double>(x);
318  }
319 
320  //============================================================================
321  // Method Description:
326  {
327  const auto numCoefficients = static_cast<uint32>(coefficients_.size());
328  if (numCoefficients == 0)
329  {
330  return {};
331  }
332 
333  NdArray<double> integralCofficients(1, numCoefficients + 1);
334  integralCofficients[0] = 0.;
335 
336  for (uint32 i = 0; i < numCoefficients; ++i)
337  {
338  integralCofficients[i + 1] = static_cast<double>(coefficients_[i]) / static_cast<double>(i + 1);
339  }
340 
341  return Poly1d<double>(integralCofficients);
342  }
343 
344  //============================================================================
345  // Method Description:
350  uint32 order() const noexcept
351  {
352  return static_cast<uint32>(coefficients_.size() - 1);
353  }
354 
355  //============================================================================
356  // Method Description:
360  void print() const
361  {
362  std::cout << *this << std::endl;
363  }
364 
365  //============================================================================
366  // 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.)))
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:
425  dtype operator()(dtype inValue) const noexcept
426  {
427  dtype polyValue = 0;
428  uint8 power = 0;
429  for (auto& coefficient : coefficients_)
430  {
431  polyValue += coefficient * utils::power(inValue, power++);
432  }
433 
434  return polyValue;
435  }
436 
437  //============================================================================
438  // Method Description:
444  Poly1d<dtype> operator+(const Poly1d<dtype>& inOtherPoly) const
445  {
446  return Poly1d<dtype>(*this) += inOtherPoly;
447  }
448 
449  //============================================================================
450  // Method Description:
457  {
458  if (this->coefficients_.size() < inOtherPoly.coefficients_.size())
459  {
460  for (size_t i = 0; i < coefficients_.size(); ++i)
461  {
462  coefficients_[i] += inOtherPoly.coefficients_[i];
463  }
464  for (size_t i = coefficients_.size(); i < inOtherPoly.coefficients_.size(); ++i)
465  {
466  coefficients_.push_back(inOtherPoly.coefficients_[i]);
467  }
468  }
469  else
470  {
471  for (size_t i = 0; i < inOtherPoly.coefficients_.size(); ++i)
472  {
473  coefficients_[i] += inOtherPoly.coefficients_[i];
474  }
475  }
476 
477  return *this;
478  }
479 
480  //============================================================================
481  // Method Description:
487  Poly1d<dtype> operator-(const Poly1d<dtype>& inOtherPoly) const
488  {
489  return Poly1d<dtype>(*this) -= inOtherPoly;
490  }
491 
492  //============================================================================
493  // Method Description:
500  {
501  if (this->coefficients_.size() < inOtherPoly.coefficients_.size())
502  {
503  for (size_t i = 0; i < coefficients_.size(); ++i)
504  {
505  coefficients_[i] -= inOtherPoly.coefficients_[i];
506  }
507  for (size_t i = coefficients_.size(); i < inOtherPoly.coefficients_.size(); ++i)
508  {
509  coefficients_.push_back(-inOtherPoly.coefficients_[i]);
510  }
511  }
512  else
513  {
514  for (size_t i = 0; i < inOtherPoly.coefficients_.size(); ++i)
515  {
516  coefficients_[i] -= inOtherPoly.coefficients_[i];
517  }
518  }
519 
520  return *this;
521  }
522 
523  //============================================================================
524  // Method Description:
530  Poly1d<dtype> operator*(const Poly1d<dtype>& inOtherPoly) const
531  {
532  return Poly1d<dtype>(*this) *= inOtherPoly;
533  }
534 
535  //============================================================================
536  // Method Description:
543  {
544  const uint32 finalCoefficientsSize = order() + inOtherPoly.order() + 1;
545  std::vector<dtype> coeffsA(finalCoefficientsSize, 0);
546  std::vector<dtype> coeffsB(finalCoefficientsSize, 0);
547  stl_algorithms::copy(coefficients_.begin(), coefficients_.end(), coeffsA.begin());
548  stl_algorithms::copy(inOtherPoly.coefficients_.cbegin(),
549  inOtherPoly.coefficients_.cend(),
550  coeffsB.begin());
551 
552  // now multiply out the coefficients
553  std::vector<dtype> finalCoefficients(finalCoefficientsSize, 0);
554  for (uint32 i = 0; i < finalCoefficientsSize; ++i)
555  {
556  for (uint32 k = 0; k <= i; ++k)
557  {
558  finalCoefficients[i] += coeffsA[k] * coeffsB[i - k];
559  }
560  }
561 
562  this->coefficients_ = finalCoefficients;
563  return *this;
564  }
565 
566  //============================================================================
567  // Method Description:
574  {
575  return Poly1d(*this) ^= inPower;
576  }
577 
578  //============================================================================
579  // Method Description:
586  {
587  if (inPower == 0)
588  {
589  coefficients_.clear();
590  coefficients_.push_back(1);
591  return *this;
592  }
593  if (inPower == 1)
594  {
595  return *this;
596  }
597 
598  auto thisPoly(*this);
599  for (uint32 power = 1; power < inPower; ++power)
600  {
601  *this *= thisPoly;
602  }
603 
604  return *this;
605  }
606 
607  //============================================================================
608  // Method Description:
615  friend std::ostream& operator<<(std::ostream& inOStream, const Poly1d<dtype>& inPoly)
616  {
617  inOStream << inPoly.str() << std::endl;
618  return inOStream;
619  }
620 
621  private:
622  std::vector<dtype> coefficients_{};
623  };
624  } // namespace polynomial
625 } // 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:3252
size_type size() const noexcept
Definition: NdArrayCore.hpp:4105
Shape shape() const noexcept
Definition: NdArrayCore.hpp:4092
bool issquare() const noexcept
Definition: NdArrayCore.hpp:2796
bool isflat() const noexcept
Definition: NdArrayCore.hpp:2744
NdArray< dtype > transpose() const
Definition: NdArrayCore.hpp:4457
NdArray< dtype > dot(const NdArray< dtype > &inOtherArray) const
Definition: NdArrayCore.hpp:2520
Definition: Poly1d.hpp:58
Poly1d< dtype > & operator*=(const Poly1d< dtype > &inOtherPoly)
Definition: Poly1d.hpp:542
Poly1d< dtype > deriv() const
Definition: Poly1d.hpp:161
Poly1d< dtypeOut > astype() const
Definition: Poly1d.hpp:130
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:456
Poly1d< double > integ() const
Definition: Poly1d.hpp:325
Poly1d< dtype > operator^(uint32 inPower) const
Definition: Poly1d.hpp:573
friend std::ostream & operator<<(std::ostream &inOStream, const Poly1d< dtype > &inPoly)
Definition: Poly1d.hpp:615
Poly1d< dtype > operator+(const Poly1d< dtype > &inOtherPoly) const
Definition: Poly1d.hpp:444
Poly1d< dtype > & operator-=(const Poly1d< dtype > &inOtherPoly)
Definition: Poly1d.hpp:499
Poly1d< dtype > & operator^=(uint32 inPower)
Definition: Poly1d.hpp:585
std::string str() const
Definition: Poly1d.hpp:371
Poly1d< dtype > operator*(const Poly1d< dtype > &inOtherPoly) const
Definition: Poly1d.hpp:530
void print() const
Definition: Poly1d.hpp:360
uint32 order() const noexcept
Definition: Poly1d.hpp:350
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:425
double area(double a, double b) const
Definition: Poly1d.hpp:112
Poly1d< dtype > operator-(const Poly1d< dtype > &inOtherPoly) const
Definition: Poly1d.hpp:487
NdArray< double > inv(const NdArray< dtype > &inArray)
Definition: inv.hpp:56
OutputIt transform(InputIt first, InputIt last, OutputIt destination, UnaryOperation unaryFunction)
Definition: StlAlgorithms.hpp:784
OutputIt copy(InputIt first, InputIt last, OutputIt destination) noexcept
Definition: StlAlgorithms.hpp:99
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:51
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