← Back to branch summary

~jrjohansson/qutip/master

« back to all changes in this revision

Viewing changes to qutip/qobj.py

  • Committer: Paul Nation
  • Date: 2011-04-21 04:46:56 UTC
  • Revision ID: git-v1:dd4c966b490aa468dfbd28cef66694df4bf235c8

Show diffs side-by-side

added added

removed removed

Lines of Context:
qutip/qobj.py
 
1
#This file is part of QuTIP.
 
2
#
 
3
#    QuTIP is free software: you can redistribute it and/or modify
 
4
#    it under the terms of the GNU General Public License as published by
 
5
#    the Free Software Foundation, either version 3 of the License, or
 
6
#   (at your option) any later version.
 
7
#
 
8
#    QuTIP is distributed in the hope that it will be useful,
 
9
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
 
10
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
11
#    GNU General Public License for more details.
 
12
#
 
13
#    You should have received a copy of the GNU General Public License
 
14
#    along with QuTIP.  If not, see <http://www.gnu.org/licenses/>.
 
15
#
 
16
# Copyright (C) 2011, Paul D. Nation & Robert J. Johansson
 
17
#
 
18
###########################################################################
 
19
from scipy import *
 
20
import scipy.sparse as sp
 
21
import scipy.linalg as la
 
22
from istests import *
 
23
 
 
24
##
 
25
# @package qobj
 
26
# Quantum Object class. Describes the composition of a quantum system and
 
27
# implements common operations.
 
28
#
 
29
 
 
30
 
 
31
class qobj():
 
32
    """
 
33
    @brief Quantum object class.
 
34
    requires: scipy, scipy.sparse.csr_matrix, scipy.linalg
 
35
    """
 
36
    ################## Define qobj class #################
 
37
    __array_priority__=100 #sets QOBJ priority above numpy arrays
 
38
    def __init__(self,inpt=array([[0]]),dims=[[],[]],shape=[]):
 
39
        """
 
40
        qobj constructor. Optionally takes the dimension array and/or
 
41
        shape array as arguments.
 
42
 
 
43
        @param array    Data for vector/matrix representation of the quantum object
 
44
        @param dims
 
45
        @param shape
 
46
        """
 
47
        if isinstance(inpt,qobj):#if input is already QOBJ then return identical copy
 
48
            self.data=sp.csr_matrix(inpt.data) #make sure matrix is sparse (safety check)
 
49
            if not any(dims):
 
50
                self.dims=inpt.dims
 
51
            else:
 
52
                self.dims=dims
 
53
            if not any(shape):
 
54
                self.shape=inpt.shape
 
55
            else:
 
56
                self.shape=shape
 
57
        else:#if input is NOT QOBJ
 
58
            #if input is int, float, or complex then convert to array
 
59
            if isinstance(inpt,(int,float,complex)):
 
60
                inpt=array([[inpt]])
 
61
            #case where input is array or sparse
 
62
            if (isinstance(inpt,ndarray)) or (isinstance(inpt,sp.csr_matrix)):
 
63
                self.data=sp.csr_matrix(inpt) #data stored as space array
 
64
                if not any(dims):
 
65
                    self.dims=[[inpt.shape[0]],[inpt.shape[1]]] #list of object dimensions
 
66
                else:
 
67
                    self.dims=dims
 
68
                if not any(shape):
 
69
                    self.shape=[inpt.shape[0],inpt.shape[1]] # list of matrix dimensions
 
70
                else:
 
71
                    self.shape=shape
 
72
            elif isinstance(inpt,list):# case where input is not array or sparse, i.e. a list
 
73
                if len(array(inpt).shape)==1:#if list has only one dimension (i.e [5,4])
 
74
                    inpt=array([inpt])
 
75
                else:#if list has two dimensions (i.e [[5,4]])
 
76
                    inpt=array(inpt)
 
77
                    self.data=sp.csr_matrix(inpt)
 
78
                    if not any(dims):
 
79
                        self.dims=[[inpt.shape[0]],[inpt.shape[1]]]
 
80
                    else:
 
81
                        self.dims=dims
 
82
                    if not any(shape):
 
83
                        self.shape=[inpt.shape[0],inpt.shape[1]]
 
84
                    else:
 
85
                        self.shape=shape
 
86
    
 
87
    ##### Definition of PLUS with qobj on LEFT (ex. qobj+4) #########                
 
88
    def __add__(self, other): #defines left addition for qobj class
 
89
        if classcheck(other)=='eseries':
 
90
            return other.__radd__(self)
 
91
        other=qobj(other)
 
92
        if prod(other.shape)==1 and prod(self.shape)!=1: #case for scalar quantum object
 
93
            dat=array(other.full())[0][0]
 
94
            out=qobj()
 
95
            out.data=self.data+dat*sp.csr_matrix(ones(self.shape))
 
