Generalizes HHL-based linear differential equation solvers (#35)

* replace `\otime` with `join`

* temp

* genearlise HHL based lin diff eq

Co-authored-by: divyanshu.gupta <[email protected]>

src/QuDiffHHL.jl

@@ -51,7 +51,7 @@ function prepare_init_state(g::Function,alg::LDEMSAlgHHL,tspan::NTuple{2, Float6
         b = bval(alg,h*(i - 1) + tspan[1],h) do t g(t) end
         init_state[Int(sz*(i - 1) + 1):Int(sz*(i))] = h*b
     end
-    return init_state
+    return init_state, N-1, N_t, sz
 end
 
 function array_qudiff(g::Function,alg::LDEMSAlgHHL,tspan::NTuple{2, Float64},h::Float64)
@@ -80,20 +80,35 @@ function array_qudiff(g::Function,alg::LDEMSAlgHHL,tspan::NTuple{2, Float64},h::
     return A_
 end
 
-function DiffEqBase.solve(prob::QuLDEProblem{uType,tType,isinplace, F, P}, alg::LDEMSAlgHHL; dt = (prob.tspan[2]-prob.tspan[1])/100, kwargs...) where {uType,tType,isinplace, F, P}
+function _array_qudiff(A::AbstractMatrix{T}, alg, tspan, dt) where {T}
+    At(t) = A
+    array_qudiff(alg, tspan, dt) do t At(t) end
+end
+
+_array_qudiff(A::Function, alg, tspan, dt) = array_qudiff(A, alg, tspan, dt)
+
+function _prepare_init_state(b::Vector{T}, alg, tspan, x, dt) where {T}
+    bt(t) = b 
+    prepare_init_state(alg, tspan, x, dt) do t bt(t) end
+end
+
+_prepare_init_state(b::Function, alg, tspan, x, dt) = prepare_init_state(b, alg, tspan, x, dt)
+
+function DiffEqBase.solve(prob::QuLDEProblem{uType,tType,isinplace, F, P}, alg::LDEMSAlgHHL; dt = (prob.tspan[2]-prob.tspan[1])/10, kwargs...) where {uType,tType,isinplace, F, P}
     A = prob.A
     b = prob.b
     tspan = prob.tspan
     x = prob.u0
     nreg = alg.nreg
-
-    matx = array_qudiff(alg, tspan, dt) do t A(t) end
-    initstate = prepare_init_state(alg, tspan, x, dt) do t b(t) end
+    matx = _array_qudiff(A, alg, tspan, dt)
+    initstate, N_p, N_t, sz  = _prepare_init_state(b, alg, tspan, x, dt) 
     λ = maximum(eigvals(matx))
     C_value = minimum(eigvals(matx) .|> abs)*0.01;
     matx = 1/(λ*2)*matx
     initstate = initstate*1/(2*λ) |> normalize!
     res = hhlsolve(matx,initstate, nreg, C_value)
     res = res/λ
-    return res
+    N = Int(log2(sz))
+    r = res[(N_p + 1)*2 + 2^N - 1: (N_p + 1)*2 + 2^N + 2*N_t - 2]# To ensure we have a power of 2 dimension for matrix
+    return r
 end;

test/QuDiffHHL_tests.jl

@@ -27,26 +27,21 @@ end
     f(u,p,t) = M*u + v;
     prob = ODEProblem(f, x, tspan)
     qprob = QuLDEProblem(A, b, x, tspan)
-    sol = solve(prob, Tsit5(), dt = h, adaptive = :false)
+    sol = solve(prob, Tsit5(), dt = h, adaptive = false)
     s = vcat(sol.u...)
 
     res = solve(qprob, QuEuler(nreg), dt = h)
-    r = res[(N_t + 1)*2 + 2^N - 1: (N_t + 1)*2 + 2^N + N_t - 3] # range of relevant values in the obtained state.
-    @test isapprox.(s, r, atol = 0.5) |> all
+    @test isapprox.(s, res, atol = 0.5) |> all
 
     res = solve(qprob, QuLeapfrog(nreg), dt = h)
-    r = res[(N_t + 1)*2 + 2^N - 1: (N_t + 1)*2 + 2^N + N_t - 3] # range of relevant values in the obtained state.
-    @test isapprox.(s, r, atol = 0.3) |> all
+    @test isapprox.(s, res, atol = 0.3) |> all
 
     res = solve(qprob, QuAB2(nreg), dt = h)
-    r = res[(N_t + 1)*2 + 2^N - 1: (N_t + 1)*2 + 2^N + N_t - 3] # range of relevant values in the obtained state.
-    @test isapprox.(s, r, atol = 0.3) |> all
+    @test isapprox.(s, res, atol = 0.3) |> all
 
     res = solve(qprob, QuAB3(nreg), dt = h)
-    r = res[(N_t + 1)*2 + 2^N - 1: (N_t + 1)*2 + 2^N + N_t - 3] # range of relevant values in the obtained state.
-    @test isapprox.(s, r, atol = 0.3) |> all
+    @test isapprox.(s, res, atol = 0.3) |> all
 
     res = solve(qprob, QuAB4(nreg),dt = h)
-    r = res[(N_t + 1)*2 + 2^N - 1: (N_t + 1)*2 + 2^N + N_t - 3] # range of relevant values in the obtained state.
-    @test isapprox.(s, r, atol = 0.3) |> all
+    @test isapprox.(s, res, atol = 0.3) |> all
 end;