cubic spline vero

cubic_spline.cpp

// to compile run: g++ -std=gnu++11 cubic_spline.cpp -lmgl
#include <iostream>
#include <eigen3/Eigen/Dense>
#include <mgl2/mgl.h>

using namespace Eigen;
using namespace std;

MatrixXd cubicSpline(const VectorXd &T, const VectorXd &Y)
{
	// returns the matrix representing the spline interpolating the data
	// with abscissae T and ordinatae Y. Each column represents the coefficients
	// of the cubic polynomial on a subinterval.
	// Assumes T is sorted, has no repeated elements and T.size() == Y.size().

	int n = T.size() - 1; // T and Y have length n+1

	// TODO: build the spline matrix with polynomials' coefficients

	// calculate h's
	VectorXd h = T.tail(n) - T.head(n);

	MatrixXd A = MatrixXd::Zero(n - 1, n - 1);
	A.diagonal() = (T.segment(2, n - 1) - T.segment(0, n - 1)) / 3;
	A.diagonal(1) = h.segment(1, n - 2) / 6;
	A.diagonal(-1) = h.segment(1, n - 2) / 6;

	VectorXd slope = (Y.tail(n) - Y.head(n)).cwiseQuotient(h);
	VectorXd r = slope.tail(n - 1) - slope.head(n - 1);

	VectorXd sigma(n + 1);
	sigma.segment(1, n - 1) = A.partialPivLu().solve(r);
	sigma(0) = 0;
	sigma(n) = 0;

	MatrixXd spline(4, n);
	spline.row(0) = T.head(n);
	spline.row(1) = slope - h.cwiseProduct(2 * sigma.head(n) + sigma.tail(n)) / 6;
	spline.row(2) = sigma.head(n) / 2;
	spline.row(3) = (sigma.tail(n) - sigma.head(n)).cwiseQuotient(6 * h);

	return spline;
	/*
	*/
	exit(1);
}

VectorXd evalCubicSpline(const MatrixXd &S, const VectorXd &T, const VectorXd &evalT)
{
	// Returns the values of the spline S calculated in the points X.
	// Assumes T is sorted, with no repetetions.

	int n = evalT.size();
	VectorXd out(n);

	for (int i = 0; i < n; i++)
	{
		// TODO: search for evaluate point corresponding piece in T
		// get polynomial for that section
		// evaluate it
		double t = evalT[i];
		int poly = -1;
		for (int j = 0; j < T.size() - 1; j++)
		{
			if (T[j] <= t && T[j + 1] >= t)
			{
				poly = j;
			}
		}

		if (poly == -1)
		{
			cout << "could not find corresponding index!" << endl;
			exit(1);
		}

		out(i) = S(0, poly) + t * S(1, poly) + std::pow(t, 2) * S(2, poly) + std::pow(t, 3) * S(3, poly);
	}

	// TODO: fill out

	return out;
}

int main()
{
	// tests
	VectorXd T(9);
	VectorXd Y(9);
	T << 0, 0.4802, 0.7634, 1, 1.232, 1.407, 1.585, 1.879, 2;
	Y << 0., 0.338, 0.7456, 0, -1.234, 0, 1.62, -2.123, 0;

	int len = 1 << 9;
	VectorXd evalT = VectorXd::LinSpaced(len, T(0), T(T.size() - 1));

	VectorXd evalSpline = evalCubicSpline(cubicSpline(T, Y), T, evalT);

	mglData datx, daty;
	datx.Link(evalT.data(), len);
	daty.Link(evalSpline.data(), len);
	mglGraph gr;
	gr.SetRanges(0, 2, -3, 3);
	gr.Plot(datx, daty, "0");
	gr.WriteFrame("spline.eps");
}