add SymPyPythonCall extension; test extensions (#398)

* add SymPyPythonCall extension; test extensions

* test SymPyPythonCall

* learn to spell

* adjust compat bounds

* punt

* oops

* another tweak

Author: jverzani

Project.toml

@@ -8,19 +8,29 @@ CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
 
+
+[weakdeps]
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+IntervalRootFinding = "d2bf35a9-74e0-55ec-b149-d360ff49b807"
+SymPy = "24249f21-da20-56a4-8eb1-6a02cf4ae2e6"
+SymPyPythonCall = "bc8888f7-b21e-4b7c-a06a-5d9c9496438c"
+
 [compat]
 ChainRulesCore = "1"
 CommonSolve = "0.1, 0.2"
 ForwardDiff = "0.10"
 IntervalRootFinding = "0.5"
 SymPy = "1"
+SymPyPythonCall = "0.1,1"
 Setfield = "0.7, 0.8, 1"
 julia = "1.0"
 
+
 [extensions]
 RootsForwardDiffExt = "ForwardDiff"
 RootsIntervalRootFindingExt = "IntervalRootFinding"
 RootsSymPyExt = "SymPy"
+RootsSymPyPythonCallExt = "SymPyPythonCall"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
@@ -33,14 +43,10 @@ Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 SymPy = "24249f21-da20-56a4-8eb1-6a02cf4ae2e6"
+SymPyPythonCall = "bc8888f7-b21e-4b7c-a06a-5d9c9496438c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [targets]
-test = ["Aqua", "ChainRulesTestUtils", "JSON", "SpecialFunctions", "Statistics", "Test", "BenchmarkTools", "ForwardDiff", "Polynomials", "SymPy", "Unitful", "Zygote"]
-
-[weakdeps]
-ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
-IntervalRootFinding = "d2bf35a9-74e0-55ec-b149-d360ff49b807"
-SymPy = "24249f21-da20-56a4-8eb1-6a02cf4ae2e6"
+test = ["Aqua", "ChainRulesTestUtils", "JSON", "SpecialFunctions", "Statistics", "Test", "BenchmarkTools", "ForwardDiff", "Polynomials", "Unitful", "Zygote", "IntervalRootFinding"]

ext/RootsSymPyPythonCallExt.jl

@@ -0,0 +1,15 @@
+module RootsSymPyPythonCallExt
+
+using Roots
+using SymPyPythonCall
+
+## Allow equations to specify a problem to solve
+function Roots.Callable_Function(M::Roots.AbstractUnivariateZeroMethod, f::SymPyPythonCall.Sym, p=nothing)
+    if f.is_Equality == true
+        f = lhs(f) - rhs(f)
+    end
+    Roots.Callable_Function(M, lambdify(f), p)
+end
+
+
+end

test/runtests.jl

@@ -32,6 +32,8 @@ include("./test_simple.jl")
 
 include("./test_composable.jl")
 VERSION >= v"1.6.0" && include("./test_allocations.jl")
+VERSION >= v"1.9.0" && include("./test_extensions.jl")
+
 
 #include("./runbenchmarks.jl")
 #include("./test_derivative_free_interactive.jl")

test/test_extensions.jl

@@ -0,0 +1,55 @@
+#=
+using SymPy
+@testset "SymPy" begin
+    SymPy.@syms x
+    @test find_zero(cos(x) ~ 1/2, (0, pi/2)) ≈ find_zero(x -> cos(x) - 1/2, (0, pi/2))
+    @test find_zero(1/2 ~ cos(x), (0, pi/2)) ≈ find_zero(x -> 1/2 - cos(x), (0, pi/2))
+    @test find_zero(cos(x) ~ x/2, (0, pi/2)) ≈ find_zero(x -> cos(x) - x/2, (0, pi/2))
+end
+=#
+
+#=
+using SymPyPythonCall
+@testset "SymPythonCall" begin
+    SymPyPythonCall.@syms x
+    @test find_zero(cos(x) ~ 1/2, (0, pi/2)) ≈ find_zero(x -> cos(x) - 1/2, (0, pi/2))
+    @test find_zero(1/2 ~ cos(x), (0, pi/2)) ≈ find_zero(x -> 1/2 - cos(x), (0, pi/2))
+    @test find_zero(cos(x) ~ x/2, (0, pi/2)) ≈ find_zero(x -> cos(x) - x/2, (0, pi/2))
+end
+=#
+
+using ForwardDiff
+@testset "ForwardDiff" begin
+    f(x, p) = x^2 - p
+    Z  = ZeroProblem(f, (0, 1000))
+    F(p) = solve(Z, Roots.Bisection(), p)
+    for p ∈ (3,5,7,11)
+        @test F(p) ≈ sqrt(p)
+        @test ForwardDiff.derivative(F, p) ≈ 1 / (2sqrt(p))
+    end
+
+    # Hessian is *broken*
+    f(x, p) = x^2 - sum(p.^2)
+    Z = ZeroProblem(f, (0, 1000))
+    F(p) = solve(Z, Roots.Bisection(), p)
+    Z = ZeroProblem(f, (0, 1000))
+    F(p) = solve(Z, Roots.Bisection(), p)
+    hess(f, p) = ForwardDiff.jacobian(p -> ForwardDiff.gradient(F, p), p)
+    for p ∈ ([1,2], [1,3], [1,4])
+        @test F(p) ≈ sqrt(sum(p.^2))
+        @test_throws DimensionMismatch ForwardDiff.hessian(F, p)
+        a, b = p
+        n = sqrt(a^2 + b^2)^3
+        @test hess(F, p) ≈ [b^2 -a*b; -a*b a^2]/n
+    end
+
+end
+
+using IntervalRootFinding
+@testset "IntervalRootFinding" begin
+    f(x) = sin(x + sin(x + sin(x)))
+    @test find_zeros(f, (-5, 5)) ≈ [-pi, 0, pi]
+    out = find_zeros(f, -5..5, Roots.Newton())
+    @test sort(out.zeros) ≈ sort([-pi,0,pi])
+    @test isempty(out.unknown)
+end