1# -*- coding: utf-8 -*- 2""" 3Create pauli spin matrices and other quantum mechanical operators 4for collections of single-spin particles. 5Manipulate quantum states and simulate quantum computer circuits. 6""" 7from __future__ import print_function 8import sys 9import numpy as np 10import copy 11import itertools 12 13__version__ = '2.3.2' 14if sys.version_info.major == 3: 15 unicode = str 16 17i = 1j 18debug = False 19sign = lambda a: '+' if a>0 else '-' if a<0 else '+' 20check = 1 21 22base_reprs = { 23 'ud':{'up':'u','down':'d'}, 24 '01':{'up':'0','down':'1'}, 25 'arrow':{'up':u'\u2191','down':u'\u2193'} 26} 27 28base_repr = base_reprs['ud'] 29 30def set_base_repr(arg): 31 '''Set the printout of the base representation: 32 33 :param str arg: 'ud', '01', or 'arrow' 34 35 These will cause state vector bases to be printed out as, 36 for example 37 38 - `arg = 'ud'` results in :math:`\\Ket{u}`, :math:`\\Ket{ud}`, etc. 39 - `arg = '01'` results in :math:`\\Ket{0}`, :math:`\\Ket{01}`, etc. 40 - `arg = 'arrow'` results in :math:`\\Ket{\\uparrow}`, :math:`\\Ket{\\uparrow\\downarrow}`, etc. 41 42 ''' 43 global base_repr 44 base_repr = base_reprs[arg] 45 state.bases = {1:['|%s>'%base_repr[x] for x in ['up','down']]} 46 state.base_reprs = arg 47 48class state(object): 49 '''Quantum state 50 51 A quantum state has the property of being either a bra or a ket vector. 52 States `a` and `b` can 53 54 - Inner product: `a.H*b` = :math:`\\Braket{a | b}` 55 - Kronecker product: (expand the number of particles in the state) `a**b` = :math:`\\Ket{a} \\otimes \\Ket{b} = \\Ket{ab}` 56 - Hermetion transpose: (convert bra to ket and vice-versa) `a.H` does either 57 :math:`\\Ket{a} \\rightarrow \\Bra{a}` or :math:`\\Bra{a} \\rightarrow \\Ket{a}` 58 - Normalize: `a.N` = :math:`\\frac{1}{\\sqrt{\\Braket{a|a}}} \\Ket{a}` 59 - Access the wave function as a column vector: `a.phi` = :math:`[\\Braket{a|e_i}]` where :math:`\\Ket{e_i}` are the basis vectors 60 61 ''' 62 63 bases = {1:['|%s>'%base_repr[x] for x in ['up','down']]} 64 base_reprs = 'ud' 65 66 def __init__(self,s, bra_flag = False): 67 '''initialize with either a basis ket, or a wave function 68 ''' 69 assert isinstance(s,(str,unicode,np.matrix,list)) 70 if isinstance(s,list): 71 s = np.matrix(s) 72 if isinstance(s,np.matrix): 73 if s.dtype not in (float,complex): 74 s = s.astype(float) 75 self.n = n = s.size 76 self.phi = np.matrix( np.ravel(s).reshape((n,1)) ) 77 self.n_particles = n_particles = int(np.log2(n)) 78 elif isinstance(s,(str,unicode)): 79 self.n_particles = n_particles = len(s.lstrip('|').rstrip('>')) 80 81 if n_particles not in self.bases: 82 self.basis = create_basis(n_particles) 83 self.bases[n_particles]=self.basis 84 else: 85 self.basis = self.bases[n_particles] 86 87 if isinstance(s,(str,unicode)): 88 if s.startswith('<') or s.endswith('|'): 89 raise Exception('%s not a ket'%s) 90 if not s.startswith('|'): 91 s = '|'+s 92 if not s.endswith('>'): 93 s = s+'>' 94 assert s in self.basis, '%s not in basis'%s 95 k = self.basis.index(s) 96 self.name = s 97 self.n = len(self.basis) 98 self.phi = np.matrix([0.]*self.n).T # the wave function 99 self.phi[k] = 1. 100 101 self.bra = False 102 self.kind = 'ket' 103 104 if bra_flag: 105 self._bra() 106 107 def _bra(self): 108 '''convert ket to bra 109 ''' 110 if self.bra: 111 return 112 self.phi = self.phi.H 113 if hasattr(self,'name') and self.name in self.basis: 114 self.name = '<'+self.name[1:-1]+'|' 115 self.basis = ['<'+s[1:-1]+'|' for s in self.basis] 116 self.bra = True 117 self.kind = 'bra' 118 return self 119 120 def _ket(self): 121 '''convert bra to ket 122 ''' 123 if not self.bra: 124 return 125 self.phi = self.phi.H 126 if hasattr(self,'name') and self.name in self.basis: 127 self.name = '|'+self.name[1:-1]+'>' 128 self.basis = ['|'+s[1:-1]+'>' for s in self.basis] 129 self.bra = False 130 self.kind = 'ket' 131 return self 132 133 def __unicode__(self): 134 phi = self.phi 135 first_one = True 136 s = '0' 137 print_options = np.get_printoptions() 138 prec = print_options['precision'] 139 if print_options['suppress']: 140 coef_format = '%%0.%df '%prec 141 else: 142 coef_format = '%%0.%dg '%prec 143 real_coef_format = coef_format 144 imag_coef_format = '%si'%coef_format 145 complex_coef_format = '( %s%%s %s )'%(real_coef_format,imag_coef_format) 146 sign_real_coef_format = '%%s %s'%real_coef_format 147 sign_imag_coef_format = '%si'%sign_real_coef_format 148 sign_complex_coef_format = '+ %s'%complex_coef_format 149 150 for k in range(self.n): 151 phi_k = phi.item(k) 152 if phi_k != 0: 153 154 if first_one: 155 s = '' 156 if isinstance(phi_k,(float,int)): 157 if phi_k == 1: 158 coef = '' 159 elif phi_k == -1: 160 coef = '- ' 161 else: 162 coef = real_coef_format%phi_k #'%g '%phi_k 163 else: # complex 164 if phi_k.imag == 0: 165 if phi_k.real == 1: 166 coef = '' 167 elif phi_k.real == -1: 168 coef = '- ' 169 else: 170 coef = real_coef_format%phi_k.real #'%g '%phi_k.real 171 elif phi_k.real == 0: 172 if phi_k.imag == 1: 173 coef = ' i ' 174 elif phi_k.imag == -1: 175 coef = ' - i ' 176 else: 177 coef = imag_coef_format%phi_k.imag #'%g i '%phi_k.imag 178 else: 179 #coef = '( %g %s %g i )'%(phi_k.real,sign(phi_k.imag),abs(phi_k.imag)) 180 coef = complex_coef_format%(phi_k.real,sign(phi_k.imag),abs(phi_k.imag)) 181 first_one = False 182 183 else: 184 if isinstance(phi_k,(float,int)): 185 if phi_k == 1: 186 coef = '+ ' 187 elif phi_k == -1: 188 coef = '- ' 189 else: 190 #coef = '%s %g '%(sign(phi_k),abs(phi_k)) 191 coef = sign_real_coef_format%(sign(phi_k),abs(phi_k)) 192 else: # complex 193 if phi_k.imag == 0: 194 if phi_k.real == 1: 195 coef = ' + ' 196 elif phi_k.real == -1: 197 coef = ' - ' 198 else: 199 #coef = '%s %g '%(sign(phi_k.real),abs(phi_k.real)) 200 coef = sign_real_coef_format%(sign(phi_k.real),abs(phi_k.real)) 201 elif phi_k.real == 0: 202 if phi_k.imag == 1: 203 coef = ' + i ' 204 elif phi_k.imag == -1: 205 coef = ' - i ' 206 else: 207 #coef = '%s %g i '%(sign(phi_k.imag),abs(phi_k.imag)) 208 coef = sign_imag_coef_format%(sign(phi_k.imag),abs(phi_k.imag)) 209 else: 210 #coef = '+ ( %g %s %g i )'%(phi_k.real,sign(phi_k.imag),abs(phi_k.imag)) 211 coef = sign_complex_coef_format%(phi_k.real,sign(phi_k.imag),abs(phi_k.imag)) 212 213 s += '%s%s '%(coef, self.basis[k]) 214 return s 215 216 def __repr__(self): 217 if sys.version_info.major == 2: 218 return self.__unicode__().encode('utf-8') 219 elif sys.version_info.major == 3: 220 return str(self.__unicode__()) 221 222 def __add__(self,another): 223 """|a> + |b> 224 """ 225 assert isinstance(another,self.__class__) 226 assert self.bra == another.bra 227 return self.__class__(self.phi + another.phi, bra_flag = self.bra) 228 229 def __sub__(self,another): 230 """|a> - |b> 231 """ 232 assert isinstance(another,self.__class__) 233 assert self.bra == another.bra 234 return self.__class__(self.phi - another.phi) 235 236 def __mul__(self,another): 237 """ 238 inner product: <a|b> 239 operation on bra: <a|A 240 scalar product: <a| alpha or <a| alpha 241 outer product: |a><b| 242 """ 243 if debug: print ('state.__mul__') 244 assert isinstance(another,(int,float,complex,self.__class__,operator,str)) 245 if isinstance(another,self.__class__): # either inner or outer product 246 if self.bra: # <a|b> 247 assert not another.bra,'cannot multiply two bras' 248 r = (self.phi*another.phi).item(0) 249 if isinstance(r,complex): 250 if r.imag == 0: 251 r = r.real 252 return r 253 else: # |a><b| 254 assert another.bra,'cannot multiply two kets' 255 return operator(self.phi*another.phi) 256 elif isinstance(another,operator): # <a|A 257 assert self.bra 258 phi = (self.phi * another.op).T 259 return bra(phi) 260 elif isinstance(another,str): # '<a|'*A 261 another = self.__class__(another) 262 return self*another 263 elif isinstance(another,(int,float,complex)): 264 if self.bra: 265 r = bra(self.phi*another) # <a| alpha 266 else: 267 r = ket(self.phi*another) # |a> alpha 268 return r 269 270 def __rmul__(self,num): 271 """scalar product alpha |a> or alpha <a| 272 """ 273 if debug: print ('state_base.__rmul__') 274 assert isinstance(num,(int,float,complex)) 275 if self.bra: 276 r = bra(self.phi*num) # <a| alpha 277 else: 278 r = ket(self.phi*num) # |a> alpha 279 return r 280 281 def __div__(self,other): 282 """scalar product |a> / alpha 283 If both arguments are states, this is a test of "divides", returning (True/False,byfactor) 284 """ 285 if isinstance(other,state): # check if one state is a multiple of another 286 return self.divides(other) 287 assert isinstance(other,(int,float,complex)) 288 r = self.copy() 289 r.phi = self.phi/other 290 return r 291 292 def divides(self,other): 293 assert isinstance(other,state) 294 a = self.phi.flat 295 b = other.phi.flat 296 q = None 297 l = [] 298 for aa,bb in zip(a,b): 299 if aa == 0 and bb == 0: 300 l.append(True) 301 elif bb == 0: 302 l.append(False) 303 else: 304 f = aa/bb 305 if q is None: 306 q = f 307 l.append(q == f) 308 r = np.array(l).all() 309 if r: 310 return(r,q) 311 else: 312 return(r,None) 313 314 def __truediv__(self,num): 315 """scalar product |a> / alpha 316 """ 317 assert isinstance(num,(int,float,complex)) 318 if self.bra: 319 r = bra(self.phi/num) # <a| alpha 320 else: 321 r = ket(self.phi/num) # |a> alpha 322 return r 323 324 def __neg__(self): 325 return (-1)*self # - |a> 326 327 @property 328 def H(self): 329 '''Hermetian conjugate 330 ''' 331 if self.bra: 332 return self.copy()._ket() 333 return self.copy()._bra() 334 335 @property 336 def N(self): 337 '''Normalized 338 ''' 339 return self.normalized() 340 341 # def __getattr__(self,dot_op): 342 # '''Unitary operators: 343 # x.H - Hermetian Transpose 344 # x.N - Normalize 345 # ''' 346 # if dot_op == 'H': 347 # if self.bra: 348 # return self.copy()._ket() 349 # return self.copy()._bra() 350 # if dot_op == 'N': 351 # return self.normalized() 352 # else: 353 # raise AttributeError(self.__repr__(),' object has no attribute ',dot_op) 354 355 def __pow__(self,other): 356 '''kroneker product 357 ''' 358 return self.kron(other) 359 360 def __eq__(self,other): 361 assert isinstance(other,self.__class__) 362 assert self.bra == other.bra 363 return (self.phi == other.phi).all() 364 365 def copy(self): 366 return copy.deepcopy(self) 367 368 def normalize(self): 369 """normalizes the wave function 370 """ 371 if self.bra: 372 norm = (self.phi*self.phi.H).item(0) 373 else: 374 norm = (self.phi.H*self.phi).item(0) 375 self.phi /= np.sqrt(norm) 376 377 def normalized(self): 378 """returns the normalized version 379 """ 380 if self.bra: 381 norm = (self.phi*self.phi.H).item(0) 382 else: 383 norm = (self.phi.H*self.phi).item(0) 384 return self/np.sqrt(norm) 385 386 def density(self): 387 """generate the density matrix 388 :math:`\\rho = \\ket{s}\\bra{s}` 389 """ 390 return (self*self.H).op 391 392 def _index(self): 393 '''makes an index from string representations of basis states to 394 interger offset within the state vector 395 ''' 396 s1 = ['u','d'] 397 s = copy.copy(s1) 398 for k in range(self.n_particles-1): 399 s = [x+y for x,y in list(itertools.product(s,s1))] 400 401 d = dict(zip(s,range(self.n))) 402 return d 403 404 def density1(self,particle_no=None): 405 '''generate the density matrix for one 406 of the particles, `particle_no` 407 Particle numbers are 0...n_particles-1 408 ''' 409 if particle_no == None: # old style call 410 return (self.density1(0),self.density1(1)) 411 assert particle_no < self.n_particles 412 d = self._index() 413 p = particle_no 414 pstr = '_' 415 rho = self.density() 416 s1 = ['u','d'] 417 keys = sorted(d.keys())[::-1] # 'uuu','uud',... 418 f = set([key[:p]+pstr+key[p+1:] for key in keys]) # 'u_u','u_d'... 419 ind_keys_u,ind_keys_d = [[ff.replace(pstr,s) for ff in f] for s in s1] # [['uuu','uud'...],['udu','udd'...]] 420 rho_p = np.matrix(np.zeros((2,2))) 421 for z in zip(ind_keys_u,ind_keys_d): 422 ind2x2 = [[(z0,z1) for z0 in z] for z1 in z] # [[('dud', 'dud'), ('ddd', 'dud')], [('dud', 'ddd'), ('ddd', 'ddd')]] 423 for i in range(2): 424 for j in range(2): 425 ii,jj = ind2x2[i][j] 426 rho_p[i,j] += rho[ d[ii],d[jj] ] 427 return rho_p 428 429 def correlation(self): 430 '''compute the :math:`n_{particles}\\times n_{particles}` correlation matrix 431 ''' 432 n_particles = self.n_particles 433 u = ket([1,0]) 434 d = ket([0,1]) 435 Sz = u*u.H - d*d.H 436 S0 = u*u.H + d*d.H # operator(np.identity(2)) 437 c = np.zeros((n_particles, n_particles)) 438 A = [] 439 for p1 in range(n_particles): 440 Ap = (Sz if p1==0 else S0) 441 for p2 in range(1,n_particles): 442 Ap = Ap**(Sz if p2==p1 else S0) 443 A.append(Ap) 444 for i in range(n_particles): 445 for j in range(n_particles): 446 c[i,j] = (A[i]*self).H*(A[j]*self) 447 return c 448 449 def prob(self,A,s): 450 """determine the probability of this state being in 451 state s after the measurement A 452 453 :param operator A: 454 :param state s: 455 456 :return: `a.prob(A,s) = a.H*A*s =` :math:`\\Braket{a | A | s}` 457 458 """ 459 return np.abs(s.H*A*self)**2 460 461 def kron(self,another): 462 assert isinstance(another,self.__class__) 463 assert self.bra == another.bra 464 r = state([0.]*self.n*another.n,self.bra) 465 r.phi = np.kron(self.phi,another.phi) 466 return r 467 468 def entangled(self): 469 '''determine if a 2-particle pure state is entangled. 470 This works only for n=2 particle states. 471 472 :return: True or False 473 474 ''' 475 assert self.n_particles == 2 476 S = entropy(ptrace(self.density(),[0])) 477 if np.isclose(S,0.): 478 return False 479 else: 480 return True 481 482class operator(object): 483 ''' 484 Quantum operator 485 486 Operators can 487 488 - Multiply :class:`state`, `A*s` = :math:`A \\Ket{s}` 489 - Multiply each other, `A*B` = :math:`A\\times B` 490 - Kronecker-product `A**B` = :math:`A\\otimes B`. 491 492 Also 493 494 - `A.H` is the Hermetian transpose: :math:`A^H` 495 - `A.N` is normalization, so that :math:`\\lvert det(\\bar A) \\rvert = 1` 496 497 :param numpy.ndmatrix A: Initialize with a square, Hermetian matrix `A` 498 ''' 499 ops = ['s0','sx','sy','sz'] 500 '''a list of operators that can be specified by a string argumant to :class:`operator`''' 501 502 sig_0 = np.matrix([[1,0],[0,1]]) 503 sig_x = np.matrix([[0,1],[1,0]]) 504 sig_y = np.matrix([[0,-i],[i,0]]) 505 sig_z = np.matrix([[1,0],[0,-1]]) 506 op_mats = [sig_0,sig_x,sig_y,sig_z] 507 508 @property 509 def matrix(self): 510 '''returns the matrix representation of the operator 511 ''' 512 return self.op 513 514 def __init__(self, A): 515 '''initialize given a square, Hermetian matrix A 516 ''' 517 if debug: print ('operator.__init__') 518 if not isinstance(A,(str,np.matrix)): 519 A = np.matrix(A) 520 assert isinstance(A,(str,np.matrix)) 521 if isinstance(A,str): 522 assert A in self.ops 523 self.name = A 524 k = self.ops.index(A) 525 self.op = A = self.__class__.op_mats[k] 526 n,m = A.shape 527 assert n == m 528 self.op = A 529 n_particles = int(np.log2(n)) 530 self.basis = create_basis(n_particles) 531 w,v = np.linalg.eig(A) 532 self.eigenstates = [v[:,k] for k in range(n)] 533 self.observables = list(w) 534 535 def __repr__(self): 536 return str(self.op) 537 538 def __mul__(self,another): 539 """ 540 multiplication on the left: operator*another 541 scalar multiplication: A*alpha 542 operation on ket: A |a> 543 operator x operator: A*B 544 operator x matrix: A*sig 545 operator x string (basis state ket): A*'|u>' 546 """ 547 if debug: print ('%r.__mul__(%r)'%(self.__class__,type(another))) 548 assert isinstance(another,(str,state,self.__class__,np.matrix,int,float,complex)) 549 if isinstance(another,(int,float,complex)): # alpha A 550 return self.__class__(self.op*another) 551 elif isinstance(another,state): # A |a> 552 assert not another.bra 553 return state(self.op*another.phi) 554 elif isinstance(another,self.__class__): # A*B 555 return self.__class__(self.op*another.op) 556 elif isinstance(another,np.matrix): # assume this operator*density_matrix 557 return self.__class__(self.op*another) 558 elif isinstance(another,str): # convert the string to a basis state 559 return self*ket(another) 560 561 def __rmul__(self,another): 562 """ 563 multiplication on the right: another*operator 564 scalar multiply : A*alpha 565 by bra: <a| A 566 by matrix: sig*A 567 by string (basis state bra): '<a|'*A 568 """ 569 if debug: print ('%r.__rmul__(%r)'%(self.__class__,type(another))) 570 assert isinstance(another,(str,state,np.matrix,int,float,complex)) 571 if isinstance(another,(int,float,complex)): 572 return self.__class__(self.op*another) # A*alpha 573 elif isinstance(another,state): 574 return bra(another.phi*self.op) # <a| A 575 elif isinstance(another,np.matrix): 576 return self.__class__(another*self.op) # sig*A 577 elif isinstance(another,str): 578 return bra(another)*self # convert the string to a basis state 579 580 def __add__(self,another): 581 """ A+B 582 """ 583 assert isinstance(another,self.__class__) 584 return self.__class__(self.op+another.op) # A + B 585 586 def __sub__(self,another): 587 """ A-B 588 """ 589 assert isinstance(another,self.__class__) 590 return self.__class__(self.op-another.op) # A - B 591 592 def __div__(self,num): 593 """ 594 divide by a scalar A/alpha 595 """ 596 assert isinstance(num,(int,float,complex)) 597 return self.__class__(self.op/num) 598 599 def __truediv__(self,num): 600 """ 601 divide by a scalar A/alpha 602 """ 603 assert isinstance(num,(int,float,complex)) 604 return self.__class__(self.op/num) 605 606 def __neg__(self): # -A 607 """ 608 scalar negation -A 609 """ 610 return (-1)*self 611 612 def __eq__(self,other): 613 assert isinstance(other,self.__class__) 614 return (self.op == other.op).all() 615 616 def __pow__(self,other): 617 return self.kron(other) 618 619 @property 620 def H(self): 621 '''Hermetian conjugate''' 622 r = self.copy() 623 r.op = self.op.conj().T 624 return r 625 626 @property 627 def N(self): 628 '''Normalized''' 629 return self.normalized() 630 631 # def __getattr__(self,dot_op): 632 # '''Unitary operator: 633 # A.N - Normalize to make unitary 634 # ''' 635 # if dot_op == 'H': 636 # r = self.copy() 637 # r.op = self.op.conj().T 638 # return r 639 # if dot_op == 'N': 640 # return self.normalized() 641 # else: 642 # raise AttributeError(self.__repr__(),' object has no attribute ',dot_op) 643 644 def copy(self): 645 return copy.deepcopy(self) 646 647 def kron(self,other): 648 assert isinstance(other,self.__class__) 649 return self.__class__(np.kron(self.op,other.op)) 650 651 def normalized(self): 652 '''return the normalized version of the operator 653 such that :math:`\\lvert det(\\bar A) \\rvert = 1`. 654 ''' 655 norm = np.sqrt(np.abs(np.linalg.det(self.matrix))) 656 if np.isclose(norm,0.): 657 raise Error('operator is close to singular, cannot normalize') 658 return self/norm 659 660 def measure(self,s): 661 '''produce the measurement using operator as an observable. 662 This produces a single number equal to one of the eigenvalues :math:`\\lambda_i` 663 of `A` according to a probability distribution 664 :math:`P({\\lambda_i}) = \\Braket{ \\Phi_i | s }^2` and 665 collapses the wave function to the eigenstate :math:`\\Ket{\\Phi_i}`. 666 667 Returns the tuple (observed value, collaped state) = :math:`(\\lambda_i, \\Ket{ \\Phi_i })`. 668 669 :param state s: 670 ''' 671 s = ket(s) 672 n = s.n 673 assert n == self.op.shape[0] 674 E = np.block(self.eigenstates) 675 c = np.array(E.H*s.phi) 676 c = [c.item(k) for k in range(n)] 677 P = [np.abs(x)**2 for x in c] 678 k = np.random.choice( range(n), p=P) 679 lam = self.observables[k] 680 #print('lam = %r'%lam) 681 k_set = np.where(np.isclose(self.observables,lam*np.ones(n)))[0] # repeated eigenvalues 682 #print('k_set = %r'%k_set) 683 ss = ket([0]*n) 684 for k in k_set: 685 #print('c[k] = %r'%c[k]) 686 ss += c[k]*ket(self.eigenstates[k]) 687 ss.normalize() 688 return(lam,ss) 689 690 def average(self,s): 691 '''calculate the average value of the observable for a given pure state, :math:`\\Ket{s}` 692 or, for (in general) a mixed state, given the density matrix :math:`\\rho` 693 694 .. math:: 695 696 \\Braket{A} = \\Braket{s | A | s}\ or \ \\Braket{A} = \\rm{tr}\\left(\\rho A\\right) 697 698 :param s: either the :class:`state` or the :class:`ndarray` density matrix 699 700 ''' 701 if isinstance(s,state): 702 return bra(s)*self*ket(s) 703 elif isinstance(s, np.ndarray): 704 rho = np.matrix(s) 705 return (rho*self.matrix).trace().item(0) 706 707 def eig(self): 708 '''calcluate the eigenvalues and eigenvectors of the operator. 709 Returns (list of eigenvalues, list of eigenvectors) = 710 :math:`([\\lambda_1,\\lambda_2,\\dots],[\\Ket{\\Phi_1},\\Ket{\\Phi_2},\\dots])` 711 712 Depricated 5/2/18 - now precomputed for all operators 713 and available as properties `observables` and `eigenstates` 714 ''' 715 ev,evec = np.linalg.eig(self.matrix) 716 evec = [ket(x) for x in evec.T] 717 return ev,evec 718 719def ket(s): 720 """generate a state vector given the state as a string or a wave function 721 or generate the ket version of a bra vector 722 723 :param s: specifier. If `s` is a: 724 725 - :class:`str`: generate one of the basis states (e.g. `u`, `d`, `uu`, `ud`, `uud`, etc.) 726 - :class:`numpy.ndarray`: generate state from the wave function vector 727 - :class:`state`: convert the state (bra or ket) to a ket 728 729 :return: :math:`\\ket{s}` 730 :rtype: :class:`state` 731 732 """ 733 assert isinstance(s,(str,unicode,np.matrix,list,state)) 734 if isinstance(s,(str,unicode)): 735 if s.startswith('<') and s.endswith('|'): 736 s = '|'+s[1:-1]+'>' # turn it into a ket string 737 return state(s) 738 elif isinstance(s,(np.matrix,list)): 739 return state(s) 740 elif isinstance(s,state): 741 if not s.bra: 742 return s 743 return s.copy()._ket() 744 raise Exception('%r is not a valid state'%s) 745 746def bra(s): 747 """generate a state vector given the state as a string or wave function 748 or generate the bra version of a ket vector 749 750 :param s: specifier. If `s` is a: 751 752 - :class:`str`: generate one of the basis states (e.g. `u`, `d`, `uu`, `ud`, `uud`, etc.) as a bra 753 - :class:`numpy.ndarray`: generate state from the wave function vector 754 - :class:`state`: convert the state (bra or ket) to a bra 755 756 :return: :math:`\\bra{s}` 757 :rtype: :class:`state` 758 759 """ 760 assert isinstance(s,(str,np.matrix,list,state)) 761 if isinstance(s,str): 762 s = ket(s) 763 return bra(s) 764 elif isinstance(s,(np.matrix,list)): 765 return bra( ket(s) ) 766 elif isinstance(s,state): 767 if s.bra: 768 return s 769 return s.copy()._bra() 770 else: 771 raise Exception('% is not a valid state'%s) 772 773def create_basis(n): 774 '''create a basis set for n particles 775 ''' 776 #s1 = ['u','d'] 777 s1 = [base_repr['up'],base_repr['down']] 778 s = copy.copy(s1) 779 for k in range(n-1): 780 s = [x+y for x,y in list(itertools.product(s,s1))] 781 s = ['|'+x+'>' for x in s] 782 return s 783 784def density(p,s): 785 '''Create the mixed state density matrix, :math:`\\rho` given 786 a set of probabilities and a list of pure states. 787 788 .. math:: 789 790 \\rho = \\sum_i p_i \\Ket{s_i} \\bra{s_i} 791 792 :param list p: the list of probabilities :math:`p_i` 793 :param list s: the list of pure states :math:`\\bra{s_i}` 794 795 The probabilities must all be between zero and one and sum to one. 796 797 :rtype: :class:`numpy.matrix` 798 799 ''' 800 p = np.array(p) 801 assert ((p>=0) &(p<=1)).all() 802 assert np.isclose(p.sum(),1.0) 803 assert len(p) == len(s) 804 rho = 0 805 for p_i,psi_i in zip(p,s): 806 rho += p_i*psi_i.density() 807 return rho 808 809def entropy(rho,frac=False,decohere=False): 810 """Calculate the Von Neumann entropy given the density matrix. 811 812 .. math:: 813 814 S = -\\rm{tr}\\left(\\rho\, \\log(\\rho)\\right) 815 816 where :math:`\\rho` is the density matrix. 817 818 :param numpy.matrix rho: the density matrix 819 :param bool frac: If `frac=True`, return :math:`S/S_{max}` where :math:`S_{max} = \\log(n)` the maximum entropy. 820 :param bool decohere: If decohere=True then assume decoherent population (off-diagonal elements of :math:`\\rho` set equal to zero) 821 822 :return: entropy, :math:`S` 823 :rtype: float 824 825 """ 826 if decohere: 827 w = np.diag(rho) 828 else: 829 w,v = np.linalg.eig(rho) 830 n = w.size 831 S = 0 832 for k in range(n): 833 if w[k] > 0 and not np.isclose(w[k],1.0): 834 S -= w[k]*np.log(w[k]) 835 if frac: 836 S_max = np.log(float(n)) 837 return S/S_max 838 else: 839 return S 840 841def set_printoptions(**karg): 842 '''feed print options to numpy. 843 help(numpy.set_printoptions) for more info 844 Typical use case: 845 846 .. code-block python 847 >>> set_printoptions(precisionn=3,suppress=True) 848 849 to set the print precision to 3 digits after the decimal point and 850 suppress scientific notation 851 ''' 852 np.set_printoptions(**karg) 853 854def _mosh(s1,s2,ps): 855 '''mosh string 1 with string 2 putting string 2 characters in posisions ps 856 ''' 857 n = len(s1) + len(s2) 858 psk = [ x for x in range(n) if x not in ps ] 859 s = np.array(['']*n) 860 s[ps] = list(s2) 861 s[psk] = list(s1) 862 return ''.join(list(s)) 863 864def ptrace(rho,ps): 865 '''**Partial Trace**: generate the density matrix with specified particles traced out. 866 867 :param np.matrix rho: the multi-particle density matrix 868 :param list ps: list of particles to trace out; given by number 0,1...n_particles-1 869 870 If the list is empty, return the full density matrix. 871 If the list has all the particles, then returns the trace of the full density matrix (always = 1) 872 The list must have unique particle numbers which are all < n_particles 873 874 :rtype: :class:`numpy.matrix` 875 876 See: https://en.wikipedia.org/wiki/Partial_trace 877 ''' 878 if len(ps) == 0: return s 879 n_particles = int(np.log2(rho.shape[0])) 880 n = n_particles 881 assert np.array([x < n for x in ps]).all(),'all elements of ps %r must be < number of particles %r'%(ps,[x < n for x in ps]) 882 assert np.unique(ps).size == len(ps),'elements of ps must be unique: %r'%ps 883 ps = sorted(ps) 884 nr = 2**len(ps) 885 nk = 2**(n - len(ps)) 886 n = 2**n 887 if nk == 1: return 1. 888 d,dr,dk = [state([1]*x)._index() for x in [n,nr,nk]] 889 di,dri,dki = [sorted(x.keys(),reverse=True) for x in [d,dr,dk]] 890 Ms = [ [[ (_mosh(x,z,ps),_mosh(y,z,ps)) for x in dki] for y in dki] for z in dri] 891 Msv = [[[rho[ d[x], d[y] ] for (x,y) in row ] for row in mat ] for mat in Ms ] 892 den = np.sum( [ np.array(x) for x in Msv], axis=0 ) 893 return np.matrix(den) 894 895def Cab(rho): 896 '''calculate the 'concurrence' of a two particle system. 897 ref: Coffman, V., Kundu, J., & Wootters, W. K. (2000). 898 Distributed Entanglement. Physical Review A, 61(5), 5–9. 899 http://doi.org/10.1103/PhysRevA.61.052306 900 901 :param np.matrix rho: 902 903 Presently this only calculates on a two particle system (4x4 density matrix) 904 ''' 905 assert rho.shape == (4,4) 906 sy = operator('sy') 907 rho_tilde = ((sy**sy)*(rho.conj())*(sy**sy)).op 908 lam = np.linalg.eig(rho*rho_tilde)[0].real 909 lam = np.clip(lam,0,np.inf) 910 lam = np.sqrt(np.sort(lam))[::-1] 911 cab = np.max([lam[0]-lam[1]-lam[2]-lam[3],0]) 912 return cab 913 914#============ Quantum Gates ============= 915from scipy.linalg import sqrtm,block_diag 916from numpy import diag,matrix 917 918class gate(operator): 919 '''Quantum gate, a subclass of :class:`operator` 920 921 Generate a gate with a constructor call: 922 923 g = **gate** (*name*) 924 925 where name is one of the names in the list `gate.gates`. These include 'Hadamard','NOT', etc. 926 One can also force any operator to be a gate if the name is one of the pre-defined operators 927 (given in the `operator.ops` list), or with: 928 929 g = **gate** (*matrix* | *operator*, [name = *name*]) 930 931 Some gates, `'Rz'` and `'XX'` need an argument, the angle :math:`\\phi`. These are 932 created by passing `phi` as the second argument: 933 934 g = **gate** (*name*, phi = *phi*) 935 936 ''' 937 gates = ['H','Hadamard','X','Y','Z','Rz','NOT','sNOT','cNOT','S','SWAP','sSWAP','cSWAP','ccNOT','XX'] 938 '''A list of the predefined gates that can be specified by string argument to :class:`gate`''' 939 940 url = 'https://en.wikipedia.org/wiki/Quantum_logic_gate' 941 '''wikipedia reference''' 942 943 # def __new__(cls,arg,**kwargs): 944 # if debug: print('gate.__new__') 945 # if isinstance(arg,operator): 946 # self = super(gate,cls).__new__(cls) 947 # self.__init__(arg) 948 # return self 949 # if not isinstance(arg,str): # i.e. matrix or list 950 # self = super(gate,cls).__new__(cls) 951 # self.__init__(arg) 952 # return self 953 # else: 954 # self = super(gate,cls).__new__(cls) 955 # self.__init__(arg) 956 # return(self) 957 958 def __init__(self,arg1, **kwargs): 959 if debug: print('gate.__init__') 960 if isinstance(arg1,operator): 961 self.op = arg1.op 962 name = '' 963 if hasattr(arg1,'name'): 964 name = arg1.name 965 self.name = '%s gate'%name 966 self.meta = kwargs 967 return 968 if not isinstance(arg1,str): # matrix or list 969 super(self.__class__,self).__init__(arg1) 970 name = '' 971 if hasattr(self,'name'): 972 name = self.name 973 self.name = '%s gate'%name 974 self.meta = kwargs 975 return 976 977 name = arg1 978 if name in operator.ops: 979 super(self.__class__,self).__init__(name) 980 self.name = '%s gate'%name 981 self.meta = kwargs 982 return 983 984 if name in ['H','Hadamard']: 985 self.op = (sx+sz).N.matrix 986 self.name = 'Hadamard (H) gate' 987 elif name in ['X','NOT']: 988 self.op = sx.matrix 989 self.name = 'NOT (X) gate' 990 elif name == 'Y': 991 self.op = sy.matrix 992 self.name = 'Y gate' 993 elif name == 'Z': 994 self.op = sz.matrix 995 self.name = 'Z gate' 996 elif name == 'Rz': 997 assert 'phi' in kwargs 998 phi = kwargs['phi'] 999 i = 1j 1000 if (isinstance(phi,str) and phi == 'inf') or (phi == np.inf): 1001 phi = np.inf 1002 self.op = matrix(diag([1.,0])) 1003 else: 1004 self.op = matrix(diag([1,np.exp(i*phi)])) 1005 self.name = 'Rz (phi = %0.3f radians)'%phi 1006 elif name == 'XX': 1007 assert 'phi' in kwargs 1008 phi = kwargs['phi'] 1009 self.name = 'Ising (XX) gate (phi = %0.3f radians)'%phi 1010 R = gate('Rz',phi=phi) 1011 X = gate('X') 1012 I = np.eye(2) 1013 i = 1j 1014 self.op = matrix(np.block([[ I , -i*(X*R).op ], 1015 [-i*(R*X).op, I ]] ))/np.sqrt(2.) 1016 elif name in ['S','SWAP']: 1017 self.op = matrix([[1,0,0,0], 1018 [0,0,1,0], 1019 [0,1,0,0], 1020 [0,0,0,1]]) 1021 self.name = 'SWAP (S) gate' 1022 1023 self.meta = kwargs 1024 1025 def __repr__(self): 1026 s = '%s\n'%self.name 1027 return s + super(self.__class__,self).__repr__() 1028 1029 def copy(self): 1030 """make a copy of the gate 1031 """ 1032 new = gate(self.op.copy()) 1033 return new 1034 1035 def sqrt(self): 1036 '''create a square root gate from a gate. 1037 ''' 1038 s = sqrtm(self.matrix) 1039 r = gate(s) 1040 r.name = 'square root %s'%self.name 1041 return r 1042 1043 def controlled(self): 1044 '''convert a gate to a controlled gate 1045 (introduces one more qubit as the control bit) 1046 ''' 1047 n = self.op.shape[0] 1048 I = np.eye(n) 1049 r = gate(matrix(block_diag(I,self.op))) 1050 r.name = 'controlled %s'%self.name 1051 return r 1052 1053def controlled(A): 1054 '''controlled gate 1055 ''' 1056 n = A.op.shape[0] 1057 I = np.eye(n) 1058 r = operator(block_diag(I,A.op)) 1059 return r 1060 1061def Rz(phi): 1062 '''phase rotation gate 1063 ''' 1064 i = 1j 1065 r = operator(diag([1,np.exp(i*phi)])) 1066 return r 1067 1068def XX(phi): 1069 '''Ising (XX) gate 1070 ''' 1071 R = Rz(phi) 1072 I = np.eye(2) 1073 i = 1j 1074 xx = operator(np.block([[I, -i*(X*R).op], 1075 [-i*(R*X).op, I]]))/np.sqrt(2.) 1076 return xx 1077 1078if 'sphinx' not in sys.modules: 1079 sx = operator('sx') 1080 sy = operator('sy') 1081 sz = operator('sz') 1082 H = gate('H') 1083 X = gate('X') 1084 Y = gate('Y') 1085 Z = gate('Z') 1086 # H = (sx+sz).N 1087 # X = sx 1088 # Y = sy 1089 # Z = sz 1090 1091 NOT = X 1092 # sNOT = NOT.sqrt() 1093 # cNOT = NOT.controlled() 1094 sNOT = operator(sqrtm(NOT.op)) 1095 cNOT = controlled(NOT) 1096 1097 S = operator([[1,0,0,0], 1098 [0,0,1,0], 1099 [0,1,0,0], 1100 [0,0,0,1]]) 1101 1102 SWAP = S 1103 sSWAP = operator(sqrtm(SWAP.op)) 1104 cSWAP = operator(block_diag(np.eye(4),SWAP.op)) 1105 cX = cNOT 1106 I = np.eye(2) 1107 ccNOT = operator(block_diag(I,I,NOT.op)) 1108 1109 1110#====================== Tests ======================= 1111import sys 1112# https://stackoverflow.com/questions/5067604/determine-function-name-from-within-that-function-without-using-traceback 1113# for current func name, specify 0 or no argument. 1114# for name of caller of current func, specify 1. 1115# for name of caller of caller of current func, specify 2. etc. 1116 1117def currentFuncName(n=0): 1118 return sys._getframe(n + 1).f_code.co_name 1119 1120def head(): 1121 bling = '-'*60 1122 print ('<%s> %s'%(currentFuncName(1),bling)) 1123 1124def tail(): 1125 blin2 = '-'*30 1126 print ('<%s> %s PASSED %s'%(currentFuncName(1),blin2,blin2)) 1127 1128def test1(): 1129 head() 1130 uuu = ket('|uuu>') 1131 print ('%r has %d particles'%(uuu,uuu.n_particles)) 1132 print (uuu) 1133 tail() 1134 1135def test2(): 1136 head() 1137 a = ket('|ud>') + ket('|du>') 1138 print (a) 1139 b = a.H 1140 print (b) 1141 a.normalize() 1142 print (a) 1143 print (a.H) 1144 tail() 1145 1146def test3(): 1147 head() 1148 c = bra('<uuduud|') 1149 print ('*'*60) 1150 print (c.basis) 1151 print ('*'*60) 1152 print (c.bases) 1153 print ('*'*60) 1154 tail() 1155 1156def test4(): 1157 global u,d,A,r,l 1158 head() 1159 u = state([1,0]) 1160 d = state([0,1]) 1161 A = u*u.H - d*d.H 1162 r = (u+d).normalized() 1163 l = (u-d).normalized() 1164 A*r 1165 print (A*r,l) 1166 assert A*r == l 1167 tail() 1168 1169def test5(): 1170 head() 1171 u = state('|u>') 1172 d = state('|d>') 1173 sz = u*u.H - d*d.H 1174 s0 = u*u.H + d*d.H 1175 sz1 = sz**s0 1176 sz2 = s0**sz 1177 s = u**(u+d).normalized() 1178 sz1.measure(s) 1179 sz2.measure(s) 1180 1181 uu = u**u 1182 ud = u**d 1183 du = d**u 1184 dd = d**d 1185 s = (uu+ud).normalized() 1186 sz1.measure(s) # (+1,uu+ud) 1187 sz2.measure(s) # { (+1,uu) or (-1,ud) } 1188 s = (du+dd).normalized() 1189 sz1.measure(s) # (-1,du+dd) 1190 sz2.measure(s) # { (+1,uu) or (-1,dd) } 1191 1192 s = uu + dd 1193 s.normalize() 1194 (sz**s0).measure(s) # { (+1,uu} or {-1,dd})} 1195 s = (u+d)**(u+d) 1196 s.normalize() 1197 (sz**s0).measure(s) # { (+1,uu+ud ) or (-1,du+dd) } 1198 1199 tail() 1200 1201def test6(): 1202 head() 1203 u = ket([1,0]) 1204 d = ket([0,1]) 1205 s = (u**d - d**u).N 1206 Sz = u*u.H - d*d.H 1207 S0 = u*u.H + d*d.H 1208 s = ket(list(np.random.normal(size=(8)))).N 1209 c = s.correlation() 1210 print(c) 1211 P = [s.prob(Sz**S0**S0,ket(x)) for x in s.basis] 1212 print(np.sum(P)) 1213 s1 = (u**d - 0*d**u).N**(u**d + d**u).N 1214 s2 = (0*u**d - d**u).N**(u**d + d**u).N 1215 rho = 0.5*s1.density() + 0.5*s2.density() 1216 print(rho) 1217 tail() 1218 1219def all_tests(): 1220 test1() 1221 test2() 1222 test3() 1223 test4() 1224 test5() 1225 test6() 1226 test_spin() 1227 test_entangled() 1228 test_ptrace() 1229 test_readme() 1230 1231s0 = operator('s0') 1232sx = operator('sx') 1233sy = operator('sy') 1234sz = operator('sz') 1235u = state('|u>') 1236d = state('|d>') 1237 1238uu = u**u 1239ud = u**d 1240du = d**u 1241dd = d**d 1242# singlet state 1243s = (ud - du).N 1244# triplet states 1245t1 = uu 1246t2 = (ud + du).N 1247t3 = dd 1248 1249# spin of an electron 1250sx_e = sx/2. 1251sy_e = sy/2. 1252sz_e = sz/2. 1253spinx = (sx_e**s0) + (s0**sx_e) 1254spiny = (sy_e**s0) + (s0**sy_e) 1255spinz = (sz_e**s0) + (s0**sz_e) 1256spinx2 = spinx*spinx 1257spiny2 = spiny*spiny 1258spinz2 = spinz*spinz 1259spin2 = spinx2 + spiny2 + spinz2 1260 1261def test_spin(): 1262 """ calculate the spin of an electron and pairs of electrons 1263 """ 1264 head() 1265 print ('spin of the electron along z axis is 1/2:') 1266 print (u.H*sz_e*u) 1267 print ('magnitude of the spin is sqrt(3)/2:') 1268 print (' quantum number is J=1/2 so') 1269 print (' spin = sqrt( J*(J+1) ) = sqrt( (1/2)*(3/2) ) = 0.866') 1270 print (np.sqrt( u.H*(sx_e*sx_e + sy_e*sy_e + sz_e*sz_e)*u )) 1271 print ('magnitude of spin of the singlet state is zero:') 1272 print (np.sqrt( s.H*spin2*s )) 1273 print ('magnitudes of spins of the triplet states are sqrt(2):') 1274 print (' quantum number of paired triplet state is J=1 so') 1275 print (' spin = sqrt( J*(J+1) ) = sqrt(2) = 1.414') 1276 print (np.sqrt( t1.H*spin2*t1 )) 1277 print (np.sqrt( t2.H*spin2*t2 )) 1278 print (np.sqrt( t3.H*spin2*t3 )) 1279 print ('spins of triplet states projected onto z axis:') 1280 print (' these are the m quantum nummbers and should be +1, 0, -1') 1281 print (t1.H*spinz*t1) 1282 print (t2.H*spinz*t2) 1283 print (t3.H*spinz*t3) 1284 tail() 1285 1286def test_entangled(s = s): 1287 """ the argument can be any mixed state of two particles. 1288 the singlet state is the default argument. 1289 """ 1290 head() 1291 print ('(perhaps) entangled state:') 1292 print (s) 1293 print ('Alice and Bob:') 1294 print (s.density()) 1295 print ("Alice's view:") 1296 rhoA = ptrace(s.density(),[1]) 1297 rhoB = ptrace(s.density(),[0]) 1298 #rhoA,rhoB = s.density1() 1299 print (rhoA) 1300 print ("Bob's view:") 1301 print (rhoB) 1302 print ("entropy of Alices's view") 1303 S = entropy(rhoA) 1304 S_frac = entropy(rhoA,frac=True) 1305 print ('S = %.3f which is %.1f%% of max entropy'%(S,S_frac*100)) 1306 print ("entropy of Bob's view") 1307 S = entropy(rhoB) 1308 S_frac = entropy(rhoB,frac=True) 1309 print ('S = %.3f which is %.1f%% of max entropy'%(S,S_frac*100)) 1310 print ('tests for entanglement') 1311 rho = s.density() 1312 rhoA = s.density1(0) 1313 print ('density matrix trace test:') 1314 print ('trace(rho_Alice^2) =',np.trace(rhoA*rhoA)) 1315 print ('(if it is < 1, the particles are entangled)') 1316 print ('correlation tests <AB> - <A><B>') 1317 cor = [] 1318 scor = [] 1319 for a,sa in zip([sx,sy,sz],['x','y','z']): 1320 cor_row = [] 1321 scor_row = [] 1322 for b,sb in zip([sx,sy,sz],['x','y','z']): 1323 scor_row.append('c'+sa+sb) 1324 sigma = (a**s0) # measures first particle without affecting 2nd 1325 tau = (s0**b) # measures second particle without affecting first 1326 c = s.H*(sigma*tau)*s - (s.H*sigma*s)*(s.H*tau*s) 1327 cor_row.append(c) 1328 scor.append(scor_row) 1329 cor.append(cor_row) 1330 cor = np.matrix(cor).real 1331 scor = np.matrix(scor) 1332 print ('correlation ') 1333 print (scor) 1334 print (cor) 1335 # 1336 # cxx = s.H*((sx**s0)*(s0**sx))*s - (s.H*((sx**s0))*s)*(s.H*((s0**sx))*s) 1337 # cyy = s.H*((sy**s0)*(s0**sy))*s - (s.H*((sy**s0))*s)*(s.H*((s0**sy))*s) 1338 # czz = s.H*((sz**s0)*(s0**sz))*s - (s.H*((sz**s0))*s)*(s.H*((s0**sz))*s) 1339 # print 'correlation x,y,z = ',cxx,cyy,czz 1340 print ('(if any correlation != 0, the particles are entangled)') 1341 tail() 1342 1343def test_ptrace(): 1344 head() 1345 s = u**d - d**u 1346 s = (s**s).N 1347 den = ptrace(s.density(),[1,3]) 1348 den = ptrace(s.density(),[]) 1349 den = ptrace(s.density(),range(s.n_particles)) 1350 tail() 1351 1352def test_readme(): 1353 head() 1354 # from the README.rst 1355 print('Spin states') 1356 #from qspin2 import bra,ket,u,d,s0,sx,sy,sz 1357 u = ket('|u>') 1358 d = ket('|d>') 1359 print(u) 1360 print(d) 1361 print(u+d) 1362 i=1j 1363 print(u+i*d) 1364 print('Operators') 1365 print(sx) 1366 print(sy) 1367 print(sz) 1368 print(sz*u) 1369 print(sz*d) 1370 print(u.H*sz*u) 1371 print(u) 1372 print(u.phi) 1373 print('eigenvalues and eigenvectors') 1374 #import numpy as np 1375 print(sz) 1376 ev,evec = np.linalg.eig(sz.matrix) 1377 print(ev) 1378 print(evec) 1379 print(sx) # spin x 1380 ev, evec = np.linalg.eig(sx.matrix) 1381 print(ev) 1382 print(evec) 1383 ev, evec = sx.eig() 1384 print(ev) 1385 print(evec) 1386 print('conditional probabilities') 1387 l = (u+d).N 1388 print(np.abs(bra(l)*ket(u))**2) 1389 print(np.abs(l.H*u)**2) 1390 print(l.prob(sz,u)) 1391 print('string representation of state') 1392 set_base_repr('01') 1393 u = ket('0') 1394 d = ket('1') 1395 s = (u**d - d**u).N 1396 print(s) 1397 set_base_repr('arrow') 1398 u = ket(u'\u2193') 1399 d = ket(u'\u2191') 1400 s = (u**d-d**u).N 1401 print(s) 1402 set_base_repr('ud') 1403 u = ket('u') 1404 d = ket('d') 1405 1406 print('partial trace') 1407 s = (u**d - d**u).N 1408 rho = s.density() 1409 rhoA = ptrace(rho,[1]) 1410 print(rhoA) 1411 1412 print('entanglement') 1413 #from qspin2 import ket 1414 u = ket('u') 1415 d = ket('d') 1416 s = (u**d - d**u).N 1417 print(s.entangled()) 1418 tail()