@@ -24,6 +24,7 @@ Base.signbit(x::SymbolicUtils.Sym{Number, Nothing}) = false
2424 verbose: print a detailed report
2525 complex_plane: generate random test points on the complex plane (if false, the points will be on real axis)
2626 homotopy: use the homotopy algorithm to generate the basis
27+ use_optim: use Optim.jl `minimize` function instead of the STLSQ algorithm (**experimental**)
2728
2829 output:
2930 -------
@@ -36,7 +37,7 @@ function integrate(eq, x = nothing; abstol = 1e-6, num_steps = 2, num_trials = 5
3637 radius = 1.0 ,
3738 show_basis = false , opt = STLSQ (exp .(- 10 : 1 : 0 )), bypass = false ,
3839 symbolic = true , max_basis = 100 , verbose = false , complex_plane = true ,
39- homotopy = true )
40+ homotopy = true , use_optim = false )
4041 eq = expand (eq)
4142 eq = apply_div_rule (eq)
4243
@@ -61,7 +62,7 @@ function integrate(eq, x = nothing; abstol = 1e-6, num_steps = 2, num_trials = 5
6162
6263 s, u, ϵ = integrate_sum (eq, x, l; bypass, abstol, num_trials, num_steps,
6364 radius, show_basis, opt, symbolic,
64- max_basis, verbose, complex_plane, homotopy)
65+ max_basis, verbose, complex_plane, homotopy, use_optim )
6566 return simplify (s), u, ϵ
6667end
6768
@@ -235,164 +236,12 @@ end
235236"""
236237function try_integrate (T, eq, x, basis, radius; kwargs... )
237238 args = Dict (kwargs)
238- abstol, opt, complex_plane, verbose = args[:abstol ], args[:opt ], args[:complex_plane ],
239- args[:verbose ]
240-
239+ use_optim = args[:use_optim ]
241240 basis = basis[2 : end ] # remove 1 from the beginning
242- n = length (basis)
243-
244- # A is an nxn matrix holding the values of the fragments at n random points
245- A = zeros (Complex{T}, (n, n))
246- X = zeros (Complex{T}, n)
247-
248- init_basis_matrix! (T, A, X, x, eq, basis, radius, complex_plane; abstol)
249-
250- y₁, ϵ₁ = sparse_fit (T, A, x, basis, opt; abstol)
251- if ϵ₁ < abstol
252- return y₁, ϵ₁
253- end
254-
255- y₂, ϵ₂ = find_singlet (T, A, basis; abstol)
256- if ϵ₂ < abstol
257- return y₂, ϵ₂
258- end
259-
260- if n < 8 # here, 8 is arbitrary and signifies a small basis
261- y₃, ϵ₃ = find_dense (T, A, basis; abstol)
262- if ϵ₃ < abstol
263- return y₃, ϵ₃
264- end
265- end
266-
267- # moving toward the poles
268- modify_basis_matrix! (T, A, X, x, eq, basis, radius; abstol)
269- y₄, ϵ₄ = sparse_fit (T, A, x, basis, opt; abstol)
270-
271- if ϵ₄ < abstol || ϵ₄ < ϵ₁
272- return y₄, ϵ₄
273- else
274- return y₁, ϵ₁
275- end
276- end
277-
278- function init_basis_matrix! (T, A, X, x, eq, basis, radius, complex_plane; abstol = 1e-6 )
279- n = size (A, 1 )
280- k = 1
281- i = 1
282-
283- eq_fun = fun! (eq, x)
284- Δbasis_fun = deriv_fun! .(basis, x)
285-
286- while k <= n
287- try
288- x₀ = test_point (complex_plane, radius)
289- X[k] = x₀ # move_toward_roots_poles(x₀, x, eq)
290- b₀ = eq_fun (X[k])
291-
292- if is_proper (b₀)
293- for j in 1 : n
294- A[k, j] = Δbasis_fun[j](X[k]) / b₀
295- end
296- if all (is_proper, A[k, :])
297- k += 1
298- end
299- end
300- catch e
301- println (" Error from init_basis_matrix!: "
302- end
303- end
304- end
305-
306- function move_toward_roots_poles (z, x, eq; n = 1 , max_r = 100.0 )
307- eq_fun = fun! (eq, x)
308- Δeq_fun = deriv_fun! (eq, x)
309- is_root = rand () < 0.5
310- z₀ = z
311- for i in 1 : n
312- dz = eq_fun (z) / Δeq_fun (z)
313- if is_root
314- z -= dz
315- else
316- z += dz
317- end
318- if abs (z) > max_r
319- return z₀
320- end
321- end
322- return z
323- end
324-
325- function modify_basis_matrix! (T, A, X, x, eq, basis, radius; abstol = 1e-6 )
326- n = size (A, 1 )
327- eq_fun = fun! (eq, x)
328- Δeq_fun = deriv_fun! (eq, x)
329- Δbasis_fun = deriv_fun! .(basis, x)
330-
331- for k in 1 : n
332- # One Newton iteration toward the poles
333- # note the + sign instead of the usual - in Newton-Raphson's method. This is
334- # because we are moving toward the poles and not zeros.
335-
336- X[k] += eq_fun (X[k]) / Δeq_fun (X[k])
337- b₀ = eq_fun (X[k])
338- for j in 1 : n
339- A[k, j] = Δbasis_fun[j](X[k]) / b₀
340- end
341- end
342- end
343-
344- function sparse_fit (T, A, x, basis, opt; abstol = 1e-6 )
345- n = length (basis)
346- # find a linearly independent subset of the basis
347- l = find_independent_subset (A; abstol)
348- A, basis, n = A[l, l], basis[l], sum (l)
349-
350- try
351- b = ones (n)
352- # q₀ = A \ b
353- q₀ = DataDrivenDiffEq. init (opt, A, b)
354- @views sparse_regression! (q₀, A, permutedims (b)' , opt, maxiter = 1000 )
355- ϵ = rms (A * q₀ - b)
356- q = nice_parameter .(q₀)
357- if sum (iscomplex .(q)) > 2
358- return nothing , Inf
359- end # eliminating complex coefficients
360- return sum (q[i] * expr (basis[i]) for i in 1 : length (basis) if q[i] != 0 ;
361- init = zero (x)),
362- abs (ϵ)
363- catch e
364- println (" Error from sparse_fit"
365- return nothing , Inf
366- end
367- end
368241
369- function find_singlet (T, A, basis; abstol)
370- σ = vec (std (A; dims = 1 ))
371- μ = vec (mean (A; dims = 1 ))
372- l = (σ .< abstol) .* (abs .(μ) .> abstol)
373- if sum (l) == 1
374- k = findfirst (l)
375- return nice_parameter (1 / μ[k]) * expr (basis[k]), σ[k]
242+ if use_optim
243+ return solve_optim (T, eq, x, basis, radius; kwargs... )
376244 else
377- return nothing , Inf
378- end
379- end
380-
381- function find_dense (T, A, basis; abstol = 1e-6 )
382- n = size (A, 1 )
383- b = ones (T, n)
384-
385- try
386- q = A \ b
387- if minimum (abs .(q)) > abstol
388- ϵ = maximum (abs .(A * q .- b))
389- if ϵ < abstol
390- y = sum (nice_parameter .(q) .* expr .(basis))
391- return y, ϵ
392- end
393- end
394- catch e
395- #
245+ return solve_sparse (T, eq, x, basis, radius; kwargs... )
396246 end
397- return nothing , Inf
398247end
0 commit comments