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
34from numpy.linalg import norm
35from numpy.testing import assert_, run_module_suite
36
37from qutip.random_objects import rand_dm, rand_unitary, rand_kraus_map
38from qutip.subsystem_apply import subsystem_apply
39from qutip.superop_reps import kraus_to_super
40from qutip.superoperator import mat2vec, vec2mat
41from qutip.tensor import tensor
42from qutip.qobj import Qobj
43
44
45class TestSubsysApply(object):
46    """
47    A test class for the QuTiP function for applying superoperators to
48    subsystems.
49    The four tests below determine whether efficient numerics, naive numerics
50    and semi-analytic results are identical.
51    """
52
53    def test_SimpleSingleApply(self):
54        """
55        Non-composite system, operator on Hilbert space.
56        """
57        tol = 1e-12
58        rho_3 = rand_dm(3)
59        single_op = rand_unitary(3)
60        analytic_result = single_op * rho_3 * single_op.dag()
61        naive_result = subsystem_apply(rho_3, single_op, [True],
62                                       reference=True)
63        naive_diff = (analytic_result - naive_result).full()
64        naive_diff_norm = norm(naive_diff)
65        assert_(naive_diff_norm < tol,
66                msg="SimpleSingle: naive_diff_norm {} "
67                    "is beyond tolerance {}".format(
68                        naive_diff_norm, tol))
69
70        efficient_result = subsystem_apply(rho_3, single_op, [True])
71        efficient_diff = (efficient_result - analytic_result).full()
72        efficient_diff_norm = norm(efficient_diff)
73        assert_(efficient_diff_norm < tol,
74                msg="SimpleSingle: efficient_diff_norm {} "
75                    "is beyond tolerance {}".format(
76                        efficient_diff_norm, tol))
77
78    def test_SimpleSuperApply(self):
79        """
80        Non-composite system, operator on Liouville space.
81        """
82        tol = 1e-12
83        rho_3 = rand_dm(3)
84        superop = kraus_to_super(rand_kraus_map(3))
85        analytic_result = vec2mat(superop.full() @ mat2vec(rho_3.full()))
86        naive_result = subsystem_apply(rho_3, superop, [True],
87                                       reference=True)
88        naive_diff = (analytic_result - naive_result).full()
89        naive_diff_norm = norm(naive_diff)
90        assert_(naive_diff_norm < tol,
91                msg="SimpleSuper: naive_diff_norm {} "
92                    "is beyond tolerance {}".format(
93                        naive_diff_norm, tol))
94
95        efficient_result = subsystem_apply(rho_3, superop, [True])
96        efficient_diff = (efficient_result - analytic_result).full()
97        efficient_diff_norm = norm(efficient_diff)
98        assert_(efficient_diff_norm < tol,
99                msg="SimpleSuper: efficient_diff_norm {} "
100                    "is beyond tolerance {}".format(
101                        efficient_diff_norm, tol))
102
103    def test_ComplexSingleApply(self):
104        """
105        Composite system, operator on Hilbert space.
106        """
107        tol = 1e-12
108        rho_list = list(map(rand_dm, [2, 3, 2, 3, 2]))
109        rho_input = tensor(rho_list)
110        single_op = rand_unitary(3)
111
112        analytic_result = rho_list
113        analytic_result[1] = single_op * analytic_result[1] * single_op.dag()
114        analytic_result[3] = single_op * analytic_result[3] * single_op.dag()
115        analytic_result = tensor(analytic_result)
116
117        naive_result = subsystem_apply(rho_input, single_op,
118                                       [False, True, False, True, False],
119                                       reference=True)
120        naive_diff = (analytic_result - naive_result).full()
121        naive_diff_norm = norm(naive_diff)
122        assert_(naive_diff_norm < tol,
123                msg="ComplexSingle: naive_diff_norm {} "
124                    "is beyond tolerance {}".format(
125                        naive_diff_norm, tol))
126
127        efficient_result = subsystem_apply(rho_input, single_op,
128                                           [False, True, False, True, False])
129        efficient_diff = (efficient_result - analytic_result).full()
130        efficient_diff_norm = norm(efficient_diff)
131        assert_(efficient_diff_norm < tol,
132                msg="ComplexSingle: efficient_diff_norm {} "
133                    "is beyond tolerance {}".format(
134                        efficient_diff_norm, tol))
135
136    def test_ComplexSuperApply(self):
137        """
138        Superoperator: Efficient numerics and reference return same result,
139        acting on non-composite system
140        """
141        tol = 1e-10
142        rho_list = list(map(rand_dm, [2, 3, 2, 3, 2]))
143        rho_input = tensor(rho_list)
144        superop = kraus_to_super(rand_kraus_map(3))
145
146        analytic_result = rho_list
147        analytic_result[1] = Qobj(
148            vec2mat(superop.full() @ mat2vec(analytic_result[1].full())))
149        analytic_result[3] = Qobj(
150            vec2mat(superop.full() @ mat2vec(analytic_result[3].full())))
151        analytic_result = tensor(analytic_result)
152
153        naive_result = subsystem_apply(rho_input, superop,
154                                       [False, True, False, True, False],
155                                       reference=True)
156        naive_diff = (analytic_result - naive_result).full()
157        naive_diff_norm = norm(naive_diff)
158        assert_(naive_diff_norm < tol,
159                msg="ComplexSuper: naive_diff_norm {} "
160                    "is beyond tolerance {}".format(
161                        naive_diff_norm, tol))
162
163        efficient_result = subsystem_apply(rho_input, superop,
164                                           [False, True, False, True, False])
165        efficient_diff = (efficient_result - analytic_result).full()
166        efficient_diff_norm = norm(efficient_diff)
167        assert_(efficient_diff_norm < tol,
168                msg="ComplexSuper: efficient_diff_norm {} "
169                    "is beyond tolerance {}".format(
170                        efficient_diff_norm, tol))
171
172
173if __name__ == "__main__":
174    run_module_suite()
175