optim and rules updated

Author: shahriariravanian

src/homotopy.jl

@@ -77,7 +77,7 @@ function generate_homotopy(eq, x)
 
         S += expand((1 + h₁) * (1 + h₂))
     end
-    
+
     S = simplify(S)
 
     unique([one(x); [equivalent(t, x) for t in terms(S)]])
@@ -91,8 +91,8 @@ function ∂(x)
 end
 
 partial_int_rules = [
-					 # trigonometric functions
-					 @rule 𝛷(sin(~x)) => (cos(~x), ∂(~x))
+                     # trigonometric functions
+                     @rule 𝛷(sin(~x)) => (cos(~x), ∂(~x))
                      @rule 𝛷(cos(~x)) => (sin(~x), ∂(~x))
                      @rule 𝛷(tan(~x)) => (log(cos(~x)), ∂(~x))
                      @rule 𝛷(csc(~x)) => (log(csc(~x) + cot(~x)), ∂(~x))
@@ -119,12 +119,12 @@ partial_int_rules = [
                      @rule 𝛷(^(csch(~x), -1)) => (cosh(~x), ∂(~x))
                      @rule 𝛷(^(sech(~x), -1)) => (sinh(~x), ∂(~x))
                      @rule 𝛷(^(coth(~x), -1)) => (log(cosh(~x)), ∂(~x))
-					 # inverse trigonometric functions
+                     # inverse trigonometric functions
                      @rule 𝛷(asin(~x)) => (~x * asin(~x) + sqrt(1 - ~x * ~x), ∂(~x))
                      @rule 𝛷(acos(~x)) => (~x * acos(~x) + sqrt(1 - ~x * ~x), ∂(~x))
                      @rule 𝛷(atan(~x)) => (~x * atan(~x) + log(~x * ~x + 1), ∂(~x))
-                     @rule 𝛷(acsc(~x)) => (~x * acsc(~x) + atanh(1 - ^(~x,-2)), ∂(~x))
-                     @rule 𝛷(asec(~x)) => (~x * asec(~x) + acosh(~x), ∂(~x))  
+                     @rule 𝛷(acsc(~x)) => (~x * acsc(~x) + atanh(1 - ^(~x, -2)), ∂(~x))
+                     @rule 𝛷(asec(~x)) => (~x * asec(~x) + acosh(~x), ∂(~x))
                      @rule 𝛷(acot(~x)) => (~x * acot(~x) + log(~x * ~x + 1), ∂(~x))
                      # inverse hyperbolic functions
                      @rule 𝛷(asinh(~x)) => (~x * asinh(~x) + sqrt(~x * ~x + 1), ∂(~x))
@@ -134,19 +134,25 @@ partial_int_rules = [
                      @rule 𝛷(asech(~x)) => (asech(~x), ∂(~x))
                      @rule 𝛷(acoth(~x)) => (~x * acot(~x) + log(~x + 1), ∂(~x))
                      # logarithmic and exponential functions
-                     @rule 𝛷(log(~x)) => (~x + ~x * log(~x) + sum(candidate_pow_minus(~x, -1); init = one(~x)), ∂(~x))
+                     @rule 𝛷(log(~x)) => (~x + ~x * log(~x) +
+                                          sum(candidate_pow_minus(~x, -1); init = one(~x)),
+                                          ∂(~x))
                      @rule 𝛷(exp(~x)) => (exp(~x), ∂(~x))
                      @rule 𝛷(^(exp(~x), ~k::is_neg)) => (^(exp(-~x), -~k), ∂(~x))
                      # square-root functions
-                     @rule 𝛷(^(~x, ~k::is_abs_half)) => (sum(candidate_sqrt(~x, ~k); init = one(~x)), 1);
+                     @rule 𝛷(^(~x, ~k::is_abs_half)) => (sum(candidate_sqrt(~x, ~k);
+                                                             init = one(~x)), 1);
                      @rule 𝛷(sqrt(~x)) => (sum(candidate_sqrt(~x, 0.5); init = one(~x)), 1);
-                     @rule 𝛷(^(sqrt(~x), -1)) => 𝛷(^(~x, -0.5))                                                      
-					 # rational functions                                                              
-                     @rule 𝛷(^(~x::is_poly, ~k::is_neg)) => (sum(candidate_pow_minus(~x, ~k); init = one(~x)), 1)
+                     @rule 𝛷(^(sqrt(~x), -1)) => 𝛷(^(~x, -0.5))
+                     # rational functions                                                              
+                     @rule 𝛷(^(~x::is_poly, ~k::is_neg)) => (sum(candidate_pow_minus(~x,
+                                                                                     ~k);
+                                                                 init = one(~x)), 1)
                      @rule 𝛷(^(~x, -1)) => (log(~x), ∂(~x))
-                     @rule 𝛷(^(~x, ~k::is_neg_int)) => (sum(^(~x, i) for i in (~k + 1):-1), ∂(~x))
+                     @rule 𝛷(^(~x, ~k::is_neg_int)) => (sum(^(~x, i) for i in (~k + 1):-1),
+                                                        ∂(~x))
                      @rule 𝛷(1 / ~x) => 𝛷(^(~x, -1))
-                     @rule 𝛷(^(~x, ~k)) => (^(~x, ~k + 1), ∂(~x))                     
+                     @rule 𝛷(^(~x, ~k)) => (^(~x, ~k + 1), ∂(~x))
                      @rule 𝛷(1) => (𝑥, 1)
                      @rule 𝛷(~x) => ((~x + ^(~x, 2)), ∂(~x))]
 

src/integral.jl

@@ -37,7 +37,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, use_optim=false)
+                   homotopy = true, use_optim = false)
     eq = expand(eq)
     eq = apply_div_rule(eq)
 
