1.0.1 formatted

Author: shahriariravanian

src/cache.jl

@@ -1,10 +1,10 @@
 const DEBUG_CACHE = true
 
 mutable struct ExprCache
-	eq		# the primary expression
-	f		# compiled eq
-	δeq		# the symbolic derivative of eq
-	δf		# compiled δeq
+    eq::Any# the primary expression
+    f::Any# compiled eq
+    δeq::Any# the symbolic derivative of eq
+    δf::Any# compiled δeq
 end
 
 cache(eq::ExprCache) = eq
@@ -14,46 +14,45 @@ expr(c::ExprCache) = c.eq
 expr(c) = c
 
 function deriv!(c::ExprCache, x)
-	if c.δeq == nothing
-		c.δeq = expand_derivatives(Differential(x)(expr(c)))		
-	end
-	return c.δeq
+    if c.δeq == nothing
+        c.δeq = expand_derivatives(Differential(x)(expr(c)))
+    end
+    return c.δeq
 end
 
 function deriv!(c, x)
-	if DEBUG_CACHE
-		error("ExprCache object expected")
-	end
-	return expand_derivatives(Differential(x)(c))
+    if DEBUG_CACHE
+        error("ExprCache object expected")
+    end
+    return expand_derivatives(Differential(x)(c))
 end
 
 function fun!(c::ExprCache, x)
-	if c.f == nothing
-		c.f = build_function(expr(c), x; expression = false)
-	end
-	return c.f
+    if c.f == nothing
+        c.f = build_function(expr(c), x; expression = false)
+    end
+    return c.f
 end
 
 function fun!(c, x)
-	if DEBUG_CACHE
-		error("ExprCache object expected")
-	end
-	return build_function(c, x; expression = false)
+    if DEBUG_CACHE
+        error("ExprCache object expected")
+    end
+    return build_function(c, x; expression = false)
 end
 
 function deriv_fun!(c::ExprCache, x)
-	if c.δf == nothing
-		c.δf = build_function(deriv!(c, x), x; expression = false)
-	end
-	return c.δf
+    if c.δf == nothing
+        c.δf = build_function(deriv!(c, x), x; expression = false)
+    end
+    return c.δf
 end
 
 function deriv_fun!(c, x)
-	if DEBUG_CACHE
-		error("ExprCache object expected")
-	end
-	return build_function(deriv!(c, x), x; expression = false)
+    if DEBUG_CACHE
+        error("ExprCache object expected")
+    end
+    return build_function(deriv!(c, x), x; expression = false)
 end
 
 Base.show(io::IO, c::ExprCache) = show(expr(c))
-

src/homotopy.jl

@@ -20,7 +20,7 @@ function transformer(eq::SymbolicUtils.Mul, f)
     return prod(transformer(t, f) for t in arguments(eq); init = 1)
 end
 function transformer(eq::SymbolicUtils.Div, f)
-    return transformer(arguments(eq)[1], f) * transformer(arguments(eq)[2] ^ -1, f)
+    return transformer(arguments(eq)[1], f) * transformer(arguments(eq)[2]^-1, f)
 end
 
 function transformer(eq::SymbolicUtils.Pow, f)
@@ -63,7 +63,7 @@ function substitute_x(eq, x, sub)
 end
 
 function generate_homotopy(eq, x)
-	eq = eq isa Num ? eq.val : eq
+    eq = eq isa Num ? eq.val : eq
     q, sub = transform(eq, x)
     S = 0
 
@@ -88,36 +88,32 @@ function ∂(x)
     return isequal(d, 0) ? 1 : d
 end
 
-partial_int_rules = [
-                     @rule 𝛷(sin(~x)) => (cos(~x), ∂(~x))
+partial_int_rules = [@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))
                      @rule 𝛷(sec(~x)) => (log(sec(~x) + tan(~x)), ∂(~x))
                      @rule 𝛷(cot(~x)) => (log(sin(~x)), ∂(~x))
-                     
                      @rule 𝛷(sinh(~x)) => (cosh(~x), ∂(~x))
                      @rule 𝛷(cosh(~x)) => (sinh(~x), ∂(~x))
                      @rule 𝛷(tanh(~x)) => (log(cosh(~x)), ∂(~x))
                      @rule 𝛷(csch(~x)) => (log(tanh(~x / 2)), ∂(~x))
                      @rule 𝛷(sech(~x)) => (atan(sinh(~x)), ∂(~x))
                      @rule 𝛷(coth(~x)) => (log(sinh(~x)), ∂(~x))