96
            out.dims=self.dims
 
97
            out.shape=self.shape
 
98
            return qobj(out)
 
99
        elif prod(self.shape)==1 and prod(other.shape)!=1:#case for scalar quantum object
 
100
            dat=array(self.full())[0][0]
 
101
            out=qobj()
 
102
            out.data=dat*sp.csr_matrix(ones(other.shape))+other.data
 
103
            out.dims=other.dims
 
104
            out.shape=other.shape
 
105
            return qobj(out)
 
106
        elif self.dims!=other.dims:
 
107
            raise TypeError('Incompatible quantum object dimensions')
 
108
        elif self.shape!=other.shape:
 
109
            raise TypeError('Matrix shapes do not match')
 
110
        else:#case for matching quantum objects
 
111
            out=qobj()
 
112
            out.data=self.data+other.data
 
113
            out.dims=self.dims
 
114
            out.shape=self.shape
 
115
            return qobj(out)
 
116
 
 
117
    ##### Definition of PLUS with qobj on RIGHT (ex. 4+qobj) ###############
 
118
    def __radd__(self,other): #defines left addition for qobj class
 
119
        return self+other
 
120
 
 
121
    ##### Definition of SUBTRACTION with qobj on LEFT (ex. qobj-4) #########
 
122
    def __sub__(self,other):
 
123
        return self+(-other)
 
124
 
 
125
    ##### Definition of SUBTRACTION with qobj on RIGHT (ex. 4-qobj) #########
 
126
    def __rsub__(self,other):
 
127
        return (-self)+other
 
128
    
 
129
    ##### Definition of Multiplication with qobj on left (ex. qobj*5) #########
 
130
    def __mul__(self,other):
 
131
        if isinstance(other,qobj): #if both are quantum objects
 
132
            if prod(other.shape)==1 and prod(self.shape)!=1:#case for scalar quantum object
 
133
                dat=array(other.data.todense())[0][0]
 
134
                out=qobj()
 
135
                out.data=dat * (sp.csr_matrix(ones(self.shape))*self.data)
 
136
                out.dims=self.dims
 
137
                out.shape=self.shape
 
138
                return qobj(out)
 
139
            elif prod(self.shape)==1 and prod(other.shape)!=1:#take care of right mul as well for scalar qobjs
 
140
                dat=array(self.data.todense())[0][0]
 
141
                out=qobj()
 
142
                out.data=dat*sp.csr_matrix(ones(other.shape))*other.data
 
143
                out.dims=other.dims
 
144
                out.shape=other.shape
 
145
                return qobj(out)
 
146
            elif self.shape[1]==other.shape[0] and self.dims[1]==other.dims[0]:
 
147
                out=qobj()
 
148
                out.data=self.data*other.data
 
149
                out.dims  = [self.dims[0],  other.dims[1]]
 
150
                out.shape = [self.shape[0], other.shape[1]]
 
151
                return qobj(out)
 
152
            else:
 
153
                raise TypeError("Incompatible qobj shapes")
 
154
 
 
155
        if isinstance(other, list): # if other is a list, do element-wise multiplication
 
156
            prod_list = []
 
157
            for item in other:
 
158
                prod_list.append(self * item)
 
159
            return prod_list
 
160
 
 
161
        if classcheck(other)=='eseries':
 
162
            return other.__rmul__(self)
 
163
 
 
164
        if isinstance(other,(int,float,complex)): #if other is int,float,complex
 
165
            out=qobj()
 
166
            out.data=self.data*other
 
167
            out.dims=self.dims
 
168
            out.shape=self.shape
 
169
            return qobj(out)
 
170
        else:
 
171
            raise TypeError("Incompatible object for multiplication")
 
172
    
 
173
    ##### Definition of Multiplication with qobj on right (ex. 5*qobj) #########
 
174
    def __rmul__(self,other):
 
175
        if isinstance(other,qobj): #if both are quantum objects
 
176
            if prod(other.shape)==1 and prod(self.shape)!=1:#case for scalar quantum object
 
177
                dat=array(other.data.todense())[0][0]
 
178
                out=qobj()
 
179
                out.data=dat * (sp.csr_matrix(ones(self.shape))*self.data)
 
180
                out.dims=self.dims
 
181
                out.shape=self.shape
 
182
                return qobj(out)
 
183
            elif prod(self.shape)==1 and prod(other.shape)!=1:#take care of right mul as well for scalar qobjs
 
184
                dat=array(self.data.todense())[0][0]
 
185
                out=qobj()
 
186
                out.data=dat*sp.csr_matrix(ones(other.shape))*other.data
 
187
                out.dims=other.dims
 
188
                out.shape=other.shape
 
189
                return qobj(out)
 
190
            elif self.shape[1]==other.shape[0] and self.dims[1]==other.dims[0]:
 
