main.cpp

//
//  main.cpp
//  SparseFunction
//
//  Created by Marco Battiato on 3/5/23.
//
//  Demonstration of sparse-grid approximation of multivariate functions.
//
//  A multidimensional basis set is built from a one-dimensional basis characterised
//  by an order. Rather than taking the full Cartesian product of the 1-D bases —
//  which costs O(N^D) coefficients — only those D-dimensional basis functions whose
//  total order (the sum of the orders across all directions) is at most N are kept.
//  This keeps the basis tractable as the dimension grows, at the cost of resolving
//  high-order behaviour in any single direction.
//
//  Monomials (TaylorBasis) are currently the only 1-D basis implemented.
//

#include "SparseFunction.hpp"
#include "BasisFunctionsSets.hpp"
#include <iostream>
#include <iomanip>

using std::cout;
using namespace SparseFunction;

// Number of coefficients a full Cartesian-product basis would require, for comparison.
constexpr long fullTensorProductSize(int order, int nVar) {
    long n = 1;
    for (int i = 0; i < nVar; ++i) n *= (order + 1);
    return n;
}

int main() {

    constexpr int order = 5, numVar = 6;

    // A polynomial in 6 variables containing every monomial of total degree <= 5.
    // TaylorBasis selects plain monomials as the underlying 1-D basis.
    sparseFunction<order, numVar, TaylorBasis> poly;

    // The point of the sparse construction: the coefficient count stays manageable
    // where the full Cartesian product would not.
    cout << "Variables:            " << numVar << "\n"
         << "Maximum total order:  " << order << "\n"
         << "Sparse basis size:    " << poly.numBasisFunction() << " coefficients\n"
         << "Full product would be: " << fullTensorProductSize(order, numVar) << " coefficients\n\n";

    // ---- Fit the polynomial to sampled data -------------------------------------

    // Sample at random points rather than on a regular grid, so the fit is not
    // relying on any special structure in the sample locations.
    constexpr int nDataPoints = 500;
    using Array = Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic>;

    const Array x = Array::Random(nDataPoints, numVar);

    // Target values drawn from a known polynomial:  3.1*x0 + 2*x1 + 3*x1*x2 + 1.2
    // Since this target lies inside the basis, a correct fit should recover these
    // coefficients almost exactly and leave a near-zero residual.
    const Array y = 3.1 * x.col(0) + 2.0 * x.col(1) + 3.0 * x.col(1) * x.col(2) + 1.2;

    // fitData solves the least-squares problem for the coefficients and returns the
    // normalised mean squared residual.
    const double residual = poly.fitData(x, y);

    cout << "Normalised residual:  " << std::scientific << std::setprecision(3)
         << residual << "\n\n";

    // Fitted coefficients, listed against the order vector of their monomial.
    cout << "Fitted coefficients:\n";
    poly.printCoeff();
    cout << "\n";

    // ---- Evaluate ----------------------------------------------------------------

    // The fitted polynomial can now be evaluated at arbitrary points. Compare the
    // first few predictions against the values they were fitted to.
    cout << std::fixed << std::setprecision(4)
         << "First 8 fitted vs. target values:\n"
         << "  fitted: " << poly(x).head(8).transpose() << "\n"
         << "  target: " << y.head(8).transpose() << "\n";

    return 0;
}