From 1f281608dc827873745b8b1081c2b4f2884bb04d Mon Sep 17 00:00:00 2001 From: Alberto Mercurio Date: Mon, 10 Feb 2025 05:39:56 +0100 Subject: [PATCH 1/6] Fix instabilities of smesolve --- src/time_evolution/smesolve.jl | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/time_evolution/smesolve.jl b/src/time_evolution/smesolve.jl index 6dab420af..fa7cfed33 100644 --- a/src/time_evolution/smesolve.jl +++ b/src/time_evolution/smesolve.jl @@ -3,7 +3,7 @@ export smesolveProblem, smesolveEnsembleProblem, smesolve _smesolve_generate_state(u, dims) = QuantumObject(vec2mat(u), type = Operator, dims = dims) function _smesolve_update_coeff(u, p, t, op_vec) - return real(dot(u, op_vec)) / 2 #this is Tr[Sn * ρ + ρ * Sn'] + return dot(op_vec, u) #this is Tr[Sn * ρ] end _smesolve_ScalarOperator(op_vec) = @@ -100,11 +100,13 @@ function smesolveProblem( K = get_data(L_evo) Id = I(prod(dims)) + Id_op = IdentityOperator(prod(dims)^2) D_l = map(sc_ops_evo_data) do op - # TODO: Implement the three-argument dot function for SciMLOperators.jl - # Currently, we are assuming a time-independent MatrixOperator + # TODO: # Currently, we are assuming a time-independent MatrixOperator + # Also, the u state may become non-hermitian, so Tr[Sn * ρ + ρ * Sn'] != real(Tr[Sn * ρ]) / 2 op_vec = mat2vec(adjoint(op.A)) - return _spre(op, Id) + _spost(op', Id) + _smesolve_ScalarOperator(op_vec) * IdentityOperator(prod(dims)^2) + op_vec_dag = mat2vec(op.A) + return _spre(op, Id) + _spost(op', Id) + _smesolve_ScalarOperator(op_vec) * Id_op + _smesolve_ScalarOperator(op_vec_dag) * Id_op end D = DiffusionOperator(D_l) From b04ddeb451c31d40a78a38ba9ca6400cbf496bc6 Mon Sep 17 00:00:00 2001 From: Alberto Mercurio Date: Mon, 10 Feb 2025 10:50:29 +0100 Subject: [PATCH 2/6] Fix instabilities of smesolve --- src/time_evolution/smesolve.jl | 3 +++ test/core-test/time_evolution.jl | 2 +- 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/src/time_evolution/smesolve.jl b/src/time_evolution/smesolve.jl index fa7cfed33..ad8abc243 100644 --- a/src/time_evolution/smesolve.jl +++ b/src/time_evolution/smesolve.jl @@ -39,6 +39,7 @@ is the Lindblad superoperator, and ```math \mathcal{H}[\hat{O}] \rho = \hat{O} \rho + \rho \hat{O}^\dagger - \mathrm{Tr}[\hat{O} \rho + \rho \hat{O}^\dagger] \rho, +``` Above, ``\hat{C}_n`` represent the operators related to pure dissipation, while ``\hat{S}_n`` are the measurement operators. The ``dW_n(t)`` term is the real Wiener increment associated to ``\hat{S}_n``. See [Wiseman2009Quantum](@cite) for more details. @@ -162,6 +163,7 @@ is the Lindblad superoperator, and ```math \mathcal{H}[\hat{O}] \rho = \hat{O} \rho + \rho \hat{O}^\dagger - \mathrm{Tr}[\hat{O} \rho + \rho \hat{O}^\dagger] \rho, +``` Above, ``\hat{C}_n`` represent the operators related to pure dissipation, while ``\hat{S}_n`` are the measurement operators. The ``dW_n(t)`` term is the real Wiener increment associated to ``\hat{S}_n``. See [Wiseman2009Quantum](@cite) for more details. @@ -273,6 +275,7 @@ is the Lindblad superoperator, and ```math \mathcal{H}[\hat{O}] \rho = \hat{O} \rho + \rho \hat{O}^\dagger - \mathrm{Tr}[\hat{O} \rho + \rho \hat{O}^\dagger] \rho, +``` Above, ``\hat{C}_n`` represent the operators related to pure dissipation, while ``\hat{S}_n`` are the measurement operators. The ``dW_n(t)`` term is the real Wiener increment associated to ``\hat{S}_n``. See [Wiseman2009Quantum](@cite) for more details. diff --git a/test/core-test/time_evolution.jl b/test/core-test/time_evolution.jl index e09690365..f2c2b6447 100644 --- a/test/core-test/time_evolution.jl +++ b/test/core-test/time_evolution.jl @@ -407,7 +407,7 @@ ntraj = 100, progress_bar = Val(false), ) - @test allocs_tot < 2710000 # TODO: Fix this high number of allocations + @test allocs_tot < 2750000 # TODO: Fix this high number of allocations allocs_tot = @allocations smesolve( H, From e7671188dbb30881547312fe8270b9d5fa765c7e Mon Sep 17 00:00:00 2001 From: Alberto Mercurio Date: Mon, 10 Feb 2025 10:52:01 +0100 Subject: [PATCH 3/6] [no ci] Format code --- src/time_evolution/smesolve.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/time_evolution/smesolve.jl b/src/time_evolution/smesolve.jl index ad8abc243..972a486ca 100644 --- a/src/time_evolution/smesolve.jl +++ b/src/time_evolution/smesolve.jl @@ -107,7 +107,10 @@ function smesolveProblem( # Also, the u state may become non-hermitian, so Tr[Sn * ρ + ρ * Sn'] != real(Tr[Sn * ρ]) / 2 op_vec = mat2vec(adjoint(op.A)) op_vec_dag = mat2vec(op.A) - return _spre(op, Id) + _spost(op', Id) + _smesolve_ScalarOperator(op_vec) * Id_op + _smesolve_ScalarOperator(op_vec_dag) * Id_op + return _spre(op, Id) + + _spost(op', Id) + + _smesolve_ScalarOperator(op_vec) * Id_op + + _smesolve_ScalarOperator(op_vec_dag) * Id_op end D = DiffusionOperator(D_l) From 6924d9ce20fdd71878f90721493f678bd513c8f1 Mon Sep 17 00:00:00 2001 From: Alberto Mercurio Date: Mon, 10 Feb 2025 11:06:00 +0100 Subject: [PATCH 4/6] Fix the definition of the stochastic term --- src/time_evolution/smesolve.jl | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/time_evolution/smesolve.jl b/src/time_evolution/smesolve.jl index 972a486ca..ca1dbc698 100644 --- a/src/time_evolution/smesolve.jl +++ b/src/time_evolution/smesolve.jl @@ -3,7 +3,7 @@ export smesolveProblem, smesolveEnsembleProblem, smesolve _smesolve_generate_state(u, dims) = QuantumObject(vec2mat(u), type = Operator, dims = dims) function _smesolve_update_coeff(u, p, t, op_vec) - return dot(op_vec, u) #this is Tr[Sn * ρ] + return 2 * real(dot(op_vec, u)) #this is Tr[Sn * ρ + ρ * Sn'] end _smesolve_ScalarOperator(op_vec) = @@ -106,11 +106,9 @@ function smesolveProblem( # TODO: # Currently, we are assuming a time-independent MatrixOperator # Also, the u state may become non-hermitian, so Tr[Sn * ρ + ρ * Sn'] != real(Tr[Sn * ρ]) / 2 op_vec = mat2vec(adjoint(op.A)) - op_vec_dag = mat2vec(op.A) return _spre(op, Id) + _spost(op', Id) + - _smesolve_ScalarOperator(op_vec) * Id_op + - _smesolve_ScalarOperator(op_vec_dag) * Id_op + _smesolve_ScalarOperator(op_vec) * Id_op end D = DiffusionOperator(D_l) From fdc7aef1dc2745482529214e0124065dcbb863de Mon Sep 17 00:00:00 2001 From: Alberto Mercurio Date: Mon, 10 Feb 2025 11:08:08 +0100 Subject: [PATCH 5/6] Add changelog --- CHANGELOG.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index e11ac2e8a..0b426e88b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased](https://github.com/qutip/QuantumToolbox.jl/tree/main) - Rename `sparse_to_dense` as `to_dense` and `dense_to_sparse` as `to_sparse`. ([#392]) +- Fix erroneous definition of the stochastic term in `smesolve`. ([#393]) ## [v0.26.0] Release date: 2025-02-09 @@ -118,3 +119,4 @@ Release date: 2024-11-13 [#388]: https://github.com/qutip/QuantumToolbox.jl/issues/388 [#389]: https://github.com/qutip/QuantumToolbox.jl/issues/389 [#392]: https://github.com/qutip/QuantumToolbox.jl/issues/392 +[#393]: https://github.com/qutip/QuantumToolbox.jl/issues/393 From 210a78e20da22675c3122ba0e43bdcc25d713bd5 Mon Sep 17 00:00:00 2001 From: Alberto Mercurio Date: Mon, 10 Feb 2025 11:22:23 +0100 Subject: [PATCH 6/6] Format code --- src/time_evolution/smesolve.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/time_evolution/smesolve.jl b/src/time_evolution/smesolve.jl index ca1dbc698..f571ac4d3 100644 --- a/src/time_evolution/smesolve.jl +++ b/src/time_evolution/smesolve.jl @@ -106,9 +106,7 @@ function smesolveProblem( # TODO: # Currently, we are assuming a time-independent MatrixOperator # Also, the u state may become non-hermitian, so Tr[Sn * ρ + ρ * Sn'] != real(Tr[Sn * ρ]) / 2 op_vec = mat2vec(adjoint(op.A)) - return _spre(op, Id) + - _spost(op', Id) + - _smesolve_ScalarOperator(op_vec) * Id_op + return _spre(op, Id) + _spost(op', Id) + _smesolve_ScalarOperator(op_vec) * Id_op end D = DiffusionOperator(D_l)