From 0da0588d324165f3cc8997d2c314d2ff1f971dd2 Mon Sep 17 00:00:00 2001 From: Daniel Ingraham Date: Wed, 22 Jun 2022 14:37:39 -0400 Subject: [PATCH 01/32] Make `brent` cs-safe, and test that we can diff through it --- src/roots.jl | 20 ++++++++++---------- test/runtests.jl | 15 +++++++++++++++ 2 files changed, 25 insertions(+), 10 deletions(-) diff --git a/src/roots.jl b/src/roots.jl index bcdcdd8..8017452 100644 --- a/src/roots.jl +++ b/src/roots.jl @@ -35,9 +35,9 @@ function brent(f, a, b; args=(), atol=2e-12, rtol=4*eps(), maxiter=100) funcalls = 2 iterations = 0 - if fpre*fcur > 0 + if real(fpre)*real(fcur) > 0 error_num = "SIGNERR" - return 0.0, (iter=iterations, fcalls=funcalls, flag=error_num) + return zero(fpre), (iter=iterations, fcalls=funcalls, flag=error_num) end if fpre == zero(fpre) error_num = "CONVERGED" @@ -50,12 +50,12 @@ function brent(f, a, b; args=(), atol=2e-12, rtol=4*eps(), maxiter=100) for i = 1:maxiter iterations += 1 - if fpre*fcur < 0 + if real(fpre)*real(fcur) < 0 xblk = xpre fblk = fpre spre = scur = xcur - xpre end - if abs(fblk) < abs(fcur) + if abs(real(fblk)) < abs(real(fcur)) xpre = xcur xcur = xblk xblk = xpre @@ -65,14 +65,14 @@ function brent(f, a, b; args=(), atol=2e-12, rtol=4*eps(), maxiter=100) fblk = fpre end - delta = (atol + rtol*abs(xcur))/2 + delta = (atol + rtol*abs_cs_safe(xcur))/2 sbis = (xblk - xcur)/2 - if fcur == zero(fcur) || abs(sbis) < delta + if real(fcur) == zero(real(fcur)) || abs(real(sbis)) < real(delta) error_num = "CONVERGED" return xcur, (iter=iterations, fcalls=funcalls, flag=error_num) end - if abs(spre) > delta && abs(fcur) < abs(fpre) + if abs(real(spre)) > real(delta) && abs(real(fcur)) < abs(real(fpre)) if xpre == xblk # interpolate stry = -fcur*(xcur - xpre)/(fcur - fpre) @@ -82,7 +82,7 @@ function brent(f, a, b; args=(), atol=2e-12, rtol=4*eps(), maxiter=100) dblk = (fblk - fcur)/(xblk - xcur) stry = -fcur*(fblk*dblk - fpre*dpre)/(dblk*dpre*(fblk - fpre)) end - if 2*abs(stry) < min(abs(spre), 3*abs(sbis) - delta) + if 2*abs(real(stry)) < min(abs(real(spre)), 3*abs(real(sbis)) - real(delta)) # good short step spre = scur scur = stry @@ -98,10 +98,10 @@ function brent(f, a, b; args=(), atol=2e-12, rtol=4*eps(), maxiter=100) end xpre = xcur; fpre = fcur - if abs(scur) > delta + if abs(real(scur)) > real(delta) xcur += scur else - xcur += (sbis > 0 ? delta : -delta) + xcur += (real(sbis) > 0 ? delta : -delta) end fcur = f(xcur, args...) diff --git a/test/runtests.jl b/test/runtests.jl index 1c7d39c..8366656 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -125,6 +125,21 @@ atol = 1e-15 xstar, _ = brent(f, 1, 4, atol=atol) @test isapprox(xstar, pi, atol=atol) +# Test that we can diff through FLOWMath.brent. +function g(x) + # xL and xR need to bracket the root of the function, so may need to adjust if `x0` below is changed. + xL, xR = 0.0, 10.0 + ystar, info = brent(y->cos(y)-x*y, xL, xR) + info.flag == "CONVERGED" || error("brent solver failed to converge: info = $info") + return sin(ystar) + cos(ystar)^2 +end +x0 = 3.0 +dfdx_fd = ForwardDiff.derivative(g, x0) +dfdx_finitediff = FiniteDiff.finite_difference_derivative(g, x0, Val{:forward}) +dfdx_cs = FiniteDiff.finite_difference_derivative(g, x0, Val{:complex}) +@test abs(dfdx_finitediff - dfdx_fd) < 1e-8 +@test abs(dfdx_cs - dfdx_fd) < 1e-12 + # ------------------------- # ------ abs_smooth --------- From d85808cdbfd743bbab5c2c470f8e5b70d621f285 Mon Sep 17 00:00:00 2001 From: Daniel Ingraham Date: Thu, 23 Jun 2022 10:20:06 -0400 Subject: [PATCH 02/32] bump version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 6c1f8d2..a96b5b9 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "FLOWMath" uuid = "6cb5d3fb-0fe8-4cc2-bd89-9fe0b19a99d3" authors = ["Andrew Ning "] -version = "0.3.3" +version = "0.3.4" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" From 9dd334d47676c512e32c7354ecde0076865db900 Mon Sep 17 00:00:00 2001 From: Judd Mehr Date: Wed, 28 Sep 2022 14:11:10 -0600 Subject: [PATCH 03/32] add second derivative function and matching test --- src/FLOWMath.jl | 1 + src/interpolate.jl | 23 +++++++++++++++++++++++ test/runtests.jl | 4 ++++ 3 files changed, 28 insertions(+) diff --git a/src/FLOWMath.jl b/src/FLOWMath.jl index b5006aa..a12804f 100644 --- a/src/FLOWMath.jl +++ b/src/FLOWMath.jl @@ -24,6 +24,7 @@ export quintic_blend include("interpolate.jl") export Akima export derivative +export derivative2 export gradient export akima export linear diff --git a/src/interpolate.jl b/src/interpolate.jl index 01b3f10..45e9497 100644 --- a/src/interpolate.jl +++ b/src/interpolate.jl @@ -169,6 +169,29 @@ function derivative(spline::Akima, x) return dydx end +""" + derivative2(spline, x) + +Computes the second derivative of an Akima spline at x. + +**Arguments** +- `spline::Akima}`: an Akima spline +- `x::Float`: the evaluation point(s) + +**Returns** +- `d2ydx2::Float`: second derivative at x using akima spline. +""" +function derivative2(spline::Akima, x) + + j = findindex(spline.xdata, x) + + # evaluate polynomial + dx = x - spline.xdata[j] + d2ydx2 = 2.0*spline.p2[j] + 6.0*spline.p3[j]*dx + + return d2ydx2 +end + """ gradient(spline, x) diff --git a/test/runtests.jl b/test/runtests.jl index 1c7d39c..5c8e80c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -586,6 +586,10 @@ dydx = derivative(spl, pi/16) @test isapprox(dydx, 1.0180260961104746, atol=1e-12) +d2ydx2 = derivative2(spl, pi/16) + +@test isapprox(d2ydx2, -0.7003004339939436, atol=1e-12) + dydx = gradient(spl, xpt) dydx_test = [1.1640128599466306, 1.0180260961104746, 0.889005522603991, 0.7769511394271804, 0.6818629465800423, 0.5473879346340718, 0.3889191062220746, 0.20645646134405082, 1.8419755006770065e-16, -0.13818978335121798, -0.33430576382780836, -0.5883479414297708, -0.900316316157106, -0.9003163161571062, -0.9003163161571062, -0.9003163161571061, -0.9003163161571061, -0.9003163161571064, -0.9003163161571064, -0.9003163161571065, -0.9003163161571061, -0.5883479414297709, -0.33430576382780913, -0.13818978335121862, -1.8419755006770062e-16, 0.20645646134405057, 0.3889191062220747, 0.5473879346340714, 0.6818629465800421, 0.7769511394271802, 0.8890055226039906, 1.0180260961104746, 1.1640128599466313] From bcca0fde2c72fff9be75bdb998929005a32756a6 Mon Sep 17 00:00:00 2001 From: Judd Mehr Date: Wed, 28 Sep 2022 17:21:59 -0600 Subject: [PATCH 04/32] change function name to --- src/interpolate.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/interpolate.jl b/src/interpolate.jl index 45e9497..5ca80b8 100644 --- a/src/interpolate.jl +++ b/src/interpolate.jl @@ -170,7 +170,7 @@ function derivative(spline::Akima, x) end """ - derivative2(spline, x) + second_derivative(spline, x) Computes the second derivative of an Akima spline at x. @@ -181,7 +181,7 @@ Computes the second derivative of an Akima spline at x. **Returns** - `d2ydx2::Float`: second derivative at x using akima spline. """ -function derivative2(spline::Akima, x) +function second_derivative(spline::Akima, x) j = findindex(spline.xdata, x) From c68936185f4df5ca1e98f5518c716d51a668cefa Mon Sep 17 00:00:00 2001 From: Judd Mehr Date: Wed, 28 Sep 2022 17:25:56 -0600 Subject: [PATCH 05/32] add second derivative to docs --- docs/src/index.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/src/index.md b/docs/src/index.md index fee7c47..94c38f4 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -87,12 +87,14 @@ You can also compute the derivative and/or gradient of the spline. ```@example akima dydx = derivative(spline, pi/2) dydx = gradient(spline, xpt) +d2ydx2 = second_derivative(spline, pi/2) nothing # hide ``` ```@docs derivative gradient +second_derivative ``` From 2e5b179d57ba115e7bbcfce716717511e6b46bbe Mon Sep 17 00:00:00 2001 From: Judd Mehr Date: Wed, 28 Sep 2022 17:29:52 -0600 Subject: [PATCH 06/32] fix new function name --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 5c8e80c..2eb934b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -586,7 +586,7 @@ dydx = derivative(spl, pi/16) @test isapprox(dydx, 1.0180260961104746, atol=1e-12) -d2ydx2 = derivative2(spl, pi/16) +d2ydx2 = second_derivative(spl, pi/16) @test isapprox(d2ydx2, -0.7003004339939436, atol=1e-12) From 0dc430e36d00b36e24db4711181037c15c97d65e Mon Sep 17 00:00:00 2001 From: Judd Mehr Date: Wed, 28 Sep 2022 17:33:38 -0600 Subject: [PATCH 07/32] update function export and run tests --- src/FLOWMath.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/FLOWMath.jl b/src/FLOWMath.jl index a12804f..076e382 100644 --- a/src/FLOWMath.jl +++ b/src/FLOWMath.jl @@ -24,7 +24,7 @@ export quintic_blend include("interpolate.jl") export Akima export derivative -export derivative2 +export second_derivative export gradient export akima export linear From 906f64901640c7bec5a93131bf8714fb9ad8a6ef Mon Sep 17 00:00:00 2001 From: Andrew Ning Date: Tue, 4 Oct 2022 10:15:37 -0600 Subject: [PATCH 08/32] minor doc update --- docs/src/index.md | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/docs/src/index.md b/docs/src/index.md index 94c38f4..67a52ec 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -47,6 +47,10 @@ An Akima spline is a 1D spline that avoids overshooting issues common with many Interpolation is perhaps clearest through plotting so we'll load a plotting package for this examples. +```@setup akima +using PyPlot +``` + ```@example akima using PyPlot using FLOWMath: akima, Akima, derivative, gradient @@ -338,7 +342,8 @@ quintic_blend ``` ### Complex-step safe functions -The [complex-step derivative approximation](https://doi.org/10.1145/838250.838251) can be used to easily and accurately approximate first derivatives. + +The [complex-step derivative approximation](https://doi.org/10.1145/838250.838251) can be used to easily and accurately approximate first derivatives. This is particularly useful to verify derivatives computed via other means like AD (in contrast to comparing against finite differencing, which suffers from inaccuracies). However, the function `f` one wishes to differentiate must be composed of functions that are compatible with the method. Most elementary functions are, but a few common ones are not: From 796f57c405a4df978036f5180c50c824a354fb2a Mon Sep 17 00:00:00 2001 From: Taylor McDonnell Date: Mon, 17 Apr 2023 17:54:20 -0600 Subject: [PATCH 09/32] add smoothing methods for more than two functions --- src/smooth.jl | 52 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 52 insertions(+) diff --git a/src/smooth.jl b/src/smooth.jl index cb44c71..d7fc9e2 100644 --- a/src/smooth.jl +++ b/src/smooth.jl @@ -137,6 +137,23 @@ function sigmoid_blend(f1x, f2x, x, xt, hardness=50) return f1x + sx*(f2x-f1x) end +""" + sigmoid_blend(fx::Tuple, xt::Tuple, x, hardness=50) + +Smoothly transitions the results of the functions in `fx` using the sigmoid function, +with the transition between the functions at the locations in `xt`. `hardness` controls +the sharpness of the transition between the functions. +""" +function sigmoid_blend(fx::Tuple, xt::Tuple, x, hardness=50) + new_fx = sigmoid_blend.(fx[1:end-1], fx[2:end], x, xt, hardness) + if length(new_fx) == 1 + return only(new_fx) + else + new_xt = (xt[1:end-1] .+ xt[2:end]) ./ 2 + return sigmoid_blend(new_fx, new_xt, x, hardness) + end +end + """ cubic_blend(f1x, f2x, x, xt, delta_x) @@ -156,6 +173,23 @@ function cubic_blend(f1x, f2x, x, xt, delta_x) end end +""" + cubic_blend(fx::Tuple, xt::Tuple, x, delta_x) + +Smoothly transitions the results of the functions in `fx` using cubic polynomials, +with the transition between the functions at the locations in `xt`. `delta_x` is the half +width of the smoothing interval. The resulting function is C1 continuous +""" +function cubic_blend(fx::Tuple, xt, x, delta_x) + new_fx = cubic_blend.(fx[1:end-1], fx[2:end], x, xt, delta_x) + if length(new_fx) == 1 + return only(new_fx) + else + new_xt = (xt[1:end-1] .+ xt[2:end]) ./ 2 + return cubic_blend(new_fx, new_xt, x, delta_x) + end +end + """ quintic_blend(f1x, f2x, x, xt, delta_x) @@ -175,4 +209,22 @@ function quintic_blend(f1x, f2x, x, xt, delta_x) end end + +""" + quintic_blend(fx::Tuple, xt::Tuple, x, delta_x) + +Smoothly transitions the results of the functions in `fx` using quintic polynomials, +with the transition between the functions at the locations in `xt`. `delta_x` is the half +width of the smoothing interval. The resulting function is C2 continuous +""" +function quintic_blend(fx::Tuple, xt, x, delta_x) + new_fx = quintic_blend.(fx[1:end-1], fx[2:end], x, xt, delta_x) + if length(new_fx) == 1 + return only(new_fx) + else + new_xt = (xt[1:end-1] .+ xt[2:end]) ./ 2 + return quintic_blend(new_fx, new_xt, x, delta_x) + end +end + # TODO AN: add smooth max/min with cubic splines From 12f18e205da764f4c717530b13f651540e1fbb22 Mon Sep 17 00:00:00 2001 From: tylercritchfield Date: Wed, 6 Sep 2023 12:00:28 -0600 Subject: [PATCH 10/32] expose eps parameter in Akima spline construction --- src/interpolate.jl | 24 +++++++++++++++--------- 1 file changed, 15 insertions(+), 9 deletions(-) diff --git a/src/interpolate.jl b/src/interpolate.jl index 5ca80b8..03169c3 100644 --- a/src/interpolate.jl +++ b/src/interpolate.jl @@ -57,18 +57,24 @@ end Creates an akima spline at node points: `xdata`, `ydata`. This is a 1D spline that avoids overshooting issues common with many other polynomial splines resulting in a more natural curve. It also only depends on local points (`i-2`...`i+2`) allow for more efficient -computation. `delta_x` is the half width of a smoothing interval used for the absolute +computation. + +`delta_x` is the half width of a smoothing interval used for the absolute value function. Set `delta_x=0` to recover the original akima spline. The smoothing is only useful if you want to differentiate xdata and ydata. In many case the nodal points -are fixed so this is not needed. Returns an akima spline object (Akima struct). -This function, only performs construction of the spline, not evaluation. -This is useful if you want to evaluate the same mesh at multiple different conditions. -A convenience method exists below to perform both in one shot. +are fixed so this is not needed. + +`eps` is a cutoff used to avoid dividing by zero in the +weighting function. Default is `1e-30` but this could be raised to machine precision in +some cases to improve derivatives. (E.g., when the denominator in line 105 is very small.) + +Returns an akima spline object (Akima struct). This function only performs construction +of the spline, not evaluation. This is useful if you want to evaluate the same mesh at +multiple different conditions. A convenience method exists below to perform both in one shot. """ -function Akima(xdata::AbstractVector{TFX}, ydata::AbstractVector{TFY}, delta_x=0.0) where {TFX, TFY} +function Akima(xdata::AbstractVector{TFX}, ydata::AbstractVector{TFY}, delta_x=0.0, eps=1e-30) where {TFX, TFY} # setup - eps = 1e-30 n = length(xdata) TCoeff = promote_type(TFX, TFY) @@ -132,7 +138,7 @@ end (spline::Akima)(x::AbstractVector) = spline.(x) """ - akima(x, y, xpt) + akima(x, y, xpt, delta=0.0, eps=1e-30) A convenience method to perform construction and evaluation of the spline in one step. See docstring for Akima for more details. @@ -144,7 +150,7 @@ See docstring for Akima for more details. **Returns** - `ypt::Vector{Float} or ::Float`: interpolated value(s) at xpt using akima spline. """ -akima(x, y, xpt) = Akima(x, y)(xpt) +akima(x, y, xpt, delta=0.0, eps=1e-30) = Akima(x, y, delta, eps)(xpt) """ derivative(spline, x) From 207b3ddb5117d7263f1242e2f09291ffc757c272 Mon Sep 17 00:00:00 2001 From: Daniel Ingraham Date: Mon, 18 Mar 2024 13:19:23 -0400 Subject: [PATCH 11/32] Small complex-step-safe change to `sigmoid` --- src/smooth.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/smooth.jl b/src/smooth.jl index d7fc9e2..277da94 100644 --- a/src/smooth.jl +++ b/src/smooth.jl @@ -116,7 +116,7 @@ end Sigmoid function, implemented with branching to avoid NaNs """ function sigmoid(x) - if x >= zero(x) + if real(x) >= zero(real(x)) z = exp(-x) return one(z) / (one(z) + z) else From 00cb54b890e446184925f0846a6c14ed32bc4764 Mon Sep 17 00:00:00 2001 From: Judd Mehr Date: Fri, 17 May 2024 15:32:20 -0600 Subject: [PATCH 12/32] add default arguments to akima docstring --- src/interpolate.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/interpolate.jl b/src/interpolate.jl index 03169c3..f9a3853 100644 --- a/src/interpolate.jl +++ b/src/interpolate.jl @@ -146,6 +146,11 @@ See docstring for Akima for more details. **Arguments** - `x, y::Vector{Float}`: the node points - `xpt::Vector{Float} or ::Float`: the evaluation point(s) +- `delta_x::Float=0.0` : the half width of a smoothing interval used for the absolute +value function. +- `eps::Float=1e-30` : a cutoff used to avoid dividing by zero in the +weighting function. Default is `1e-30` but this could be raised to machine precision in +some cases to improve derivatives. **Returns** - `ypt::Vector{Float} or ::Float`: interpolated value(s) at xpt using akima spline. From 31ebbc14ecb42a2a92469224ecd294d95c5e4518 Mon Sep 17 00:00:00 2001 From: Judd Mehr Date: Tue, 25 Jun 2024 09:04:34 -0600 Subject: [PATCH 13/32] Update test.yaml update julia-runtest according to deprecation warning and errors. --- .github/workflows/test.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/test.yaml b/.github/workflows/test.yaml index 8e7f77d..7262c15 100644 --- a/.github/workflows/test.yaml +++ b/.github/workflows/test.yaml @@ -21,4 +21,4 @@ jobs: - uses: julia-actions/setup-julia@v1 with: version: ${{ matrix.julia-version }} - - uses: julia-actions/julia-runtest@master + - uses: julia-actions/julia-runtest@v1 From f3728780fc5129bc37f094dcca642a2697bf7b1d Mon Sep 17 00:00:00 2001 From: juddmehr Date: Tue, 25 Jun 2024 10:40:23 -0600 Subject: [PATCH 14/32] minor updates to interpolation docs --- src/interpolate.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/interpolate.jl b/src/interpolate.jl index 03169c3..cd55203 100644 --- a/src/interpolate.jl +++ b/src/interpolate.jl @@ -25,7 +25,7 @@ function findindex(xvec, x) # this version prevents extrapolation # if i == 0 # throw(DomainError(x, "x falls below range of provided spline data")) - # elseif i == n + # elseif i == n # if x == xvec[n] # i = n - 1 # else @@ -57,19 +57,19 @@ end Creates an akima spline at node points: `xdata`, `ydata`. This is a 1D spline that avoids overshooting issues common with many other polynomial splines resulting in a more natural curve. It also only depends on local points (`i-2`...`i+2`) allow for more efficient -computation. +computation. `delta_x` is the half width of a smoothing interval used for the absolute value function. Set `delta_x=0` to recover the original akima spline. The smoothing is only useful if you want to differentiate xdata and ydata. In many case the nodal points -are fixed so this is not needed. +are fixed so this is not needed. -`eps` is a cutoff used to avoid dividing by zero in the -weighting function. Default is `1e-30` but this could be raised to machine precision in +`eps` is a cutoff used to avoid dividing by zero in the +weighting function. Default is `1e-30` but this could be raised to machine precision in some cases to improve derivatives. (E.g., when the denominator in line 105 is very small.) Returns an akima spline object (Akima struct). This function only performs construction -of the spline, not evaluation. This is useful if you want to evaluate the same mesh at +of the spline, not evaluation. This is useful if you want to evaluate the same mesh at multiple different conditions. A convenience method exists below to perform both in one shot. """ function Akima(xdata::AbstractVector{TFX}, ydata::AbstractVector{TFY}, delta_x=0.0, eps=1e-30) where {TFX, TFY} @@ -241,7 +241,7 @@ end """ linear(xdata, ydata, x::AbstractVector) - + Convenience function to perform linear interpolation at multiple points. """ linear(xdata, ydata, x::AbstractVector) = linear.(Ref(xdata), Ref(ydata), x) From f3bf371903b50e3e4349dea52a93b652904a9f0c Mon Sep 17 00:00:00 2001 From: juddmehr Date: Tue, 25 Jun 2024 10:44:36 -0600 Subject: [PATCH 15/32] put using PyPlot in setup blocks --- docs/src/index.md | 22 +++++++++++++++------- 1 file changed, 15 insertions(+), 7 deletions(-) diff --git a/docs/src/index.md b/docs/src/index.md index 67a52ec..efdc389 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -52,7 +52,6 @@ using PyPlot ``` ```@example akima -using PyPlot using FLOWMath: akima, Akima, derivative, gradient x = 0:pi/4:2*pi @@ -104,7 +103,7 @@ second_derivative ### Linear Spline -Linear interpolation is straightforward. +Linear interpolation is straightforward. ```@example linear using FLOWMath: linear, derivative, gradient @@ -213,7 +212,11 @@ interp4d The absolute value function is not differentiable at x = 0. The below function smoothly adds a small quadratic function in place of the cusp with a half-width given by `delta_x`. This small rounding at the bottom can prevent numerical issues with gradient-based optimization. -```@example +```@setup abbs +using PyPlot +``` + +```@example abss using FLOWMath: abs_smooth x = range(-2.0, 2.0, length=100) @@ -221,7 +224,6 @@ delta_x = 0.1 y = abs_smooth.(x, delta_x) -using PyPlot figure() plot(x, y) savefig("abs.svg"); nothing # hide @@ -236,7 +238,7 @@ abs_smooth ### Kreisselmeier-Steinhauser Constraint Aggregation Function -The Kreisselmeier-Steinhauser (KS) function is often used with constrained gradient-based optimization problems to smoothly aggregate an arbitrary number of constraints into a single constraint. It may also be used as a smooth approximation of the maximum function (or minimum function). A salient feature of this function is that it is guaranteed to overestimate the maximum function (or underestimate the minimum function). This feature of the function can be used to ensure that the resulting constraint is conservative. +The Kreisselmeier-Steinhauser (KS) function is often used with constrained gradient-based optimization problems to smoothly aggregate an arbitrary number of constraints into a single constraint. It may also be used as a smooth approximation of the maximum function (or minimum function). A salient feature of this function is that it is guaranteed to overestimate the maximum function (or underestimate the minimum function). This feature of the function can be used to ensure that the resulting constraint is conservative. We provide two implementations of this function: `ksmax` and `ksmin`. `ksmax` and `ksmin` may be used to smoothly approximate the maximum and minimum functions, respectively. Both functions take the optional parameter `hardness` that controls the smoothness of the resulting function. As `hardness` increases the function more and more closely approximates the maximum (or minimum) function. @@ -261,6 +263,10 @@ ksmin The sigmoid function may be used to smoothly blend the results of two continuous one-dimensional functions. The method implemented in this package uses a user-specified transition location (`xt`) and scales the input of the sigmoid function using the input `hardness` in order to adjust the smoothness of the transition between the two functions. +```setup sb +using PyPlot +``` + ```@example sb using FLOWMath: sigmoid_blend @@ -282,7 +288,6 @@ xt = 0.0 hardness = 25 y = sigmoid_blend.(f1x, f2x, x, xt, hardness) -using PyPlot figure() plot(x, f1x) plot(x, f2x) @@ -301,6 +306,10 @@ sigmoid_blend Cubic or quintic polynomials can also be used to construct a piecewise function that smoothly blends two functions. The advantage of this approach compared to `sigmoid_blend` is that the blending can be restricted to a small interval defined by the half-width `delta_x`. The disadvantage of this approach is that the resulting function is only C1 continuous when `cubic_blend` is used, and C2 continuous when `quintic_blend` is used. The method implemented in this package uses a user-specified transition location (`xt`). The smoothness of the transition between the two functions can be adjusted by modifying `delta_x`, which is the half-width of the transition interval. +```@setup poly +using PyPlot +``` + ```@example poly using FLOWMath: cubic_blend, quintic_blend @@ -324,7 +333,6 @@ delta_x = 0.1 y1 = cubic_blend.(f1x, f2x, x, xt, delta_x) y2 = quintic_blend.(f1x, f2x, x, xt, delta_x) -using PyPlot figure() plot(x, f1x) plot(x, f2x) From ef8801b2091657837c6e446de87dd17b1041d9fe Mon Sep 17 00:00:00 2001 From: juddmehr Date: Tue, 25 Jun 2024 11:16:05 -0600 Subject: [PATCH 16/32] cleanup docs and remove references to external packages that docs can't find. --- docs/Manifest.toml | 386 ++++++++++++++++++++++++++++++++------------- docs/Project.toml | 2 +- docs/make.jl | 12 +- docs/src/index.md | 27 ++-- src/cs_safe.jl | 10 -- src/interpolate.jl | 106 ++++++------- src/quadrature.jl | 10 +- src/roots.jl | 41 ++--- src/smooth.jl | 75 +++++---- 9 files changed, 406 insertions(+), 263 deletions(-) diff --git a/docs/Manifest.toml b/docs/Manifest.toml index a11b446..1e9c69a 100644 --- a/docs/Manifest.toml +++ b/docs/Manifest.toml @@ -1,188 +1,360 @@ # This file is machine-generated - editing it directly is not advised -[[Base64]] +julia_version = "1.10.0" +manifest_format = "2.0" +project_hash = "e3a22c471fb9a9ace33800905d36686a71c62915" + +[[deps.ANSIColoredPrinters]] +git-tree-sha1 = "574baf8110975760d391c710b6341da1afa48d8c" +uuid = "a4c015fc-c6ff-483c-b24f-f7ea428134e9" +version = "0.0.1" + +[[deps.AbstractTrees]] +git-tree-sha1 = "2d9c9a55f9c93e8887ad391fbae72f8ef55e1177" +uuid = "1520ce14-60c1-5f80-bbc7-55ef81b5835c" +version = "0.4.5" + +[[deps.ArgTools]] +uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f" +version = "1.1.1" + +[[deps.Artifacts]] +uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" + +[[deps.Base64]] uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" -[[ColorTypes]] +[[deps.CodecZlib]] +deps = ["TranscodingStreams", "Zlib_jll"] +git-tree-sha1 = "59939d8a997469ee05c4b4944560a820f9ba0d73" +uuid = "944b1d66-785c-5afd-91f1-9de20f533193" +version = "0.7.4" + +[[deps.ColorTypes]] deps = ["FixedPointNumbers", "Random"] -git-tree-sha1 = "b9de8dc6106e09c79f3f776c27c62360d30e5eb8" +git-tree-sha1 = "b10d0b65641d57b8b4d5e234446582de5047050d" uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f" -version = "0.9.1" +version = "0.11.5" -[[Colors]] -deps = ["ColorTypes", "FixedPointNumbers", "InteractiveUtils", "Printf", "Reexport"] -git-tree-sha1 = "177d8b959d3c103a6d57574c38ee79c81059c31b" +[[deps.Colors]] +deps = ["ColorTypes", "FixedPointNumbers", "Reexport"] +git-tree-sha1 = "362a287c3aa50601b0bc359053d5c2468f0e7ce0" uuid = "5ae59095-9a9b-59fe-a467-6f913c188581" -version = "0.11.2" +version = "0.12.11" -[[Compat]] -deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"] -git-tree-sha1 = "ed2c4abadf84c53d9e58510b5fc48912c2336fbb" -uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" -version = "2.2.0" +[[deps.CompilerSupportLibraries_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" +version = "1.0.5+1" -[[Conda]] -deps = ["JSON", "VersionParsing"] -git-tree-sha1 = "9a11d428dcdc425072af4aea19ab1e8c3e01c032" +[[deps.Conda]] +deps = ["Downloads", "JSON", "VersionParsing"] +git-tree-sha1 = "51cab8e982c5b598eea9c8ceaced4b58d9dd37c9" uuid = "8f4d0f93-b110-5947-807f-2305c1781a2d" -version = "1.3.0" - -[[DataStructures]] -deps = ["InteractiveUtils", "OrderedCollections"] -git-tree-sha1 = "b7720de347734f4716d1815b00ce5664ed6bbfd4" -uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" -version = "0.17.9" +version = "1.10.0" -[[Dates]] +[[deps.Dates]] deps = ["Printf"] uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" -[[DelimitedFiles]] -deps = ["Mmap"] -uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab" - -[[Distributed]] -deps = ["Random", "Serialization", "Sockets"] -uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b" - -[[DocStringExtensions]] -deps = ["LibGit2", "Markdown", "Pkg", "Test"] -git-tree-sha1 = "88bb0edb352b16608036faadcc071adda068582a" +[[deps.DocStringExtensions]] +deps = ["LibGit2"] +git-tree-sha1 = "2fb1e02f2b635d0845df5d7c167fec4dd739b00d" uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" -version = "0.8.1" +version = "0.9.3" -[[Documenter]] -deps = ["Base64", "Dates", "DocStringExtensions", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "REPL", "Test", "Unicode"] -git-tree-sha1 = "d497bcc45bb98a1fbe19445a774cfafeabc6c6df" +[[deps.Documenter]] +deps = ["ANSIColoredPrinters", "AbstractTrees", "Base64", "CodecZlib", "Dates", "DocStringExtensions", "Downloads", "Git", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "MarkdownAST", "Pkg", "PrecompileTools", "REPL", "RegistryInstances", "SHA", "TOML", "Test", "Unicode"] +git-tree-sha1 = "5461b2a67beb9089980e2f8f25145186b6d34f91" uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4" -version = "0.24.5" +version = "1.4.1" -[[FixedPointNumbers]] -git-tree-sha1 = "4aaea64dd0c30ad79037084f8ca2b94348e65eaa" -uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93" -version = "0.7.1" +[[deps.Downloads]] +deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"] +uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6" +version = "1.6.0" -[[InteractiveUtils]] +[[deps.Expat_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "1c6317308b9dc757616f0b5cb379db10494443a7" +uuid = "2e619515-83b5-522b-bb60-26c02a35a201" +version = "2.6.2+0" + +[[deps.FileWatching]] +uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee" + +[[deps.FixedPointNumbers]] +deps = ["Statistics"] +git-tree-sha1 = "05882d6995ae5c12bb5f36dd2ed3f61c98cbb172" +uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93" +version = "0.8.5" + +[[deps.Git]] +deps = ["Git_jll"] +git-tree-sha1 = "04eff47b1354d702c3a85e8ab23d539bb7d5957e" +uuid = "d7ba0133-e1db-5d97-8f8c-041e4b3a1eb2" +version = "1.3.1" + +[[deps.Git_jll]] +deps = ["Artifacts", "Expat_jll", "JLLWrappers", "LibCURL_jll", "Libdl", "Libiconv_jll", "OpenSSL_jll", "PCRE2_jll", "Zlib_jll"] +git-tree-sha1 = "d18fb8a1f3609361ebda9bf029b60fd0f120c809" +uuid = "f8c6e375-362e-5223-8a59-34ff63f689eb" +version = "2.44.0+2" + +[[deps.IOCapture]] +deps = ["Logging", "Random"] +git-tree-sha1 = "b6d6bfdd7ce25b0f9b2f6b3dd56b2673a66c8770" +uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89" +version = "0.2.5" + +[[deps.InteractiveUtils]] deps = ["Markdown"] uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" -[[JSON]] +[[deps.JLLWrappers]] +deps = ["Artifacts", "Preferences"] +git-tree-sha1 = "7e5d6779a1e09a36db2a7b6cff50942a0a7d0fca" +uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210" +version = "1.5.0" + +[[deps.JSON]] deps = ["Dates", "Mmap", "Parsers", "Unicode"] -git-tree-sha1 = "b34d7cef7b337321e97d22242c3c2b91f476748e" +git-tree-sha1 = "31e996f0a15c7b280ba9f76636b3ff9e2ae58c9a" uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" -version = "0.21.0" +version = "0.21.4" -[[LaTeXStrings]] -deps = ["Compat"] -git-tree-sha1 = "7ab9b8788cfab2bdde22adf9004bda7ad9954b6c" +[[deps.LaTeXStrings]] +git-tree-sha1 = "50901ebc375ed41dbf8058da26f9de442febbbec" uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" -version = "1.0.3" +version = "1.3.1" -[[LibGit2]] +[[deps.LazilyInitializedFields]] +git-tree-sha1 = "8f7f3cabab0fd1800699663533b6d5cb3fc0e612" +uuid = "0e77f7df-68c5-4e49-93ce-4cd80f5598bf" +version = "1.2.2" + +[[deps.LibCURL]] +deps = ["LibCURL_jll", "MozillaCACerts_jll"] +uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21" +version = "0.6.4" + +[[deps.LibCURL_jll]] +deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"] +uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0" +version = "8.4.0+0" + +[[deps.LibGit2]] +deps = ["Base64", "LibGit2_jll", "NetworkOptions", "Printf", "SHA"] uuid = "76f85450-5226-5b5a-8eaa-529ad045b433" -[[Libdl]] +[[deps.LibGit2_jll]] +deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll"] +uuid = "e37daf67-58a4-590a-8e99-b0245dd2ffc5" +version = "1.6.4+0" + +[[deps.LibSSH2_jll]] +deps = ["Artifacts", "Libdl", "MbedTLS_jll"] +uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8" +version = "1.11.0+1" + +[[deps.Libdl]] uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" -[[LinearAlgebra]] -deps = ["Libdl"] +[[deps.Libiconv_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "f9557a255370125b405568f9767d6d195822a175" +uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531" +version = "1.17.0+0" + +[[deps.LinearAlgebra]] +deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"] uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -[[Logging]] +[[deps.Logging]] uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" -[[MacroTools]] -deps = ["DataStructures", "Markdown", "Random"] -git-tree-sha1 = "e2fc7a55bb2224e203bbd8b59f72b91323233458" +[[deps.MacroTools]] +deps = ["Markdown", "Random"] +git-tree-sha1 = "2fa9ee3e63fd3a4f7a9a4f4744a52f4856de82df" uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" -version = "0.5.3" +version = "0.5.13" -[[Markdown]] +[[deps.Markdown]] deps = ["Base64"] uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" -[[Mmap]] +[[deps.MarkdownAST]] +deps = ["AbstractTrees", "Markdown"] +git-tree-sha1 = "465a70f0fc7d443a00dcdc3267a497397b8a3899" +uuid = "d0879d2d-cac2-40c8-9cee-1863dc0c7391" +version = "0.1.2" + +[[deps.MbedTLS_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1" +version = "2.28.2+1" + +[[deps.Mmap]] uuid = "a63ad114-7e13-5084-954f-fe012c677804" -[[OrderedCollections]] -deps = ["Random", "Serialization", "Test"] -git-tree-sha1 = "c4c13474d23c60d20a67b217f1d7f22a40edf8f1" -uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" -version = "1.1.0" +[[deps.MozillaCACerts_jll]] +uuid = "14a3606d-f60d-562e-9121-12d972cd8159" +version = "2023.1.10" + +[[deps.NetworkOptions]] +uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908" +version = "1.2.0" -[[Parsers]] -deps = ["Dates", "Test"] -git-tree-sha1 = "0139ba59ce9bc680e2925aec5b7db79065d60556" +[[deps.OpenBLAS_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"] +uuid = "4536629a-c528-5b80-bd46-f80d51c5b363" +version = "0.3.23+2" + +[[deps.OpenSSL_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "a028ee3cb5641cccc4c24e90c36b0a4f7707bdf5" +uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95" +version = "3.0.14+0" + +[[deps.PCRE2_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "efcefdf7-47ab-520b-bdef-62a2eaa19f15" +version = "10.42.0+1" + +[[deps.Parsers]] +deps = ["Dates", "PrecompileTools", "UUIDs"] +git-tree-sha1 = "8489905bcdbcfac64d1daa51ca07c0d8f0283821" uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" -version = "0.3.10" +version = "2.8.1" -[[Pkg]] -deps = ["Dates", "LibGit2", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"] +[[deps.Pkg]] +deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"] uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" +version = "1.10.0" + +[[deps.PrecompileTools]] +deps = ["Preferences"] +git-tree-sha1 = "5aa36f7049a63a1528fe8f7c3f2113413ffd4e1f" +uuid = "aea7be01-6a6a-4083-8856-8a6e6704d82a" +version = "1.2.1" + +[[deps.Preferences]] +deps = ["TOML"] +git-tree-sha1 = "9306f6085165d270f7e3db02af26a400d580f5c6" +uuid = "21216c6a-2e73-6563-6e65-726566657250" +version = "1.4.3" -[[Printf]] +[[deps.Printf]] deps = ["Unicode"] uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" -[[PyCall]] -deps = ["Conda", "Dates", "Libdl", "LinearAlgebra", "MacroTools", "Pkg", "Serialization", "Statistics", "Test", "VersionParsing"] -git-tree-sha1 = "6e5bac1b1faf3575731a6a5b76f638f2389561d3" +[[deps.PyCall]] +deps = ["Conda", "Dates", "Libdl", "LinearAlgebra", "MacroTools", "Serialization", "VersionParsing"] +git-tree-sha1 = "9816a3826b0ebf49ab4926e2b18842ad8b5c8f04" uuid = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0" -version = "1.91.2" +version = "1.96.4" -[[PyPlot]] +[[deps.PyPlot]] deps = ["Colors", "LaTeXStrings", "PyCall", "Sockets", "Test", "VersionParsing"] -git-tree-sha1 = "ccecc72cf5b41a5de686bd76999040050a8a3472" +git-tree-sha1 = "9220a9dae0369f431168d60adab635f88aca7857" uuid = "d330b81b-6aea-500a-939a-2ce795aea3ee" -version = "2.8.2" +version = "2.11.2" -[[REPL]] -deps = ["InteractiveUtils", "Markdown", "Sockets"] +[[deps.REPL]] +deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"] uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" -[[Random]] -deps = ["Serialization"] +[[deps.Random]] +deps = ["SHA"] uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" -[[Reexport]] -deps = ["Pkg"] -git-tree-sha1 = "7b1d07f411bc8ddb7977ec7f377b97b158514fe0" +[[deps.Reexport]] +git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b" uuid = "189a3867-3050-52da-a836-e630ba90ab69" -version = "0.2.0" +version = "1.2.2" -[[SHA]] +[[deps.RegistryInstances]] +deps = ["LazilyInitializedFields", "Pkg", "TOML", "Tar"] +git-tree-sha1 = "ffd19052caf598b8653b99404058fce14828be51" +uuid = "2792f1a3-b283-48e8-9a74-f99dce5104f3" +version = "0.1.0" + +[[deps.SHA]] uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" +version = "0.7.0" -[[Serialization]] +[[deps.Serialization]] uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" -[[SharedArrays]] -deps = ["Distributed", "Mmap", "Random", "Serialization"] -uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383" - -[[Sockets]] +[[deps.Sockets]] uuid = "6462fe0b-24de-5631-8697-dd941f90decc" -[[SparseArrays]] -deps = ["LinearAlgebra", "Random"] +[[deps.SparseArrays]] +deps = ["Libdl", "LinearAlgebra", "Random", "Serialization", "SuiteSparse_jll"] uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +version = "1.10.0" -[[Statistics]] +[[deps.Statistics]] deps = ["LinearAlgebra", "SparseArrays"] uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +version = "1.10.0" -[[Test]] -deps = ["Distributed", "InteractiveUtils", "Logging", "Random"] +[[deps.SuiteSparse_jll]] +deps = ["Artifacts", "Libdl", "libblastrampoline_jll"] +uuid = "bea87d4a-7f5b-5778-9afe-8cc45184846c" +version = "7.2.1+1" + +[[deps.TOML]] +deps = ["Dates"] +uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76" +version = "1.0.3" + +[[deps.Tar]] +deps = ["ArgTools", "SHA"] +uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e" +version = "1.10.0" + +[[deps.Test]] +deps = ["InteractiveUtils", "Logging", "Random", "Serialization"] uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" -[[UUIDs]] +[[deps.TranscodingStreams]] +git-tree-sha1 = "d73336d81cafdc277ff45558bb7eaa2b04a8e472" +uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa" +version = "0.10.10" +weakdeps = ["Random", "Test"] + + [deps.TranscodingStreams.extensions] + TestExt = ["Test", "Random"] + +[[deps.UUIDs]] deps = ["Random", "SHA"] uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" -[[Unicode]] +[[deps.Unicode]] uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" -[[VersionParsing]] -git-tree-sha1 = "80229be1f670524750d905f8fc8148e5a8c4537f" +[[deps.VersionParsing]] +git-tree-sha1 = "58d6e80b4ee071f5efd07fda82cb9fbe17200868" uuid = "81def892-9a0e-5fdd-b105-ffc91e053289" -version = "1.2.0" +version = "1.3.0" + +[[deps.Zlib_jll]] +deps = ["Libdl"] +uuid = "83775a58-1f1d-513f-b197-d71354ab007a" +version = "1.2.13+1" + +[[deps.libblastrampoline_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "8e850b90-86db-534c-a0d3-1478176c7d93" +version = "5.8.0+1" + +[[deps.nghttp2_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d" +version = "1.52.0+1" + +[[deps.p7zip_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0" +version = "17.4.0+2" diff --git a/docs/Project.toml b/docs/Project.toml index 7c05cc0..61e0fa9 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -3,4 +3,4 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee" [compat] -Documenter = "0.24" +Documenter = "1.4.1" diff --git a/docs/make.jl b/docs/make.jl index c175794..ddec09b 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,17 +1,13 @@ -using Documenter, FLOWMath +using Documenter#, FLOWMath makedocs(; modules=[FLOWMath], format=Documenter.HTML(), - pages=[ - "Home" => "index.md", - ], + pages=["Home" => "index.md"], repo="https://github.com/byuflowlab/FLOWMath.jl/blob/{commit}{path}#L{line}", sitename="FLOWMath.jl", authors="Andrew Ning ", - # assets=String[], + checkdocs=:exports, ) -deploydocs( - repo="github.com/byuflowlab/FLOWMath.jl.git", -) +deploydocs(; repo="github.com/byuflowlab/FLOWMath.jl.git") diff --git a/docs/src/index.md b/docs/src/index.md index efdc389..3aa4db6 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -9,7 +9,7 @@ Examples of the available methods are shown below. More examples are available This is just simple trapezoidal integration using vectors. [Gaussian quadrature](https://github.com/JuliaMath/QuadGK.jl) is much better, but for times when we need to define a mesh for other purposes and cannot use an adaptive method a simple trapezoidal integration fits the bill. ```@example trapz -using FLOWMath: trapz +using ..FLOWMath: trapz x = range(0.0, stop=pi+1e-15, step=pi/100) y = sin.(x) @@ -27,7 +27,7 @@ trapz Brent's method is an effective 1D root finding method as it combines bracketing methods (bisection) with fast quadratic interpolation. Thus, you can get near quadratic convergence but with safeguarding. ```@example -using FLOWMath: brent +using ..FLOWMath: brent f(x) = x^2 - 1.0 xstar, outputs = brent(f, -2.0, 0) @@ -52,7 +52,7 @@ using PyPlot ``` ```@example akima -using FLOWMath: akima, Akima, derivative, gradient +using ..FLOWMath: akima, Akima, derivative, second_derivative, gradient x = 0:pi/4:2*pi y = sin.(x) @@ -106,7 +106,7 @@ second_derivative Linear interpolation is straightforward. ```@example linear -using FLOWMath: linear, derivative, gradient +using ..FLOWMath: linear, derivative, gradient xvec = [1.0, 2.0, 4.0, 5.0] yvec = [2.0, 3.0, 5.0, 8.0] @@ -140,7 +140,7 @@ The functions `interp2d`, `interp3d`, and `interp4d` are generic and will accept 2D: ```@example akima -using FLOWMath: interp2d +using ..FLOWMath: interp2d x = -5.0:5.0 y = -5.0:5.0 @@ -166,7 +166,7 @@ savefig("contour.svg"); nothing # hide 4D: ```@example akima -using FLOWMath: interp4d +using ..FLOWMath: interp4d x = -1:0.2:1 y = -1:0.2:1 @@ -212,12 +212,12 @@ interp4d The absolute value function is not differentiable at x = 0. The below function smoothly adds a small quadratic function in place of the cusp with a half-width given by `delta_x`. This small rounding at the bottom can prevent numerical issues with gradient-based optimization. -```@setup abbs +```@setup abss using PyPlot ``` ```@example abss -using FLOWMath: abs_smooth +using ..FLOWMath: abs_smooth x = range(-2.0, 2.0, length=100) delta_x = 0.1 @@ -243,7 +243,7 @@ The Kreisselmeier-Steinhauser (KS) function is often used with constrained gradi We provide two implementations of this function: `ksmax` and `ksmin`. `ksmax` and `ksmin` may be used to smoothly approximate the maximum and minimum functions, respectively. Both functions take the optional parameter `hardness` that controls the smoothness of the resulting function. As `hardness` increases the function more and more closely approximates the maximum (or minimum) function. ```@example ks -using FLOWMath: ksmax, ksmin +using ..FLOWMath: ksmax, ksmin x = [1.2, 0.0, 0.5] hardness = 100 @@ -257,18 +257,20 @@ min_x = ksmin(x, hardness) ```@docs ksmax ksmin +ksmax_adaptive +ksmin_adaptive ``` ### Blending functions using the sigmoid function The sigmoid function may be used to smoothly blend the results of two continuous one-dimensional functions. The method implemented in this package uses a user-specified transition location (`xt`) and scales the input of the sigmoid function using the input `hardness` in order to adjust the smoothness of the transition between the two functions. -```setup sb +```@setup sb using PyPlot ``` ```@example sb -using FLOWMath: sigmoid_blend +using ..FLOWMath: sigmoid_blend x = 0.1 f1x = x @@ -300,6 +302,7 @@ savefig("sigmoid.svg"); nothing # hide ```@docs sigmoid_blend +sigmoid ``` ### Blending functions using cubic or quintic polynomials @@ -311,7 +314,7 @@ using PyPlot ``` ```@example poly -using FLOWMath: cubic_blend, quintic_blend +using ..FLOWMath: cubic_blend, quintic_blend x = 0.05 f1x = x diff --git a/src/cs_safe.jl b/src/cs_safe.jl index 300c124..8c3b8b5 100644 --- a/src/cs_safe.jl +++ b/src/cs_safe.jl @@ -5,8 +5,6 @@ using LinearAlgebra: norm, dot abs_cs_safe(x) Calculate the absolute value of `x` in a manner compatible with the complex-step derivative approximation. - -See also: [`abs`](@ref). """ abs_cs_safe @@ -22,8 +20,6 @@ end abs2_cs_safe(x) Calculate the squared absolute value of `x` in a manner compatible with the complex-step derivative approximation. - -See also: [`abs2`](@ref). """ abs2_cs_safe @@ -39,8 +35,6 @@ end norm_cs_safe(x, p) Calculate the `p`-norm value of iterable `x` in a manner compatible with the complex-step derivative approximation. - -See also: [`norm`](@ref). """ norm_cs_safe @@ -56,8 +50,6 @@ end dot_cs_safe(a, b) Calculate the dot product of vectors `a` and `b` in a manner compatible with the complex-step derivative approximation. - -See also: [`norm`](@ref). """ dot_cs_safe @@ -76,8 +68,6 @@ end atan_cs_safe(y, x) Calculate the two-argument arctangent function in a manner compatible with the complex-step derivative approximation. - -See also: [`atan`](@ref). """ atan_cs_safe diff --git a/src/interpolate.jl b/src/interpolate.jl index dd03847..ef4008f 100644 --- a/src/interpolate.jl +++ b/src/interpolate.jl @@ -5,13 +5,11 @@ Interpolation Methods using OffsetArrays: OffsetVector - """ private function, find index in array where x would be inserted for interpolation between xvec[i] and xvec[i+1] """ function findindex(xvec, x) - n = length(xvec) i = searchsortedlast(real(xvec), real(x)) @@ -39,18 +37,15 @@ end ## ------ Akima Interpolation --------- # struct used internally -struct Akima{TX, TY, TCoeff} - +struct Akima{TX,TY,TCoeff} xdata::TX ydata::TY p0::Vector{TCoeff} p1::Vector{TCoeff} p2::Vector{TCoeff} p3::Vector{TCoeff} - end - """ Akima(xdata, ydata, delta_x=0.0) @@ -72,65 +67,66 @@ Returns an akima spline object (Akima struct). This function only performs const of the spline, not evaluation. This is useful if you want to evaluate the same mesh at multiple different conditions. A convenience method exists below to perform both in one shot. """ -function Akima(xdata::AbstractVector{TFX}, ydata::AbstractVector{TFY}, delta_x=0.0, eps=1e-30) where {TFX, TFY} +function Akima( + xdata::AbstractVector{TFX}, ydata::AbstractVector{TFY}, delta_x=0.0, eps=1e-30 +) where {TFX,TFY} # setup n = length(xdata) TCoeff = promote_type(TFX, TFY) # compute segment slopes - m = OffsetVector(zeros(TCoeff, n+3), -1:n+1) - for i = 1:n-1 - m[i] = (ydata[i+1] - ydata[i]) / (xdata[i+1] - xdata[i]) + m = OffsetVector(zeros(TCoeff, n + 3), -1:(n + 1)) + for i in 1:(n - 1) + m[i] = (ydata[i + 1] - ydata[i]) / (xdata[i + 1] - xdata[i]) end # estimation for end points - m[0] = 2.0*m[1] - m[2] - m[-1] = 2.0*m[0] - m[1] - m[n] = 2.0*m[n-1] - m[n-2] - m[n+1] = 2.0*m[n] - m[n-1] + m[0] = 2.0 * m[1] - m[2] + m[-1] = 2.0 * m[0] - m[1] + m[n] = 2.0 * m[n - 1] - m[n - 2] + m[n + 1] = 2.0 * m[n] - m[n - 1] # slope at points t = zeros(TCoeff, n) - for i = 1:n - m1 = m[i-2] - m2 = m[i-1] + for i in 1:n + m1 = m[i - 2] + m2 = m[i - 1] m3 = m[i] - m4 = m[i+1] + m4 = m[i + 1] w1 = abs_smooth(m4 - m3, delta_x) w2 = abs_smooth(m2 - m1, delta_x) if (real(w1) < eps && real(w2) < eps) - t[i] = 0.5*(m2 + m3) # special case to avoid divide by zero + t[i] = 0.5 * (m2 + m3) # special case to avoid divide by zero else - t[i] = (w1*m2 + w2*m3) / (w1 + w2) + t[i] = (w1 * m2 + w2 * m3) / (w1 + w2) end end # polynomial cofficients - p0 = zeros(TCoeff, n-1) - p1 = zeros(TCoeff, n-1) - p2 = zeros(TCoeff, n-1) - p3 = zeros(TCoeff, n-1) - for i = 1:n-1 - dx = xdata[i+1] - xdata[i] + p0 = zeros(TCoeff, n - 1) + p1 = zeros(TCoeff, n - 1) + p2 = zeros(TCoeff, n - 1) + p3 = zeros(TCoeff, n - 1) + for i in 1:(n - 1) + dx = xdata[i + 1] - xdata[i] t1 = t[i] - t2 = t[i+1] + t2 = t[i + 1] p0[i] = ydata[i] p1[i] = t1 - p2[i] = (3.0*m[i] - 2.0*t1 - t2)/dx - p3[i] = (t1 + t2 - 2.0*m[i])/dx^2 + p2[i] = (3.0 * m[i] - 2.0 * t1 - t2) / dx + p3[i] = (t1 + t2 - 2.0 * m[i]) / dx^2 end return Akima(xdata, ydata, p0, p1, p2, p3) end function (spline::Akima)(x::Number) - j = findindex(spline.xdata, x) # evaluate polynomial dx = x - spline.xdata[j] - y = spline.p0[j] + spline.p1[j]*dx + spline.p2[j]*dx^2 + spline.p3[j]*dx^3 + y = spline.p0[j] + spline.p1[j] * dx + spline.p2[j] * dx^2 + spline.p3[j] * dx^3 return y end @@ -148,8 +144,8 @@ See docstring for Akima for more details. - `xpt::Vector{Float} or ::Float`: the evaluation point(s) - `delta_x::Float=0.0` : the half width of a smoothing interval used for the absolute value function. -- `eps::Float=1e-30` : a cutoff used to avoid dividing by zero in the -weighting function. Default is `1e-30` but this could be raised to machine precision in +- `eps::Float=1e-30` : a cutoff used to avoid dividing by zero in the +weighting function. Default is `1e-30` but this could be raised to machine precision in some cases to improve derivatives. **Returns** @@ -170,12 +166,11 @@ Computes the derivative of an Akima spline at x. - `dydx::Float`: derivative at x using akima spline. """ function derivative(spline::Akima, x) - j = findindex(spline.xdata, x) # evaluate polynomial dx = x - spline.xdata[j] - dydx = spline.p1[j] + 2*spline.p2[j]*dx + 3*spline.p3[j]*dx^2 + dydx = spline.p1[j] + 2 * spline.p2[j] * dx + 3 * spline.p3[j] * dx^2 return dydx end @@ -193,12 +188,11 @@ Computes the second derivative of an Akima spline at x. - `d2ydx2::Float`: second derivative at x using akima spline. """ function second_derivative(spline::Akima, x) - j = findindex(spline.xdata, x) # evaluate polynomial dx = x - spline.xdata[j] - d2ydx2 = 2.0*spline.p2[j] + 6.0*spline.p3[j]*dx + d2ydx2 = 2.0 * spline.p2[j] + 6.0 * spline.p3[j] * dx return d2ydx2 end @@ -217,7 +211,6 @@ Computes the gradient of a Akima spline at x. """ gradient(spline::Akima, x) = derivative.(Ref(spline), x) - # ------ Linear Interpolation --------- """ @@ -234,12 +227,11 @@ Linear interpolation. - `y::Float64`: value at x using linear interpolation """ function linear(xdata, ydata, x::Number) - i = findindex(xdata, x) # lienar interpolation - eta = (x - xdata[i]) / (xdata[i+1] - xdata[i]) - y = ydata[i] + eta*(ydata[i+1] - ydata[i]) + eta = (x - xdata[i]) / (xdata[i + 1] - xdata[i]) + y = ydata[i] + eta * (ydata[i + 1] - ydata[i]) return y end @@ -251,14 +243,12 @@ Convenience function to perform linear interpolation at multiple points. """ linear(xdata, ydata, x::AbstractVector) = linear.(Ref(xdata), Ref(ydata), x) - """ derivative of linear interpolation at `x::Number` """ function derivative(xdata, ydata, x) - i = findindex(xdata, x) - dydx = (ydata[i+1] - ydata[i]) / (xdata[i+1] - xdata[i]) + dydx = (ydata[i + 1] - ydata[i]) / (xdata[i + 1] - xdata[i]) return dydx end @@ -268,7 +258,6 @@ gradient of linear interpolation at `x::Vector` """ gradient(xdata, ydata, x) = derivative.(Ref(xdata), Ref(ydata), x) - # ------- higher order recursive 1D interpolation ------ """ @@ -289,7 +278,6 @@ but it is generalizable to any spline approach and any dimension. - `fhat::Matrix{Float}`: where fhat[i, j] is the estimate function value at xpt[i], ypt[j] """ function interp2d(interp1d, xdata, ydata, fdata, xpt, ypt) - ny = length(ydata) nxpt = length(xpt) nypt = length(ypt) @@ -299,14 +287,13 @@ function interp2d(interp1d, xdata, ydata, fdata, xpt, ypt) yinterp = Array{R}(undef, ny, nxpt) output = Array{R}(undef, nxpt, nypt) - for i = 1:ny + for i in 1:ny yinterp[i, :] .= interp1d(xdata, fdata[:, i], xpt) end - for i = 1:nxpt + for i in 1:nxpt output[i, :] .= interp1d(ydata, yinterp[:, i], ypt) end - return output end @@ -316,7 +303,6 @@ end Same as interp2d, except in three dimension. """ function interp3d(interp1d, xdata, ydata, zdata, fdata, xpt, ypt, zpt) - nz = length(zdata) nxpt = length(xpt) nypt = length(ypt) @@ -326,11 +312,11 @@ function interp3d(interp1d, xdata, ydata, zdata, fdata, xpt, ypt, zpt) zinterp = Array{R}(undef, nz, nxpt, nypt) output = Array{R}(undef, nxpt, nypt, nzpt) - for i = 1:nz + for i in 1:nz zinterp[i, :, :] .= interp2d(interp1d, xdata, ydata, fdata[:, :, i], xpt, ypt) end - for j = 1:nypt - for i = 1:nxpt + for j in 1:nypt + for i in 1:nxpt output[i, j, :] .= interp1d(zdata, zinterp[:, i, j], zpt) end end @@ -338,14 +324,12 @@ function interp3d(interp1d, xdata, ydata, zdata, fdata, xpt, ypt, zpt) return output end - """ interp4d(interp1d, xdata, ydata, zdata, fdata, xpt, ypt, zpt) Same as interp3d, except in four dimensions. """ function interp4d(interp1d, xdata, ydata, zdata, tdata, fdata, xpt, ypt, zpt, tpt) - nt = length(tdata) nxpt = length(xpt) nypt = length(ypt) @@ -356,12 +340,14 @@ function interp4d(interp1d, xdata, ydata, zdata, tdata, fdata, xpt, ypt, zpt, tp tinterp = Array{R}(undef, nt, nxpt, nypt, nzpt) output = Array{R}(undef, nxpt, nypt, nzpt, ntpt) - for i = 1:nt - tinterp[i, :, :, :] .= interp3d(interp1d, xdata, ydata, zdata, fdata[:, :, :, i], xpt, ypt, zpt) + for i in 1:nt + tinterp[i, :, :, :] .= interp3d( + interp1d, xdata, ydata, zdata, fdata[:, :, :, i], xpt, ypt, zpt + ) end - for k = 1:nzpt - for j = 1:nypt - for i = 1:nxpt + for k in 1:nzpt + for j in 1:nypt + for i in 1:nxpt output[i, j, k, :] .= interp1d(tdata, tinterp[:, i, j, k], tpt) end end diff --git a/src/quadrature.jl b/src/quadrature.jl index 022fd79..aff29a9 100644 --- a/src/quadrature.jl +++ b/src/quadrature.jl @@ -1,18 +1,14 @@ # ---- Quadrature Methods --------- - - - """ trapz(x, y) Integrate y w.r.t. x using the trapezoidal method. """ function trapz(x, y) - integral = 0.0 - for i = 1:length(x)-1 - integral += (x[i+1]-x[i])*0.5*(y[i] + y[i+1]) + for i in 1:(length(x) - 1) + integral += (x[i + 1] - x[i]) * 0.5 * (y[i] + y[i + 1]) end return integral -end \ No newline at end of file +end diff --git a/src/roots.jl b/src/roots.jl index 8017452..968ed62 100644 --- a/src/roots.jl +++ b/src/roots.jl @@ -1,8 +1,5 @@ # ------ Root Finding Methods ----- - - - """ brent(f, a, b; args=(), atol=2e-12, rtol=4*eps(), maxiter=100) @@ -23,19 +20,22 @@ - 'fcalls::Int`: number of function calls - 'flag::String`: a convergence/error message. """ -function brent(f, a, b; args=(), atol=2e-12, rtol=4*eps(), maxiter=100) - - xpre = a; xcur = b +function brent(f, a, b; args=(), atol=2e-12, rtol=4 * eps(), maxiter=100) + xpre = a + xcur = b # xblk = 0.0; fblk = 0.0; spre = 0.0; scur = 0.0 error_num = "INPROGRESS" fpre = f(xpre, args...) fcur = f(xcur, args...) - xblk = zero(fpre); fblk = zero(fpre); spre = zero(fpre); scur = zero(fpre) + xblk = zero(fpre) + fblk = zero(fpre) + spre = zero(fpre) + scur = zero(fpre) funcalls = 2 iterations = 0 - - if real(fpre)*real(fcur) > 0 + + if real(fpre) * real(fcur) > 0 error_num = "SIGNERR" return zero(fpre), (iter=iterations, fcalls=funcalls, flag=error_num) end @@ -48,9 +48,9 @@ function brent(f, a, b; args=(), atol=2e-12, rtol=4*eps(), maxiter=100) return xcur, (iter=iterations, fcalls=funcalls, flag=error_num) end - for i = 1:maxiter + for i in 1:maxiter iterations += 1 - if real(fpre)*real(fcur) < 0 + if real(fpre) * real(fcur) < 0 xblk = xpre fblk = fpre spre = scur = xcur - xpre @@ -65,8 +65,8 @@ function brent(f, a, b; args=(), atol=2e-12, rtol=4*eps(), maxiter=100) fblk = fpre end - delta = (atol + rtol*abs_cs_safe(xcur))/2 - sbis = (xblk - xcur)/2 + delta = (atol + rtol * abs_cs_safe(xcur)) / 2 + sbis = (xblk - xcur) / 2 if real(fcur) == zero(real(fcur)) || abs(real(sbis)) < real(delta) error_num = "CONVERGED" return xcur, (iter=iterations, fcalls=funcalls, flag=error_num) @@ -75,14 +75,14 @@ function brent(f, a, b; args=(), atol=2e-12, rtol=4*eps(), maxiter=100) if abs(real(spre)) > real(delta) && abs(real(fcur)) < abs(real(fpre)) if xpre == xblk # interpolate - stry = -fcur*(xcur - xpre)/(fcur - fpre) + stry = -fcur * (xcur - xpre) / (fcur - fpre) else # extrapolate - dpre = (fpre - fcur)/(xpre - xcur) - dblk = (fblk - fcur)/(xblk - xcur) - stry = -fcur*(fblk*dblk - fpre*dpre)/(dblk*dpre*(fblk - fpre)) + dpre = (fpre - fcur) / (xpre - xcur) + dblk = (fblk - fcur) / (xblk - xcur) + stry = -fcur * (fblk * dblk - fpre * dpre) / (dblk * dpre * (fblk - fpre)) end - if 2*abs(real(stry)) < min(abs(real(spre)), 3*abs(real(sbis)) - real(delta)) + if 2 * abs(real(stry)) < min(abs(real(spre)), 3 * abs(real(sbis)) - real(delta)) # good short step spre = scur scur = stry @@ -91,13 +91,14 @@ function brent(f, a, b; args=(), atol=2e-12, rtol=4*eps(), maxiter=100) spre = sbis scur = sbis end - else + else # bisect spre = sbis scur = sbis end - xpre = xcur; fpre = fcur + xpre = xcur + fpre = fcur if abs(real(scur)) > real(delta) xcur += scur else diff --git a/src/smooth.jl b/src/smooth.jl index 277da94..d50af55 100644 --- a/src/smooth.jl +++ b/src/smooth.jl @@ -1,6 +1,5 @@ # ------ Smoothing functions often used to permit continuous differentiability ----- - """ abs_smooth(x, delta_x) @@ -9,13 +8,11 @@ delta_x is the half width of the smoothing interval. Typically usage is with gradient-based optimization. """ function abs_smooth(x, delta_x) - if real(x) < delta_x && real(x) > -delta_x - return x^2/(2*delta_x) + delta_x/2 + return x^2 / (2 * delta_x) + delta_x / 2 else return abs_cs_safe(x) end - end """ @@ -27,7 +24,7 @@ overestimate the maximum function, i.e. `maximum(x) <= ksmax(x, hardness)`. """ function ksmax(x, hardness=50) k = maximum(x) - return 1.0/hardness*log(sum(exp.(hardness*(x.-k)))) .+ k + return 1.0 / hardness * log(sum(exp.(hardness * (x .- k)))) .+ k end """ @@ -53,11 +50,13 @@ function ksmax_adaptive(x, hardness=50; tol=1e-6, smoothing_fraction=0.1) # check if derivative of the KS function wrt hardness is within the tolerance deriv = ksmax_h(x, hardness) target_deriv = -tol - if abs(deriv) > tol * (1-smoothing_fraction) + if abs(deriv) > tol * (1 - smoothing_fraction) # increase hardness by applying Newton's method once in log-log space - dderiv = ksmax_hh(x, hardness)*hardness/deriv - new_hardness = 10^(log10(hardness) + log10(target_deriv/deriv)/dderiv) - hardness = quintic_blend(hardness, new_hardness, -deriv, -target_deriv, smoothing_fraction*tol) + dderiv = ksmax_hh(x, hardness) * hardness / deriv + new_hardness = 10^(log10(hardness) + log10(target_deriv / deriv) / dderiv) + hardness = quintic_blend( + hardness, new_hardness, -deriv, -target_deriv, smoothing_fraction * tol + ) end # use the new hardness to compute ksmax return ksmax(x, hardness) @@ -73,8 +72,9 @@ Newton's method rather than the secant method for increasing hardness values. Some blending is also used to ensure that result is C1 continuous. `smoothing_fraction` controls the smoothness of this blending. """ -ksmin_adaptive(x, hardness=50; tol=1e-6, smoothing_fraction=0.1) = - -ksmax_adaptive(-x, hardness, tol=tol, smoothing_fraction=smoothing_fraction) +function ksmin_adaptive(x, hardness=50; tol=1e-6, smoothing_fraction=0.1) + return -ksmax_adaptive(-x, hardness; tol=tol, smoothing_fraction=smoothing_fraction) +end """ ksmax_h(x, hardness) @@ -84,11 +84,11 @@ function with respect to `hardness`. """ function ksmax_h(x, hardness) k = maximum(x) - tmp1 = exp.(hardness*(x.-k)) - tmp2 = sum((x.-k).*tmp1) + tmp1 = exp.(hardness * (x .- k)) + tmp2 = sum((x .- k) .* tmp1) tmp3 = sum(tmp1) - tmp4 = 1.0/hardness*log(tmp3) - return 1.0/hardness*(tmp2/tmp3 - tmp4) + tmp4 = 1.0 / hardness * log(tmp3) + return 1.0 / hardness * (tmp2 / tmp3 - tmp4) end """ @@ -99,15 +99,15 @@ function with respect to `hardness`. """ function ksmax_hh(x, hardness) k = maximum(x) - tmp1 = exp.(hardness*(x.-k)) - tmp2 = sum((x.-k).*tmp1) - tmp2_h = sum((x.-k).^2 .* tmp1) + tmp1 = exp.(hardness * (x .- k)) + tmp2 = sum((x .- k) .* tmp1) + tmp2_h = sum((x .- k) .^ 2 .* tmp1) tmp3 = sum(tmp1) - tmp3_h = sum((x.-k).*tmp1) - tmp4 = 1.0/hardness*log(tmp3) - tmp4_h = 1.0/hardness*(tmp2/tmp3 - tmp4) - return -1.0/hardness^2*(tmp2/tmp3 - tmp4) + - 1.0/hardness*(tmp2_h/tmp3 - tmp2*tmp3_h/tmp3^2 - tmp4_h) + tmp3_h = sum((x .- k) .* tmp1) + tmp4 = 1.0 / hardness * log(tmp3) + tmp4_h = 1.0 / hardness * (tmp2 / tmp3 - tmp4) + return -1.0 / hardness^2 * (tmp2 / tmp3 - tmp4) + + 1.0 / hardness * (tmp2_h / tmp3 - tmp2 * tmp3_h / tmp3^2 - tmp4_h) end """ @@ -133,8 +133,8 @@ with the transition between the functions located at `xt`. `hardness` controls t sharpness of the transition between the two functions. """ function sigmoid_blend(f1x, f2x, x, xt, hardness=50) - sx = sigmoid(hardness*(x-xt)) - return f1x + sx*(f2x-f1x) + sx = sigmoid(hardness * (x - xt)) + return f1x + sx * (f2x - f1x) end """ @@ -145,11 +145,11 @@ with the transition between the functions at the locations in `xt`. `hardness` c the sharpness of the transition between the functions. """ function sigmoid_blend(fx::Tuple, xt::Tuple, x, hardness=50) - new_fx = sigmoid_blend.(fx[1:end-1], fx[2:end], x, xt, hardness) + new_fx = sigmoid_blend.(fx[1:(end - 1)], fx[2:end], x, xt, hardness) if length(new_fx) == 1 return only(new_fx) else - new_xt = (xt[1:end-1] .+ xt[2:end]) ./ 2 + new_xt = (xt[1:(end - 1)] .+ xt[2:end]) ./ 2 return sigmoid_blend(new_fx, new_xt, x, hardness) end end @@ -167,9 +167,9 @@ function cubic_blend(f1x, f2x, x, xt, delta_x) elseif x >= xt + delta_x return f2x else - xp = (x-xt)/(2*delta_x) + 1/2 - sx = -2*xp^3 + 3*xp^2 - return f1x + sx*(f2x-f1x) + xp = (x - xt) / (2 * delta_x) + 1 / 2 + sx = -2 * xp^3 + 3 * xp^2 + return f1x + sx * (f2x - f1x) end end @@ -181,11 +181,11 @@ with the transition between the functions at the locations in `xt`. `delta_x` is width of the smoothing interval. The resulting function is C1 continuous """ function cubic_blend(fx::Tuple, xt, x, delta_x) - new_fx = cubic_blend.(fx[1:end-1], fx[2:end], x, xt, delta_x) + new_fx = cubic_blend.(fx[1:(end - 1)], fx[2:end], x, xt, delta_x) if length(new_fx) == 1 return only(new_fx) else - new_xt = (xt[1:end-1] .+ xt[2:end]) ./ 2 + new_xt = (xt[1:(end - 1)] .+ xt[2:end]) ./ 2 return cubic_blend(new_fx, new_xt, x, delta_x) end end @@ -203,13 +203,12 @@ function quintic_blend(f1x, f2x, x, xt, delta_x) elseif x >= xt + delta_x return f2x else - xp = (x-xt)/(2*delta_x) + 1/2 - sx = 6*xp^5 - 15*xp^4 + 10*xp^3 - return f1x + sx*(f2x-f1x) + xp = (x - xt) / (2 * delta_x) + 1 / 2 + sx = 6 * xp^5 - 15 * xp^4 + 10 * xp^3 + return f1x + sx * (f2x - f1x) end end - """ quintic_blend(fx::Tuple, xt::Tuple, x, delta_x) @@ -218,11 +217,11 @@ with the transition between the functions at the locations in `xt`. `delta_x` is width of the smoothing interval. The resulting function is C2 continuous """ function quintic_blend(fx::Tuple, xt, x, delta_x) - new_fx = quintic_blend.(fx[1:end-1], fx[2:end], x, xt, delta_x) + new_fx = quintic_blend.(fx[1:(end - 1)], fx[2:end], x, xt, delta_x) if length(new_fx) == 1 return only(new_fx) else - new_xt = (xt[1:end-1] .+ xt[2:end]) ./ 2 + new_xt = (xt[1:(end - 1)] .+ xt[2:end]) ./ 2 return quintic_blend(new_fx, new_xt, x, delta_x) end end From 250bed257e145b2909c8f225532b39db9108de51 Mon Sep 17 00:00:00 2001 From: juddmehr Date: Tue, 25 Jun 2024 11:22:03 -0600 Subject: [PATCH 17/32] update makedocs call to fix warnings --- docs/make.jl | 8 +++++--- src/cs_safe.jl | 16 +++++++--------- 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index ddec09b..bd00706 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,10 +1,12 @@ -using Documenter#, FLOWMath +using Documenter, FLOWMath makedocs(; modules=[FLOWMath], - format=Documenter.HTML(), + format=Documenter.HTML(; + repolink="https://github.com/byuflowlab/FLOWMath.jl/blob/{commit}{path}#L{line}", + edit_link="master", + ), pages=["Home" => "index.md"], - repo="https://github.com/byuflowlab/FLOWMath.jl/blob/{commit}{path}#L{line}", sitename="FLOWMath.jl", authors="Andrew Ning ", checkdocs=:exports, diff --git a/src/cs_safe.jl b/src/cs_safe.jl index 8c3b8b5..e9f85a0 100644 --- a/src/cs_safe.jl +++ b/src/cs_safe.jl @@ -1,6 +1,5 @@ using LinearAlgebra: norm, dot - """ abs_cs_safe(x) @@ -9,7 +8,7 @@ Calculate the absolute value of `x` in a manner compatible with the complex-step abs_cs_safe @inline function abs_cs_safe(x::T) where {T<:Complex} - return x*sign(real(x)) + return x * sign(real(x)) end @inline function abs_cs_safe(x) @@ -39,7 +38,7 @@ Calculate the `p`-norm value of iterable `x` in a manner compatible with the com norm_cs_safe @inline function norm_cs_safe(x::AbstractArray{T}, p::Real=2) where {T<:Complex} - return sum(x.^p)^(1/p) + return sum(x .^ p)^(1 / p) end @inline function norm_cs_safe(x, p::Real=2) @@ -63,7 +62,6 @@ end return dot(a, b) end - """ atan_cs_safe(y, x) @@ -77,11 +75,11 @@ end @inline function atan_cs_safe(y::T, x::T) where {T<:Complex} # Stolen from openmdao/utils/cs_safe.py - # a = np.real(y) - # b = np.imag(y) - # c = np.real(x) - # d = np.imag(x) - # return np.arctan2(a, c) + 1j * (c * b - a * d) / (a**2 + c**2) + # a = np.real(y) + # b = np.imag(y) + # c = np.real(x) + # d = np.imag(x) + # return np.arctan2(a, c) + 1j * (c * b - a * d) / (a**2 + c**2) a = real(y) b = imag(y) c = real(x) From 61527935101ce69eb84e6062069f0f01a639b0ae Mon Sep 17 00:00:00 2001 From: juddmehr Date: Tue, 25 Jun 2024 11:26:18 -0600 Subject: [PATCH 18/32] update pacakges --- docs/Manifest.toml | 2 +- docs/Project.toml | 3 --- 2 files changed, 1 insertion(+), 4 deletions(-) diff --git a/docs/Manifest.toml b/docs/Manifest.toml index 1e9c69a..03f8da3 100644 --- a/docs/Manifest.toml +++ b/docs/Manifest.toml @@ -2,7 +2,7 @@ julia_version = "1.10.0" manifest_format = "2.0" -project_hash = "e3a22c471fb9a9ace33800905d36686a71c62915" +project_hash = "c13332b07c14bfd72e34d786090157c7220cfb4b" [[deps.ANSIColoredPrinters]] git-tree-sha1 = "574baf8110975760d391c710b6341da1afa48d8c" diff --git a/docs/Project.toml b/docs/Project.toml index 61e0fa9..b2a24d9 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,6 +1,3 @@ [deps] Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee" - -[compat] -Documenter = "1.4.1" From 6c8d5faa8b007eea037a644ac80fc0aa22ddf091 Mon Sep 17 00:00:00 2001 From: juddmehr Date: Tue, 25 Jun 2024 11:28:23 -0600 Subject: [PATCH 19/32] update compat --- Project.toml | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index a96b5b9..fed3cc4 100644 --- a/Project.toml +++ b/Project.toml @@ -8,8 +8,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" [compat] -julia = "1.2" -OffsetArrays = "1" +julia = "1.6" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" From fee9e918629d0df33722673eaa0513806a3e4a9c Mon Sep 17 00:00:00 2001 From: juddmehr Date: Tue, 25 Jun 2024 11:29:53 -0600 Subject: [PATCH 20/32] update test setup to current versions --- .github/workflows/test.yaml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/test.yaml b/.github/workflows/test.yaml index 7262c15..ab580ba 100644 --- a/.github/workflows/test.yaml +++ b/.github/workflows/test.yaml @@ -12,13 +12,13 @@ jobs: runs-on: ${{ matrix.os }} strategy: matrix: - julia-version: ['1.6'] + julia-version: ['1.6', '1.10'] julia-arch: [x64] os: [ubuntu-latest, windows-latest, macOS-latest] steps: - - uses: actions/checkout@v1.0.0 - - uses: julia-actions/setup-julia@v1 + - uses: actions/checkout@v4 + - uses: julia-actions/setup-julia@v2 with: version: ${{ matrix.julia-version }} - uses: julia-actions/julia-runtest@v1 From a0e061632ef2735455c15ca55befdfb042bbb354 Mon Sep 17 00:00:00 2001 From: juddmehr Date: Tue, 25 Jun 2024 11:36:15 -0600 Subject: [PATCH 21/32] update workflows --- .github/workflows/docs.yaml | 4 ++-- .github/workflows/test.yaml | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/docs.yaml b/.github/workflows/docs.yaml index 4b26bd7..adcb311 100644 --- a/.github/workflows/docs.yaml +++ b/.github/workflows/docs.yaml @@ -14,7 +14,7 @@ jobs: julia-arch: [x86] os: [ubuntu-latest] steps: - - uses: actions/checkout@v1.0.0 + - uses: actions/checkout@v4 - uses: julia-actions/setup-julia@latest with: version: ${{ matrix.julia-version }} @@ -26,4 +26,4 @@ jobs: env: GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # For authentication with GitHub Actions token DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # For authentication with SSH deploy key - run: julia --project=docs/ docs/make.jl \ No newline at end of file + run: julia --project=docs/ docs/make.jl diff --git a/.github/workflows/test.yaml b/.github/workflows/test.yaml index ab580ba..cec7089 100644 --- a/.github/workflows/test.yaml +++ b/.github/workflows/test.yaml @@ -18,7 +18,7 @@ jobs: steps: - uses: actions/checkout@v4 - - uses: julia-actions/setup-julia@v2 + - uses: julia-actions/setup-julia@latest with: version: ${{ matrix.julia-version }} - uses: julia-actions/julia-runtest@v1 From 494a25854134dd964a0717d9aed827f402a9704c Mon Sep 17 00:00:00 2001 From: juddmehr Date: Tue, 25 Jun 2024 11:39:48 -0600 Subject: [PATCH 22/32] turn off failfast to see if things fail everywhere --- .github/workflows/test.yaml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/test.yaml b/.github/workflows/test.yaml index cec7089..8a2d422 100644 --- a/.github/workflows/test.yaml +++ b/.github/workflows/test.yaml @@ -11,6 +11,7 @@ jobs: test: runs-on: ${{ matrix.os }} strategy: + fail-fast: false matrix: julia-version: ['1.6', '1.10'] julia-arch: [x64] From f96c089932696e05adf30347038bfb1b3d7ef7b4 Mon Sep 17 00:00:00 2001 From: juddmehr Date: Tue, 25 Jun 2024 11:45:30 -0600 Subject: [PATCH 23/32] update docs julia version --- .github/workflows/docs.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/docs.yaml b/.github/workflows/docs.yaml index adcb311..2e02b51 100644 --- a/.github/workflows/docs.yaml +++ b/.github/workflows/docs.yaml @@ -10,7 +10,7 @@ jobs: runs-on: ${{ matrix.os }} strategy: matrix: - julia-version: [1.3.0] + julia-version: ['1.10'] julia-arch: [x86] os: [ubuntu-latest] steps: From a573aca32cdb095a78282f44ff9f04d851babf00 Mon Sep 17 00:00:00 2001 From: juddmehr Date: Tue, 25 Jun 2024 11:50:23 -0600 Subject: [PATCH 24/32] add dependabot stuff for actions --- .github/dependabot.yml | 7 +++++++ 1 file changed, 7 insertions(+) create mode 100644 .github/dependabot.yml diff --git a/.github/dependabot.yml b/.github/dependabot.yml new file mode 100644 index 0000000..b75fee9 --- /dev/null +++ b/.github/dependabot.yml @@ -0,0 +1,7 @@ +# https://docs.github.com/github/administering-a-repository/configuration-options-for-dependency-updates +version: 2 +updates: + - package-ecosystem: "github-actions" + directory: "/" + schedule: + interval: "monthly" From a74813bdbb84947358316eec8312d12a00f28795 Mon Sep 17 00:00:00 2001 From: juddmehr Date: Tue, 25 Jun 2024 11:52:11 -0600 Subject: [PATCH 25/32] bump version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index fed3cc4..e23b7ee 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "FLOWMath" uuid = "6cb5d3fb-0fe8-4cc2-bd89-9fe0b19a99d3" authors = ["Andrew Ning "] -version = "0.3.4" +version = "0.4.0" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" From eb0f9767c7e017d3f5ab9061ebf32e0f20ea5899 Mon Sep 17 00:00:00 2001 From: Judd Mehr Date: Wed, 26 Jun 2024 07:57:02 -0600 Subject: [PATCH 26/32] Update docs.yaml add permissions according to documenter docs --- .github/workflows/docs.yaml | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/.github/workflows/docs.yaml b/.github/workflows/docs.yaml index 2e02b51..f08c1bf 100644 --- a/.github/workflows/docs.yaml +++ b/.github/workflows/docs.yaml @@ -7,6 +7,10 @@ on: jobs: build: + permissions: + contents: write # Required when authenticating with `GITHUB_TOKEN`, not needed when authenticating with SSH deploy keys + pull-requests: read # Required when using `push_preview=true` + statuses: write # Optional, used to report documentation build statuses runs-on: ${{ matrix.os }} strategy: matrix: From 2cd5ec540ee66dabbc33e9c8d1d6f340f794af4e Mon Sep 17 00:00:00 2001 From: Judd Mehr Date: Wed, 26 Jun 2024 08:00:31 -0600 Subject: [PATCH 27/32] Update docs.yaml allow manual workflow run --- .github/workflows/docs.yaml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/docs.yaml b/.github/workflows/docs.yaml index f08c1bf..a6d78c4 100644 --- a/.github/workflows/docs.yaml +++ b/.github/workflows/docs.yaml @@ -4,6 +4,7 @@ on: push: paths: - 'docs/**' + workflow_dispatch: jobs: build: From 8d818ac97bacda2c5a54011387b49bf1fa50160f Mon Sep 17 00:00:00 2001 From: juddmehr Date: Wed, 26 Jun 2024 08:43:35 -0600 Subject: [PATCH 28/32] bump version @JuliaRegistrator register --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index e23b7ee..e9008fe 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "FLOWMath" uuid = "6cb5d3fb-0fe8-4cc2-bd89-9fe0b19a99d3" authors = ["Andrew Ning "] -version = "0.4.0" +version = "0.4.1" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" From 1727c2a44c216e270e2b46131782d593ce35b50f Mon Sep 17 00:00:00 2001 From: Judd Mehr Date: Thu, 27 Jun 2024 10:56:53 -0600 Subject: [PATCH 29/32] @JuliaRegistrator register --- .github/workflows/test.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/test.yaml b/.github/workflows/test.yaml index 8a2d422..a40f4e7 100644 --- a/.github/workflows/test.yaml +++ b/.github/workflows/test.yaml @@ -13,7 +13,7 @@ jobs: strategy: fail-fast: false matrix: - julia-version: ['1.6', '1.10'] + julia-version: ['1.10'] julia-arch: [x64] os: [ubuntu-latest, windows-latest, macOS-latest] From c8021827963466589f2d653ea385525946262a68 Mon Sep 17 00:00:00 2001 From: Judd Mehr Date: Fri, 28 Jun 2024 14:57:08 -0600 Subject: [PATCH 30/32] @JuliaRegistrator register change to run tests on pushes to master to capture changes in workflows. Also allow manually run tests. --- .github/workflows/test.yaml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/test.yaml b/.github/workflows/test.yaml index a40f4e7..7ed0483 100644 --- a/.github/workflows/test.yaml +++ b/.github/workflows/test.yaml @@ -2,10 +2,10 @@ name: Run tests on: push: - paths: - - '**.jl' + branches: + - master pull_request: - + workflow_dispatch: jobs: test: From fb8f6e964f36828a6aca960aa0433af735f6a516 Mon Sep 17 00:00:00 2001 From: juddmehr Date: Mon, 1 Jul 2024 14:36:15 -0600 Subject: [PATCH 31/32] don't skip version numbers --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index e9008fe..e23b7ee 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "FLOWMath" uuid = "6cb5d3fb-0fe8-4cc2-bd89-9fe0b19a99d3" authors = ["Andrew Ning "] -version = "0.4.1" +version = "0.4.0" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" From f041a93bd7b5e16575ca1931552929dc1b3d34d1 Mon Sep 17 00:00:00 2001 From: Judd Mehr Date: Mon, 1 Jul 2024 14:46:53 -0600 Subject: [PATCH 32/32] Update Project.toml Add a compat for OffsetArrays --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index e23b7ee..052515d 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,7 @@ OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" [compat] julia = "1.6" +OffsetArrays="1" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"