@@ -236,13 +236,12 @@ end
 """
 function try_integrate(T, eq, x, basis, radius; kwargs...)
     args = Dict(kwargs)
-	use_optim = args[:use_optim]
+    use_optim = args[:use_optim]
     basis = basis[2:end]    # remove 1 from the beginning
-	
-	if use_optim
-		return solve_optim(T, eq, x, basis, radius; kwargs...)
-	else
-		return solve_sparse(T, eq, x, basis, radius; kwargs...)
-	end
-end
 
+    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

@@ -1,96 +1,95 @@
 using Optim
 
-function solve_optim(T, eq, x, basis, radius, rounds=2; kwargs...)
+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)    
+    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 = find_independent_subset(A; abstol)
     A, basis, n = A[:, l], basis[l], sum(l)
     p = rank_basis(A, basis)
-    
-   	# q, ϵ = find_minimizer(A, λ)
-   	
-   	# if ϵ > abstol
-   	# 	return 0, ϵ
-   	# end
-   	
-   	# qa = q
-   	# μ = maximum(abs.(qa))
-    
+
+    # q, ϵ = find_minimizer(A, λ)
+
+    # if ϵ > abstol
+    # 	return 0, ϵ
+    # end
+
+    # qa = q
+    # μ = maximum(abs.(qa))
+
     # for ρ in exp10.(-1:-1:-8)    
-	# 	l = abs.(qa) .> ρ * μ	    
-
-	qm = zeros(n)
-	ϵm = 1e6
-	lm = qm .> 0		
-
-	for i = 1:n
-		l = (1:n .<= i)
-	 	q, ϵ = find_minimizer(A[:, l], λ)
-	 	nz = sum(abs.(q) .> abstol)
-	 	
-	 	# println(i, '\t', ϵ, '\t', nz)
-		
-		if ϵ*nz < ϵm
-			ϵm = ϵ*nz
-			qm = q
-			lm = l			
-		end
-    end  
-   	
-   	if ϵm < abstol
-   		return reconstruct(qm, basis[lm]), ϵm
-   	else
-   		return 0, ϵm
-   	end
+    # 	l = abs.(qa) .> ρ * μ	    
+
+    qm = zeros(n)
+    ϵm = 1e6
+    lm = qm .> 0
+
+    for i in 1:n
+        l = (1:n .<= i)
+        q, ϵ = find_minimizer(A[:, l], λ)
+        nz = sum(abs.(q) .> abstol)
+
+        # println(i, '\t', ϵ, '\t', nz)
+
+        if ϵ * nz < ϵm
+            ϵm = ϵ * nz
+            qm = q
+            lm = l
+        end
+    end
+
+    if ϵm < abstol
+        return reconstruct(qm, basis[lm]), ϵm
+    else
+        return 0, ϵm
+    end
 end
 
 function reconstruct(q, basis)
-	q = nice_parameter.(q)    
-	y = sum(q[i] * expr(basis[i]) for i=1:length(basis))
-	return y
+    q = nice_parameter.(q)
+    y = sum(q[i] * expr(basis[i]) for i in 1:length(basis))
+    return y
 end
 
 # returns a vector of indices of basis elems from the most important to the least
 function rank_basis(A, basis)
-	n, m = size(A)
-	q = A \ ones(n)
-	w = [abs(q[i])*norm(A[:,i]) for i=1:m] 
-	p = sortperm(w; rev=true)
-	return p
+    n, m = size(A)
+    q = A \ ones(n)
+    w = [abs(q[i]) * norm(A[:, i]) for i in 1:m]
+    p = sortperm(w; rev = true)
+    return p
 end
 
 clamp_down(x, η) = abs(x) < η ? 0 : x
 
 function find_minimizer(A, λ)
-	n, m = size(A)
-	
-	f = function(q)
-		q .= clamp_down.(q, maximum(abs.(q))*1e-6)
-		t = A*q .- 1 
-		l = sum(t' * t)
-		l += λ * norm(q, 1)
-		return l
-	end
-	
-	g! = function(G, q)			
-		t = A*q .- 1
-		G .= 2 * real(A' * t)		
-	end
-	
-	q0 = A \ ones(n) # randn(m)	
-	res = Optim.optimize(f, g!, q0, LBFGS())	
-	q = Optim.minimizer(res)
-	ϵ = Optim.minimum(res)
-	
-	return q, ϵ
-end
+    n, m = size(A)
 
+    f = function (q)
+        q .= clamp_down.(q, maximum(abs.(q)) * 1e-6)
+        t = A * q .- 1
+        l = sum(t' * t)
+        l += λ * norm(q, 1)
+        return l
+    end
+
+    g! = function (G, q)
+        t = A * q .- 1
+        G .= 2 * real(A' * t)
+    end
+
+    q0 = A \ ones(n) # randn(m)	
+    res = Optim.optimize(f, g!, q0, LBFGS())
+    q = Optim.minimizer(res)
+    ϵ = Optim.minimum(res)
+
+    return q, ϵ
+end

src/sparse.jl

@@ -1,5 +1,4 @@
 
-
 function solve_sparse(T, eq, x, basis, radius; kwargs...)
     args = Dict(kwargs)
     abstol, opt, complex_plane, verbose = args[:abstol], args[:opt], args[:complex_plane],
@@ -107,19 +106,18 @@ function modify_basis_matrix!(T, A, X, x, eq, basis, radius; abstol = 1e-6)
     end
 end
 
-
 function sparse_fit(T, A, x, basis, opt; abstol = 1e-6)
     # find a linearly independent subset of the basis
     l = find_independent_subset(A; abstol)
     A, basis = A[l, l], basis[l]
     n, m = size(A)
 
     try
-        b = ones((n,1))        
+        b = ones((n, 1))
         q₀ = DataDrivenDiffEq.init(opt, A, b)
         # @views sparse_regression!(q₀, A, permutedims(b)', opt, maxiter = 1000)
         @views sparse_regression!(q₀, A, b, opt, maxiter = 1000)
-        	ϵ = rms(A * q₀ .- b)
+        ϵ = rms(A * q₀ .- b)
         q = nice_parameter.(q₀)
         if sum(iscomplex.(q)) > 2
             return nothing, Inf