Merge pull request #3 from dgan181/HHL_linear_diff

WIP : Forward Euler implementation using HHL

.travis.yml

@@ -1,38 +1,18 @@
-## Documentation: http://docs.travis-ci.com/user/languages/julia/
+# Documentation: http://docs.travis-ci.com/user/languages/julia/
 language: julia
 os:
   - linux
   - osx
 julia:
-  - 0.7
   - 1.0
+  - 1.1
   - nightly
+matrix:
+  allow_failures:
+    - julia: nightly
+  fast_finish: true
 notifications:
   email: false
-git:
-  depth: 99999999
-
-## uncomment the following lines to allow failures on nightly julia
-## (tests will run but not make your overall status red)
-#matrix:
-#  allow_failures:
-#  - julia: nightly
-
-## uncomment and modify the following lines to manually install system packages
-#addons:
-#  apt: # apt-get for linux
-#    packages:
-#    - gfortran
-#before_script: # homebrew for mac
-#  - if [ $TRAVIS_OS_NAME = osx ]; then brew install gcc; fi
-
-## uncomment the following lines to override the default test script
-#script:
-#  - julia -e 'Pkg.clone(pwd()); Pkg.build("QuAlgorithmZoo"); Pkg.test("QuAlgorithmZoo"; coverage=true)'
 after_success:
-  # push coverage results to Coveralls
-  - julia -e 'using Pkg; cd(Pkg.dir("QuAlgorithmZoo")); Pkg.add("Coverage"); using Coverage; Coveralls.submit(Coveralls.process_folder())'
-  # push coverage results to Codecov
-  - julia -e 'using Pkg; cd(Pkg.dir("QuAlgorithmZoo")); Pkg.add("Coverage"); using Coverage; Codecov.submit(Codecov.process_folder())'
-  # add Documenter
-  - julia -e 'using Pkg; cd(Pkg.dir("QuAlgorithmZoo")); include(joinpath("docs", "make.jl"))'
+  - julia -e 'using Pkg; Pkg.add("Coverage"); using Coverage; Codecov.submit(process_folder())'
+  - julia -e 'using Pkg; Pkg.add("Coverage"); using Coverage; Coveralls.submit(process_folder())'

Project.toml

@@ -13,3 +13,7 @@ Yao = "5872b779-8223-5990-8dd0-5abbb0748c8c"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+
+[targets]
+test = ["Test", "OrdinaryDiffEq"]

src/QuAlgorithmZoo.jl

@@ -29,6 +29,7 @@ include("PhaseEstimation.jl")
 include("HHL.jl")
 include("hamiltonian_solvers.jl")
 include("HadamardTest.jl")
+include("lin_diffEq_HHL.jl")
 
 
 end # module

src/lin_diffEq_HHL.jl

@@ -0,0 +1,86 @@
+export Array_QuEuler, prepare_init_state, solve_QuEuler
+
+"""
+    Based on : arxiv.org/abs/1010.2745v2
+
+    QuArray_Euler(N_t,N,h,A)
+    prepare_init_state(b,x,h,N_t)
+
+    x' = Ax + b
+
+    * A - input matrix.
+    * b - input vector.
+    * x - inital vector
+    * N - dimension of b (as a power of 2).
+    * h - step size.
+    * tspan - time span.
+"""
+function prepare_init_state(tspan::Tuple,x::Vector,h::Float64,g::Function)
+  init_state = x;
+  N_t = round(2*(tspan[2] - tspan[1])/h + 3)
+  b = similar(g(1))
+  for i = 1:N_t
+   if i < (N_t+1)/2
+      b = g(h*i + tspan[1])
+      init_state = [init_state;h*b]
+    else
+      init_state = [init_state;zero(b)]
+   end
+
+  end
+  init_state = [init_state;zero(init_state)]
+  init_state
+end
+
+function Array_QuEuler(tspan::Tuple,h::Float64,g::Function)
+  N_t = round(2*(tspan[2] - tspan[1])/h + 3)
+  I_mat = Matrix{Float64}(I, size(g(1)));
+  A_ = I_mat;
+  zero_mat = zero(I_mat);
+  tmp_A = Array{Float64,2}(undef,size(g(1)))
+
+  for j = 2: N_t +1
+    A_ = [A_ zero_mat]
+  end
+
+  for i = 2 : N_t+1
+    tmp_A =  g(i*h + tspan[1])
+    tmp_A = -1*(I_mat + h*tmp_A)
+    tA_ = copy(zero_mat)
+    if i<3
+      tA_ = [tmp_A I_mat]
+    else
+      for j = 1:i-3
+        tA_ = [tA_ zero_mat]
+      end
+      if i < (N_t + 1)/2 + 1
+        tA_ = [tA_ tmp_A I_mat]
+      else
+        tA_ = [tA_ -1*I_mat I_mat]
+      end
+    end
+    if i< N_t+1
+      for j = i+1: N_t+1
+        tA_ = [tA_ zero_mat]
+      end
+    end
+    A_ = [A_;tA_]
+
+  end
+    A_ = [zero(A_) A_;A_' zero(A_)]
+    A_
+end
+
+function solve_QuEuler(A::Function, b::Function, x::Vector,tspan::Tuple, h::Float64)
+
+  mat = Array_QuEuler(tspan,h,A)
+  state = prepare_init_state(tspan,x,h,b)
+  λ = maximum(eigvals(mat))
+  C_value = minimum(eigvals(mat) .|> abs)*0.1;
+  n_reg = 12;
+  mat = 1/(λ*2)*mat
+  state = state*1/(2*λ) |> normalize!
+  res = hhlsolve(mat,state, n_reg, C_value)
+  res = res/λ
+  res
+end

test/lin_diffEq_test.jl

@@ -0,0 +1,32 @@
+using Yao
+using Yao.Intrinsics
+using QuAlgorithmZoo
+using Test, LinearAlgebra
+using OrdinaryDiffEq
+
+function diffEq_problem(nbit::Int)
+    siz = 1<<nbit
+    A = (rand(ComplexF64, siz,siz))
+    A = (A+A')/2
+    b = normalize!(rand(ComplexF64, siz))
+    x = normalize!(rand(ComplexF64, siz))
+    A, b, x
+end
+
+@testset "Linear_differential_equation_Euler_HHL" begin
+N = 1
+h = 0.02
+tspan = (0.0,0.6)
+M, v, x = diffEq_problem(N)
+A(t) = M
+b(t) = v
+res = solve_QuEuler(A, b, x, tspan, h)
+
+f(u,p,t) = M*u + v;
+prob = ODEProblem(f,x,tspan)
+sol = solve(prob,Euler(),dt = 0.02,adaptive = :false)
+s = vcat(sol.u...)
+N_t = round(2*(tspan[2] - tspan[1])/h + 3);
+r = res[Int((N_t+1)*2+2^N+1): Int((N_t+1)*2+2^N + N_t - 1)] # range of relevant values in the obtained state.
+@test isapprox.(s,r,atol = 0.05) |> all
+end