1# This file is part of QuTiP: Quantum Toolbox in Python.
2#
3#    Copyright (c) 2011 and later, Paul D. Nation and Robert J. Johansson.
4#    All rights reserved.
5#
6#    Redistribution and use in source and binary forms, with or without
7#    modification, are permitted provided that the following conditions are
8#    met:
9#
10#    1. Redistributions of source code must retain the above copyright notice,
11#       this list of conditions and the following disclaimer.
12#
13#    2. Redistributions in binary form must reproduce the above copyright
14#       notice, this list of conditions and the following disclaimer in the
15#       documentation and/or other materials provided with the distribution.
16#
17#    3. Neither the name of the QuTiP: Quantum Toolbox in Python nor the names
18#       of its contributors may be used to endorse or promote products derived
19#       from this software without specific prior written permission.
20#
21#    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
22#    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
23#    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
24#    PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
25#    HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
26#    SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
27#    LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
28#    DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
29#    THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
30#    (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
31#    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
32###############################################################################
33
34import numpy as np
35from numpy.random import rand
36import scipy.linalg as la
37from numpy.testing import assert_, assert_equal, run_module_suite
38import scipy
39
40from qutip import (rand_dm, rand_unitary, spre, spost, vector_to_operator,
41                   operator_to_vector, mat2vec, vec2mat, vec2mat_index,
42                   mat2vec_index, tensor, sprepost, to_super, reshuffle,
43                   identity, destroy, create, qeye, QobjEvo, Qobj)
44from qutip.superoperator import liouvillian, liouvillian_ref, \
45                                lindblad_dissipator
46
47
48def f(t, args):
49    return t*(1-0.5j)
50
51
52class TestMatVec:
53    """
54    A test class for the QuTiP function for matrix/vector conversion.
55    """
56
57    def testOperatorVector(self):
58        """
59        Superoperator: Operator - vector - operator conversion.
60        """
61        N = 3
62        rho1 = rand_dm(N)
63        rho2 = vector_to_operator(operator_to_vector(rho1))
64
65        assert_((rho1 - rho2).norm() < 1e-8)
66
67    def testOperatorVectorTensor(self):
68        """
69        Superoperator: Operator - vector - operator conversion with a tensor product state.
70        """
71        Na = 3
72        Nb = 2
73        rhoa = rand_dm(Na)
74        rhob = rand_dm(Nb)
75        rho1 = tensor(rhoa, rhob)
76        rho2 = vector_to_operator(operator_to_vector(rho1))
77
78        assert_((rho1 - rho2).norm() < 1e-8)
79
80    def testOperatorVectorNotSquare(self):
81        """
82        Superoperator: Operator - vector - operator conversion for non-square matrix.
83        """
84        op1 = Qobj(np.random.rand(6).reshape((3, 2)))
85        op2 = vector_to_operator(operator_to_vector(op1))
86
87        assert_((op1 - op2).norm() < 1e-8)
88
89    def testOperatorSpreAppl(self):
90        """
91        Superoperator: apply operator and superoperator from left (spre)
92        """
93        N = 3
94        rho = rand_dm(N)
95        U = rand_unitary(N)
96
97        rho1 = U * rho
98        rho2_vec = spre(U) * operator_to_vector(rho)
99        rho2 = vector_to_operator(rho2_vec)
100
101        assert_((rho1 - rho2).norm() < 1e-8)
102
103    def testOperatorSpostAppl(self):
104        """
105        Superoperator: apply operator and superoperator from right (spost)
106        """
107        N = 3
108        rho = rand_dm(N)
109        U = rand_unitary(N)
110
111        rho1 = rho * U
112        rho2_vec = spost(U) * operator_to_vector(rho)
113        rho2 = vector_to_operator(rho2_vec)
114
115        assert_((rho1 - rho2).norm() < 1e-8)
116
117    def testOperatorUnitaryTransform(self):
118        """
119        Superoperator: Unitary transformation with operators and superoperators
120        """
121        N = 3
122        rho = rand_dm(N)
123        U = rand_unitary(N)
124
125        rho1 = U * rho * U.dag()
126        rho2_vec = spre(U) * spost(U.dag()) * operator_to_vector(rho)
127        rho2 = vector_to_operator(rho2_vec)
128
129        assert_((rho1 - rho2).norm() < 1e-8)
130
131    def testMatrixVecMat(self):
132        """
133        Superoperator: Conversion matrix to vector to matrix
134        """
135        M = rand(10, 10)
136        V = mat2vec(M)
137        M2 = vec2mat(V)
138        assert_(la.norm(M - M2) == 0.0)
139
140    def testVecMatVec(self):
141        """
142        Superoperator: Conversion vector to matrix to vector
143        """
144        V = rand(100)     # a row vector
145        M = vec2mat(V)
146        V2 = mat2vec(M).T  # mat2vec returns a column vector
147        assert_(la.norm(V - V2) == 0.0)
148
149    def testVecMatIndexConversion(self):
150        """
151        Superoperator: Conversion between matrix and vector indices
152        """
153        N = 10
154        for I in range(N * N):
155            i, j = vec2mat_index(N, I)
156            I2 = mat2vec_index(N, i, j)
157            assert_(I == I2)
158
159    def testVecMatIndexCompability(self):
160        """
161        Superoperator: Compatibility between matrix/vector and
162        corresponding index conversions.
163        """
164        N = 10
165        M = rand(N, N)
166        V = mat2vec(M)
167        for I in range(N * N):
168            i, j = vec2mat_index(N, I)
169            assert_(V[I][0] == M[i, j])
170
171    def test_reshuffle(self):
172        U1 = rand_unitary(2)
173        U2 = rand_unitary(3)
174        U3 = rand_unitary(4)
175
176        U = tensor(U1, U2, U3)
177        S = to_super(U)
178        S_col = reshuffle(S)
179
180        assert_equal(S_col.dims[0], [[2, 2], [3, 3], [4, 4]])
181
182        assert_(reshuffle(S_col) == S)
183
184    def test_sprepost(self):
185        U1 = rand_unitary(3)
186        U2 = rand_unitary(3)
187
188        S1 = spre(U1) * spost(U2)
189        S2 = sprepost(U1, U2)
190
191        assert_(S1 == S2)
192
193    def testLiouvillianImplem(self):
194        """
195        Superoperator: Randomized comparison of standard and reference
196        Liouvillian functions.
197        """
198        N1 = 3
199        N2 = 4
200        N3 = 5
201
202        a1 = tensor(rand_dm(N1, density=0.75), identity(N2), identity(N3))
203        a2 = tensor(identity(N1), rand_dm(N2, density=0.75), identity(N3))
204        a3 = tensor(identity(N1), identity(N2), rand_dm(N3, density=0.75))
205        H = a1.dag() * a1 + a2.dag() * a2 + a3.dag() * a3
206
207        c_ops = [np.sqrt(0.01) * a1, np.sqrt(0.025) * a2, np.sqrt(0.05) * a3]
208
209        L1 = liouvillian(H, c_ops)
210        L2 = liouvillian_ref(H, c_ops)
211
212        assert_((L1 - L2).norm('max') < 1e-8)
213
214
215
216class TestSuper_td:
217    """
218    A test class for the QuTiP superoperator functions.
219    """
220
221    def setup_method(self):
222        N = 3
223        self.t1 = QobjEvo([qeye(N)*(1.+0.1j),[create(N)*(1.-0.1j),f]])
224        self.t2 = QobjEvo([destroy(N)*(1.-0.2j)])
225        self.t3 = QobjEvo([[destroy(N)*create(N)*(1.+0.2j),f]])
226        self.q1 = qeye(N)*(1.+0.3j)
227        self.q2 = destroy(N)*(1.-0.3j)
228        self.q3 = destroy(N)*create(N)*(1.+0.4j)
229
230    def test_spre_td(self):
231        "Superoperator: spre, time-dependent"
232        assert_(spre(self.t1)(.5) == spre(self.t1(.5)))
233
234    def test_spost_td(self):
235        "Superoperator: spre, time-dependent"
236        assert_(spost(self.t1)(.5) == spost(self.t1(.5)))
237
238    def test_sprepost_td(self):
239        "Superoperator: sprepost, time-dependent"
240        # left QobjEvo
241        assert_(sprepost(self.t1, self.q2)(.5) ==
242                sprepost(self.t1(.5), self.q2))
243        # left QobjEvo
244        assert_(sprepost(self.q2, self.t1)(.5) ==
245                sprepost(self.q2, self.t1(.5)))
246        # left 2 QobjEvo, one cte
247        assert_(sprepost(self.t1, self.t2)(.5) ==
248                sprepost(self.t1(.5), self.t2(.5)))
249
250    def test_operator_vector_td(self):
251        "Superoperator: operator_to_vector, time-dependent"
252        assert_(operator_to_vector(self.t1)(.5) ==
253                operator_to_vector(self.t1(.5)))
254        vec = operator_to_vector(self.t1)
255        assert_(vector_to_operator(vec)(.5) == vector_to_operator(vec(.5)))
256
257    def test_liouvillian_td(self):
258        "Superoperator: liouvillian, time-dependent"
259        assert_(liouvillian(self.t1)(0.5) == liouvillian(self.t1(0.5)))
260        assert_(liouvillian(None, [self.t2])(0.5) ==
261                liouvillian(None, [self.t2(0.5)]))
262        assert_(liouvillian(self.t1, [self.t2, self.q1, self.t3],
263                            chi=[1,2,3])(0.5) ==
264                liouvillian(self.t1(0.5), [self.t2(0.5), self.q1, self.t3(0.5)],
265                            chi=[1,2,3]))
266
267    def test_lindblad_dissipator_td(self):
268        "Superoperator: lindblad_dissipator, time-dependent"
269        assert_(lindblad_dissipator(self.t2)(.5) ==
270                lindblad_dissipator(self.t2(.5)))
271        assert_(lindblad_dissipator(self.t2, self.q1)(.5) ==
272                lindblad_dissipator(self.t2(.5), self.q1))
273        assert_(lindblad_dissipator(self.q1, self.t2)(.5) ==
274                lindblad_dissipator(self.q1, self.t2(.5)))
275
276
277
278if __name__ == "__main__":
279    run_module_suite()
280