Polynomial Regression (Quadratic Fit) in C++

PolynomialRegression.h

#pragma once
/**
 * https://gist.github.com/chrisengelsma/108f7ab0a746323beaaf7d6634cf4add
 * 
 * PURPOSE:
 *
 *  Polynomial Regression aims to fit a non-linear relationship to a set of
 *  points. It approximates this by solving a series of linear equations using 
 *  a least-squares approach.
 *
 *  We can model the expected value y as an nth degree polynomial, yielding
 *  the general polynomial regression model:
 *
 *  y = a0 + a1 * x + a2 * x^2 + ... + an * x^n
 *
 * LICENSE:
 *
 * MIT License
 * 
 * Copyright (c) 2020 Chris Engelsma
 * 
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 * 
 * The above copyright notice and this permission notice shall be included in all
 * copies or substantial portions of the Software.
 * 
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 *
 * @author Chris Engelsma
 */
#include <vector>
#include <stdlib.h>
#include <cmath>
#include <stdexcept>

namespace PolynomialRegression
{

  using size_s = long long int;

  template <typename TYPE>
  bool fit_impl(const TYPE *x, const TYPE *y, size_s N, const int &order, std::vector<TYPE> &coeffs)
  {
    int n = order;
    int np1 = n + 1;
    int np2 = n + 2;
    int tnp1 = 2 * n + 1;
    TYPE tmp;

    // X = vector that stores values of sigma(xi^2n)
    std::vector<TYPE> X(tnp1);
    for (int i = 0; i < tnp1; ++i)
    {
      X[i] = 0;
      for (int j = 0; j < N; ++j)
        X[i] += (TYPE)pow(x[j], i);
    }

    // a = vector to store final coefficients.
    std::vector<TYPE> a(np1, 0);

    // B = normal augmented matrix that stores the equations.
    std::vector<std::vector<TYPE>> B(np1, std::vector<TYPE>(np2, 0));

    for (int i = 0; i <= n; ++i)
      for (int j = 0; j <= n; ++j)
        B[i][j] = X[i + j];

    // Y = vector to store values of sigma(xi^n * yi)
    std::vector<TYPE> Y(np1);
    for (int i = 0; i < np1; ++i)
    {
      Y[i] = (TYPE)0;
      for (int j = 0; j < N; ++j)
      {
        Y[i] += (TYPE)pow(x[j], i) * y[j];
      }
    }

    // Load values of Y as last column of B
    for (int i = 0; i <= n; ++i)
      B[i][np1] = Y[i];

    n += 1;
    int nm1 = n - 1;

    // Pivotisation of the B matrix.
    for (int i = 0; i < n; ++i)
      for (int k = i + 1; k < n; ++k)
        if (B[i][i] < B[k][i])
          for (int j = 0; j <= n; ++j)
          {
            tmp = B[i][j];
            B[i][j] = B[k][j];
            B[k][j] = tmp;
          }

    // Performs the Gaussian elimination.
    // (1) Make all elements below the pivot equals to zero
    //     or eliminate the variable.
    for (int i = 0; i < nm1; ++i)
      for (int k = i + 1; k < n; ++k)
      {
        TYPE t = B[k][i] / B[i][i];
        for (int j = 0; j <= n; ++j)
          B[k][j] -= t * B[i][j]; // (1)
      }

    // Back substitution.
    // (1) Set the variable as the rhs of last equation
    // (2) Subtract all lhs values except the target coefficient.
    // (3) Divide rhs by coefficient of variable being calculated.
    for (int i = nm1; i >= 0; --i)
    {
      a[i] = B[i][n]; // (1)
      for (int j = 0; j < n; ++j)
        if (j != i)
          a[i] -= B[i][j] * a[j]; // (2)
      a[i] /= B[i][i];            // (3)
    }

    coeffs.resize(a.size());
    for (size_t i = 0; i < a.size(); ++i)
      coeffs[i] = a[i];

    return true;
  }

  inline bool fit(const double *x, const double *y, size_s N, const int &order,
                  std::vector<double> &coeffs)
  {
    return fit_impl(x, y, N, order, coeffs);
  }

  inline void eval(const double *x, double *y, size_s count, std::vector<double> &coef)
  {
    const int np = coef.size();

    for (size_s i = 0; i < count; ++i)
    {
      const double x_i = x[i];

      if (2 == np)
      {
        y[i] = coef[0] + x_i * coef[1];
      }
      else if (3 == np)
      {
        y[i] = coef[0] + x_i * coef[1] + x_i * x_i * coef[2];
      }
      else

      {
        double v = coef[0];
        for (int p = 1; p < np; ++p)
        {
          double x_p = x_i;

          // Assume p is small and therefore `pow` is poor
          for (int _ = 1; _ < p; ++_)
          {
            x_p *= x_i;
          }

          v += x_p * coef[p];
        }
        y[i] = v;
      }
    }
  }
};