Introducing Dependency Injection into Interpolation Code

Greetings!

In one of my previous blogs I have touched upon the idea of programming to interfaces. I did not expect that the blog would be getting so many hits as it does. As far as writing maintainable loosely coupled code, I have a long way to go to become good at. Thus I’ve decided to investigate the most useful design patters to learn and introduce into the hobby and professional coding practices. And I am not going to start small…

One of the seriously misunderstood software design patterns is dependency injection. Many people like writing about it on their sites and blogs, but not many people do a decent job at describing it correctly. I am not claiming that I will now fix this for my blog readers. Instead, I will tel you that I am reading a fantastic book about it: Dependency Injection in .NET by Mark Seemann. So, if you are looking for a good material on this topic, please check out this recommendation.

My last post was on linear and cubic spline interpolation. The implementation can be improved, since I’ve hard-coded at least one inflexible dependency: outputting the results to the console screen. Anyone who wishes to use my example code without using the console screen (e.g. writing to a file instead) would need to rewrite at least twelve lines of code!  We can definitely improve upon this by decoupling the console writer and substituting a generic class that is ‘injected‘ with the user-specified writer through a constructor (aka constructor injection).

To begin with, I made a small change to the Tridiagonal class and removed console-based output. The Solve method now throws ArgumentException whenever passed matrix is not square or something else went wrong. I then made a similar change to the extensions. After all, it is not a good idea to tie extensions to a specific output class like Console. The LinearInterpolation and CubicSplineInterpolation classes now have private constructors. The basic checks on user input is now done by the initializing methods, which could have been defined as properties. Also, the interpolating classes do not display any messages when something goes wrong. Instead, they throw appropriate exceptions, which the calling method will deal with. Finally, I have added a new UserOutput namespace which defines an interface and a very generic class that will work with any concrete implementation of Write method declared by the interface:

namespace UserOutput
{
#region GENERIC_OUTPUT
internal interface IWriter
{
void Write(string message);
}

internal class Output
{

internal Output(IWriter writer)
{
if (writer == null)
{
throw new ArgumentNullException("Null writer interface.");
}

this.writer = writer;
}

internal void Display(string message)
{
writer.Write(message);
}
}
#endregion
}


The concrete implementation is another class defined as:

internal class ConsoleWriter : UserOutput.IWriter
{
public void Write(string message)
{
Console.WriteLine(message);
}
}


The Main method creates an instance of the IWriter interface and sets it to the ConsoleWriter. The linear and cubic spline interpolations are called through the IInterpolate interface which they both define. We no longer need two separate instances of each class. The full code can be downloaded here as a pdf file.

Linear and Cubic Spline Interpolation

In this post on numerical methods I will share with you the theoretical background and the implementation of the two types of interpolations: linear and natural cubic spline. These interpolations are often used within the financial industry. For example, discount rates for some maturity for which there is no point on an yield curve is often interpolated using cubic spline interpolation. Another example can be in risk management, like a position’s Value-at-Risk calculation obtained from a partial or a 2-dimensional revaluation grid.

Linear Interpolation

Let’s begin with the simplest method – linear interpolation. The idea is that we are given a set of numerical points and function values at these points. The task is to use the given set and approximate the function’s value at some different points. That is, given $x_i$ where $i=0,...,n-1$ our task is to estimate $f(x)$ for $x_0 \geq x\leq x_n$. Of course we may require to go outside of the range of our set of points, which would require extrapolation (or projection outside the known function values).

Almost all interpolation techniques are based around the concept of function approximation. Mathematically, the backbone is simple:

$f(x)=y_{i}+\frac{y_{i+1}-y_i}{x_{i+1}-x_i} (x-x_i)$

The above expression tells us that the value of the function we are approximating at point x will be around $y_i$ and that $x_{i} \geq x\leq x_{i+1}$. We can re-write this expression as follows:

$s_{k}(x)=a_{k}+b_{k}(x-x_{k})$

where we have substituted $y_{i}$ for $a_{k}$ and $\frac{y_{i+1}-y_i}{x_{i+1}-x_i}$ for $b_{k}$. Our linear interpolation is now taking a form of linear regression around $a_{k}$.

Linear interpolation is the most basic type of interpolations. It works remarkably well for smooth functions with sufficient number of points. However, because it is such a basic method, interpolating more complex functions requires a little bit more work.

