Introduction
Microsoft Excel has a lot of useful functionality and is widely used for data analysis. However it does not offer any built-in functions for data interpolation. This article describes how to implement interpolation routines using C# and use Excel-DNA to call them from Excel.
Background
Interpolation means estimating a value of the function based on its tabulated values. There are many different interpolation methods suitable for particular modeling needs. The purpose of the article is to show how interpolation routines can be implemented in C# and exposed to Excel with a help of Excel-DNA library. Peculiarities of the different interpolation methods are not discussed. The code examples below focus on 1D linear interpolation, but C# code and corresponding add-in support 1D interpolation, 2D interpolation on regular rectangular grid and scattered data in arbitrary number of dimensions.
Using the code
Developing functions for Excel is not easy. There are several wrapper libraries that greatly simplify interaction with Excel. For example Excel-DNA library allows to write Excel functions in C# with almost no boilerplate code https://exceldna.codeplex.com/.
Here is example of C# function exported by Excel-DNA to Excel. The function computes a sum of two numbers. Adding ExcelFunction attribute to the static public method instructs Excel-DNA to export it to Excel.
using ExcelDna.Integration;
public static class Examples
{
[ExcelFunction(Description = "adds two numbers")]
public static double acq_sum(double a, double b)
{
return a + b;
}
}
Assume that we have implementation of linear interpolation in the the class ACQ.Math.Interpolation.LinearInterpolation, which can be used as shown in the code below.
double[] x = new double[] { 0, 1 };
double[] y = new double[] { 0, 1 };
ACQ.Math.Interpolation.InterpolationInterface interpolator = new ACQ.Math.Interpolation.LinearInterpolation(x, y);
double xi = 0.5;
double yi = interpolator.Eval(xi);
Using this class, linear interpolation in Excel can be implemented with just a few lines of code.
using ExcelDna.Integration;
public static class Examples
{
[ExcelFunction(Description = "Linear Interpolation", IsThreadSafe = true)]
public static object acq_linear_interpolation(double xi, double[] x, double[] y, object bounds)
{
if (ExcelDnaUtil.IsInFunctionWizard())
return ExcelError.ExcelErrorRef;
else
{
ACQ.Math.Interpolation.InterpolationInterface interpolator = null;
try
{
interpolator = new ACQ.Math.Interpolation.LinearInterpolation(x, y);
interpolator.Bounds = ExcelHelper.CheckValue(bounds, true);
}
catch (Exception ex)
{
}
if (interpolator != null)
return ExcelHelper.CheckNan(interpolator.Eval(xi));
else
return ExcelError.ExcelErrorNA;
}
}
}
}
Here we used Excel-DNA function to check if function is being called from Excel Wizard ( ExcelDnaUtil.IsInFunctionWizard() ).
Interpolation constructor checks that all of the following conditions are satisfied and throws an exception otherwise.
- Array of nodes (x) should have the same size as array of function values (y)
- Array of nodes (x) should be sorted (from small to large values)
- Array of nodes (x) should not have repeated values (i.e. all nodes should be different, x[i] != x[i+1])
Interplator Eval method does not throw exceptions but returns a NaN value to indicate that something is wrong. Excel does not support NaN or Infinity, therefore we return ExcelError.ExcelErrorNA value in case when interpolation results in error. ExcelError.ExcelErrorNA will appear in Excel as #N/A.
static class ExcelHelper
{
internal static object CheckNan(double value)
{
object result;
if (Double.IsNaN(value) || Double.IsInfinity(value))
{
result = ExcelError.ExcelErrorNA;
}
else
{
result = value;
}
return result;
}
internal static T CheckValue<T>(object value, T defaultValue) where T : struct
{
T result = value is T ? (T)value : defaultValue;
return result;
}
}
Interpolator Eval function returns NaN when specified point xi is located outside tabulated function values (i.e. x> x[n-1] or xi < x[0]) and bounds is false. When bounds is true and xi is outside interpolation range, interpolator return y[0] when xi < x[0] and y[n-1] when xi > x[n-1]. Bounds argument is optional therefore is as declared as object instead of boolean. When bounds argument is not specified in Excel, Excel-DNA will pass the following value ExcelDna.Integration.ExcelMissing, and bounds will be set to default value of true using CheckValue function.
If linear interpolation is all you need then you are done. Unfortunately interpolator is constructed every time this function is called. Constructing linear interpolator is cheap, since it only checks validity of the input. However other interpolation algorithms require a lot more work and therefore it is not efficient to call constructor for each interpolation point. For example natural cubic spline is constructed by solving tridiagonal system of equations. Is it possible to call constructor once, store the interpolator object somewhere and then use it to compute interpolated values? The solution presented below does exactly that.
As far as I know Excel does not offer any built-in way to store temporary objects. Therefore we need to develop out own object cache where we can store objects.
The class shown below based on various ideas from Excel-DNA discussion boards that revolve around use of ExcelAsyncUtil.Observe to manage object lifetime. To the best of my knowledge it was first suggested by Govert van Drimmelen (creator of Excel-DNA). Observe function allows to link the object to the input, so when input changes new object is automatically constructed and old object is deleted. This works seamlessly with Excel automatic calculation.
The HandleStorage
class stores object handles in the dictionary and provides method for creating handles:CreateHandle
, retrieving objects TryGetObject and calling methods of stored objects TryReadObject. Notice the use of reader-writer lock to synchronize the access to the underlying dictionary of handles.
class HandleStorage
{
private ReaderWriterLockSlim m_lock = new ReaderWriterLockSlim();
private Dictionary<string, Handle> m_storage = new Dictionary<string, Handle>();
internal object CreateHandle(string objectType, object[] parameters, Func<string, object[], object> maker)
{
return ExcelAsyncUtil.Observe(objectType, parameters, () =>
{
var value = maker(objectType, parameters);
var handle = new Handle(this, objectType, value);
m_lock.EnterWriteLock();
try
{
m_storage.Add(handle.Name, handle);
}
finally
{
m_lock.ExitWriteLock();
}
return handle;
});
}
internal bool TryGetObject<T>(string name, out T value) where T : class
{
bool found = false;
value = default(T);
m_lock.EnterReadLock();
try
{
Handle handle;
if (m_storage.TryGetValue(name, out handle))
{
if (handle.Value is T)
{
value = handle.Value as T;
found = true;
}
}
}
finally
{
m_lock.ExitReadLock();
}
return found;
}
internal Tuple<bool, TResult> TryReadObject<T, TResult, TArg>(string name, Func<T, TArg, TResult> reader, TArg argument) where T : class
{
bool valid = false;
TResult result = default(TResult);
T obj = default(T);
m_lock.EnterReadLock();
try
{
Handle handle;
if (m_storage.TryGetValue(name, out handle))
{
if (handle.Value is T)
{
obj = handle.Value as T;
if (obj != null)
{
result = reader(obj, argument);
valid = true;
}
}
}
}
finally
{
m_lock.ExitReadLock();
}
return new Tuple<bool, TResult>(valid, result);
}
internal void Remove(Handle handle)
{
object value;
if (TryGetObject(handle.Name, out value))
{
m_lock.EnterWriteLock();
try
{
m_storage.Remove(handle.Name);
}
finally
{
m_lock.ExitWriteLock();
}
var disp = value as IDisposable;
if (disp != null)
{
disp.Dispose();
}
}
}
}
Since we only need one object storage we create the GlobalCache object that internally has static instance of HandleStorage
class GlobalCache
{
private static HandleStorage m_storage = new HandleStorage();
internal static object CreateHandle(string objectType, object[] parameters, Func<string, object[], object> maker)
{
return m_storage.CreateHandle(objectType, parameters, maker);
}
internal static bool TryGetObject<T>(string name, out T value) where T : class
{
return m_storage.TryGetObject<T>(name, out value);
}
}
The final piece of the design that allows us to use object handles in excel is Handle class. This class implements IExcelObservable interface and creates a unique string id for each object that will be passed to Excel. Notice that constructor also takes a reference to the handle storage, so that handle can be removed from storage when handle is disposed.
class Handle : IExcelObservable, IDisposable
{
private static readonly object m_lock = new object();
private static int m_index;
private readonly HandleStorage m_storage;
private IExcelObserver m_observer;
private readonly string m_name;
private readonly object m_value;
public Handle(HandleStorage storage, string objectType, object value)
{
m_storage = storage;
m_value = value;
lock (m_lock)
{
m_name = String.Format("{0}:{1}", objectType, m_index++);
}
}
public IDisposable Subscribe(IExcelObserver observer)
{
m_observer = observer;
m_observer.OnNext(m_name);
return this;
}
public void Dispose()
{
m_storage.Remove(this);
}
public string Name
{
get
{
return m_name;
}
}
public object Value
{
get
{
return m_value;
}
}
}
Now the interpolator objects can be easily created and stored in cache for further use. The string id of object handle is returned to Excel.
public class ExcelInterpolator
{
private static readonly object m_sync = new object();
private static readonly string m_tag = "#acqInterpolator";
private static readonly string m_defaultInterpolator = "Linear";
[ExcelFunction(Description = "Create Interpolator object")]
public static object acq_interpolator_create(
[ExcelArgument(Description = "Array of nodes")] double[] x,
[ExcelArgument(Description = "Array of values")] double[] y,
[ExcelArgument(Description = "linear, quadratic, cubic, hermite, akima, steffen etc")] object method,
[ExcelArgument(Description = "Out of range value: false (num error), true (closest)")] object bounds)
{
if (ExcelDnaUtil.IsInFunctionWizard())
return ExcelError.ExcelErrorRef;
else
{
return ACQ.Excel.Handles.GlobalCache.CreateHandle(m_tag, new object[] { x, y, method, bounds, "acq_interpolator_create" },
(objectType, parameters) =>
{
ACQ.Math.Interpolation.InterpolationInterface interpolator = null;
try
{
string interpolation_method = ExcelHelper.Check(method, m_defaultInterpolator);
interpolator = ACQ.Math.Interpolation.InterpolationFactory.GetInterpolator(interpolation_method, x, y);
interpolator.Bounds = ExcelHelper.CheckValue(bounds, true);
}
catch (Exception ex)
{
}
if (interpolator == null)
return ExcelError.ExcelErrorNull;
else
return interpolator;
});
}
}
}
Interpolation objects can be retrieved from the handle storage and used for interpolation
[ExcelFunction(Description = "Evaluate interpolation at specified point")]
public static object acq_interpolator_eval(
[ExcelArgument(Description = "Interpolator object")] string handle,
[ExcelArgument(Description = "Interpolation point")] double x)
{
ACQ.Math.Interpolation.InterpolationInterface interpolator;
if (ACQ.Excel.Handles.GlobalCache.TryGetObject<ACQ.Math.Interpolation.InterpolationInterface>(handle, out interpolator))
{
if (interpolator != null)
{
return ExcelHelper.CheckNan(interpolator.Eval(x));
}
}
return ExcelError.ExcelErrorRef;
}
Alternatively one can use read method ( TryReadObject ) to compute interpolation in the thread safe way.
[ExcelFunction(Description = "Evaluate interpolation at specified point (thread safe version)", IsThreadSafe = true)]
public static object acq_interpolator_eval_tsafe(
[ExcelArgument(Description = "Interpolator object")] string handle,
[ExcelArgument(Description = "Interpolation point")] double x)
{
Tuple<bool, double> results = ACQ.Excel.Handles.GlobalCache.TryReadObject<ACQ.Math.Interpolation.InterpolationInterface, double, double>(handle, (interpolator, point) =>
{
return interpolator.Eval(point);
}, x);
if (results.Item1)
{
return ExcelHelper.CheckNan(results.Item2);
}
else
{
return ExcelError.ExcelErrorRef;
}
}
Interpolation method is specified as string (not case sensitive), and corresponding interpolation class is found using reflection. Therefore new interpolation methods automatically become available in Excel when interpolation class is added to the corresponding namespace. This is done in two steps. First dictionary of all interpolation types is created in static constructor of InterpolationFactory class. The interpolation type can be easily found from the name of the method provided that interpolation class is located in ACQ.Math.Interpolation namespace and it implements InterpolationInterface. Activator.CreateInstance is used to create an interpolation object.
public class InterpolationFactory
{
private static Dictionary<string, Type> m_interpolation_types = new Dictionary<string, Type>();
static InterpolationFactory()
{
Type base_type = typeof(InterpolationInterface);
Type[] types = GetClassTypes(System.Reflection.Assembly.GetExecutingAssembly(), base_type.Namespace);
foreach(Type t in types)
{
if (!t.IsAbstract && base_type.IsAssignableFrom(t))
{
m_interpolation_types[t.FullName.ToLower()] = t;
}
}
}
private static Type[] GetClassTypes(Assembly assembly, string nameSpace)
{
return assembly.GetTypes().Where(t => t.IsClass && String.Equals(t.Namespace, nameSpace)).ToArray();
}
public static Type GetInterpolationType(string method)
{
string name = String.Format("ACQ.Math.Interpolation.{0}Interpolation", method).ToLower();
Type result;
m_interpolation_types.TryGetValue(name, out result); return result;
}
public static InterpolationInterface GetInterpolator(Type type, params object[] arguments)
{
InterpolationInterface interpolator = Activator.CreateInstance(type, arguments) as InterpolationInterface;
return interpolator;
}
public static InterpolationInterface GetInterpolator(string method, params object[] arguments)
{
InterpolationInterface interpolator = null;
Type interpolator_type = GetInterpolationType(method);
if (interpolator_type != null)
{
interpolator = Activator.CreateInstance(interpolator_type, arguments) as InterpolationInterface;
}
return interpolator;
}
}
Adding new interpolation method
All interpolation methods implement InterpolationInterface
public interface InterpolationInterface
{
double Eval(double x);
double EvalDeriv(double x); bool Bounds { get; set; }
}
There is also an abstract base class ( InterpolationBase ) that implements methods that are common to all interpolation algorithms such as input check and index lookup. Abbreviated version is shown below
public abstract class InterpolationBase : InterpolationInterface
{
protected readonly double[] m_x;
protected readonly double[] m_y;
protected bool m_bounds = true;
public InterpolationBase(double[] x, double[] y, bool bCopyData = true)
{
if (x == null || y == null)
throw new ArgumentNullException("interpolation arrays can not be null");
if (x.Length != y.Length)
throw new ArgumentException("interpolation x and y arrays should have the same length");
if (x.Length < 2)
{
throw new ArgumentException("interpolation array should have at least 2 nodes");
}
for (int i = 0; i < x.Length - 1; i++)
{
if (x[i + 1] <= x[i])
{
throw new ArgumentException("interpolation nodes should be ordered");
}
}
if (bCopyData)
{
m_x = (double[])x.Clone();
m_y = (double[])y.Clone();
}
else
{
m_x = x;
m_y = y;
}
}
public bool Bounds
{
set
{
m_bounds = value;
}
get
{
return m_bounds;
}
}
public abstract double Eval(double x);
public virtual double EvalDeriv(double x)
{
return Double.NaN; }
protected int FindIndex(double x, out double value)
{
int index = 0;
value = Double.NaN;
if (x < m_x[0]) {
if (m_bounds)
{
value = m_y[0];
}
}
else if (x > m_x[m_x.Length - 1]) {
if (m_bounds)
{
value = m_y[m_y.Length - 1];
}
}
else
{
index = Array.BinarySearch<double>(m_x, x);
if (index == 0) {
index = 1; }
else if (index < 0)
{
index = ~index; }
}
return index;
}
protected void Bracket(int index1, out double x0, out double x1, out double y0, out double y1)
{
int index0 = index1 - 1;
x0 = m_x[index0];
x1 = m_x[index1];
y0 = m_y[index0];
y1 = m_y[index1];
}
Here is the full implementation of linear interpolation. The index look up, input check and the range check are done in base class.
public class LinearInterpolation : InterpolationBase
{
public LinearInterpolation(double[] x, double[] y)
: base(x, y)
{ }
public override double Eval(double x)
{
double value;
int index = FindIndex(x, out value);
if(index > 0 )
{
double x0, x1, y0, y1;
Bracket(index, out x0, out x1, out y0, out y1);
double dx = x1 - x0;
double b = (x - x0) / dx;
value = y0 + b * (y1 - y0);
}
return value;
}
}
2D interpolation
2D interpolation on rectangular grid can be easily implemented by doing 1D interpolation in each dimension. The code below shows an example.
public class BiInterpolation<T> : InterpolationBase2D where T : InterpolationInterface
{
public BiInterpolation(double[] x1, double[] x2, double[,] y)
: base(x1, x2, y, true)
{
}
public BiInterpolation(double[] x1, double[] x2, double[,] y, bool copyData)
: base(x1, x2, y, copyData)
{
}
public override double Eval(double x1, double x2)
{
double value = Double.NaN;
int index1, index2;
FindIndex(x1, x2, out index1, out index2);
if (index1 > 0 && index2 > 0)
{
int n1 = m_x1.Length;
int n2 = m_x2.Length;
double[] yt = new double[n1];
double[] y2 = new double[n2];
InterpolationInterface interpolator;
Type interpolator_type = typeof(T);
for(int i=0; i<n2; i++)
{
for (int j = 0; j < n1; j++)
{
yt[j] = m_y[i, j];
}
interpolator = Activator.CreateInstance(interpolator_type, m_x1, yt) as InterpolationInterface;
interpolator.Bounds = false;
y2[i] = interpolator.Eval(x1);
}
interpolator = Activator.CreateInstance(interpolator_type, m_x2, y2) as InterpolationInterface;
return interpolator.Eval(x2);
}
return value;
}
}
Notice that class interpolates over complete grind and makes more assumptions about interpolation method. This is very inefficient because most interpolation algorithms are local and only need function values in the adjacent nodes. But it allows to turn any 1D interpolation method into 2D. The code below shows an example of 2D Steffen interpolation.
public class BisteffenInterpolation : BiInterpolation<SteffenInterpolation>
{
public BisteffenInterpolation(double[] x1, double[] x2, double[,] y)
: base(x1, x2, y, true)
{
}
public BisteffenInterpolation(double[] x1, double[] x2, double[,] y, bool copyData)
: base(x1, x2, y, copyData)
{
}
}
This deficiency can be easily rectified by adding the size of the support to the 2D interpolation class, this will allow interpolation algorithm to tell 2D interpolator how many nodes it needs to perform interpolation. The code below shows how support size is used to adjust interpolation window. Setting support size to zero, tells interpolator to use full grid (for example for natural cubic spline interpolation).
public class BiInterpolation<T> : InterpolationBase2D where T : InterpolationInterface
{
public BiInterpolation(double[] x1, double[] x2, double[,] y)
: base(x1, x2, y, true)
{
}
public BiInterpolation(double[] x1, double[] x2, double[,] y, bool copyData)
: base(x1, x2, y, copyData)
{
}
protected virtual int SupportSize
{
get
{
return 0;
}
}
public override double Eval(double x1, double x2)
{
double value = Double.NaN;
int index1, index2;
FindIndex(x1, x2, out index1, out index2);
if (index1 > 0 && index2 > 0)
{
int sn = this.SupportSize;
int n1 = m_x1.Length;
int n2 = m_x2.Length;
int j0, j1, i0, i1;
if (sn <= 0)
{
j0 = i0 = 0;
j1 = n1;
i1 = n2;
}
else
{
j0 = System.Math.Max(0, index1 - sn);
j1 = System.Math.Min(n1, index1 + sn);
i0 = System.Math.Max(0, index2 - sn);
i1 = System.Math.Min(n2, index2 + sn);
}
double[] x1t = new double[j1 - j0];
double[] yt = new double[j1 - j0];
double[] x2t = new double[i1 - i0];
double[] y2 = new double[i1 - i0];
InterpolationInterface interpolator;
Type interpolator_type = typeof(T);
for (int j = j0; j < j1; j++)
{
x1t[j - j0] = m_x1[j];
}
for (int i = i0; i < i1; i++)
{
for (int j = j0; j < j1; j++)
{
yt[j - j0] = m_y[i, j];
}
interpolator = Activator.CreateInstance(interpolator_type, x1t, yt) as InterpolationInterface;
interpolator.Bounds = false;
y2[i - i0] = interpolator.Eval(x1);
x2t[i - i0] = m_x2[i];
}
interpolator = Activator.CreateInstance(interpolator_type, x2t, y2) as InterpolationInterface;
return interpolator.Eval(x2);
}
return value;
}
}
Interpolation methods
The following interpolation methods are currently implemented in 1D:
- Nearest, Backward and Forward - table lookup
- Linear - linear interpolation
- Quadratic - quadratic spline interpolation
- Cubic - natural cubic spline
- Hermite and HermiteQS - local cubic spline (aka Catmull-Rom spline)
- Steffen and SteffenSoft - monotonic cubic spline
- Akima - Akima spline (cubic spline with special condition for derivatives)
- AkimaPeriodic - Akima spline with periodic boundary conditions
- Multiquadrics - Interpolation based on multiquadrics radial basis function
- ExpTension - exponential tension spline
Interpolation on 2D regular rectangular grid can be easily done by using 1D interpolation in each dimension (note that current implementation is not very efficient because it does not consider locality of 1D interpolators).
- BiLinear - linear in each dimension
- BiCubic - natural cubic spline in each dimension
- BiSteffen, BiSteffenSoft
- BiAkima
- BiHermite, BiHermiteQS
Scattered data interpolation is based on radial basis functions and currently limited to 512 interpolation nodes. The following radial basis functions are implemented (first three are the most common). Optional scale factors can be provided for each dimension.
- Linear
- Cubic
- Multiquadrics
- Gaussian
- Thinplate,
- InverseQuadratic,
- InverseMultiquadric
Here is a complete list of interpolation functions. All functions have prefix acq_ and located in ACQ category.
- acq_interpolator_create - creates interpolator object (returns a handle)
- acq_interpolator_tension_create - creates interpolator object for exponential tension spline
- acq_interpolator_eval - evaluates interpolation at specified point
- acq_interpolator_eval_deriv - evaluates interpolation derivative at specified point
- acq_interpolation - evaluates interpolation at specified point(in-situ without constructing interpolator object)
- acq_interpolator2d_create - creates 2D interpolator object
- acq_interpolator2d_eval - evaluates interpolation at specified point
- acq_interpolation2d - evaluates interpolation at specified point(in-situ without constructing interpolator object)
- acq_interpolator_scattered_create - creates RBF interpolator
- acq_interpolator_scattered_eval - evaluates RBF interpolator at specified point
- acq_interpolator_scattered_eval_x5 - evaluates RBF interpolator at specified point, coordinates are specified individually (up to 5)
Add-in installation
Add-in binary is located ACQ\Distribution folder. There are 32 and 64bit versions. There is also an Excel Demo.xlsx file, and some instructions in primer.interpolation.pdf
- Click the File tab, click Options, and then click the Add-Ins category
- In the Manage box, click Excel Add-ins, and then click Go. The Add-Ins dialog box appears
- Browse to the location of the ACQ32.xll/ACQ64.xll files, and pick the xll file based on Excel bitness.
- Make sure your Excel security settings allow you to run Add-ins
- All ACQ functions have prefix "acq_" and located in ACQ category.
The latest version is available at https://github.com/ratesquant/ACQ
Acknowledgment
This ACQ add-in is based on Excel-DNA https://exceldna.codeplex.com/. Thanks to Govert van Drimmelen (creator of Excel-DNA) for making it happen.
History
2016-05-01: First version
2016-05-02: Added example of interpolation class and references
References
- Hiroshi Akima. 1970. A New Method of Interpolation and Smooth Curve Fitting Based on Local Procedures. J. ACM 17, 4 (October 1970), 589-602
- A Simple Method for Monotonic Interpolation in One Dimension, Steffen, M., Astronomy and Astrophysics, Vol. 239, NO. NOV(II), P. 443, 1990
- P. Rentrop, "An Algorithm for the Computation of the Exponential Tension Spline", Numerische Mathematik 35(1):81-93 · February 1980
