NumCpp  1.0
A C++ implementation of the Python Numpy library
Poly1d.hpp
Go to the documentation of this file.
1 #pragma once
30 
35 #include "NumCpp/Core/Types.hpp"
37 #include "NumCpp/Linalg/inv.hpp"
38 #include "NumCpp/NdArray.hpp"
40 #include "NumCpp/Utils/num2str.hpp"
41 #include "NumCpp/Utils/power.hpp"
42 
43 #include <iostream>
44 #include <string>
45 #include <type_traits>
46 #include <utility>
47 #include <vector>
48 
49 namespace nc
50 {
51  namespace polynomial
52  {
53  //================================================================================
57  template<typename dtype>
58  class Poly1d
59  {
60  private:
61  STATIC_ASSERT_ARITHMETIC(dtype);
62 
63  public:
64  //============================================================================
65  // Method Description:
69  Poly1d() = default;
70 
71  //============================================================================
72  // Method Description:
79  Poly1d(const NdArray<dtype>& inValues, bool isRoots = false)
80  {
81  if (inValues.size() > DtypeInfo<uint8>::max())
82  {
83  THROW_INVALID_ARGUMENT_ERROR("can only make a polynomial of order " + utils::num2str(DtypeInfo<uint8>::max()));
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
135  {
136  return static_cast<dtypeOut>(value);
137  };
138 
139  stl_algorithms::transform(coefficients_.begin(), coefficients_.end(), newCoefficients.begin(), function);
140 
141  return Poly1d<dtypeOut>(newCoefficients);
142  }
143 
144  //============================================================================
145  // Method Description:
152  {
153  auto coefficientsCopy = coefficients_;
154  return NdArray<dtype>(coefficientsCopy);
155  }
156 
157  //============================================================================
158  // Method Description:
163  {
164  const uint32 numCoefficients = static_cast<uint32>(coefficients_.size());
165  if (numCoefficients == 0)
166  {
167  return {};
168  }
169  else if (numCoefficients == 1)
170  {
171  return Poly1d<dtype>({ 0 });
172  }
173 
174  NdArray<dtype> derivativeCofficients(1, numCoefficients - 1);
175 
176  uint32 counter = 0;
177  for (uint32 i = 1; i < numCoefficients; ++i)
178  {
179  derivativeCofficients[counter++] = coefficients_[i] * i;
180  }
181 
182  return Poly1d<dtype>(derivativeCofficients);
183  }
184 
185  //============================================================================
186  // Method Description:
193  static Poly1d<double> fit(const NdArray<dtype>& xValues, const NdArray<dtype>& yValues, uint8 polyOrder)
194  {
195  const auto numMeasurements = xValues.size();
196 
197  if (yValues.size() != numMeasurements)
198  {
199  THROW_INVALID_ARGUMENT_ERROR("Input x and y arrays must be of equal size.");
200  }
201 
202  if (!xValues.isflat())
203  {
204  THROW_INVALID_ARGUMENT_ERROR("Input x must be a flattened [1, n] or [n, 1] array.");
205  }
206 
207  if (!yValues.isflat())
208  {
209  THROW_INVALID_ARGUMENT_ERROR("Input y must be a flattened [n, 1] array.");
210  }
211 
212  NdArray<double> a(numMeasurements, polyOrder + 1);
213  for (uint32 measIdx = 0; measIdx < numMeasurements; ++measIdx)
214  {
215  const auto xDouble = static_cast<double>(xValues[measIdx]);
216  for (uint8 order = 0; order < a.numCols(); ++order)
217  {
218  a(measIdx, order) = utils::power(xDouble, order);
219  }
220  }
221 
222  NdArray<double> aInv;
223  if (a.issquare())
224  {
225  aInv = linalg::inv(a);
226  }
227  else
228  {
229  // psuedo-inverse
230  auto aT = a.transpose();
231  auto aTaInv = linalg::inv(aT.dot(a));
232  aInv = aTaInv.dot(aT);
233  }
234 
235  auto x = aInv.dot(yValues.template astype<double>());
236  return Poly1d<double>(x);
237  }
238 
239  //============================================================================
240  // Method Description:
248  static Poly1d<double> fit(const NdArray<dtype>& xValues, const NdArray<dtype>& yValues,
249  const NdArray<dtype>& weights, uint8 polyOrder)
250  {
251  const auto numMeasurements = xValues.size();
252 
253  if (yValues.size() != numMeasurements)
254  {
255  THROW_INVALID_ARGUMENT_ERROR("Input x and y arrays must be of equal size.");
256  }
257 
258  if (weights.size() != numMeasurements)
259  {
260  THROW_INVALID_ARGUMENT_ERROR("Input x and weights arrays must be of equal size.");
261  }
262 
263  if (!xValues.isflat())
264  {
265  THROW_INVALID_ARGUMENT_ERROR("Input x must be a flattened [1, n] or [n, 1] array.");
266  }
267 
268  if (!yValues.isflat())
269  {
270  THROW_INVALID_ARGUMENT_ERROR("Input y must be a flattened [n, 1] array.");
271  }
272 
273  if (!weights.isflat())
274  {
275  THROW_INVALID_ARGUMENT_ERROR("Input weights must be a flattened [1, n] or [n, 1] array.");
276  }
277 
278  NdArray<double> a(numMeasurements, polyOrder + 1);
279  for (uint32 measIdx = 0; measIdx < numMeasurements; ++measIdx)
280  {
281  const auto xDouble = static_cast<double>(xValues[measIdx]);
282  for (uint8 order = 0; order < a.numCols(); ++order)
283  {
284  a(measIdx, order) = utils::power(xDouble, order);
285  }
286  }
287 
288  NdArray<double> aWeighted(a.shape());
289  NdArray<double> yWeighted(yValues.shape());
290 
291  for (uint32 measIdx = 0; measIdx < numMeasurements; ++measIdx)
292  {
293  const auto weight = static_cast<double>(weights[measIdx]);
294 
295  yWeighted[measIdx] = yValues[measIdx] * weight;
296  for (uint8 order = 0; order < a.numCols(); ++order)
297  {
298  aWeighted(measIdx, order) = a(measIdx, order) * weight;
299  }
300  }
301 
302  NdArray<double> aInv;
303  if (aWeighted.issquare())
304  {
305  aInv = linalg::inv(aWeighted);
306  }
307  else
308  {
309  // psuedo-inverse
310  auto aT = a.transpose();
311  auto aTaInv = linalg::inv(aT.dot(aWeighted));
312  aInv = aTaInv.dot(aT);
313  }
314 
315  auto x = aInv.dot(yWeighted);
316  return Poly1d<double>(x);
317  }
318 
319  //============================================================================
320  // Method Description:
325  {
326  const uint32 numCoefficients = static_cast<uint32>(coefficients_.size());
327  if (numCoefficients == 0)
328  {
329  return {};
330  }
331 
332  NdArray<double> integralCofficients(1, numCoefficients + 1);
333  integralCofficients[0] = 0.0;
334 
335  for (uint32 i = 0; i < numCoefficients; ++i)
336  {
337  integralCofficients[i + 1] = static_cast<double>(coefficients_[i]) / static_cast<double>(i + 1);
338  }
339 
340  return Poly1d<double>(integralCofficients);
341  }
342 
343  //============================================================================
344  // 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:
372  std::string str() const
373  {
374  const uint32 numCoeffients = static_cast<uint32>(coefficients_.size());
375 
376  std::string repr = "Poly1d<";
377  uint32 power = 0;
378  for (auto& coefficient : coefficients_)
379  {
381  {
382  if (coefficient == 0)
383  {
384  ++power;
385  continue;
386  }
387  }
388  else
389  {
390  if (utils::essentiallyEqual(coefficient, static_cast<dtype>(0.0)))
391  {
392  ++power;
393  continue;
394  }
395  }
396 
397  repr += utils::num2str(coefficient);
398 
399  if (power > 1)
400  {
401  repr += "x^" + utils::num2str(power);
402  }
403  else if (power == 1)
404  {
405  repr += "x";
406  }
407 
408  ++power;
409 
410  if (power < numCoeffients)
411  {
412  repr += " + ";
413  }
414  }
415 
416  return repr + ">";
417  }
418 
419  //============================================================================
420  // Method Description:
428  dtype operator()(dtype inValue) const noexcept
429  {
430  dtype polyValue = 0;
431  uint8 power = 0;
432  for (auto& coefficient : coefficients_)
433  {
434  polyValue += coefficient * utils::power(inValue, power++);
435  }
436 
437  return polyValue;
438  }
439 
440  //============================================================================
441  // Method Description:
449  Poly1d<dtype> operator+(const Poly1d<dtype>& inOtherPoly) const
450  {
451  return Poly1d<dtype>(*this) += inOtherPoly;
452  }
453 
454  //============================================================================
455  // Method Description:
463  Poly1d<dtype>& operator+=(const Poly1d<dtype>& inOtherPoly)
464  {
465  if (this->coefficients_.size() < inOtherPoly.coefficients_.size())
466  {
467  for (size_t i = 0; i < coefficients_.size(); ++i)
468  {
469  coefficients_[i] += inOtherPoly.coefficients_[i];
470  }
471  for (size_t i = coefficients_.size(); i < inOtherPoly.coefficients_.size(); ++i)
472  {
473  coefficients_.push_back(inOtherPoly.coefficients_[i]);
474  }
475  }
476  else
477  {
478  for (size_t i = 0; i < inOtherPoly.coefficients_.size(); ++i)
479  {
480  coefficients_[i] += inOtherPoly.coefficients_[i];
481  }
482  }
483 
484  return *this;
485  }
486 
487  //============================================================================
488  // Method Description:
496  Poly1d<dtype> operator-(const Poly1d<dtype>& inOtherPoly) const
497  {
498  return Poly1d<dtype>(*this) -= inOtherPoly;
499  }
500 
501  //============================================================================
502  // Method Description:
510  Poly1d<dtype>& operator-=(const Poly1d<dtype>& inOtherPoly)
511  {
512  if (this->coefficients_.size() < inOtherPoly.coefficients_.size())
513  {
514  for (size_t i = 0; i < coefficients_.size(); ++i)
515  {
516  coefficients_[i] -= inOtherPoly.coefficients_[i];
517  }
518  for (size_t i = coefficients_.size(); i < inOtherPoly.coefficients_.size(); ++i)
519  {
520  coefficients_.push_back(-inOtherPoly.coefficients_[i]);
521  }
522  }
523  else
524  {
525  for (size_t i = 0; i < inOtherPoly.coefficients_.size(); ++i)
526  {
527  coefficients_[i] -= inOtherPoly.coefficients_[i];
528  }
529  }
530 
531  return *this;
532  }
533 
534  //============================================================================
535  // Method Description:
543  Poly1d<dtype> operator*(const Poly1d<dtype>& inOtherPoly) const
544  {
545  return Poly1d<dtype>(*this) *= inOtherPoly;
546  }
547 
548  //============================================================================
549  // Method Description:
557  Poly1d<dtype>& operator*=(const Poly1d<dtype>& inOtherPoly)
558  {
559  const uint32 finalCoefficientsSize = order() + inOtherPoly.order() + 1;
560  std::vector<dtype> coeffsA(finalCoefficientsSize, 0);
561  std::vector<dtype> coeffsB(finalCoefficientsSize, 0);
562  stl_algorithms::copy(coefficients_.begin(), coefficients_.end(), coeffsA.begin());
563  stl_algorithms::copy(inOtherPoly.coefficients_.cbegin(), inOtherPoly.coefficients_.cend(), coeffsB.begin());
564 
565  // now multiply out the coefficients
566  std::vector<dtype> finalCoefficients(finalCoefficientsSize, 0);
567  for (uint32 i = 0; i < finalCoefficientsSize; ++i)
568  {
569  for (uint32 k = 0; k <= i; ++k)
570  {
571  finalCoefficients[i] += coeffsA[k] * coeffsB[i - k];
572  }
573  }
574 
575  this->coefficients_ = finalCoefficients;
576  return *this;
577  }
578 
579  //============================================================================
580  // Method Description:
589  {
590  return Poly1d(*this) ^= inPower;
591  }
592 
593  //============================================================================
594  // Method Description:
603  {
604  if (inPower == 0)
605  {
606  coefficients_.clear();
607  coefficients_.push_back(1);
608  return *this;
609  }
610  else if (inPower == 1)
611  {
612  return *this;
613  }
614 
615  auto thisPoly(*this);
616  for (uint32 power = 1; power < inPower; ++power)
617  {
618  *this *= thisPoly;
619  }
620 
621  return *this;
622  }
623 
624  //============================================================================
625  // Method Description:
633  friend std::ostream& operator<<(std::ostream& inOStream, const Poly1d<dtype>& inPoly)
634  {
635  inOStream << inPoly.str() << std::endl;
636  return inOStream;
637  }
638 
639  private:
640  std::vector<dtype> coefficients_{};
641  };
642  }
643 }
StaticAsserts.hpp
nc::polynomial::Poly1d::operator*=
Poly1d< dtype > & operator*=(const Poly1d< dtype > &inOtherPoly)
Definition: Poly1d.hpp:557
nc::linalg::inv
NdArray< double > inv(const NdArray< dtype > &inArray)
Definition: inv.hpp:55
nc::NdArray::shape
Shape shape() const noexcept
Definition: NdArrayCore.hpp:4296
nc::polynomial::Poly1d::area
double area(double a, double b) const
Definition: Poly1d.hpp:112
Error.hpp
nc::NdArray::dot
NdArray< dtype > dot(const NdArray< dtype > &inOtherArray) const
Definition: NdArrayCore.hpp:2635
nc::polynomial::Poly1d::integ
Poly1d< double > integ() const
Definition: Poly1d.hpp:324
nc::utils::essentiallyEqual
bool essentiallyEqual(dtype inValue1, dtype inValue2) noexcept
Definition: essentiallyEqual.hpp:53
nc::uint8
std::uint8_t uint8
Definition: Types.hpp:43
nc::utils::num2str
std::string num2str(dtype inNumber)
Definition: num2str.hpp:47
nc::NdArray::issquare
bool issquare() const noexcept
Definition: NdArrayCore.hpp:2936
nc::polynomial::Poly1d::operator^
Poly1d< dtype > operator^(uint32 inPower) const
Definition: Poly1d.hpp:588
nc::NdArray::transpose
NdArray< dtype > transpose() const
Definition: NdArrayCore.hpp:4591
nc::polynomial::Poly1d::operator+
Poly1d< dtype > operator+(const Poly1d< dtype > &inOtherPoly) const
Definition: Poly1d.hpp:449
nc::polynomial::Poly1d
Definition: Poly1d.hpp:58
nc::NdArray< dtype >
nc::polynomial::Poly1d::operator^=
Poly1d< dtype > & operator^=(uint32 inPower)
Definition: Poly1d.hpp:602
nc::uint32
std::uint32_t uint32
Definition: Types.hpp:41
nc::polynomial::Poly1d::order
uint32 order() const noexcept
Definition: Poly1d.hpp:350
nc::NdArray::numCols
uint32 numCols() const noexcept
Definition: NdArrayCore.hpp:3415
nc::polynomial::Poly1d::fit
static Poly1d< double > fit(const NdArray< dtype > &xValues, const NdArray< dtype > &yValues, uint8 polyOrder)
Definition: Poly1d.hpp:193
nc::polynomial::Poly1d::operator()
dtype operator()(dtype inValue) const noexcept
Definition: Poly1d.hpp:428
NdArray.hpp
num2str.hpp
nc::polynomial::Poly1d::operator-=
Poly1d< dtype > & operator-=(const Poly1d< dtype > &inOtherPoly)
Definition: Poly1d.hpp:510
nc::polynomial::Poly1d::fit
static Poly1d< double > fit(const NdArray< dtype > &xValues, const NdArray< dtype > &yValues, const NdArray< dtype > &weights, uint8 polyOrder)
Definition: Poly1d.hpp:248
nc::polynomial::Poly1d::Poly1d
Poly1d(const NdArray< dtype > &inValues, bool isRoots=false)
Definition: Poly1d.hpp:79
nc::polynomial::Poly1d::deriv
Poly1d< dtype > deriv() const
Definition: Poly1d.hpp:162
nc::stl_algorithms::copy
OutputIt copy(InputIt first, InputIt last, OutputIt destination) noexcept
Definition: StlAlgorithms.hpp:96
nc::NdArray::isflat
bool isflat() const noexcept
Definition: NdArrayCore.hpp:2874
nc::NdArray::size
size_type size() const noexcept
Definition: NdArrayCore.hpp:4310
power.hpp
nc::polynomial::Poly1d::operator<<
friend std::ostream & operator<<(std::ostream &inOStream, const Poly1d< dtype > &inPoly)
Definition: Poly1d.hpp:633
nc
Definition: Coordinate.hpp:45
nc::swap
void swap(NdArray< dtype > &inArray1, NdArray< dtype > &inArray2) noexcept
Definition: swap.hpp:43
nc::power
constexpr dtype power(dtype inValue, uint8 inExponent) noexcept
Definition: Functions/power.hpp:54
nc::polynomial::Poly1d::str
std::string str() const
Definition: Poly1d.hpp:372
nc::DtypeInfo
Holds info about the dtype.
Definition: DtypeInfo.hpp:41
DtypeInfo.hpp
essentiallyEqual.hpp
nc::polynomial::Poly1d::coefficients
NdArray< dtype > coefficients() const
Definition: Poly1d.hpp:151
THROW_INVALID_ARGUMENT_ERROR
#define THROW_INVALID_ARGUMENT_ERROR(msg)
Definition: Error.hpp:37
nc::polynomial::Poly1d::operator*
Poly1d< dtype > operator*(const Poly1d< dtype > &inOtherPoly) const
Definition: Poly1d.hpp:543
nc::utils::power
dtype power(dtype inValue, uint8 inPower) noexcept
Definition: Utils/power.hpp:49
nc::stl_algorithms::transform
OutputIt transform(InputIt first, InputIt last, OutputIt destination, UnaryOperation unaryFunction) noexcept
Definition: StlAlgorithms.hpp:703
nc::polynomial::Poly1d::Poly1d
Poly1d()=default
StlAlgorithms.hpp
Types.hpp
diagflat.hpp
nc::polynomial::Poly1d::operator-
Poly1d< dtype > operator-(const Poly1d< dtype > &inOtherPoly) const
Definition: Poly1d.hpp:496
nc::polynomial::Poly1d::astype
Poly1d< dtypeOut > astype() const
Definition: Poly1d.hpp:130
nc::polynomial::Poly1d::operator+=
Poly1d< dtype > & operator+=(const Poly1d< dtype > &inOtherPoly)
Definition: Poly1d.hpp:463
nc::polynomial::Poly1d::print
void print() const
Definition: Poly1d.hpp:360
inv.hpp