SupportValues.h

Go to the documentation of this file.

 
#pragma once
 
#include "../../Matrix.h"
#include "../lin_alg/gaussSeidel.h"
 
 
class PolynomialBase
{
public:
  PolynomialBase([[maybe_unused]] const Matrix<double>& X, [[maybe_unused]] const Matrix<double>& Y) { }
 
  virtual Matrix<double> Evaluate(const Matrix<double>&) const = 0;
 
  virtual std::string Function() const = 0;
 
  friend std::ostream& operator<<(std::ostream& ostr, const PolynomialBase& base) {
    ostr << base.Function() << std::flush;
    return ostr;
  }
};
 
class MonomBase : public PolynomialBase
{
public:
  Matrix<double> A;
 
  MonomBase(const Matrix<double>& X, const Matrix<double>& Y)
    : PolynomialBase(X, Y) {
    auto x = X;
    // ensure required orientation
    if(x.columns() > x.rows()) { x = x.Transpose(); }
    Matrix<double> out(0, x.rows(), x.rows());
    for(size_t i = 0; i < out.rows(); ++i) {
      for(size_t j = 0; j < out.rows(); ++j) { out(i, j) = pow(x(i, 0), j); }
    }
    if(out.Determinant() == 0.0) {
      std::cerr << "Error computing Vandermonde-Matrix! Most likely passed equal support values" << std::endl;
      return;
    }
    A = gaussSeidel(out, Y);
  }
  virtual Matrix<double> Evaluate(const Matrix<double>& X) const {
    Matrix<double> Y(A(0, 0), X.rows(), 1);
    for(size_t i = 0; i < A.rows(); ++i) {
      Y += A(i, 0) * X.Apply([i](const double& in) { return pow(in, i); });
    }
    return Y;
  };
  virtual std::string Function() const {
    std::string out;
    for(size_t row = 0; row < A.rows(); ++row) {
      if(row != 0) { out += " + "; }
      out += format("%d", int(A(row, 0)));
      if(row > 0) { out += format("x^%d", row); }
    }
    return out;
  }
};
 
class LagrangeBase : public PolynomialBase
{
  Matrix<double> x;
  Matrix<double> y;
 
public:
  LagrangeBase(const Matrix<double>& X, const Matrix<double>& Y)
    : PolynomialBase(X, Y)
    , x(X)
    , y(Y) { }
 
  [[nodiscard]] double GetCoefficient(const double& xk, const size_t i) const {
    double out = 1.0;
    for(size_t j = 0; j < x.rows(); ++j) {
      if(i == j) continue;
 
      out *= (xk - x(j, 0)) / (x(i, 0) - x(j, 0));
    }
    return out;
  }
 
  virtual Matrix<double> Evaluate(const Matrix<double>& in) const {
    Matrix<double> out(0, in.rows(), 1);
    for(size_t k = 0; k < in.rows(); ++k) {
      for(size_t i = 0; i < y.rows(); ++i) { out(k, 0) += y(i, 0) * GetCoefficient(in(k, 0), i); }
    }
    return out;
  }
  [[nodiscard]] std::string buildLx(size_t i) const {
    std::string Lx;
 
    bool needsOperator = false;
    for(size_t j = 0; j < x.rows(); ++j) {
      if(i == j) {
        needsOperator = false;
        continue;
      }
 
      // we need to exit, if one element of the product is 0
      if(x(i, 0) - x(j, 0) == 0) return "";
 
      if(needsOperator) { Lx += " * "; }
      if(x(i, 0) - x(j, 0) < 0) { Lx += "-"; }
 
      // Dividend
      if(x(j, 0) != 0) {
        if(x(j, 0) < 0) {
          Lx += format("(x + %.3f)", std::abs(x(j, 0)));
        } else {
          Lx += format("(x - %.3f)", x(j, 0));
        }
      } else {
        Lx += "x";
      }
 
      // Divisor
      Lx += format("/%.3f", std::abs(x(i, 0) - x(j, 0)));
      needsOperator = true;
    }
    return Lx;
  }
  virtual std::string Function() const {
    std::string out;
    for(size_t i = 0; i < y.rows(); ++i) {
      if(y(i, 0) == 0) continue;
 
      if(i > 0) { out += " + "; }
      out += format("%d * (", int(y(i, 0)));
      out += buildLx(i);
      out += ")";
    }
    return out;
  }
};
 
class NewtonBase : public PolynomialBase
{
  Matrix<double> X;
 
public:
  Matrix<double> b;
 
  NewtonBase(const Matrix<double>& x, const Matrix<double>& y)
    : PolynomialBase(x, y)
    , X(x) {
    auto m = x.rows();
    Matrix<double> F(0, m, m);
    F.SetColumn(0, y);
    for(size_t j = 1; j < m; ++j) {
      for(size_t i = j; i < m; ++i) { F(i, j) = (F(i, j - 1) - F(i - 1, j - 1)) / (x(i, 0) - x(i - j, 0)); }
    }
    b = diag_elements(F);
  }
  virtual std::string Function() const {
    std::string out;
    for(size_t i = 0; i < b.rows(); ++i) {
      if(i > 0) { out += " + "; }
      out += format("%.3f", b(i, 0));
      if(i > 0) {
        for(size_t j = 0; j < i; ++j) {
          if(X(j, 0) != 0) {
            if(X(j, 0) > 0) out += format("* (x - %.3f)", X(j, 0));
            else
              out += format("* (x + %.3f)", std::abs(X(j, 0)));
          } else
            out += "* x";
        }
      }
    }
    return out;
  }
 
  double GetCoefficient(size_t index, const double& x) const {
    double out = 1;
    for(size_t i = 0; i < index - 1; ++i) { out *= (x - X(i, 0)); }
    return out;
  }
  virtual Matrix<double> Evaluate(const Matrix<double>& xIn) const {
    Matrix<double> y(b(b.rows() - 1, 0), xIn.rows(), 1);
    for(int j = b.rows() - 2; j >= 0; --j) {
      auto bAct = b(j, 0);
      auto xAct = X(j, 0);
      y         = xIn
          .Apply([xAct](const double& in) { return in - xAct; })  // (xIn - X[j])
          .HadamardMulti(y)                                       // (xIn - X[j]) * y
          .Apply([bAct](const double& in) { return bAct + in; }); // y = b[j] + (xIn - X[j]) * y
    }
    return y;
  }
};