diff --git a/changelog.md b/changelog.md index 2265348..5762232 100644 --- a/changelog.md +++ b/changelog.md @@ -1,3 +1,37 @@ +# v0.8.7 +The 1D interpolation methods now support extrapolation using these methods: +- `Constant`: Set all points outside the range of the interpolator to `extrapValue`. +- `Edge`: Use the value of the left/right edge. +- `Linear`: Uses linear extrapolation using the two points closest to the edge. +- `Native` (default): Uses the native method of the interpolator to extrapolate. For Linear1D it will be a linear extrapolation, and for Cubic and Hermite splines it will be cubic extrapolation. +- `Error`: Raises an `ValueError` if `x` is outside the range. + +These are passed in as an argument to `eval` and `derivEval`: +```nim +let valEdge = interp.eval(x, Edge) +let valConstant = interp.eval(x, Constant, NaN) +``` + +# v0.8.6 +- `levmarq` now accepts `yError`. +- `paramUncertainties` allows you to calculate the uncertainties of fitted parameters. +- `chi2` test added + +# v0.8.5 +Fix rbf bug. + +# v0.8.4 +With radial basis function interpolation, `numericalnim` finally gets an interpolation method which works on scattered data in arbitrary dimensions! + +Basic usage: +``` +let interp = newRbf(points, values) +let result = interp.eval(evalPoints) +``` + +# v0.8.1-v0.8.3 +CI-related bug fixes. + # v0.8.0 - 09.05.2022 ## Optimization has joined the chat Multi-variate optimization and differentiation has been introduced. diff --git a/src/numericalnim/interpolate.nim b/src/numericalnim/interpolate.nim index bfb634f..61ce78d 100644 --- a/src/numericalnim/interpolate.nim +++ b/src/numericalnim/interpolate.nim @@ -1,4 +1,4 @@ -import strformat, math, tables, algorithm +import std/[strformat, math, tables, algorithm] import arraymancer, cdt/[dt, vectors, edges, types] import ./utils, @@ -17,6 +17,12 @@ export rbf ## - Hermite spline (recommended): cubic spline that works with many types of values. Accepts derivatives if available. ## - Cubic spline: cubic spline that only works with `float`s. ## - Linear spline: Linear spline that works with many types of values. +## +## ### Extrapolation +## Extrapolation is supported for all 1D interpolators by passing the type of extrapolation as an argument of `eval`. +## The default is to use the interpolator's native method to extrapolate. This means that Linear does linear extrapolation, +## while Hermite and Cubic performs cubic extrapolation. Other options are using a constant value, using the values of the edges, +## linear extrapolation and raising an error if the x-value is outside the domain. runnableExamples: import numericalnim, std/[math, sequtils] @@ -28,6 +34,8 @@ runnableExamples: let val = interp.eval(0.314) + let valExtrapolate = interp.eval(1.1, Edge) + ## ## 2D interpolation ## - Bicubic: Works on gridded data. ## - Bilinear: Works on gridded data. @@ -70,6 +78,7 @@ runnableExamples: type InterpolatorType*[T] = ref object X*: seq[float] + Y*: seq[T] coeffs_f*: Tensor[float] coeffs_T*: Tensor[T] high*: int @@ -77,6 +86,8 @@ type eval_handler*: EvalHandler[T] deriveval_handler*: EvalHandler[T] EvalHandler*[T] = proc(self: InterpolatorType[T], x: float): T {.nimcall.} + ExtrapolateKind* = enum + Constant, Edge, Linear, Native, Error Eval2DHandler*[T] = proc (self: Interpolator2DType[T], x, y: float): T {.nimcall.} Interpolator2DType*[T] = ref object z*, xGrad*, yGrad*, xyGrad*: Tensor[T] # 2D tensor @@ -101,32 +112,7 @@ type proc findInterval*(list: openArray[float], x: float): int {.inline.} = - let highIndex = list.high - if x < list[0] or list[highIndex] < x: - raise newException(ValueError, &"x = {x} isn't in the interval [{list[0]}, {list[highIndex]}]") - result = max(0, lowerbound(list, x) - 1) - #[ - ## Finds the index of the element to the left of x in list using binary search. list must be ordered. - let highIndex = list.high - if x < list[0] or list[highIndex] < x: - raise newException(ValueError, &"x = {x} isn't in the interval [{list[0]}, {list[highIndex]}]") - var upper = highIndex - var lower = 0 - var n = floorDiv(upper + lower, 2) - # find interval using binary search - for i in 0 .. highIndex: - if x < list[n]: - # x is below current interval - upper = n - n = floorDiv(upper + lower, 2) - continue - if list[n+1] < x: - # x is above current interval - lower = n + 1 - n = floorDiv(upper + lower, 2) - continue - # x is in the interval - return n]# + result = clamp(lowerbound(list, x) - 1, 0, list.high-1) # CubicSpline @@ -190,7 +176,7 @@ proc newCubicSpline*[T: SomeFloat](X: openArray[float], Y: openArray[ ## Returns a cubic spline. let (xSorted, ySorted) = sortAndTrimDataset(@X, @Y) let coeffs = constructCubicSpline(xSorted, ySorted) - result = InterpolatorType[T](X: xSorted, coeffs_T: coeffs, high: xSorted.high, + result = InterpolatorType[T](X: xSorted, Y: ySorted, coeffs_T: coeffs, high: xSorted.high, len: xSorted.len, eval_handler: eval_cubicspline, deriveval_handler: derivEval_cubicspline) @@ -249,7 +235,7 @@ proc newHermiteSpline*[T](X: openArray[float], Y, dY: openArray[ var coeffs = newTensorUninit[T](ySorted.len, 2) for i in 0 .. ySorted.high: coeffs[i, _] = @[ySorted[i], dySorted[i]].toTensor.reshape(1, 2) - result = InterpolatorType[T](X: xSorted, coeffs_T: coeffs, high: xSorted.high, + result = InterpolatorType[T](X: xSorted, Y: ySorted, coeffs_T: coeffs, high: xSorted.high, len: xSorted.len, eval_handler: eval_hermitespline, deriveval_handler: derivEval_hermitespline) @@ -269,7 +255,7 @@ proc newHermiteSpline*[T](X: openArray[float], Y: openArray[ var coeffs = newTensorUninit[T](Y.len, 2) for i in 0 .. ySorted.high: coeffs[i, _] = @[ySorted[i], dySorted[i]].toTensor.reshape(1, 2) - result = InterpolatorType[T](X: xSorted, coeffs_T: coeffs, high: xSorted.high, + result = InterpolatorType[T](X: xSorted, Y: ySorted, coeffs_T: coeffs, high: xSorted.high, len: xSorted.len, eval_handler: eval_hermitespline, deriveval_handler: derivEval_hermitespline) @@ -301,25 +287,113 @@ proc newLinear1D*[T](X: openArray[float], Y: openArray[ var coeffs = newTensor[T](ySorted.len, 1) for i in 0 .. ySorted.high: coeffs[i, 0] = ySorted[i] - result = InterpolatorType[T](X: xSorted, coeffs_T: coeffs, high: xSorted.high, + result = InterpolatorType[T](X: xSorted, Y: ySorted, coeffs_T: coeffs, high: xSorted.high, len: xSorted.len, eval_handler: eval_linear1d, deriveval_handler: derivEval_linear1d) # General Spline stuff -template eval*[T](interpolator: InterpolatorType[T], x: float): untyped = +type Missing = object +proc missing(): Missing = discard + +proc eval*[T; U](interpolator: InterpolatorType[T], x: float, extrap: ExtrapolateKind = Native, extrapValue: U = missing()): T = ## Evaluates an interpolator. - interpolator.eval_handler(interpolator, x) + ## - `x`: The point to evaluate the interpolator at. + ## - `extrap`: The extrapolation method to use. Available options are: + ## - `Constant`: Set all points outside the range of the interpolator to `extrapValue`. + ## - `Edge`: Use the value of the left/right edge. + ## - `Linear`: Uses linear extrapolation using the two points closest to the edge. + ## - `Native` (default): Uses the native method of the interpolator to extrapolate. For Linear1D it will be a linear extrapolation, and for Cubic and Hermite splines it will be cubic extrapolation. + ## - `Error`: Raises an `ValueError` if `x` is outside the range. + ## - `extrapValue`: The extrapolation value to use when `extrap = Constant`. + ## + ## > Beware: `Native` extrapolation for the cubic splines can very quickly diverge if the extrapolation is too far away from the interpolation points. + when U is Missing: + assert extrap != Constant, "When using `extrap = Constant`, a value `extrapValue` must be supplied!" + else: + when not T is U: + {.error: &"Type of `extrap` ({U}) is not the same as the type of the interpolator ({T})!".} + + let xLeft = x < interpolator.X[0] + let xRight = x > interpolator.X[^1] + if xLeft or xRight: + case extrap + of Constant: + when U is Missing: + discard + else: + return extrapValue + of Native: + discard + of Edge: + return + if xLeft: interpolator.Y[0] + else: interpolator.Y[^1] + of Linear: + let (xs, ys) = + if xLeft: + ((interpolator.X[0], interpolator.X[1]), (interpolator.Y[0], interpolator.Y[1])) + else: + ((interpolator.X[^2], interpolator.X[^1]), (interpolator.Y[^2], interpolator.Y[^1])) + let k = (x - xs[0]) / (xs[1] - xs[0]) + return ys[0] + k * (ys[1] - ys[0]) + of Error: + raise newException(ValueError, &"x = {x} isn't in the interval [{interpolator.X[0]}, {interpolator.X[^1]}]") + + result = interpolator.eval_handler(interpolator, x) + -template derivEval*[T](interpolator: InterpolatorType[T], x: float): untyped = +proc derivEval*[T; U](interpolator: InterpolatorType[T], x: float, extrap: ExtrapolateKind = Native, extrapValue: U = missing()): T = ## Evaluates the derivative of an interpolator. - interpolator.deriveval_handler(interpolator, x) - -proc eval*[T](spline: InterpolatorType[T], x: openArray[float]): seq[T] = + ## - `x`: The point to evaluate the interpolator at. + ## - `extrap`: The extrapolation method to use. Available options are: + ## - `Constant`: Set all points outside the range of the interpolator to `extrapValue`. + ## - `Edge`: Use the value of the left/right edge. + ## - `Linear`: Uses linear extrapolation using the two points closest to the edge. + ## - `Native` (default): Uses the native method of the interpolator to extrapolate. For Linear1D it will be a linear extrapolation, and for Cubic and Hermite splines it will be cubic extrapolation. + ## - `Error`: Raises an `ValueError` if `x` is outside the range. + ## - `extrapValue`: The extrapolation value to use when `extrap = Constant`. + ## + ## > Beware: `Native` extrapolation for the cubic splines can very quickly diverge if the extrapolation is too far away from the interpolation points. + when U is Missing: + assert extrap != Constant, "When using `extrap = Constant`, a value `extrapValue` must be supplied!" + else: + when not T is U: + {.error: &"Type of `extrap` ({U}) is not the same as the type of the interpolator ({T})!".} + + let xLeft = x < interpolator.X[0] + let xRight = x > interpolator.X[^1] + if xLeft or xRight: + case extrap + of Constant: + when U is Missing: + discard + else: + return extrapValue + of Native: + discard + of Edge: + return + if xLeft: interpolator.deriveval_handler(interpolator, interpolator.X[0]) + else: interpolator.deriveval_handler(interpolator, interpolator.X[^1]) + of Linear: + let (xs, ys) = + if xLeft: + ((interpolator.X[0], interpolator.X[1]), (interpolator.deriveval_handler(interpolator, interpolator.X[0]), interpolator.deriveval_handler(interpolator, interpolator.X[1]))) + else: + ((interpolator.X[^2], interpolator.X[^1]), (interpolator.deriveval_handler(interpolator, interpolator.X[^2]), interpolator.deriveval_handler(interpolator, interpolator.X[^1]))) + let k = (x - xs[0]) / (xs[1] - xs[0]) + return ys[0] + k * (ys[1] - ys[0]) + of Error: + raise newException(ValueError, &"x = {x} isn't in the interval [{interpolator.X[0]}, {interpolator.X[^1]}]") + + result = interpolator.deriveval_handler(interpolator, x) + +proc eval*[T; U](spline: InterpolatorType[T], x: openArray[float], extrap: ExtrapolateKind = Native, extrapValue: U = missing()): seq[T] = ## Evaluates an interpolator at all points in `x`. result = newSeq[T](x.len) for i, xi in x: - result[i] = eval(spline, xi) + result[i] = eval(spline, xi, extrap, extrapValue) proc toProc*[T](spline: InterpolatorType[T]): InterpolatorProc[T] = ## Returns a proc to evaluate the interpolator. @@ -329,11 +403,11 @@ converter toNumContextProc*[T](spline: InterpolatorType[T]): NumContextProc[T, f ## Convert interpolator to `NumContextProc`. result = proc(x: float, ctx: NumContext[T, float]): T = eval(spline, x) -proc derivEval*[T](spline: InterpolatorType[T], x: openArray[float]): seq[T] = +proc derivEval*[T; U](spline: InterpolatorType[T], x: openArray[float], extrap: ExtrapolateKind = Native, extrapValue: U = missing()): seq[T] = ## Evaluates the derivative of an interpolator at all points in `x`. result = newSeq[T](x.len) for i, xi in x: - result[i] = derivEval(spline, xi) + result[i] = derivEval(spline, xi, extrap, extrapValue) proc toDerivProc*[T](spline: InterpolatorType[T]): InterpolatorProc[T] = ## Returns a proc to evaluate the derivative of the interpolator. diff --git a/tests/test_interpolate.nim b/tests/test_interpolate.nim index 8f10a9a..576efe5 100644 --- a/tests/test_interpolate.nim +++ b/tests/test_interpolate.nim @@ -184,6 +184,75 @@ test "linear1DSpline derivEval, seq input": for i, val in res: check isClose(val, cos(tTest[i]), 5e-2) +let tOutside = @[-0.1, 10.1] +let yOutside = tOutside.map(f) +let dyOutside = tOutside.map(df) + +test "Extrapolate Constant": + for interp in [linear1DSpline, cubicSpline, hermiteSpline]: + let res = interp.eval(tOutside, Constant, NaN) + for x in res: + check classify(x) == fcNan + +test "Extrapolate Edge": + for interp in [linear1DSpline, cubicSpline, hermiteSpline]: + let res = interp.eval(tOutside, Edge) + check res == @[y[0], y[^1]] + +test "Extrapolate Error": + for interp in [linear1DSpline, cubicSpline, hermiteSpline]: + expect ValueError: + let res = interp.eval(tOutside, Error) + +test "Extrapolate Linear": + let exact = linear1DSpline.eval(tOutside, Native) + for interp in [linear1DSpline, cubicSpline, hermiteSpline]: + let res = interp.eval(tOutside, ExtrapolateKind.Linear) + check res == exact + +test "Extrapolate Native (Cubic)": + let res = cubicSpline.eval(tOutside, Native) + for (y1, y2) in zip(res, yOutside): + check abs(y1 - y2) < 1e-2 + +test "Extrapolate Native (Hermit)": + let res = hermiteSpline.eval(tOutside, Native) + for (y1, y2) in zip(res, yOutside): + check abs(y1 - y2) < 1.1e-2 + +test "Extrapolate Deriv Constant": + for interp in [linear1DSpline, cubicSpline, hermiteSpline]: + let res = interp.derivEval(tOutside, Constant, NaN) + for x in res: + check classify(x) == fcNan + +test "Extrapolate Deriv Edge": + for interp in [linear1DSpline, cubicSpline, hermiteSpline]: + let res = interp.derivEval(tOutside, Edge) + check abs(res[0] - dy[0]) < 1e-2 + check abs(res[1] - dy[^1]) < 4e-2 + +test "Extrapolate Deriv Error": + for interp in [linear1DSpline, cubicSpline, hermiteSpline]: + expect ValueError: + let res = interp.derivEval(tOutside, Error) + +test "Extrapolate Deriv Linear": + let exact = linear1DSpline.derivEval(tOutside, Native) + for interp in [linear1DSpline, cubicSpline, hermiteSpline]: + let res = interp.derivEval(tOutside, ExtrapolateKind.Linear) + check abs(res[0] - exact[0]) < 1e-2 + check abs(res[1] - exact[1]) < 5e-2 + +test "Extrapolate Deriv Native (Cubic)": + let res = cubicSpline.derivEval(tOutside, Native) + check abs(res[0] - dyOutside[0]) < 1e-2 + check abs(res[1] - dyOutside[^1]) < 1.05e-1 + +test "Extrapolate Deriv Native (Hermit)": + let res = hermiteSpline.derivEval(tOutside, Native) + check abs(res[0] - dyOutside[0]) < 3e-2 + check abs(res[1] - dyOutside[^1]) < 2e-1 # 2D Interpolation