#ifndef STIM_FINITE_DIFFERENCE_COEFFICIENTS
#define STIM_FINITE_DIFFERENCE_COEFFICIENTS

#include <vector>
#include <algorithm>

/// Calculate the coefficients used for finite differences on arbitrary grids. This code is based on:
///	Fornberg, "Generation of Finite Difference Formulas on Arbitrarily Spaced Grids", Mathematics of Computation, 51(184) 699-706, 1988

/// @param Z is the location where the approximation is to be accurate
/// @param X is a vector of coordinates for grid points
/// @param M is the highest derivative to be computed
/// This function returns a matrix C[a][b] where a is the derivative and b is the coefficient for x = b

template <typename T>
std::vector< std::vector<T> > diff_coefficient_matrix(T Z, std::vector<T> X, size_t M){
	size_t n = X.size();								//calculate the number of data points available
	std::vector< std::vector<T> > C(M+1);							//create an array of m arrays
	for(size_t m = 0; m < M+1; m++)						//for each derivative
		C[m].resize(n);									//allocate space for n coefficients
	
	T c1, c2, c3, c4, c5;
	c1 = 1;												//initialize the matrix
	c4 = X[0] - Z;
	C[0][0] = 1;
	for(size_t i = 2; i <= n; i++){
		size_t mn = (std::min)(i, M+1);
		c2 = 1;
		c5 = c4;
		c4 = X[i-1] - Z;
		for(size_t j = 1; j <= i-1; j++){
			c3 = X[i-1] - X[j-1];
			c2 = c2 * c3;
			if(j == i-1){
				for(size_t k = 2; k <= mn; k++){
					C[k-1][i-1] = c1 * ((k-1) * C[k-2][i-2] - c5 * C[k-1][i-2]) / c2;
				}
				C[0][i-1] = -c1 * c5 * C[0][i-2]/c2;
			}
			T c_p0 = C[0][j-1];
			T c_p1 = C[1][j-1];
			for(size_t k = 2; k <= mn; k++){
				c_p1 = C[k-1][j-1];
				C[k-1][j-1] = (c4 * c_p1 - (k-1) * c_p0) / c3;
				c_p0 = c_p1;
			}
			C[0][j-1] = c4 * C[0][j-1] / c3;
		}
		c1 = c2;
	}
	return C;
}

/// Returns the coefficients used to estimate the nth derivative at x0 given values at points in x

/// @param x0 is the point where the derivative f'(x0) is to be estimated
/// @param x is the set of points where the values of f(x) are known
/// @param n specifies the derivative, ex. 2 estimates the second derivative f''(x0)
template <typename T>
std::vector<T> diff_coefficients(T x0, std::vector<T> x, size_t n){
	std::vector< std::vector<T> > C = diff_coefficient_matrix(x0, x, n);		//calculate all derivative coefficients
	return C.back();
}



#endif