-
-					 @rule 𝛷(^(sin(~x), -1)) => (log(csc(~x) + cot(~x)), ∂(~x))
+                     @rule 𝛷(^(sin(~x), -1)) => (log(csc(~x) + cot(~x)), ∂(~x))
                      @rule 𝛷(^(cos(~x), -1)) => (log(sec(~x) + tan(~x)), ∂(~x))
                      @rule 𝛷(^(tan(~x), -1)) => (log(sin(~x)), ∂(~x))
                      @rule 𝛷(^(csc(~x), -1)) => (cos(~x), ∂(~x))
                      @rule 𝛷(^(sec(~x), -1)) => (sin(~x), ∂(~x))
                      @rule 𝛷(^(cot(~x), -1)) => (log(cos(~x)), ∂(~x))
-                     
                      @rule 𝛷(^(sinh(~x), -1)) => (log(tanh(~x / 2)), ∂(~x))
                      @rule 𝛷(^(cosh(~x), -1)) => (atan(sinh(~x)), ∂(~x))
                      @rule 𝛷(^(tanh(~x), -1)) => (log(sinh(~x)), ∂(~x))
                      @rule 𝛷(^(csch(~x), -1)) => (cosh(~x), ∂(~x))
                      @rule 𝛷(^(sech(~x), -1)) => (sinh(~x), ∂(~x))
                      @rule 𝛷(^(coth(~x), -1)) => (log(cosh(~x)), ∂(~x))
-                     
-					 # @rule 𝛷(^(sin(~x), ~k::is_neg)) => 𝛷(^(csc(~x), -~k))
+
+                     # @rule 𝛷(^(sin(~x), ~k::is_neg)) => 𝛷(^(csc(~x), -~k))
                      # @rule 𝛷(^(cos(~x), ~k::is_neg)) => 𝛷(^(sec(~x), -~k))                     
                      # @rule 𝛷(^(tan(~x), ~k::is_neg)) => 𝛷(^(cot(~x), -~k))
                      # @rule 𝛷(^(csc(~x), ~k::is_neg)) => 𝛷(^(sin(~x), -~k))
@@ -129,22 +125,19 @@ partial_int_rules = [
                      # @rule 𝛷(^(csch(~x), ~k::is_neg)) => 𝛷(^(sinh(~x), -~k))
                      # @rule 𝛷(^(sech(~x), ~k::is_neg)) => 𝛷(^(cosh(~x), -~k))
                      # @rule 𝛷(^(coth(~x), ~k::is_neg)) => 𝛷(^(tanh(~x), -~k))
-                     
-                     
+
                      @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) + acosh(~x), ∂(~x))     # needs an abs inside acosh
                      @rule 𝛷(asec(~x)) => (~x * asec(~x) + acosh(~x), ∂(~x))     # needs an abs inside acosh
                      @rule 𝛷(acot(~x)) => (~x * acot(~x) + log(~x * ~x + 1), ∂(~x))
-                     
                      @rule 𝛷(asinh(~x)) => (~x * asinh(~x) + sqrt(~x * ~x + 1), ∂(~x))
                      @rule 𝛷(acosh(~x)) => (~x * acosh(~x) + sqrt(~x * ~x - 1), ∂(~x))
                      @rule 𝛷(atanh(~x)) => (~x * atanh(~x) + log(~x + 1), ∂(~x))
                      @rule 𝛷(acsch(~x)) => (acsch(~x), ∂(~x))
                      @rule 𝛷(asech(~x)) => (asech(~x), ∂(~x))
                      @rule 𝛷(acoth(~x)) => (~x * acot(~x) + log(~x + 1), ∂(~x))
-                     
                      @rule 𝛷(log(~x)) => (~x + ~x * log(~x), ∂(~x))
                      @rule 𝛷(^(~x, ~k::is_abs_half)) => (sum(candidate_sqrt(~x, ~k);
                                                              init = one(~x)), 1);