191
                out=qobj()
 
192
                out.data=other.data * self.data
 
193
                out.dims=self.dims
 
194
                out.shape=[self.shape[0],other.shape[1]]
 
195
                return qobj(out)
 
196
            else:
 
197
                raise TypeError("Incompatible qobj shapes")
 
198
 
 
199
        if isinstance(other, list): # if other is a list, do element-wise multiplication
 
200
            prod_list = []
 
201
            for item in other:
 
202
                prod_list.append(item * self)
 
203
            return prod_list
 
204
 
 
205
        if classcheck(other)=='eseries':
 
206
            return other.__mul__(self)
 
207
 
 
208
        if isinstance(other,(int,float,complex)): #if other is int,float,complex
 
209
            out=qobj()
 
210
            out.data=other*self.data
 
211
            out.dims=self.dims
 
212
            out.shape=self.shape
 
213
            return qobj(out)
 
214
        else:
 
215
            raise TypeError("Incompatible object for multiplication")
 
216
 
 
217
    ##### Definition of Division with qobj on left (ex. qobj / sqrt(2)) #########
 
218
    def __div__(self,other):
 
219
        if isinstance(other,qobj): #if both are quantum objects
 
220
            raise TypeError("Incompatible qobj shapes [division with qobj not implemented]")
 
221
        #if isinstance(other,qobj): #if both are quantum objects
 
222
        #    raise TypeError("Incompatible qobj shapes [division with qobj not implemented]")
 
223
        if isinstance(other,(int,float,complex)): #if other is int,float,complex
 
224
            out=qobj()
 
225
            out.data=self.data/other
 
226
            out.dims=self.dims
 
227
            out.shape=self.shape
 
228
            return qobj(out)
 
229
        else:
 
230
            raise TypeError("Incompatible object for division")
 
231
    def __neg__(self):
 
232
        out=qobj()
 
233
        out.data=-self.data
 
234
        out.dims=self.dims
 
235
        out.shape=self.shape
 
236
        return qobj(out)
 
237
 
 
238
            
 
239
    def __str__(self):
 
240
        #return "Quantum object: ", dimensions = " + str(self.shape) + "\n" + str(self.data)
 
241
        print "Quantum object: " + "dims = " + str(self.dims) + ", shape = " + str(self.shape)
 
242
        print "QOBJ data = "
 
243
        print self.full()
 
244
        return ""
 
245
    
 
246
    #####---functions acting on quantum objects---######################
 
247
    def dag(self):
 
248
        #returns Adjont (dagger) of operator
 
249
        out=qobj()
 
250
        out.data=self.data.T.conj()
 
251
        out.dims=[self.dims[1],self.dims[0]]
 
252
        out.shape=[self.shape[1],self.shape[0]]
 
253
        return qobj(out)
 
254
    def norm(self):
 
255
        #returns norm of quantum object
 
256
        return la.norm(self.full())
 
257
    def tr(self):
 
258
        #returns trace of quantum object
 
259
        return sum(diag(self.full()))   
 
260
    def full(self):
 
261
        #returns dense matrix form of quantum object data
 
262
        if isinstance(self.data, ndarray):
 
263
            return self.data
 
264
        return array(self.data.todense())
 
265
    def expm(self):
 
266
        if self.dims[0][0]==self.dims[1][0]:
 
267
            return sp_expm(self)
 
268
        else:
 
269
            raise TypeError('Invalid operand for matrix exponential');
 
270
        
 
271
 
 
272
##############################################################################
 
273
#
 
274
#
 
275
# functions acting on QOBJ class
 
276
#
 
277
#
 
278
def dag(inqobj):
 
279
        if not isinstance(inqobj,qobj): #checks for qobj
 
280
                raise TypeError("Input is not a quantum object")
 
281
        outqobj=qobj()
 
282
        outqobj.data=inqobj.data.T.conj()
 
283
        outqobj.dims=[inqobj.dims[1],inqobj.dims[0]]
 
284
        outqobj.shape=[inqobj.shape[1],inqobj.shape[0]]
 
285
        return qobj(outqobj)
 
286
 
 
287
def trans(A):
 
288
    out=qobj()
 
289
    out.data=A.data.T
 
290
    out.dims=[A.dims[1],A.dims[0]]
 
291
    out.shape=[A.shape[1],A.shape[0]]
 
292
    return qobj(out)
 
293
 
 
294
 
 
295
def isherm(ops):
 
296
    if isinstance(ops,qobj):
 
297
        ops=array([ops])
 
298
    ops=array(ops)
 
299
    out=zeros(len(ops))
 
300
    for k in range(len(ops)):
 
301
        if (ops[k].dag()-ops[k]).norm()/ops[k].norm()>1e-12:#any non-zero elements
 
302
            out[k]=0
 
303
        else:
 
304
            out[k]=1
 
