From 74625774e2ff997da5407908e0392cf48279a81b Mon Sep 17 00:00:00 2001
From: Yuning <neuromancer43@gmail.com>
Date: Fri, 8 Nov 2024 16:06:22 +0100
Subject: [PATCH] add partial initialize to `RandomFunction`

---
 src/SpinShuttling.jl    | 10 +++++-----
 src/stochastics.jl      | 19 +++++++++++++++----
 test/teststochastics.jl |  3 ++-
 3 files changed, 22 insertions(+), 10 deletions(-)

diff --git a/src/SpinShuttling.jl b/src/SpinShuttling.jl
index cbfc378..99b5617 100644
--- a/src/SpinShuttling.jl
+++ b/src/SpinShuttling.jl
@@ -17,7 +17,7 @@ export ShuttlingModel, OneSpinModel, TwoSpinModel,
     TwoSpinSequentialModel, TwoSpinParallelModel, 
     RandomFunction, CompositeRandomFunction,
     OrnsteinUhlenbeckField, PinkLorentzianField
-export statefidelity, sampling, characteristicfunction, characteristicvalue
+export statefidelity, sampling, initialize!, characteristicfunction, characteristicvalue
 export dephasingmatrix, covariance, covariancematrix
 export W
 
@@ -73,12 +73,12 @@ with arbitrary shuttling path x(t).
 - `x::Function`: Shuttling path
 """
 function OneSpinModel(Ψ::Vector{<:Complex}, T::Real, N::Int,
-    B::GaussianRandomField, x::Function)
+    B::GaussianRandomField, x::Function; initialize::Bool=true)
 
     t = range(0, T, N)
     f(x::Function, t::Real) = (t, x(t)...)
     P = f.(x, t)
-    R = RandomFunction(P, B)
+    R = RandomFunction(P, B, initialize=initialize)
     model = ShuttlingModel(1, Ψ, T, N, B, [x], R)
     return model
 end
@@ -128,13 +128,13 @@ with arbitrary shuttling paths x₁(t), x₂(t).
 - `x₂::Function`: Shuttling path for the second spin
 """
 function TwoSpinModel(Ψ::Vector{<:Complex}, T::Real, N::Int,
-    B::GaussianRandomField, x₁::Function, x₂::Function)
+    B::GaussianRandomField, x₁::Function, x₂::Function; initialize::Bool=true)
 
     X = [x₁, x₂]
     t = range(0, T, N)
     f(x::Function, t::Real) = (t, x(t)...)
     P = vcat(f.(x₁, t), f.(x₂, t))
-    R = RandomFunction(P, B)
+    R = RandomFunction(P, B, initialize=initialize)
     model = ShuttlingModel(2, Ψ, T, N, B, X, R)
     return model
 end
diff --git a/src/stochastics.jl b/src/stochastics.jl
index 72a5266..83a6abe 100644
--- a/src/stochastics.jl
+++ b/src/stochastics.jl
@@ -38,15 +38,20 @@ a Gaussian random process traced from a Gaussian random field.
 - `Σ::Symmetric{<:Real}`: covariance matrices
 - `L::Matrix{<:Real}`: lower triangle matrix of Cholesky decomposition
 """
-struct RandomFunction
+mutable struct RandomFunction
     μ::Vector{<:Real}
     P::Vector{<:Point} # sample trace
     Σ::Symmetric{<:Real} # covariance matrices
-    L::Matrix{<:Real} # Lower triangle matrix of Cholesky decomposition
-    function RandomFunction(P::Vector{<:Point}, process::GaussianRandomField)
+    L::Union{Matrix{<:Real}, Nothing} # Lower triangle matrix of Cholesky decomposition
+    function RandomFunction(P::Vector{<:Point}, process::GaussianRandomField; initialize::Bool=true)
         μ = process.μ isa Function ? process.μ(P) : repeat([process.μ], length(P))
         Σ = covariancematrix(P, process)
-        L = collect(cholesky(Σ).L)
+        if initialize
+            L = collect(cholesky(Σ).L)
+        else
+            L = nothing
+        end
+        
         return new(μ, P, Σ, L)
     end
 
@@ -55,6 +60,12 @@ struct RandomFunction
     end
 end
 
+"""
+Initialize the Cholesky decomposition of the covariance matrix of a random function.
+"""
+function initialize!(R::RandomFunction)
+    R.L = collect(cholesky(R.Σ).L)
+end
 
 """
 Divide the covariance matrix of a direct summed random function into partitions. 
diff --git a/test/teststochastics.jl b/test/teststochastics.jl
index a69923f..60e1298 100644
--- a/test/teststochastics.jl
+++ b/test/teststochastics.jl
@@ -11,7 +11,8 @@ visualize=true
     t=range(0, T, N)
     P=collect(zip(t, v.*t))
     B=OrnsteinUhlenbeckField(0,[κₜ,κₓ],σ)
-    R=RandomFunction(P , B)
+    R=RandomFunction(P , B, initialize=false)
+    initialize!(R)
     @test R() isa Vector{<:Real}
     @test R.Σ isa Symmetric