Merge pull request #16 from dingraha/master

Make Akima complex stepable

Project.toml

@@ -1,7 +1,7 @@
 name = "FLOWMath"
 uuid = "6cb5d3fb-0fe8-4cc2-bd89-9fe0b19a99d3"
 authors = ["Andrew Ning <[email protected]>"]
-version = "0.3.1"
+version = "0.3.2"
 
 [deps]
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"

src/FLOWMath.jl

@@ -1,5 +1,8 @@
 module FLOWMath
 
+include("cs_safe.jl")
+export abs_cs_safe
+
 include("quadrature.jl")
 export trapz
 

src/cs_safe.jl

@@ -0,0 +1,7 @@
+function abs_cs_safe(x::T) where {T<:Complex}
+    return x*sign(real(x))
+end
+
+function abs_cs_safe(x)
+    return abs(x)
+end

src/interpolate.jl

@@ -13,7 +13,7 @@ between xvec[i] and xvec[i+1]
 function findindex(xvec, x)
 
     n = length(xvec)
-    i = searchsortedlast(xvec, x)
+    i = searchsortedlast(real(xvec), real(x))
 
     # this version allows extrapolation
     if i == 0
@@ -93,7 +93,7 @@ function Akima(xdata::AbstractVector{TFX}, ydata::AbstractVector{TFY}, delta_x=0
         m4 = m[i+1]
         w1 = abs_smooth(m4 - m3, delta_x)
         w2 = abs_smooth(m2 - m1, delta_x)
-        if (w1 < eps && w2 < eps)
+        if (real(w1) < eps && real(w2) < eps)
             t[i] = 0.5*(m2 + m3)  # special case to avoid divide by zero
         else
             t[i] = (w1*m2 + w2*m3) / (w1 + w2)

src/smooth.jl

@@ -10,10 +10,10 @@ Typically usage is with gradient-based optimization.
 """
 function abs_smooth(x, delta_x)
 
-    if x < delta_x && x > -delta_x
+    if real(x) < delta_x && real(x) > -delta_x
         return x^2/(2*delta_x) + delta_x/2
     else
-        return abs(x)
+        return abs_cs_safe(x)
     end
 
 end

test/runtests.jl

@@ -526,11 +526,17 @@ J = ForwardDiff.jacobian(wrapper, y)
 J2 = FiniteDiff.finite_difference_jacobian(wrapper, y)
 @test maximum(abs.(J - J2)) < 1e-6
 
+J2 = FiniteDiff.finite_difference_jacobian(wrapper, y, Val(:complex))
+@test maximum(abs.(J - J2)) < 1e-12
+
 wrapper2(y) = Akima(x, y, 0.1)(xpt)
 J = ForwardDiff.jacobian(wrapper2, x)
 J2 = FiniteDiff.finite_difference_jacobian(wrapper2, x)
 @test maximum(abs.(J - J2)) < 1e-6
 
+J2 = FiniteDiff.finite_difference_jacobian(wrapper2, x, Val(:complex))
+@test maximum(abs.(J - J2)) < 1e-12
+
 # ---------------------------
 
 # ----- 1D linear interpolation ------