@@ -154,7 +147,8 @@ partial_int_rules = [
                      @rule 𝛷(sqrt(~x)) => (sum(candidate_sqrt(~x, 0.5); init = one(~x)), 1);
                      @rule 𝛷(^(sqrt(~x), -1)) => 𝛷(^(~x, -0.5))
                      @rule 𝛷(^(~x, -1)) => (log(~x), ∂(~x))
-                     @rule 𝛷(^(~x, ~k::is_neg_int)) => (sum(^(~x, i) for i=~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 𝛷(exp(~x)) => (exp(~x), ∂(~x))

src/integral.jl

@@ -60,8 +60,8 @@ function integrate(eq, x = nothing; abstol = 1e-6, num_steps = 2, num_trials = 5
     # homotopy = homotopy && !bypass
 
     s, u, ϵ = integrate_sum(eq, x, l; bypass, abstol, num_trials, num_steps,
-                         radius, show_basis, opt, symbolic,
-                         max_basis, verbose, complex_plane, homotopy)
+                            radius, show_basis, opt, symbolic,
+                            max_basis, verbose, complex_plane, homotopy)
     return simplify(s), u, ϵ
 end
 
@@ -77,7 +77,7 @@ end
 function integrate_sum(eq, x, l; bypass = false, kwargs...)
     solved = 0
     unsolved = 0
-    	ϵ₀ = 0
+    ϵ₀ = 0
     ts = bypass ? [eq] : terms(eq)
 
     if length(ts) > 1
@@ -88,16 +88,16 @@ function integrate_sum(eq, x, l; bypass = false, kwargs...)
         s, u, ϵ = integrate_term(p, x, l; kwargs...)
         solved += s
         unsolved += u
-        	ϵ₀ = max(ϵ₀, ϵ)
+        ϵ₀ = max(ϵ₀, ϵ)
     end
 
     if !isequal(unsolved, 0)
-		eq = factor_rational(simplify_trigs(unsolved))		
+        eq = factor_rational(simplify_trigs(unsolved))
 
         if !isequal(eq, unsolved)
             eq = expand(eq)
             unsolved = 0
-            	ϵ₀ = 0
+            ϵ₀ = 0
             ts = bypass ? [eq] : terms(eq)
 
             if length(ts) > 1
@@ -110,7 +110,7 @@ function integrate_sum(eq, x, l; bypass = false, kwargs...)
                 s, u, ϵ = integrate_term(p, x, l; kwargs...)
                 solved += s
                 unsolved += u
-				ϵ₀ = max(ϵ₀, ϵ)
+                ϵ₀ = max(ϵ₀, ϵ)
 
                 if !isequal(u, 0)   # premature termination on the first failure
                     return 0, eq, ϵ₀
@@ -147,7 +147,7 @@ function integrate_term(eq, x, l; kwargs...)
         result(l, "Successful", y)
         return y, 0, 0
     end
-    
+
     eq = cache(eq)
     basis = generate_basis(eq, x; homotopy)
 
@@ -161,20 +161,21 @@ function integrate_term(eq, x, l; kwargs...)
     end
 
     # D = Differential(x)
-    	ϵ₀ = Inf
+    ϵ₀ = Inf
     y₀ = 0
 
     # rescue
-    	ϵᵣ = Inf
-    yᵣ = 0      
+    ϵᵣ = Inf
+    yᵣ = 0
 
     for i in 1:num_steps
         if length(basis) > max_basis
             break
         end
 
         if symbolic
-            y, ϵ = try_symbolic(Float64, expr(eq), x, expr.(basis), deriv!.(basis, x); kwargs...)
+            y, ϵ = try_symbolic(Float64, expr(eq), x, expr.(basis), deriv!.(basis, x);
+                                kwargs...)
 
             if !isequal(y, 0) && accept_solution(eq, x, y, radius) < abstol
                 result(l, "Successful symbolic", y)
@@ -188,12 +189,12 @@ function integrate_term(eq, x, l; kwargs...)
             r = radius #*sqrt(2)^j
             y, ϵ = try_integrate(Float64, eq, x, basis, r; kwargs...)
 
-            	ϵ = accept_solution(eq, x, y, r)
+            ϵ = accept_solution(eq, x, y, r)
             if ϵ < abstol
                 result(l, "Successful numeric (attempt $j out of $num_trials)", y)
                 return y, 0, ϵ
             elseif ϵ < ϵᵣ
-                	ϵᵣ = ϵ
+                ϵᵣ = ϵ
                 yᵣ = y
             end
         end
@@ -202,7 +203,7 @@ function integrate_term(eq, x, l; kwargs...)
 
         if i < num_steps
             basis = expand_basis(basis, x)
-            
+
             if show_basis
                 inform(l, "Expanding the basis (|β| = $(length(basis)))", basis)
             elseif verbose
@@ -302,27 +303,27 @@ function init_basis_matrix!(T, A, X, x, eq, basis, radius, complex_plane; abstol
     end
 end
 
-function move_toward_roots_poles(z, x, eq; n=1, max_r=100.0)
-	eq_fun = fun!(eq, x)
+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 = 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	
+    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)    
+    n = size(A, 1)
     eq_fun = fun!(eq, x)
     Δeq_fun = deriv_fun!(eq, x)
     Δbasis_fun = deriv_fun!.(basis, x)
@@ -331,7 +332,7 @@ function modify_basis_matrix!(T, A, X, x, eq, basis, radius; abstol = 1e-6)
         # 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
@@ -351,12 +352,13 @@ function sparse_fit(T, A, x, basis, opt; abstol = 1e-6)
         # q₀ = A \ b
         q₀ = DataDrivenDiffEq.init(opt, A, b)
         @views sparse_regression!(q₀, A, permutedims(b)', opt, maxiter = 1000)
-        	ϵ = rms(A * q₀ - b)
+        ϵ = 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)),
+        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)
@@ -394,4 +396,3 @@ function find_dense(T, A, basis; abstol = 1e-6)
     end
     return nothing, Inf
 end
-

src/numeric_utils.jl

@@ -66,4 +66,3 @@ function fibonacci_spiral(n, i)
     θ, r = mod((i - 1) / ϕ, 1), (i - 1) / n
     sqrt(r) * cis(2π * θ)
 end
-