305
    return out
 
306
 
 
307
##############################################################################
 
308
#      
 
309
#
 
310
# some functions for increased compatibility with quantum optics toolbox:
 
311
#
 
312
#
 
313
 
 
314
def dims(obj):
 
315
    if isinstance(obj,qobj):
 
316
        return qobj.dims
 
317
    else:
 
318
        raise TypeError("Incompatible object for dims (not a qobj)")
 
319
 
 
320
 
 
321
def shape(inpt):
 
322
    from scipy import shape as shp
 
323
    if isinstance(inpt,qobj):
 
324
        return qobj.shape
 
325
    else:
 
326
        return shp(inpt)
 
327
 
 
328
##############################################################################
 
329
#
 
330
#
 
331
# check for class type (ESERIES,FSERIES)
 
332
#
 
333
#
 
334
#
 
335
def classcheck(inpt):
 
336
    '''
 
337
    Checks for ESERIES and FSERIES class types
 
338
    '''
 
339
    from eseries import eseries
 
340
    if isinstance(inpt,eseries):
 
341
        return 'eseries'
 
342
    else:
 
343
        pass
 
344
 
 
345
 
 
346
########################################################################
 
347
def sp_expm(qo):
 
348
    #############################
 
349
    def pade(m):
 
350
        n=shape(A)[0]
 
351
        c=padecoeff(m)
 
352
        if m!=13:
 
353
            apows= [[] for jj in xrange(int(ceil((m+1)/2)))]
 
354
            apows[0]=sp.eye(n,n).tocsr()
 
355
            apows[1]=A*A
 
356
            for jj in xrange(2,int(ceil((m+1)/2))):
 
357
                apows[jj]=apows[jj-1]*apows[1]
 
358
            U=sp.lil_matrix(zeros((n,n))).tocsr(); V=sp.lil_matrix(zeros((n,n))).tocsr()
 
359
            for jj in xrange(m,0,-2):
 
360
                U=U+c[jj]*apows[jj/2]
 
361
            U=A*U
 
362
            for jj in xrange(m-1,-1,-2):
 
363
                V=V+c[jj]*apows[(jj+1)/2]
 
364
            F=la.solve((-U+V).todense(),(U+V).todense())
 
365
            return sp.lil_matrix(F).tocsr()
 
366
        elif m==13:
 
367
            A2=A*A
 
368
            A4=A2*A2
 
369
            A6=A2*A4
 
370
            U = A*(A6*(c[13]*A6+c[11]*A4+c[9]*A2)+c[7]*A6+c[5]*A4+c[3]*A2+c[1]*sp.eye(n,n).tocsr())
 
371
            V = A6*(c[12]*A6 + c[10]*A4 + c[8]*A2)+ c[6]*A6 + c[4]*A4 + c[2]*A2 + c[0]*sp.eye(n,n).tocsr()
 
372
            F=la.solve((-U+V).todense(),(U+V).todense()) 
 
373
            return sp.lil_matrix(F).tocsr()
 
374
    #################################
 
375
    A=qo.data #extract Qobj data (sparse matrix)
 
376
    m_vals=array([3,5,7,9,13])
 
377
    theta=array([0.01495585217958292,0.2539398330063230,0.9504178996162932,2.097847961257068,5.371920351148152],dtype=float)
 
378
    normA=la.norm(A.todense(),1)
 
379
    if normA<=theta[-1]:
 
380
        for ii in xrange(len(m_vals)):
 
381
            if normA<=theta[ii]:
 
382
                F=pade(m_vals[ii])
 
383
                break
 
384
    else:
 
385
        t,s=frexp(normA/theta[-1])
 
386
        s=s-(t==0.5)
 
387
        A=A/2.0**s
 
388
        F=pade(m_vals[-1])
 
389
        for i in xrange(s):
 
390
            F=F*F
 
391
    return qobj(F,dims=qo.dims,shape=qo.shape)
 
392
 
 
393
 
 
394
def padecoeff(m):
 
395
    if m==3:
 
396
        return array([120, 60, 12, 1])
 
397
    elif m==5:
 
398
        return array([30240, 15120, 3360, 420, 30, 1])
 
399
    elif m==7:
 
400
        return array([17297280, 8648640, 1995840, 277200, 25200, 1512, 56, 1])
 
401
    elif m==9:
 
402
        return array([17643225600, 8821612800, 2075673600, 302702400, 30270240,2162160, 110880, 3960, 90, 1])
 
403
    elif m==13:
 
404
        return array([64764752532480000, 32382376266240000, 7771770303897600,1187353796428800, 129060195264000, 10559470521600,670442572800, 33522128640, 1323241920,40840800, 960960, 16380, 182, 1])
 
405
 
 
406
 
 
407
 
 
408
 
 
409
 
 
410
 
 
411
 
 
412
 
 
413