Cubic Spline Interpolation

In what follows I am going to borrow a few formulations from a very good book: Guide to Scientific Computing by P.R.Turner. The idea of a spline interpolation is to extend the single polynomial of linear interpolation to higher degrees. But unlike with piece-wise polynomials, the higher degree polynomials constructed by the cubic spline are different in each interval $[x_i, x_{i+1}]$. The cubic spline can be defined as:

$s_{k}(x)=a_{k}+b_{k}(x-x_{k})+c_{k}(x-x_{k})^{2}+d_{k}(x-x_{k})^{3}$

If you compare the above expression with the linear interpolation, you will see that we have added two terms to the rhs. It is these extra terms that should allow us achieve greater precision. There are certain requirements placed on s. These are: it must be continuous, and its first and second derivatives must exist. We only need to determine $b_{k}, c_{k}$ and $d_{k}$ coefficients, since $a_{k}=f_{k}$. Take $h_{k}=x_{k+1}-x_{k}$. The continuity condition (see P.R.Turner pg. 104) gives the following equalities:

$d_{k}=\frac{c_{k+1}-c_{k}}{3h_{k}}$
$b_{k}=\frac{f_{k+1}-f{k}}{h_{k}}-\frac{h_k}{3}(c_{k+1}+2c_k)$

Thus, if we had the values for $c_{k}$, we would be able to determine all the coefficients. After some re-arranging, it is possible to obtain:

$h_{k}c_{k}+2(h_{k}+h_{k+1})c_{k+1}+h_{k+1}c_{k+2}=3(\frac{f_{k+2}-f_{k+1}}{h_{k+1}}-\frac{f_{k+1}-f_{k}}{h_k})$

for $k=0,1,...,n-3$. This is a tridiagonal system of linear equations, which can be solved in a number of ways.

Let’s now look at how we can implement the linear and cubic spline interpolation in C#.

Implementing Linear and Cubic Spline Interpolation in C#

The code is broken into five regions. Tridiagonal Matrix region defines a Tridiagonal class to solve a system of linear equations. The Extensions regions defines a few extensions to allows for matrix manipulations. The Foundation region is where the parent Interpolation class is defined. Both, Linear and Cubic Spline Interpolations classes inherit from it.

using System;
using System.Collections;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace TridiagonalMatrix
{
#region TRIDIAGONAL_MATRIX

// solves a Tridiagonal Matrix using Thomas algorithm
public sealed class Tridiagonal
{
public double[] Solve(double[,] matrix, double[] d)
{
int rows = matrix.GetLength(0);
int cols = matrix.GetLength(1);
int len = d.Length;
// this must be a square matrix
if (rows == cols && rows == len)
{

double[] b = new double[rows];
double[] a = new double[rows];
double[] c = new double[rows];

// decompose the matrix into its tri-diagonal components
for (int i = 0; i < rows; i++)
{
for (int j=0; j<cols; j++)
{
if (i == j)
b[i] = matrix[i, j];
else if (i == (j - 1))
c[i] = matrix[i, j];
else if (i == (j + 1))
a[i] = matrix[i, j];
}
}
try
{
c[0] = c[0] / b[0];
d[0] = d[0] / b[0];

for (int i = 1; i < len - 1; i++)
{
c[i] = c[i] / (b[i] - a[i] * c[i - 1]);
d[i] = (d[i] - a[i] * d[i - 1]) / (b[i] - a[i] * c[i - 1]);
}
d[len - 1] = (d[len - 1] - a[len - 1] * d[len - 2]) / (b[len - 1] - a[len - 1] * c[len - 2]);

// back-substitution step
for (int i = (len - 1); i-- > 0; )
{
d[i] = d[i] - c[i] * d[i + 1];
}

return d;
}
catch(DivideByZeroException)
{
Console.WriteLine("Division by zero was attempted. Most likely reason is that ");
Console.WriteLine("the tridiagonal matrix condition ||(b*i)||>||(a*i)+(c*i)|| is not satisified.");
return null;
}
}
else
{
Console.WriteLine("Error: the matrix must be square. The vector must be the same size as the matrix length.");
return null;
}
}
}

#endregion
}

