optim working!

Author: shahriariravanian

Project.toml

@@ -7,6 +7,7 @@ version = "1.0.1"
 DataDrivenDiffEq = "2445eb08-9709-466a-b3fc-47e12bd697a2"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
@@ -27,4 +28,3 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
 test = ["PyCall", "SymPy", "Test"]
-

src/SymbolicNumericIntegration.jl

@@ -18,6 +18,7 @@ include("homotopy.jl")
 
 include("numeric_utils.jl")
 include("sparse.jl")
+include("optim.jl")
 include("integral.jl")
 
 export integrate, generate_basis

src/integral.jl

@@ -36,7 +36,7 @@ function integrate(eq, x = nothing; abstol = 1e-6, num_steps = 2, num_trials = 5
                    radius = 1.0,
                    show_basis = false, opt = STLSQ(exp.(-10:1:0)), bypass = false,
                    symbolic = true, max_basis = 100, verbose = false, complex_plane = true,
-                   homotopy = true)
+                   homotopy = true, use_optim=false)
     eq = expand(eq)
     eq = apply_div_rule(eq)
 
@@ -61,7 +61,7 @@ function integrate(eq, x = nothing; abstol = 1e-6, num_steps = 2, num_trials = 5
 
     s, u, ϵ = integrate_sum(eq, x, l; bypass, abstol, num_trials, num_steps,
                             radius, show_basis, opt, symbolic,
-                            max_basis, verbose, complex_plane, homotopy)
+                            max_basis, verbose, complex_plane, homotopy, use_optim)
     return simplify(s), u, ϵ
 end
 
@@ -234,7 +234,14 @@ end
     integral, error
 """
 function try_integrate(T, eq, x, basis, radius; kwargs...)
+    args = Dict(kwargs)
+	use_optim = args[:use_optim]
     basis = basis[2:end]    # remove 1 from the beginning
-	return solve_sparse(T, eq, x, basis, radius; kwargs...)
+	
+	if use_optim
+		return solve_optim(T, eq, x, basis, radius; kwargs...)
+	else
+		return solve_sparse(T, eq, x, basis, radius; kwargs...)
+	end
 end
 

src/optim.jl

@@ -0,0 +1,64 @@
+using Optim
+
+function solve_optim(T, eq, x, basis, radius, rounds=2; kwargs...)
+    args = Dict(kwargs)
+    abstol, complex_plane, verbose = args[:abstol], args[:complex_plane], args[:verbose]
+    n = length(basis)
+    λ = 1e-9
+    
+    A = zeros(Complex{T}, (2n, n))
+    X = zeros(Complex{T}, 2n)
+    init_basis_matrix!(T, A, X, x, eq, basis, radius, complex_plane; abstol)    
+    # modify_basis_matrix!(T, A, X, x, eq, basis, radius; abstol)
+    
+    l = find_independent_subset(A; abstol) 
+    # l .&= rand(length(l)) .> 0.5   
+    A, basis = A[:, l], basis[l]
+   	q, ϵ = find_minimizer(A, λ)
+   	
+   	if ϵ > abstol
+   		return 0, ϵ
+   	end
+   	
+   	qa = q
+   	μ = maximum(abs.(qa))
+    
+    for ρ in exp10.(-1:-1:-8)    
+		l = abs.(qa) .> ρ * μ
+		q, ϵ = find_minimizer(A[:, l], λ)
+		if ϵ < abstol
+			q = nice_parameter.(q)    
+			basis = basis[l]
+			y = sum(q[i] * expr(basis[i]) for i=1:length(basis))
+			return y, ϵ
+		end
+    end
+   	
+   	return 0, ϵ
+end
+
+
+function find_minimizer(A, λ)
+	n, m = size(A)
+    B = real.(A' * A)
+    b = real.(A' * ones(n))
+	
+	f = function(q)
+		l = λ * sum(abs.(q))		# L1 norm
+		# l = λ * sqrt(sum(q .^ 2))	# L2 norm
+		l += sum(abs2.(A * q .- 1))
+		return l
+	end
+	
+	g! = function(G, q)		
+		G .= 2 * (B * q .- b) .+ λ * sign.(q)
+	end
+	
+	q0 = randn(m)	
+	res = Optim.optimize(f, g!, q0, LBFGS())	
+	q = Optim.minimizer(res)
+	ϵ = Optim.minimum(res)
+	
+	return q, ϵ
+end
+