From 5034fa3826739e04b804373eb8418e3d09e2609f Mon Sep 17 00:00:00 2001 From: HugoGranstrom <5092565+HugoGranstrom@users.noreply.github.com> Date: Thu, 6 Jul 2023 16:24:43 +0200 Subject: [PATCH 1/8] implement extrapolation for eval --- src/numericalnim/interpolate.nim | 85 +++++++++++++++++++------------- 1 file changed, 50 insertions(+), 35 deletions(-) diff --git a/src/numericalnim/interpolate.nim b/src/numericalnim/interpolate.nim index bfb634f..868d4a9 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, @@ -70,6 +70,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 +78,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 +104,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 +168,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 +227,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 +247,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 +279,62 @@ 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) + # check for extrapolation + 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 = ## Evaluates the derivative of an interpolator. interpolator.deriveval_handler(interpolator, x) -proc eval*[T](spline: InterpolatorType[T], x: openArray[float]): seq[T] = +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. From aebe990aaa2e851eb44273dd874f4228a97ea44e Mon Sep 17 00:00:00 2001 From: HugoGranstrom <5092565+HugoGranstrom@users.noreply.github.com> Date: Thu, 6 Jul 2023 16:24:55 +0200 Subject: [PATCH 2/8] test extrapolation --- tests/test_interpolate.nim | 34 ++++++++++++++++++++++++++++++++++ 1 file changed, 34 insertions(+) diff --git a/tests/test_interpolate.nim b/tests/test_interpolate.nim index 8f10a9a..b70ee9f 100644 --- a/tests/test_interpolate.nim +++ b/tests/test_interpolate.nim @@ -184,6 +184,40 @@ 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) + +test "Extrapolate Constant": + for interp in [linear1DSpline, cubicSpline, hermiteSpline]: + let res = interp.eval(tOutside, Constant, NaN) + for x in res: + check x.isNan() + +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 # 2D Interpolation From ed074e7be2cc17b86f452e57c2250c6025dca96f Mon Sep 17 00:00:00 2001 From: HugoGranstrom <5092565+HugoGranstrom@users.noreply.github.com> Date: Thu, 6 Jul 2023 16:35:04 +0200 Subject: [PATCH 3/8] use old method of checking nan --- tests/test_interpolate.nim | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_interpolate.nim b/tests/test_interpolate.nim index b70ee9f..363c72a 100644 --- a/tests/test_interpolate.nim +++ b/tests/test_interpolate.nim @@ -191,7 +191,7 @@ test "Extrapolate Constant": for interp in [linear1DSpline, cubicSpline, hermiteSpline]: let res = interp.eval(tOutside, Constant, NaN) for x in res: - check x.isNan() + check classify(x) == fcNan test "Extrapolate Edge": for interp in [linear1DSpline, cubicSpline, hermiteSpline]: From 445bf0a3931c467b048715c0abeba7b5d69e5749 Mon Sep 17 00:00:00 2001 From: HugoGranstrom <5092565+HugoGranstrom@users.noreply.github.com> Date: Fri, 7 Jul 2023 11:10:45 +0200 Subject: [PATCH 4/8] extrapolation for derivEval --- src/numericalnim/interpolate.nim | 40 ++++++++++++++++++++++++++++---- 1 file changed, 36 insertions(+), 4 deletions(-) diff --git a/src/numericalnim/interpolate.nim b/src/numericalnim/interpolate.nim index 868d4a9..43e5fe3 100644 --- a/src/numericalnim/interpolate.nim +++ b/src/numericalnim/interpolate.nim @@ -326,9 +326,41 @@ proc eval*[T; U](interpolator: InterpolatorType[T], x: float, extrap: Extrapolat 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) + 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`. @@ -344,11 +376,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. From 58d006ad7bd3475c760e6352cf12a556cd4dc9b9 Mon Sep 17 00:00:00 2001 From: HugoGranstrom <5092565+HugoGranstrom@users.noreply.github.com> Date: Fri, 7 Jul 2023 11:10:57 +0200 Subject: [PATCH 5/8] test extrapolation deriv --- tests/test_interpolate.nim | 35 +++++++++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) diff --git a/tests/test_interpolate.nim b/tests/test_interpolate.nim index 363c72a..576efe5 100644 --- a/tests/test_interpolate.nim +++ b/tests/test_interpolate.nim @@ -186,6 +186,7 @@ test "linear1DSpline derivEval, seq input": 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]: @@ -219,6 +220,40 @@ test "Extrapolate Native (Hermit)": 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 From 595d86777c0547babec8b98fe43b9fb2d38554b6 Mon Sep 17 00:00:00 2001 From: HugoGranstrom <5092565+HugoGranstrom@users.noreply.github.com> Date: Fri, 7 Jul 2023 11:23:04 +0200 Subject: [PATCH 6/8] add doc comments --- src/numericalnim/interpolate.nim | 21 ++++++++++++++++++++- 1 file changed, 20 insertions(+), 1 deletion(-) diff --git a/src/numericalnim/interpolate.nim b/src/numericalnim/interpolate.nim index 43e5fe3..6cc8dd2 100644 --- a/src/numericalnim/interpolate.nim +++ b/src/numericalnim/interpolate.nim @@ -290,7 +290,16 @@ proc missing(): Missing = discard proc eval*[T; U](interpolator: InterpolatorType[T], x: float, extrap: ExtrapolateKind = Native, extrapValue: U = missing()): T = ## Evaluates an interpolator. - # check for extrapolation + ## - `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: @@ -328,6 +337,16 @@ proc eval*[T; U](interpolator: InterpolatorType[T], x: float, extrap: Extrapolat proc derivEval*[T; U](interpolator: InterpolatorType[T], x: float, extrap: ExtrapolateKind = Native, extrapValue: U = missing()): T = ## Evaluates the derivative of an interpolator. + ## - `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: From 8e235e89a63a8eb32fbe727cd4647a6f1dc602c5 Mon Sep 17 00:00:00 2001 From: HugoGranstrom <5092565+HugoGranstrom@users.noreply.github.com> Date: Fri, 7 Jul 2023 11:28:54 +0200 Subject: [PATCH 7/8] add example --- src/numericalnim/interpolate.nim | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/numericalnim/interpolate.nim b/src/numericalnim/interpolate.nim index 6cc8dd2..61ce78d 100644 --- a/src/numericalnim/interpolate.nim +++ b/src/numericalnim/interpolate.nim @@ -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. From 5d2608bf2ef38697f37845e6e42e5abfda896e73 Mon Sep 17 00:00:00 2001 From: HugoGranstrom <5092565+HugoGranstrom@users.noreply.github.com> Date: Fri, 7 Jul 2023 11:36:14 +0200 Subject: [PATCH 8/8] update changelog --- changelog.md | 34 ++++++++++++++++++++++++++++++++++ 1 file changed, 34 insertions(+) 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.