namespace BasicInterpolation
{
#region EXTENSIONS
public static class Extensions
{
// returns a list sorted in ascending order
public static List<T> SortedList<T>(this T[] array)
{
List<T> l = array.ToList();
l.Sort();
return l;

}

// returns a difference between consecutive elements of an array
public static double[] Diff(this double[] array)
{
int len = array.Length - 1;
double[] diffsArray = new double[len];
for (int i = 1; i <= len;i++)
{
diffsArray[i - 1] = array[i] - array[i - 1];
}
return diffsArray;
}

// scaled an array by another array of doubles
public static double[] Scale(this double[] array, double[] scalor, bool mult=true)
{
int len = array.Length;
double[] scaledArray = new double[len];

if (mult)
{
for (int i = 0; i < len; i++)
{
scaledArray[i] = array[i] * scalor[i];
}
}
else
{
for (int i = 0; i < len; i++)
{
if (scalor[i] != 0)
{
scaledArray[i] = array[i] / scalor[i];
}
else
{
// basic fix to prevent division by zero
scalor[i] = 0.00001;
scaledArray[i] = array[i] / scalor[i];

}
}
}

return scaledArray;
}

public static double[,] Diag(this double[,] matrix, double[] diagVals)
{
int rows = matrix.GetLength(0);
int cols = matrix.GetLength(1);
// the matrix has to be scare
if (rows==cols)
{
double[,] diagMatrix = new double[rows, cols];
int k = 0;
for (int i = 0; i < rows; i++)
for (int j = 0; j < cols; j++ )
{
if(i==j)
{
diagMatrix[i, j] = diagVals[k];
k++;
}
else
{
diagMatrix[i, j] = 0;
}
}
return diagMatrix;
}
else
{
Console.WriteLine("Diag should be used on square matrix only.");
return null;
}

}
}

#endregion
#region FOUNDATION
internal interface IInterpolate
{
double? Interpolate( double p);
}

internal  abstract class Interpolation: IInterpolate
{

public Interpolation(double[] _x, double[] _y)
{
int xLength = _x.Length;
if (xLength == _y.Length && xLength > 1 && _x.Distinct().Count() == xLength)
{
x = _x;
y = _y;
}
}

// cubic spline relies on the abscissa values to be sorted
public Interpolation(double[] _x, double[] _y, bool checkSorted=true)
{
int xLength = _x.Length;
if (checkSorted)
{
if (xLength == _y.Length && xLength > 1 && _x.Distinct().Count() == xLength && Enumerable.SequenceEqual(_x.SortedList(),_x.ToList()))
{
x = _x;
y = _y;
}
}
else
{
if (xLength == _y.Length && xLength > 1 && _x.Distinct().Count() == xLength)
{
x = _x;
y = _y;
}
}
}

public double[] X
{
get
{
return x;
}
}

public double[] Y
{
get
{
return y;
}
}

public abstract double? Interpolate(double p);

private double[] x;
private double[] y;
}

#endregion

#region LINEAR
internal sealed class LinearInterpolation: Interpolation
{

public LinearInterpolation(double[] _x, double [] _y):base(_x,_y)
{
len = X.Length;
if (len > 1)
{
Console.WriteLine("Successfully set abscissa and ordinate.");
baseset = true;
// make a copy of X as a list for later use
lX = X.ToList();
}
else
Console.WriteLine("Ensure x and y are the same length and have at least 2 elements. All x values must be unique.");
}

public override double? Interpolate(double p)
{
if (baseset)
{
double? result = null;
double Rx;

try
{
// point p may be outside abscissa's range
// if it is, we return null
Rx = X.First(s => s >= p);
}
catch(ArgumentNullException)
{
return null;
}

// at this stage we know that Rx contains a valid value
// find the index of the value close to the point required to be interpolated for
int i = lX.IndexOf(Rx);

// provide for index not found and lower and upper tabulated bounds
if (i==-1)
return null;

if (i== len-1 && X[i]==p)
return Y[len - 1];

if (i == 0)
return Y[0];

// linearly interpolate between two adjacent points
double h = (X[i] - X[i - 1]);
double A = (X[i] - p) / h;
double B = (p - X[i - 1]) / h;

result = Y[i - 1] * A + Y[i] * B;

return result;

}
else
{
return null;
}
}

private bool baseset = false;
private int len;
private List<double> lX;
}

#endregion

#region  CUBIC_SPLINE
internal sealed class CubicSplineInterpolation:Interpolation
{
public CubicSplineInterpolation(double[] _x, double[] _y): base(_x, _y, true)
{
len = X.Length;
if (len > 1)
{

Console.WriteLine("Successfully set abscissa and ordinate.");
baseset = true;
// make a copy of X as a list for later use
lX = X.ToList();
}
else
Console.WriteLine("Ensure x and y are the same length and have at least 2 elements. All x values must be unique. X-values must be in ascending order.");
}

public override double? Interpolate(double p)
{
if (baseset)
{
double? result = 0;
int N = len - 1;

double[] h = X.Diff();
double[] D = Y.Diff().Scale(h, false);
double[] s = Enumerable.Repeat(3.0, D.Length).ToArray();
double[] dD3 = D.Diff().Scale(s);
double[] a = Y;

// generate tridiagonal system
double[,] H = new double[N-1, N-1];
double[] diagVals = new double[N-1];
for (int i=1; i<N;i++)
{
diagVals[i-1] = 2 * (h[i-1] + h[i]);
}

H = H.Diag(diagVals);

// H can be null if non-square matrix is passed
if (H != null)
{
for (int i = 0; i < N - 2; i++)
{
H[i, i + 1] = h[i + 1];
H[i + 1, i] = h[i + 1];
}

double[] c = new double[N + 2];
c = Enumerable.Repeat(0.0, N + 1).ToArray();

// solve tridiagonal matrix
TridiagonalMatrix.Tridiagonal le = new TridiagonalMatrix.Tridiagonal();
double[] solution = le.Solve(H, dD3);

for (int i = 1; i < N;i++ )
{
c[i] = solution[i - 1];
}

double[] b = new double[N];
double[] d = new double[N];
for (int i = 0; i < N; i++ )
{
b[i] = D[i] - (h[i] * (c[i + 1] + 2.0 * c[i])) / 3.0;
d[i] = (c[i + 1] - c[i]) / (3.0 * h[i]);
}

double Rx;

try
{
// point p may be outside abscissa's range
// if it is, we return null
Rx = X.First(m => m >= p);
}
catch
{
return null;
}

// at this stage we know that Rx contains a valid value
// find the index of the value close to the point required to be interpolated for
int iRx = lX.IndexOf(Rx);

if (iRx == -1)
return null;

if (iRx == len - 1 && X[iRx] == p)
return Y[len - 1];

if (iRx == 0)
return Y[0];

iRx = lX.IndexOf(Rx)-1;
Rx = p - X[iRx];
result = a[iRx] + Rx * (b[iRx] + Rx * (c[iRx] + Rx * d[iRx]));

return result;
}
else
{
return null;
}

}
else
return null;
}

private bool baseset = false;
private int len;
private List<double> lX;
}
#endregion
class EntryPoint
{
static void Main(string[] args)
{
// code relies on abscissa values to be sorted
// there is a check for this condition, but no fix
// f(x) = 1/(1+x^2)*sin(x)
double[] ex = { -5, -3, -2, -1, 0, 1, 2, 3, 5};
double[] ey = { 0.036882, -0.01411, -0.18186, -0.42074, 0, 0.420735, 0.181859, 0.014112, -0.03688};
double[] p = {-5,-4.5,-4,-3.5,-3.2,-3,-2.8,-2.1,-1.5,-1,-0.8,0,0.65,1,1.3,1.9,2.6,3.1,3.5,4,4.5,5 };

LinearInterpolation LI = new LinearInterpolation(ex, ey);
CubicSplineInterpolation CS = new CubicSplineInterpolation(ex, ey);

Console.WriteLine("Linear Interpolation:");

foreach (double pp in p)
{
Console.Write(LI.Interpolate(pp).ToString()+", ");

}
Console.WriteLine("Cubic Spline Interpolation:");
foreach (double pp in p)
{

Console.WriteLine(CS.Interpolate(pp).ToString());
}

Note that I wrote a console application, and this is the reason for throwing as few exceptions as possible. Like with all the code on my blog – please feel free to take it and modify as you wish. It is now time to compare performance of linear and cubic spline interpolations on a test function: $f(x)=\frac{1}{1+x^{2}}\sin x$. Using the above code we obtain the following results: