Ttykuu

http://en.wikipedia.org/wiki/Monotone_cubic_interpolation

We have implemented it using the formula from Wikipedia :

public class MonotoneCubicSplineInterpolation
 
 public static double[] Calc(double[] xs, double[] ys, double[] x_interp)
 
 var length = xs.Length;

 // Deal with length issues
 if (length != ys.Length)
 
 IPDevLoggerWrapper.Error("Need an equal count of xs and ys");
 throw new Exception("Need an equal count of xs and ys");
 
 if (length == 0)
 
 return null;
 
 if (length == 1)
 
 return new double[] ys[0];
 


 // Get consecutive differences and slopes
 var delta = new double[length - 1];
 var m = new double[length];

 for (int i = 0; i < length - 1; i++)
 
 delta[i] = (ys[i + 1] - ys[i]) / (xs[i + 1] - xs[i]);
 if (i > 0)
 
 m[i] = (delta[i - 1] + delta[i]) / 2;
 
 
 var toFix = new List<int>();
 for (int i = 1; i < length - 1; i++)
 (delta[i] < 0 && delta[i - 1] > 0))
 
 toFix.Add(i);
 
 
 foreach (var val in toFix)
 
 m[val] = 0;
 

 m[0] = delta[0];
 m[length - 1] = delta[length - 2];

 toFix.Clear();
 for (int i = 0; i < length - 1; i++)
 
 if (delta[i] == 0)
 
 toFix.Add(i);
 
 
 foreach (var val in toFix)
 
 m[val] = 0;
 m[val + 1] = 0;
 

 var alpha = new double[length - 1];
 var beta = new double[length - 1];
 var dist = new double[length - 1];
 var tau = new double[length - 1];
 for (int i = 0; i < length - 1; i++)
 
 alpha[i] = m[i] / delta[i];
 beta[i] = m[i + 1] / delta[i];
 dist[i] = Math.Pow(alpha[i], 2) + Math.Pow(beta[i], 2);
 tau[i] = 3/Math.Sqrt(dist[i]);
 

 toFix.Clear();
 for (int i = 0; i < length - 1; i++)
 
 if (dist[i] > 9)
 
 toFix.Add(i);
 
 

 foreach (var val in toFix)
 
 m[val] = tau[val] * alpha[val] * delta[val];
 m[val + 1] = tau[val] * beta[val] * delta[val];
 

 var y_interp = new double[x_interp.Length];
 int ind = 0;

 foreach (var x in x_interp)
 
 int i;
 for (i = xs.Length - 2; i >= 0; --i)
 
 if (xs[i] <= x)
 
 break;
 
 
 var h = xs[i + 1] - xs[i];
 var t = (x - xs[i])/h;
 var t2 = Math.Pow(t, 2);
 var t3 = Math.Pow(t, 3);
 var h00 = 2*t3 - 3*t2 + 1;
 var h10 = t3 - 2*t2 + t;
 var h01 = -2*t3 + 3*t2;
 var h11 = t3 - t2;
 y_interp[ind++] = h00*ys[i] + h10*h*m[i] + h01*ys[i + 1] + h11*h*m[i + 1];

 continue;
 

 return y_interp;
 
 

Please comment about style, correctness and complexity.

Don't throw Exception

throw new Exception("Need an equal count of xs and ys");

It forces client code to catch any subclass of Exception. In this case I would throw an ArgumentException.

The continue at the end of the last loop is redundant.

Here you're using double.Equals

if (delta[i] == 0)

From MSDN

The Equals method should be used with caution, because two apparently
equivalent values can be unequal due to the differing precision of the
two values.

That link covers two techniques for dealing with this.

As far as I can tell, toFix can be removed. For example,

var toFix = new List<int>();
for (int i = 1; i < length - 1; i++)

foreach (var val in toFix)

 m[val] = 0;

Can be rewritten as

for (int i = 1; i < length - 1; i++)

I've just had to write this algorithm for some work i have been trying to complete. Couple with the wiki page on monotone cubic interpolation and your example i have optimised (slightly) and re-written the example. This can be done O(n) but chosen to use two for loops and the interpolation loop - can be optimised further.

public class MonotoneCubicSpline

 private readonly double[] _x;
 private readonly double[] _y;
 private readonly int _n;
 private readonly double[] _m;

 public MonotoneCubicSpline(IList<Point> points)
 
 _x = points.OrderBy(p => p.X).Select(p => p.X).ToArray();
 _y = points.OrderBy(p => p.X).Select(p => p.Y).ToArray();
 _n = points.Count;
 _m = new double[_n];

 var delta = new double[_n - 1];

 for (var k = 0; k < _n - 1; k++)
 
 // 1. Find the slope of the section.
 delta[k] = (_y[k + 1] - _y[k]) / (_x[k + 1] - _x[k]);

 // 2. Obtain an initial estimate of the control point slope.
 if (k == 0)
 
 _m[k] = delta[k];
 
 else
 (delta[k] < 0 && delta[k - 1] > 0))
 
 _m[k] = 0;
 
 else
 
 _m[k] = (delta[k - 1] + delta[k]) / 2;
 
 
 

 _m[0] = delta[0];
 _m[_n - 1] = delta[_n - 2];

 // Improve control point slope.
 for (var k = 0; k < _n -1; k++)
 
 

 public double Interpolate(double x)
 
 if (x >= _x[_n - 1])
 
 return (float) _y[_n - 1];
 

 if (x <= _x[0])
 
 return (float) _y[0];
 

 var i = 0;
 for (i = _x.Length - 1; i >= 0; i--)
 
 if (_x[i] <= x)
 
 break;
 
 

 var h = _x[i + 1] - _x[i];
 var t = (x - _x[i]) / h;

 var t2 = Math.Pow(t, 2);
 var t3 = Math.Pow(t, 3);

 var h00 = 2 * t3 - 3 * t2 + 1;
 var h10 = t3 - 2 * t2 + t;
 var h01 = -2 * t3 + 3 * t2;
 var h11 = t3 - t2;

 var v = h00 * _y[i] + h10 * h * _m[i] + h01 * _y[i + 1] + h11 * h * _m[i + 1];

 return v;