sparse.jl added

Author: shahriariravanian

src/SymbolicNumericIntegration.jl

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

src/integral.jl

@@ -234,165 +234,7 @@ end
     integral, error
 """
 function try_integrate(T, eq, x, basis, radius; kwargs...)
-    args = Dict(kwargs)
-    abstol, opt, complex_plane, verbose = args[:abstol], args[:opt], args[:complex_plane],
-                                          args[:verbose]
-
     basis = basis[2:end]    # remove 1 from the beginning
-    n = length(basis)
-
-    # A is an nxn matrix holding the values of the fragments at n random points
-    A = zeros(Complex{T}, (n, n))
-    X = zeros(Complex{T}, n)
-
-    init_basis_matrix!(T, A, X, x, eq, basis, radius, complex_plane; abstol)
-
-    y₁, ϵ₁ = sparse_fit(T, A, x, basis, opt; abstol)
-    if ϵ₁ < abstol
-        return y₁, ϵ₁
-    end
-
-    y₂, ϵ₂ = find_singlet(T, A, basis; abstol)
-    if ϵ₂ < abstol
-        return y₂, ϵ₂
-    end
-
-    if n < 8    # here, 8 is arbitrary and signifies a small basis
-        y₃, ϵ₃ = find_dense(T, A, basis; abstol)
-        if ϵ₃ < abstol
-            return y₃, ϵ₃
-        end
-    end
-
-    # moving toward the poles
-    modify_basis_matrix!(T, A, X, x, eq, basis, radius; abstol)
-    y₄, ϵ₄ = sparse_fit(T, A, x, basis, opt; abstol)
-
-    if ϵ₄ < abstol || ϵ₄ < ϵ₁
-        return y₄, ϵ₄
-    else
-        return y₁, ϵ₁
-    end
+	return solve_sparse(T, eq, x, basis, radius; kwargs...)
 end
 
-function init_basis_matrix!(T, A, X, x, eq, basis, radius, complex_plane; abstol = 1e-6)
-    n = size(A, 1)
-    k = 1
-    i = 1
-
-    eq_fun = fun!(eq, x)
-    Δbasis_fun = deriv_fun!.(basis, x)
-
-    while k <= n
-        try
-            x₀ = test_point(complex_plane, radius)
-            X[k] = x₀ # move_toward_roots_poles(x₀, x, eq)
-            b₀ = eq_fun(X[k])
-
-            if is_proper(b₀)
-                for j in 1:n
-                    A[k, j] = Δbasis_fun[j](X[k]) / b₀
-                end
-                if all(is_proper, A[k, :])
-                    k += 1
-                end
-            end
-        catch e
-            println("Error from init_basis_matrix!: ", e)
-        end
-    end
-end
-
-function move_toward_roots_poles(z, x, eq; n = 1, max_r = 100.0)
-    eq_fun = fun!(eq, x)
-    Δeq_fun = deriv_fun!(eq, x)
-    is_root = rand() < 0.5
-    z₀ = z
-    for i in 1:n
-        dz = eq_fun(z) / Δeq_fun(z)
-        if is_root
-            z -= dz
-        else
-            z += dz
-        end
-        if abs(z) > max_r
-            return z₀
-        end
-    end
-    return z
-end
-
-function modify_basis_matrix!(T, A, X, x, eq, basis, radius; abstol = 1e-6)
-    n = size(A, 1)
-    eq_fun = fun!(eq, x)
-    Δeq_fun = deriv_fun!(eq, x)
-    Δbasis_fun = deriv_fun!.(basis, x)
-
-    for k in 1:n
-        # One Newton iteration toward the poles
-        # note the + sign instead of the usual - in Newton-Raphson's method. This is
-        # because we are moving toward the poles and not zeros.
-
-        X[k] += eq_fun(X[k]) / Δeq_fun(X[k])
-        b₀ = eq_fun(X[k])
-        for j in 1:n
-            A[k, j] = Δbasis_fun[j](X[k]) / b₀
-        end
-    end
-end
-
-function sparse_fit(T, A, x, basis, opt; abstol = 1e-6)
-    n = length(basis)
-    # find a linearly independent subset of the basis
-    l = find_independent_subset(A; abstol)
-    A, basis, n = A[l, l], basis[l], sum(l)
-
-    try
-        b = ones(n)
-        # q₀ = A \ b
-        q₀ = DataDrivenDiffEq.init(opt, A, b)
-        @views sparse_regression!(q₀, A, permutedims(b)', opt, maxiter = 1000)
-        ϵ = rms(A * q₀ - b)
-        q = nice_parameter.(q₀)
-        if sum(iscomplex.(q)) > 2
-            return nothing, Inf
-        end   # eliminating complex coefficients
-        return sum(q[i] * expr(basis[i]) for i in 1:length(basis) if q[i] != 0;
-                   init = zero(x)),
-               abs(ϵ)
-    catch e
-        println("Error from sparse_fit", e)
-        return nothing, Inf
-    end
-end
-
-function find_singlet(T, A, basis; abstol)
-    σ = vec(std(A; dims = 1))
-    μ = vec(mean(A; dims = 1))
-    l = (σ .< abstol) .* (abs.(μ) .> abstol)
-    if sum(l) == 1
-        k = findfirst(l)
-        return nice_parameter(1 / μ[k]) * expr(basis[k]), σ[k]
-    else
-        return nothing, Inf
-    end
-end
-
-function find_dense(T, A, basis; abstol = 1e-6)
-    n = size(A, 1)
-    b = ones(T, n)
-
-    try
-        q = A \ b
-        if minimum(abs.(q)) > abstol
-            ϵ = maximum(abs.(A * q .- b))
-            if ϵ < abstol
-                y = sum(nice_parameter.(q) .* expr.(basis))
-                return y, ϵ
-            end
-        end
-    catch e
-        #
-    end
-    return nothing, Inf
-end