在 Excel 中使用 Excel-DNA 进行插值






4.57/5 (7投票s)
使用 C# 和 Excel-DNA 为 Excel 添加插值函数
引言
Microsoft Excel 拥有许多有用的功能,并广泛用于数据分析。然而,它不提供任何内置的数据插值函数。本文介绍了如何使用 C# 实现插值例程,并使用 Excel-DNA 从 Excel 调用它们。
背景
插值是指根据函数表值估计函数值。有许多不同的插值方法适用于特定的建模需求。本文的目的是展示如何在 C# 中实现插值例程,并通过 Excel-DNA 库将其公开到 Excel。文章不讨论不同插值方法的特点。下面的代码示例侧重于一维线性插值,但 C# 代码和相应的插件支持一维插值、规则矩形网格上的二维插值以及任意维度的散点数据。
使用代码
为 Excel 开发函数并不容易。有几个包装器库极大地简化了与 Excel 的交互。例如,Excel-DNA 库允许使用几乎没有样板代码的 C# 编写 Excel 函数 https://exceldna.codeplex.com/。
这是 Excel-DNA 导出到 Excel 的 C# 函数示例。该函数计算两个数字的总和。将 ExcelFunction
属性添加到静态公共方法会指示 Excel-DNA 将其导出到 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; } }
假设我们在 ACQ.Math.Interpolation.LinearInterpolation
类中实现了线性插值,该类可按以下代码所示使用。
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); //should produce 0.5
使用这个类,只需几行代码就可以在 Excel 中实现线性插值。
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 { //create linear interpolator interpolator = new ACQ.Math.Interpolation.LinearInterpolation(x, y); interpolator.Bounds = ExcelHelper.CheckValue(bounds, true); } catch (Exception ex) { //Log Exception } //interpolate if (interpolator != null) return ExcelHelper.CheckNan(interpolator.Eval(xi)); else return ExcelError.ExcelErrorNA; } } } }
这里我们使用了 Excel-DNA 函数来检查函数是否从 Excel 向导调用 (ExcelDnaUtil.IsInFunctionWizard())。
插值构造函数会检查所有以下条件是否满足,否则会抛出异常。
- 节点数组 (x) 的大小应与函数值数组 (y) 的大小相同。
- 节点数组 (x) 应已排序(从小到大)。
- 节点数组 (x) 不应有重复值(即所有节点都应不同,x[i] != x[i+1])。
Interplator Eval 方法不抛出异常,而是返回 NaN 值以指示有问题。Excel 不支持 NaN 或 Infinity,因此在插值出错的情况下,我们返回 ExcelError.ExcelErrorNA
值。ExcelError.ExcelErrorNA
将在 Excel 中显示为 #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;
}
}
当指定的点 xi 位于表格函数值之外(即 x > x[n-1] 或 xi < x[0])且 bounds 为 false 时,Interpolator Eval 函数返回 NaN。当 bounds 为 true 且 xi 在插值范围之外时,插值器会返回 y[0](当 xi < x[0] 时)和 y[n-1](当 xi > x[n-1] 时)。Bounds 参数是可选的,因此声明为 object 而不是 boolean。当在 Excel 中未指定 bounds 参数时,Excel-DNA 将传递值 ExcelDna.Integration.ExcelMissing
,并且 bounds 将使用 CheckValue
函数设置为默认值 true。
如果您只需要线性插值,那么您已经完成了。不幸的是,插值器每次调用此函数时都会被构造。构造线性插值器成本很低,因为它只检查输入的有效性。然而,其他插值算法需要更多的计算量,因此每次插值点都调用构造函数效率不高。例如,自然三次样条的构造是通过求解三对角线方程组来完成的。是否可以调用一次构造函数,将插值器对象存储在某处,然后使用它来计算插值?下面的解决方案正是这样做的。
据我所知,Excel 没有提供任何内置方法来存储临时对象。因此,我们需要开发自己的对象缓存来存储对象。
下面的类基于 Excel-DNA 讨论版上的各种想法,这些想法围绕使用 ExcelAsyncUtil.Observe
来管理对象生命周期。据我所知,这是由 Govert van Drimmelen(Excel-DNA 的创建者)首次提出的。Observe 函数允许将对象链接到输入,因此当输入更改时,会自动构造新对象并删除旧对象。这与 Excel 的自动计算无缝集成。
HandleStorage
类在字典中存储对象句柄,并提供创建句柄的方法:CreateHandle
,检索对象 TryGetObject
以及调用存储对象的 TryReadObject
方法。注意使用了读写锁来同步对句柄底层字典的访问。
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();
}
}
}
}
由于我们只需要一个对象存储,因此我们创建了 GlobalCache
对象,该对象内部有一个 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);
}
}
允许我们在 Excel 中使用对象句柄的设计的最后一块是 Handle 类。该类实现 IExcelObservable 接口,并为将传递给 Excel 的每个对象创建一个唯一的字符串 ID。请注意,构造函数还引用句柄存储,以便在句柄被处置时可以从存储中删除句柄。
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;
}
}
}
现在可以轻松创建插值对象并将其存储在缓存中以供将来使用。对象句柄的字符串 ID 将返回给 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)
{
//LogError
}
if (interpolator == null)
return ExcelError.ExcelErrorNull;
else
return interpolator;
});
}
}
}
可以从句柄存储中检索插值对象并用于插值。
[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;
}
或者,可以使用 read 方法 (TryReadObject) 以线程安全的方式计算插值。
[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;
}
}
插值方法被指定为字符串(不区分大小写),并且使用反射来查找相应的插值类。因此,当插值类添加到相应的命名空间时,新的插值方法会自动在 Excel 中可用。这分两步完成。首先,在 InterpolationFactory
类的静态构造函数中创建所有插值类型的字典。只要插值类位于 ACQ.Math.Interpolation
命名空间中并且实现了 InterpolationInterface
,就可以轻松地从方法名称中找到插值类型。使用 Activator.CreateInstance 来创建插值对象。
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); //returns null if not found
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;
}
}
添加新的插值方法
所有插值方法都实现 InterpolationInterface
。
public interface InterpolationInterface
{
double Eval(double x);
double EvalDeriv(double x); //compute first derivative
bool Bounds { get; set; }
}
还有一个抽象基类 (InterpolationBase),它实现了所有插值算法共有的方法,例如输入检查和索引查找。下面显示了缩减版本。
public abstract class InterpolationBase : InterpolationInterface
{
protected readonly double[] m_x;
protected readonly double[] m_y;
protected bool m_bounds = true; //this is not needed in constructor and can be switched on or off
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");
}
//check that data is ordered
for (int i = 0; i < x.Length - 1; i++)
{
if (x[i + 1] <= x[i])
{
throw new ArgumentException("interpolation nodes should be ordered");
}
}
//Data is copied to save memory
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; //not all interpolation classes have to provide this method
}
protected int FindIndex(double x, out double value)
{
int index = 0;
value = Double.NaN;
if (x < m_x[0]) //check left boundary
{
if (m_bounds)
{
value = m_y[0];
}
}
else if (x > m_x[m_x.Length - 1]) //check right boundary
{
if (m_bounds)
{
value = m_y[m_y.Length - 1];
}
}
else
{
index = Array.BinarySearch<double>(m_x, x);
if (index == 0) // x = m_x[0]
{
index = 1; //we can always do this because we check in constructor that there are at least two nodes
}
else if (index < 0)
{
index = ~index; //binary search returns compliment of the node larger than a value
}
}
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];
}
这是线性插值的完整实现。索引查找、输入检查和范围检查是在基类中完成的。
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;
}
}
二维插值
二维矩形网格上的二维插值可以通过在每个维度上进行一维插值轻松实现。下面的代码显示了一个示例。
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)
{
//slow reference method
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;
}
}
请注意,该类在完整的网格上进行插值,并对插值方法做出更多假设。这是非常低效的,因为大多数插值算法是局部的,并且只需要相邻节点中的函数值。但它可以将任何一维插值方法转换为二维。下面的代码显示了二维 Steffen 插值的示例。
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)
{
}
}
通过向二维插值类添加支持的大小,可以轻松纠正此缺陷,这将允许插值算法告知二维插值器执行插值所需的节点数量。下面的代码显示了支持大小如何用于调整插值窗口。将支持大小设置为零,告诉插值器使用整个网格(例如,对于自然三次样条插值)。
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;
//determine interpolation window
//use all grind
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;
}
}
插值方法
目前一维中实现的插值方法如下:
- Nearest, Backward 和 Forward - 表查找
- Linear - 线性插值
- Quadratic - 二次样条插值
- Cubic - 自然三次样条
- Hermite 和 HermiteQS - 局部三次样条(也称为 Catmull-Rom 样条)
- Steffen 和 SteffenSoft - 单调三次样条
- Akima - Akima 样条(具有导数特殊条件的样条)
- AkimaPeriodic - Akima 样条,具有周期性边界条件
- Multiquadrics - 基于多二次径向基函数的插值
- ExpTension - 指数张力样条
二维规则矩形网格上的插值可以通过在每个维度上使用一维插值轻松完成(请注意,当前实现效率不高,因为它没有考虑一维插值器的局部性)。
- BiLinear - 每个维度上的线性插值
- BiCubic - 每个维度上的自然三次样条
- BiSteffen, BiSteffenSoft
- BiAkima
- BiHermite, BiHermiteQS
散点数据插值基于径向基函数,目前限制为 512 个插值节点。实现了以下径向基函数(前三个最常见)。可以为每个维度提供可选的比例因子。
- Linear
- Cubic
- Multiquadrics
- Gaussian
- Thinplate,
- InverseQuadratic,
- InverseMultiquadric
以下是插值函数的完整列表。所有函数都带有前缀 acq_,并位于 ACQ 类别中。
- acq_interpolator_create - 创建插值器对象(返回句柄)
- acq_interpolator_tension_create - 为指数张力样条创建插值器对象
- acq_interpolator_eval - 在指定点评估插值
- acq_interpolator_eval_deriv - 在指定点评估插值导数
- acq_interpolation - 在指定点评估插值(就地执行,无需构造插值器对象)
- acq_interpolator2d_create - 创建二维插值器对象
- acq_interpolator2d_eval - 在指定点评估插值
- acq_interpolation2d - 在指定点评估插值(就地执行,无需构造插值器对象)
- acq_interpolator_scattered_create - 创建 RBF 插值器
- acq_interpolator_scattered_eval - 在指定点评估 RBF 插值器
- acq_interpolator_scattered_eval_x5 - 在指定点评估 RBF 插值器,坐标单独指定(最多 5 个)
插件安装
插件二进制文件位于 ACQ\Distribution 文件夹。有 32 位和 64 位版本。还有一个 Excel Demo.xlsx 文件,以及 primer.interpolation.pdf 中的一些说明。
- 单击“文件”选项卡,然后单击“选项”,再单击“加载项”类别。
- 在“管理”框中,单击“Excel 加载项”,然后单击“转到”。此时将显示“加载项”对话框。
- 浏览到 ACQ32.xll/ACQ64.xll 文件的位置,然后根据 Excel 的位选择 xll 文件。
- 确保您的 Excel 安全设置允许您运行加载项。
- 所有 ACQ 函数都带有前缀“acq_”,并位于 ACQ 类别中。
最新版本可在 https://github.com/ratesquant/ACQ 获取。
致谢
此 ACQ 插件基于 Excel-DNA https://exceldna.codeplex.com/。感谢 Govert van Drimmelen(Excel-DNA 的创建者)的辛勤付出。
历史
2016-05-01:第一个版本
2016-05-02:添加了插值类示例和参考文献
参考文献
- 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