add symmetric integration

src/integration.jl

@@ -1,3 +1,5 @@
+const ArrayOrSubArray{T,N} = Union{Array{T,N}, SubArray{T,N}}
+
 """
 1D Simpson integral of function `f(x)` on a given array of `y=f(x₁), f(x₂)...` with constant
 increment `h`
@@ -8,7 +10,7 @@ increment `h`
 - `method::Symbol=:simpson`: the method of integration 
 ...
 """
-function integrate(y::AbstractVector{<:Real}, h::Real; method::Symbol=:simpson)::Real
+function integrate(y::ArrayOrSubArray{<:Real,1}, h::Real; method::Symbol=:simpson)::Real
     n = length(y)-1
     if method==:simpson
         n % 2 == 0 || error("`y` length (number of intervals) must be odd")
@@ -47,6 +49,10 @@ function integrate(f::Function, a::Real, b::Real, n::Int;
     return integrate(y, h, method=method)
 end
 
+function integrate(f::Function, x_range::Tuple{Real,Real,Int}; method::Symbol=:simpson)::Real
+    return integrate(f, x_range..., method=method)
+end
+
 """
 2D integral of a function `f(x, y)` over the ranges [x_min, x_max] and [y_min, y_max] using Simpson's rule or the Trapezoidal rule applied successively in each dimension.
 
@@ -66,8 +72,8 @@ result = integrate((x, y) -> x * y, (0, 1, 50), (0, 1, 50), method=:simpson)
 """
 function integrate(f::Function, x_range::Tuple{Real,Real,Int}, y_range::Tuple{Real,Real,Int}; 
     method::Symbol=:simpson)::Real
-    g = y-> integrate(x->f(x,y), x_range..., method = method)
-    return integrate(g, y_range..., method = method)
+    g = y-> integrate(x->f(x,y), x_range, method = method)
+    return integrate(g, y_range, method = method)
 end
 
 """
@@ -81,10 +87,29 @@ increments `h_x` and `h_y`.
 - `method::Symbol=:simpson`: the method of integration.
 ...
 """
-function integrate(z::AbstractMatrix{<:Real}, h_x::Real, h_y::Real; method::Symbol=:simpson)::Real
+function integrate(z::ArrayOrSubArray{<:Real,2}, h_x::Real, h_y::Real; method::Symbol=:simpson)::Real
     nrows, ncols = size(z)
     # Integrate along x direction for each y
     integral_x_direction = [integrate((@view z[:, j]), h_x, method=method) for j = 1:ncols]
     # Integrate the result in y direction
     return integrate(integral_x_direction, h_y, method=method)
+end
+
+integrate(z::ArrayOrSubArray{<:Real,2}, h; method::Symbol=:simpson)=integrate(z, h, h; method=method)
+
+"""
+Special methods for the double integral on symmetric matrix with singularity on diagonal entries.
+"""
+function integrate(z::Symmetric, h::Real)::Real
+    n, _ = size(z)
+    @assert n%2 == 1 "The matrix must be of odd size"
+    m=(n+1)÷2
+    _integrate(x::ArrayOrSubArray) = integrate(x, h; method = :trapezoid)
+    # Integrate along upper half trapzoid
+    int_upper = [_integrate((@view z[j:n, j])) for j = 1:m]|> _integrate
+    # Integrate along lower half trapzoid
+    int_lower = [_integrate((@view z[1:j, j])) for j = m:n]|> _integrate
+    # Integrate the duplicated box
+    int_box = _integrate((@view z[m:n, 1:m]))
+    return 2*(int_upper + int_lower - int_box)
 end

src/stochastics.jl

@@ -167,7 +167,7 @@ function characteristicfunction(R::RandomFunction)::Tuple{Vector{<:Real},Vector{
     dt=R.P[2,1]-R.P[1,1]
     N=size(R.Σ,1)
     @assert N%2==1
-    χ(j::Int)=exp.(1im*integrate(R.μ[1:j], dt))*exp.(-integrate((@view R.Σ[1:j,1:j]), dt, dt)/2)
+    χ(j::Int)=exp.(1im*integrate(view(R.μ, 1:j), dt))*exp.(-integrate(view(R.Σ, 1:j,1:j), dt, dt)/2)
     t=R.P[2:2:N-1,1]
     f=[χ(j) for j in 3:2:N] # only for simpson's rule
     return (t,f)
@@ -180,5 +180,5 @@ Using Simpson's rule by default.
 """
 function characteristicvalue(R::RandomFunction)::Number
     dt=R.P[2,1]-R.P[1,1]
-    return exp.(1im*integrate(R.μ, dt))*exp.(-integrate(R.Σ, dt, dt)/2)
+    return exp.(1im*integrate(R.μ, dt))*exp.(-integrate((@view R.Σ[:,:]), dt, dt)/2)
 end

test/Manifest.toml

@@ -20,9 +20,9 @@ version = "1.5.0"
 
 [[deps.Adapt]]
 deps = ["LinearAlgebra", "Requires"]
-git-tree-sha1 = "0fb305e0253fd4e833d486914367a2ee2c2e78d0"
+git-tree-sha1 = "cea4ac3f5b4bc4b3000aa55afb6e5626518948fa"
 uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
-version = "4.0.1"
+version = "4.0.3"
 
     [deps.Adapt.extensions]
     AdaptStaticArraysExt = "StaticArrays"
@@ -35,15 +35,16 @@ uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
 version = "1.1.1"
 
 [[deps.ArrayInterface]]
-deps = ["Adapt", "LinearAlgebra", "Requires", "SparseArrays", "SuiteSparse"]
-git-tree-sha1 = "881e43f1aa014a6f75c8fc0847860e00a1500846"
+deps = ["Adapt", "LinearAlgebra", "SparseArrays", "SuiteSparse"]
+git-tree-sha1 = "44691067188f6bd1b2289552a23e4b7572f4528d"
 uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
-version = "7.8.0"
+version = "7.9.0"
 
     [deps.ArrayInterface.extensions]
     ArrayInterfaceBandedMatricesExt = "BandedMatrices"
     ArrayInterfaceBlockBandedMatricesExt = "BlockBandedMatrices"
     ArrayInterfaceCUDAExt = "CUDA"
+    ArrayInterfaceChainRulesExt = "ChainRules"
     ArrayInterfaceGPUArraysCoreExt = "GPUArraysCore"
     ArrayInterfaceReverseDiffExt = "ReverseDiff"
     ArrayInterfaceStaticArraysCoreExt = "StaticArraysCore"
@@ -53,6 +54,7 @@ version = "7.8.0"
     BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
     BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0"
     CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
+    ChainRules = "082447d4-558c-5d27-93f4-14fc19e9eca2"
     GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
     ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
     StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
@@ -76,10 +78,10 @@ uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0"
 version = "1.0.8+1"
 
 [[deps.Cairo_jll]]
-deps = ["Artifacts", "Bzip2_jll", "CompilerSupportLibraries_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"]
-git-tree-sha1 = "4b859a208b2397a7a623a03449e4636bdb17bcf2"
+deps = ["Artifacts", "Bzip2_jll", "CompilerSupportLibraries_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"]
+git-tree-sha1 = "a4c43f59baa34011e303e76f5c8c91bf58415aaf"
 uuid = "83423d85-b0ee-5818-9007-b63ccbeb887a"
-version = "1.16.1+1"
+version = "1.18.0+1"
 
 [[deps.Calculus]]
 deps = ["LinearAlgebra"]
@@ -144,9 +146,9 @@ version = "1.1.0+0"
 
 [[deps.ConcurrentUtilities]]
 deps = ["Serialization", "Sockets"]
-git-tree-sha1 = "9c4708e3ed2b799e6124b5673a712dda0b596a9b"
+git-tree-sha1 = "6cbbd4d241d7e6579ab354737f4dd95ca43946e1"
 uuid = "f0e56b4a-5159-44fe-b623-3e5288b988bb"
-version = "2.3.1"
+version = "2.4.1"
 
 [[deps.ConstructionBase]]
 deps = ["LinearAlgebra"]
@@ -174,9 +176,9 @@ version = "1.16.0"
 
 [[deps.DataStructures]]
 deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
-git-tree-sha1 = "1fb174f0d48fe7d142e1109a10636bc1d14f5ac2"
+git-tree-sha1 = "0f4b5d62a88d8f59003e43c25a8a90de9eb76317"
 uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
-version = "0.18.17"
+version = "0.18.18"
 
 [[deps.Dates]]
 deps = ["Printf"]
@@ -296,9 +298,9 @@ weakdeps = ["PDMats", "SparseArrays", "Statistics"]
 
 [[deps.FiniteDiff]]
 deps = ["ArrayInterface", "LinearAlgebra", "Requires", "Setfield", "SparseArrays"]
-git-tree-sha1 = "73d1214fec245096717847c62d389a5d2ac86504"
+git-tree-sha1 = "bc0c5092d6caaea112d3c8e3b238d61563c58d5f"
 uuid = "6a86dc24-6348-571c-b903-95158fe2bd41"
-version = "2.22.0"
+version = "2.23.0"
 
     [deps.FiniteDiff.extensions]
     FiniteDiffBandedMatricesExt = "BandedMatrices"
@@ -363,15 +365,15 @@ version = "3.3.9+0"
 
 [[deps.GR]]
 deps = ["Artifacts", "Base64", "DelimitedFiles", "Downloads", "GR_jll", "HTTP", "JSON", "Libdl", "LinearAlgebra", "Pkg", "Preferences", "Printf", "Random", "Serialization", "Sockets", "TOML", "Tar", "Test", "UUIDs", "p7zip_jll"]
-git-tree-sha1 = "3458564589be207fa6a77dbbf8b97674c9836aab"
+git-tree-sha1 = "3437ade7073682993e092ca570ad68a2aba26983"
 uuid = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
-version = "0.73.2"
+version = "0.73.3"
 
 [[deps.GR_jll]]
 deps = ["Artifacts", "Bzip2_jll", "Cairo_jll", "FFMPEG_jll", "Fontconfig_jll", "FreeType2_jll", "GLFW_jll", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Libtiff_jll", "Pixman_jll", "Qt6Base_jll", "Zlib_jll", "libpng_jll"]
-git-tree-sha1 = "77f81da2964cc9fa7c0127f941e8bce37f7f1d70"
+git-tree-sha1 = "a96d5c713e6aa28c242b0d25c1347e258d6541ab"
 uuid = "d2c73de3-f751-5644-a686-071e5b155ba9"
-version = "0.73.2+0"
+version = "0.73.3+0"
 
 [[deps.Gettext_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "XML2_jll"]
@@ -381,9 +383,9 @@ version = "0.21.0+0"
 
 [[deps.Glib_jll]]
 deps = ["Artifacts", "Gettext_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Libiconv_jll", "Libmount_jll", "PCRE2_jll", "Zlib_jll"]
-git-tree-sha1 = "e94c92c7bf4819685eb80186d51c43e71d4afa17"
+git-tree-sha1 = "359a1ba2e320790ddbe4ee8b4d54a305c0ea2aff"
 uuid = "7746bdde-850d-59dc-9ae8-88ece973131d"
-version = "2.76.5+0"
+version = "2.80.0+0"
 
 [[deps.Graphite2_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
@@ -558,10 +560,10 @@ uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531"
 version = "1.17.0+0"
 
 [[deps.Libmount_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "9c30530bf0effd46e15e0fdcf2b8636e78cbbd73"
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "dae976433497a2f841baadea93d27e68f1a12a97"
 uuid = "4b2f31a3-9ecc-558c-b454-b3730dcb73e9"
-version = "2.35.0+0"
+version = "2.39.3+0"
 
 [[deps.Libtiff_jll]]
 deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "LERC_jll", "Libdl", "XZ_jll", "Zlib_jll", "Zstd_jll"]
@@ -571,9 +573,9 @@ version = "4.5.1+1"
 
 [[deps.Libuuid_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "e5edc369a598dfde567269dc6add5812cfa10cd5"
+git-tree-sha1 = "0a04a1318df1bf510beb2562cf90fb0c386f58c4"
 uuid = "38a345b3-de98-5d2b-a5d3-14cd9215e700"
-version = "2.39.3+0"
+version = "2.39.3+1"
 
 [[deps.LinearAlgebra]]
 deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"]
@@ -757,15 +759,15 @@ version = "3.1.0"
 
 [[deps.PlotUtils]]
 deps = ["ColorSchemes", "Colors", "Dates", "PrecompileTools", "Printf", "Random", "Reexport", "Statistics"]
-git-tree-sha1 = "862942baf5663da528f66d24996eb6da85218e76"
+git-tree-sha1 = "7b1a9df27f072ac4c9c7cbe5efb198489258d1f5"
 uuid = "995b91a9-d308-5afd-9ec6-746e21dbc043"
-version = "1.4.0"
+version = "1.4.1"
 
 [[deps.Plots]]
 deps = ["Base64", "Contour", "Dates", "Downloads", "FFMPEG", "FixedPointNumbers", "GR", "JLFzf", "JSON", "LaTeXStrings", "Latexify", "LinearAlgebra", "Measures", "NaNMath", "Pkg", "PlotThemes", "PlotUtils", "PrecompileTools", "Printf", "REPL", "Random", "RecipesBase", "RecipesPipeline", "Reexport", "RelocatableFolders", "Requires", "Scratch", "Showoff", "SparseArrays", "Statistics", "StatsBase", "UUIDs", "UnicodeFun", "UnitfulLatexify", "Unzip"]
-git-tree-sha1 = "c4fa93d7d66acad8f6f4ff439576da9d2e890ee0"
+git-tree-sha1 = "3c403c6590dd93b36752634115e20137e79ab4df"
 uuid = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
-version = "1.40.1"
+version = "1.40.2"
 
     [deps.Plots.extensions]
     FileIOExt = "FileIO"
@@ -783,9 +785,9 @@ version = "1.40.1"
 
 [[deps.PrecompileTools]]
 deps = ["Preferences"]
-git-tree-sha1 = "03b4c25b43cb84cee5c90aa9b5ea0a78fd848d2f"
+git-tree-sha1 = "5aa36f7049a63a1528fe8f7c3f2113413ffd4e1f"
 uuid = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
-version = "1.2.0"
+version = "1.2.1"
 
 [[deps.Preferences]]
 deps = ["TOML"]
@@ -1063,9 +1065,9 @@ version = "1.1.34+0"
 
 [[deps.XZ_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "37195dcb94a5970397ad425b95a9a26d0befce3a"
+git-tree-sha1 = "31c421e5516a6248dfb22c194519e37effbf1f30"
 uuid = "ffd25f8a-64ca-5728-b0f7-c24cf3aae800"
-version = "5.6.0+0"
+version = "5.6.1+0"
 
 [[deps.Xorg_libICE_jll]]
 deps = ["Libdl", "Pkg"]
@@ -1277,9 +1279,9 @@ version = "1.18.0+0"
 
 [[deps.libpng_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Zlib_jll"]
-git-tree-sha1 = "1ea2ebe8ffa31f9c324e8c1d6e86b4165b84a024"
+git-tree-sha1 = "d7015d2e18a5fd9a4f47de711837e980519781a4"
 uuid = "b53b4c65-9356-5827-b1ea-8c7a1a84506f"
-version = "1.6.43+0"
+version = "1.6.43+1"
 
 [[deps.libvorbis_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Ogg_jll", "Pkg"]

test/runtests.jl

@@ -1,6 +1,7 @@
 using Test
 using SpinShuttling
 
-include("testquadrature.jl")
-include("testfidelity.jl")
-include("teststochastics.jl")
+# include("testquadrature.jl")
+# include("testfidelity.jl")
+# include("teststochastics.jl")
+include("testsymmetric.jl")

test/teststochastics.jl

@@ -1,6 +1,5 @@
-using SpinShuttling: covariancematrix, CompositeRandomFunction, 
-Symmetric, Cholesky, covariancepartition, ishermitian, issymmetric, integrate, std,
-characteristicfunction
+using SpinShuttling: covariancepartition
+using LinearAlgebra: Symmetric, Cholesky, ishermitian, issymmetric, integrate
 using Plots
 using FFTW
 using LsqFit
@@ -49,36 +48,75 @@ visualize=true
     end
 end
 
-## 
+# 
 @testset "trapezoid vs simpson for covariance matrix" begin
-    T=400; L=10; σ = sqrt(2) / 20; N=1001; κₜ=1/20;κₓ=1/0.1;
-    v=L/T;
-    t=range(0, T, N)
-    P=hcat(t, v.*t)
-    M=20
+    L=10; σ = sqrt(2) / 20; N=501; κₜ=1/20;κₓ=1/0.1;
+    v=20;
+    B=OrnsteinUhlenbeckField(0,[κₜ,κₓ],σ)
+
+    M=30
     err1=zeros(M)
     err2=zeros(M)
     i=1
-    for T in T/2 .*(1 .+rand(M))
-        B=OrnsteinUhlenbeckField(0,[κₜ,κₓ],σ)
+    for T in range(1, 50, length=M)
+        t=range(0, T, N)
+        P=hcat(t, v.*t)
         R=RandomFunction(P , B)
         dt=T/N
         f1=exp(-integrate(R.Σ[:,:], dt, dt, method=:trapezoid)/2) 
         f2=exp(-integrate(R.Σ[:,:], dt, dt, method=:simpson)/2)
-        f3=1/2(1+φ(T,L,B))
+        @test isapprox(f1,f2,rtol=1e-2) 
+        f3=φ(T,L,B)
         err1[i]=abs(f1-f3)
         err2[i]=abs(f2-f3)
         i+=1
     end
-    println("std 1st order:", std(err1))
-    println("std 2st order:", std(err2))
+    println("mean 1st order:", mean(err1))
+    println("mean 2nd order:", mean(err2))
     if visualize
-        fig=plot(err1, xlabel="sample", ylabel="error", label="trapezoid")
+        fig=plot(err1, xlabel="T", ylabel="error", label="trapezoid")
         plot!(err2, label="simpson")
         display(fig)
     end
 end
 
+##
+
+## 
+@testset "symmetric integration for covariance matrix" begin
+    σ = sqrt(2) / 20; N=21; κₜ=1;κₓ=0.01;
+    v=2;
+    B=OrnsteinUhlenbeckField(0,[κₜ,κₓ],σ)
+
+    M=5
+    err1=zeros(M)
+    err1_sym=zeros(M)
+    err2=zeros(M)
+    i=1
+    T_range=range(10, 20, length=M)
+    for T in T_range
+        t=range(0, T, N)
+        P=hcat(t, v.*t)
+        R=RandomFunction(P , B)
+        dt=T/N
+        f1=exp(-integrate(R.Σ[:,:], dt, dt, method=:trapezoid)/2) 
+        f1_sym=exp(-integrate(R.Σ, dt)/2)
+        f2=exp(-integrate(R.Σ[:,:], dt, dt, method=:simpson)/2)
+        f3=φ(T,v*T,B)
+        err1[i]=abs(f1-f3)/f3
+        err1_sym[i]=abs(f1_sym-f3)/f3
+        err2[i]=abs(f2-f3)/f3
+        i+=1
+    end
+    println("mean 1st order:", mean(err1))
+    println("mean 1st order symmetric:", mean(err1_sym))
+    println("mean 2nd order:", mean(err2))
+    @test mean(err1) < 1e-2
+    @test mean(err1_sym) < 1e-2
+    @test mean(err2) < 1e-2
+end
+
+
 ##
 @testset "test 1/f noise" begin
     σ = sqrt(2) / 20; M = 1000; N=1001; κₜ=1/20;κₓ=0.1;