diff --git a/src/integration.jl b/src/integration.jl index ebdca25..a5a08ad 100644 --- a/src/integration.jl +++ b/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 \ No newline at end of file diff --git a/src/stochastics.jl b/src/stochastics.jl index 7f5dac0..0183e7c 100644 --- a/src/stochastics.jl +++ b/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 diff --git a/test/Manifest.toml b/test/Manifest.toml index 422e735..ebc40e3 100644 --- a/test/Manifest.toml +++ b/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"] diff --git a/test/runtests.jl b/test/runtests.jl index 8c35a17..8778927 100644 --- a/test/runtests.jl +++ b/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") diff --git a/test/teststochastics.jl b/test/teststochastics.jl index 4981e1f..bf708ed 100644 --- a/test/teststochastics.jl +++ b/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;