"""
Module ``Ops.py`` defines operators, transform operators into their
symmetry-adapted forms in the presence of abelian quantum numbers, and write
operators to Fortran-readable files.
See the manual for more information about the handling of Operators in
MPSPyLib.
"""
import numpy as np
from math import sqrt
from copy import deepcopy
from collections import defaultdict
import random
import string
import hashlib
from sys import version_info
if(version_info[0] < 3):
from itertools import izip as myzip
else:
myzip = zip
localDimension=0
[docs]class OperatorList():
"""
Defines and handles action on operators for any simulation. Handles
transformation to symmetry-adapted operators and writing out the
fortran-compatible files. This class still provides the dictionary
operators on the operators (``__setitem__``, ``__getitem__``,
``__delitem__``, ``__iter__``, ``__contains__``, ``__len__``). If you
need more built-in dictionary methods, file a bug report.
**Arguments**
Operators : dict
The operators passed in a dictionary to the init method are saved
in the corresponding internal dictionary ``__oper__``.
**Variables**
has_complex : bool
flag if list contains any complex operators.
gabelian : list of strings
containing the names of the generators for the abelian symmetry. This
can (and is intended) to be updated before writing out operators.
gdiscrete : list of strings
containts the names of the generators for the discrete symmetry. This
can (and is intended) to be updated before writing out operators.
``__oper__`` : dict (interal use only)
dictionary with operators on the complete local Hilbert space.
``__qoper__`` : dict (internal use only)
dictionary with the symmetry-adapted operators. Deleted when updating
generators, generated before writing out operators.
**Details**
The following keys are blocked for naming operators, since they are all
used for the qOperators: ``q_d``, ``nqns``, ``npns``, ``q_numbers``,
``q_shifts``, ``q_scales``, ``transformation``. Using those keys will
return the value of the qOperator dictionary.
"""
def __init__(self, Operators):
self.__oper__ = {}
self.__qoper__ = None
self.has_complex = False
self.gabelian = []
self.gdiscrete = []
# Set existing operators
for key in Operators:
self.__oper__[key] = Operators[key]
[docs] def __setitem__(self, key, item):
"""
Setting matrix in usual dictionary fashion.
**Arguments**
key : str (or hashable)
identifier for matrix, intended to be a string with operator
name
item : 2d np.array
Matrix representing the operator in the local Hilbert space.
"""
if(key in ['q_d', 'nqns', 'npns', 'q_numbers', 'q_shifts',
'q_scales', 'transformation']):
raise KeyError("Key '" + str(key) + "' is reserved and cannot be" +
"used.")
if(np.any(np.iscomplex(item))): self.has_complex = True
self.__oper__[key] = item
[docs] def __getitem__(self, key):
"""
Return the matrix with the corresponding key in usual dictionary
fashion. Returns as well qOperator specific values for blocked
keys ``q_d``, ``nqns``, ``npns``, ``q_numbers``, ``q_shifts``,
``q_scales``, ``transformation``.
**Arguments**
key : str (or hashable)
identifier for matrix, intended to be a string with operator
name
"""
if(key in ['q_d', 'nqns', 'npns', 'q_numbers', 'q_shifts',
'q_scales', 'transformation']):
return self.__qoper__[key]
return self.__oper__[key]
[docs] def __delitem__(self, key):
"""
Delete operator in usual dictionary fashion. Does check if operator
is at present generator of symmtry. Does not update the flag
``has_complex``.
**Arguments**
key : str (or hashable)
identifier for matrix, intended to be a string with operator
name
"""
if(not key in self.__oper__):
print("Operator " + str(key) + " to be removed is not " +
"in Operators!")
return
if((key in self.gabelian) or (key in self.gdiscrete)):
raise Exception("Operator " + str(key) + " to be removed is " +
"generator of a symmetry an cannot be removed!")
del self.__oper__[key]
[docs] def __iter__(self):
"""
Defining iteration in the usual dictionary fashion.
"""
return iter(self.__oper__)
[docs] def __len__(self):
"""
Returning the length in the usual dictionary fashion returning the
number of elements in the dictionary.
"""
return len(self.__oper__)
[docs] def keys(self):
"""
Built-in methods of dictionaries ``keys()``.
"""
return self.__oper__.keys()
[docs] def __contains__(self, key):
"""
Checking if certain key exists in the keys of the dictionary of
operators.
**Arguments**
key : str (or hashable)
identifier for matrix, intended to be a string with operator
name
"""
return (key in self.__oper__)
[docs] def has_key(self, key):
"""
Built-in method of dictionaries ``has_key()``.
**Arguments**
key
passed to built-in method ``has_key()``.
"""
return self.__oper__.has_key(key)
[docs] def remove_operator(self, key):
"""
Remove an operator from the dictionary of operators.
**Arguments**
key : str
string identifier of the operator.
"""
self.__delitem__(key)
[docs] def add_abelian(self, generators):
"""
Add the abelian symmetries. If list differs from previous list
of abelian generators, the qOperators are deleted.
**Arguments**
generators : list of strings
contains the names of the diagonal generators of the symmetries.
"""
# Ensure that it's a list
if(not isinstance(generators, list)): generators = [generators]
# Check for fast return
if(generators == self.gabelian): return
# Reset qdict
self.__qoper__ = None
for key in generators:
if(not key in self.__oper__):
raise Exception("Not all Abelian generators are found " +
"in Operators.")
self.gabelian = deepcopy(generators)
return
[docs] def add_discrete(self, generators):
"""
Add the discrete symmetries. If list differs from previous list
of discrete generators, the qOperators are deleted.
**Arguments**
generators : list of strings
contains the names of the diagonal generators of the symmetries.
"""
# Ensure generators are a list
if(not isinstance(generators, list)): generators = [generators]
# Check for fast return
if(generators == self.gdiscrete): return
# Reset qdict
self.__qoper__ = None
for key in generators:
if(not key in self.__oper__):
raise Exception("Not all discrete generators are found " +
"in Operators.")
self.gdiscrete = deepcopy(generators)
return
[docs] def get_nqns_npns(self):
"""
The number of abelian and discrete symmetries is defined over the
length of the list (assuming no duplicates).
"""
return len(self.gabelian), len(self.gdiscrete)
[docs] def get_hash(self):
"""
Hashing the symmetries of the system.
"""
hh = hashlib.sha512()
hh.update(str.encode(str(self.gabelian) + str(self.gdiscrete)))
return hh.hexdigest()
[docs] def check_symmetry(self, HamMPO, counter, eps=1e-8):
"""
Check for every generator that it obeys the symmetry within the
MPO given.
**Arguments**
HamMPO : MPO or list of MPOs
Check with built-in method the symmetry with corresponding
generators.
counter : int
if HamMPO is list, check element counter, otherwise not
referenced.
eps : float, optional
Accept non-zero commutators up to `eps`, default to
:math:`10^{-8}`.
"""
# Assume symmetries are true
has_sym = True
# Check for the abelian symmetries
for key in self.gabelian:
if(hasattr(HamMPO, '__len__')):
has_sym = (has_sym and
HamMPO[counter].has_symmetry(self.__oper__[key],
self.has_complex,
False, eps=eps))
else:
has_sym = (has_sym and
HamMPO.has_symmetry(self.__oper__[key],
self.has_complex,
False, eps=eps))
# Check for the discrete symmetries
for key in self.gdiscrete:
if(hasattr(HamMPO, '__len__')):
has_sym = (has_sym and
HamMPO[counter].has_symmetry(self.__oper__[key],
self.has_complex,
True, eps=eps))
else:
has_sym = (has_sym and
HamMPO.has_symmetry(self.__oper__[key],
self.has_complex,
True, eps=eps))
if(not has_sym):
print("")
print("Warning: Check your symmetry!")
print("")
return
[docs] def check(self, Op, Add=True, tol=1e-6):
"""
Check if an operator op is proportional to one in the operator list
Operators. If it is not, Add=True specifies to add it to the list with
a random string sequence identifier. Function returns key for operator
and the scaling in a list.
**Arguments**
Op : numpy 2d array
Check if this operator has an proportional operator in the list
of operators.
Add : boolean, optional
Add the operator if (True), otherwise not.
Default to True.
tol : float, optional
Tolerance for accepting operators as proportional.
Default to 1e-6
"""
# Reshape into vector, get non-zero indices, select non-zero entries
tmp1 = np.reshape(Op, (np.prod(Op.shape)))
idx1 = np.nonzero(tmp1)[0]
tmp1 = tmp1[idx1]
for key in sorted(self.__oper__.keys()):
# Do same steps for existing operators
tmp2 = np.reshape(self.__oper__[key], (np.prod(Op.shape)))
idx2 = np.nonzero(tmp2)[0]
tmp2 = tmp2[idx2]
# Catch operators with zeros only
if((idx1.shape[0] == 0) and (idx2.shape[0] == 0)):
return [key, 0.0]
# Check their positions
if(idx1.shape[0] != idx2.shape[0]): continue
if(np.any(idx1 != idx2)): continue
# Get scaling factor if possible
tmp = tmp1 / tmp2
if(np.any(np.abs(tmp - tmp[0]) > tol)): continue
# Matching, but no complex weight allowed:
if(np.imag(tmp[0]) != 0.0): continue
return [key, np.real(tmp[0])]
# No match found in dictionary (hash entries for name)
hstring = ''
for elem in np.reshape(Op, (np.prod(Op.shape))):
hstring += '%30.15e%30.15e'%(np.real(elem), np.imag(elem))
hh = hashlib.sha512()
hh.update(str.encode(hstring))
name = hh.hexdigest()[:8]
# Just in case 8 digits are not enough, 32 should be
if(name in self.__oper__.keys()): name = hh.hexdigest()[:32]
if(Add):
self.__oper__[name] = Op
return [name, 1.0]
[docs] def write(self, filestub, PostProcess=False):
"""
Write out a file defining the operators for use by the Fortran
routines It passes back a dict of operators which are keyed by an
increasing integer instead of the human-readable keys and the
concatenated file name.
**Arguments**
fileStub : string
base for the filename concatenated with `Ops.dat`.
PostProcess : logical, optional
If `PostProcess=True`, the operators are not written to file
and only the mapping between keys and integers is done.
Default to False
**Details**
The file (Operators) written has the following format:
1) Number of operators, number of rows and columns in the local
Hilbert space.
2) Matrix elements of each operator in one line
3) Position of the identity operator.
The file (qOperators) written has the following format:
1) Number of operators, number of symmetries (nqs)
2) For each operator, first number of tuples pairs
3) Indices for mapping of each symmetry [2 * nqs integer]
4) Indices for the possible degenerate subspace [2 integer]
5) Value (real or complex entry, either way in one line)
6) ... untill all tuples of this operator
7) ... untill all operators
"""
has_symm = (len(self.gabelian) + len(self.gdiscrete) > 0)
if(has_symm and (self.__qoper__ is None)):
self.OperatorstoQOperators()
if(self.has_complex and has_symm):
# Complex operators and symmetries
return self.write_q_complex(filestub, PostProcess=PostProcess)
elif(self.has_complex):
# Complex operators and no symmetries
return self.write_complex(filestub, PostProcess=PostProcess)
elif(has_symm):
# Real operators and symmetries
return self.write_q_real(filestub, PostProcess=PostProcess)
else:
# Real operators and no symmetries
return self.write_real(filestub, PostProcess=PostProcess)
[docs] def OperatorstoQOperators(self):
"""
Return the non-trivial elements of the operators Operators acting on the
degeneracy spaces based on the Abelian generators OutGenerators. The
generators are assumed to be diagonal. The functions returns a dictionary
for the `QOperators` (see Details for description of the dictionary).
**Details**
1) Check for each generator that it is diagonal
2) Shift values on the diagonal such that the minimum is zero
3) Scale according to the minimal difference (only considering
the neighboring diagonal elements; for zero=1e-10 scaling by one).
4) Build up unique identifiers: tuple of the (shfted scaled) quantum number
of the same i-th element on the diagonal for different generators.
5) Based on the previous step, we can calculate the degeneracy dimension
of each tuple
6) Generating a hash table: Hash all the tuples and then sort them
descending. Use the Square Root of Primers hash defined analogue in
fortran. Then create a list of desceding tuples eleminating
multiple equal entries.
7) Connecting bases (transfer matrix :math:`T`). The transformation into
the new basis for any operator :math:`O` is :math:`O' = T^t O T`.
8) When transforming each operator into the new basis, the operator is
block-shaped. Each block is extracted for the corresponding quantum
numbers. The size of the block is determined by the degeneracy of
the two quantum numbers.
There dictionary has the following keys:
* `nqns` : number of generators (integer)
* `q_shifts` : list of shifts for each generator
* `q_scales` : list of scaling factors for each generator
* `q_d` : dictionary of degeneracy dimension for tuples
* `q_numbers` : the sorted tuples a set (each entry appears only once
* `transformation` : this is the transfer matrix to be applied to the
operators.
* keys of `Operators` : dictionary containing a) the key `q` for the
quantum number in a list. Each entry in the list is a list of the
two involved quantum numbers. b) the key 'Blocks' for the corresponding
blocks in a list. Each entry in the list is a matrix.
See the manual for more information about symmetry-adapted operators.
"""
self.__qoper__ = {}
# In previous versions there were checks ...
# 1) No abelian and discrete symmetries at the same time
# 2) No discrete symmetries at all
# 3) Not more than one discrete symmetry
# Get list `Generators` with the generating operator for each
# symmetry defined in `OutGenerators`
Generators = []
# Add abelian symmetries
self.__qoper__['nqns'] = len(self.gabelian)
for key in self.gabelian:
Generators.append(deepcopy(self.__oper__[key]))
# Add discrete symmetries
self.__qoper__['npns'] = len(self.gdiscrete)
for key in self.gdiscrete:
Generators.append(deepcopy(self.__oper__[key]))
# Define auxiliary lists
shifts = []
scales = []
for Generator in Generators:
a = deepcopy(Generator)
# Check that generators are diagonal and real-valued
for i in range(len(a)):
if(np.iscomplex(a[i, i])):
# Check that diagonal is real (checks entry, not type)
raise ValueError("Complex entry on diagonal of generator!")
a[i,i]=0.0
if np.any(a!=0.0):
raise Exception("Generator passed in to OperatorstoQOperators" +
" is not diagonal!")
# Find the minimum value for that generator, define a shift so
# that this minimum value is zero
dgen = np.diagonal(Generator)
mval = np.min(dgen)
shifts.append(mval)
dgen = dgen - mval
# Find the differences between values, use to define a scaling
# factor such that all quantum numbers are integers
diffs = []
for i in range(len(dgen)-1):
diffs.append(abs(dgen[i] - dgen[i+1]))
diffsnozero = [x for x in diffs if abs(x)>1E-10]
if not diffsnozero:
# True for empty list
scales.append(1.0)
else:
scales.append(min(diffsnozero))
for j in range(len(Generators)):
for i in range(len(Generators[j])):
Generators[j][i,i]=(Generators[j][i,i]-shifts[j])/scales[j]
# LOCAL Generators are now diagonal lists of integers increasing from 0.
# (generators in the QOperators list are not shifted and scaled)
# save shifts and scales for fixing quantum numbers later.
self.__qoper__['q_shifts'] = shifts
self.__qoper__['q_scales'] = scales
# compute unique identifiers for kronecker products of quantum numbers
indices = []
for i in np.diag(Generators[0]):
indices.append([int(i)])
for j in range(1,len(Generators)):
count=0
for i in np.diag(Generators[j]):
indices[count].append(int(i))
count=count+1
# store as tuples so that we can use them as dict keys
for count in range(Generators[0].shape[0]):
indices[count] = tuple(indices[count])
#compute degeneracy dimensions
degdim = defaultdict(int)
for elt in indices:
degdim[elt] += 1
self.__qoper__['q_d'] = degdim
# sort basis according to the hashes of quantum number tuples. This
# ensures all degenerate irreps are adjacent and simplifies sorting in
# fortran. Note that we sort in descending order here, which translates
# to ascending order in fortran, as we add new operator elements to the
# head of the list. get indices of the permutation leading to the
# sorted vector
hashes = []
for index in indices:
hashes.append(SROPHash(index))
# Sort by hashes
permindices = sorted(range(len(indices)), key=hashes.__getitem__,
reverse=True)
# get only a set of the indices-we are not interested in repeated tuples
littlehashes = []
sortedindices = []
for i in range(Generators[0].shape[1]):
if not littlehashes:
# Initializing values for first iterations
curhash=SROPHash(indices[permindices[i]])
littlehashes.append(curhash)
sortedindices.append(indices[permindices[i]])
lasthash=curhash
else:
curhash=SROPHash(indices[permindices[i]])
if curhash!=lasthash:
littlehashes.append(curhash)
sortedindices.append(indices[permindices[i]])
lasthash=curhash
self.__qoper__['q_numbers']=sortedindices
# because the abelian generators are already diagonal, the unitary matrix
# connecting the bases is just a column permutation of the identity matrix
transmatrix=np.zeros((Generators[0].shape[0],Generators[0].shape[1]))
for i in range(Generators[0].shape[1]):
transmatrix[permindices[i],i]=1.0
self.__qoper__['transformation']=transmatrix
tuplesub = lambda xs, ys: tuple(x - y for x, y in myzip(xs, ys))
# transform ALL operators to the new basis and write data to qOperators
for key in self.__oper__:
self.__qoper__[key]={}
self.__qoper__[key]['q']=[]
self.__qoper__[key]['Blocks']=[]
thisOperators=np.dot(np.dot(transmatrix.transpose(),
self.__oper__[key]), transmatrix)
checkeri=[]
checkerj=[]
runningj=0
for jindex in sortedindices:
dj=degdim[jindex]
runningi=0
for iindex in sortedindices:
di=degdim[iindex]
# more pythonic comparison?
if np.any(np.abs(thisOperators[runningi:runningi+di,
runningj:runningj+dj])>1E-10):
self.__qoper__[key]['q'].append([iindex,jindex])
self.__qoper__[key]['Blocks'].append(
thisOperators[runningi:runningi+di,runningj:runningj+dj])
checkeri.append(jindex)
checkerj.append(iindex)
runningi=runningi+di
runningj=runningj+dj
# ensure that all operators are injective, that is, take tuples
# to unique tuples. This prevents us from having to traverse lists
# multiple times in Fortran.
if((len(set(checkeri)) != len(checkeri)) or
(len(set(checkerj)) != len(checkerj))):
print("Quantum number transfers:", checkeri, checkerj)
raise Exception("Operator " + key + " is not injective! See" +
" the manual for more detail")
return
[docs] def write_real(self, fileStub, PostProcess=False):
"""
Write out a file defining the operators for use by the Fortran routines
It passes back a dict of operators which are keyed by an increasing
integer instead of the human-readable keys and the complete filename.
**Arguments**
fileStub : string
base for the filename concatenated with `Ops.dat`.
PostProcess : logical, optional
If `PostProcess=True`, the operators are not written to file
and only the mapping between keys and integers is done.
Default to False
"""
IntOperators = {}
count = 0
for key in sorted(self.__oper__.keys()):
IntOperators[count] = self.__oper__[key]
IntOperators[key] = count+1
count = count + 1
dl = int(IntOperators[0].shape[0])
dr = int(IntOperators[0].shape[1])
Opsfilename = fileStub + 'Ops.dat'
if(not PostProcess):
Opsfile = open(Opsfilename, 'w')
Opsfile.write('%16i'%(count) + '%16i'%(dl) + '%16i'%(dr) + '\n')
for i in range(count):
ostring=''
for j in np.reshape(IntOperators[i],dl*dr):
ostring+='%30.15E'%(j)
Opsfile.write(ostring+'\n')
Opsfile.write('%16i'%(IntOperators['I'])+'\n')
Opsfile.close()
return IntOperators, Opsfilename
[docs] def write_complex(self, fileStub, PostProcess=False):
"""
Write out a file defining the complex operators for use by the Fortran
routines. It passes back a dict of operators which are keyed by an
increasing integer instead of the human-readable keys and the complete
filename.
**Arguments**
fileStub : string
base for the filename concatenated with `Ops.dat`.
PostProcess : logical, optional
If `PostProcess=True`, the operators are not written to file
and only the mapping between keys and integers is done.
Default to False
"""
IntOperators = {}
count = 0
for key in sorted(self.__oper__.keys()):
IntOperators[count] = self.__oper__[key]
IntOperators[key] = count + 1
count = count + 1
dl = int(IntOperators[0].shape[0])
dr = int(IntOperators[0].shape[1])
Opsfilename = fileStub + 'Ops.dat'
if(not PostProcess):
Opsfile = open(Opsfilename, 'w')
Opsfile.write('%16i'%(count) + '%16i'%(dl) + '%16i'%(dr) + '\n')
for i in range(count):
ostring = ''
for j in np.reshape(IntOperators[i], dl*dr):
ostring+='%30.15E'%(np.real(j)) + '%30.15E'%(np.imag(j))
Opsfile.write(ostring+'\n')
Opsfile.write('%16i'%(IntOperators['I'])+'\n')
Opsfile.close()
return IntOperators, Opsfilename
[docs] def write_q_real(self, fileStub, PostProcess=False):
"""
Write out a file defining the qoperators for use by the Fortran routines
Also passes back a dict of operators which are keyed by an increasing integer
instead of the human-readable keys.
**Arguments**
fileStub : string
base for the filename concatenated with 'Ops.dat'
PostProcess: logical, optional
If `PostProcess=True`, the operators are not written to file
and only the mapping between keys and integers is done.
Default to False
"""
forbiddenkeys = ['q_d', 'nqns', 'npns', 'q_numbers', 'q_shifts',
'q_scales', 'transformation']
IntOperators = {}
count = 0
for key in sorted(self.__qoper__.keys()):
if key not in forbiddenkeys:
IntOperators[count]=self.__qoper__[key]
IntOperators[key]=count+1
count=count+1
else:
# do not raise error: forbidden keys are only no matrices
pass
Opsfilename=fileStub+'Ops.dat'
if not PostProcess:
Opsfile=open(Opsfilename,'w')
# number of operators, number of quantum numbers
Opsfile.write('%16i'%(count) + '\n')
for i in range(count):
Operator=IntOperators[i]
# number of index tuple pairs
Opsfile.write('%16i'%(len(Operator['q']))+'\n')
for q in range(len(Operator['q'])):
#write out index tuple pair
tstr=''
for j in Operator['q'][q][0]:
# the j^{th} element of the q^{th} i-tuple
tstr+='%16i'%(j)
for j in Operator['q'][q][1]:
# the j^{th} element of the q^{th} j-tuple
tstr+='%16i'%(j)
Opsfile.write(tstr+'\n')
ostring=''
#write out degeneracy dimensions
Opsfile.write('%16i'%(self.__qoper__['q_d'][Operator['q'][q][0]]) +
'%16i'%(self.__qoper__['q_d'][Operator['q'][q][1]])+'\n')
#write out block hami
for j in np.reshape(Operator['Blocks'][q],
self.__qoper__['q_d'][Operator['q'][q][0]] *
self.__qoper__['q_d'][Operator['q'][q][1]]):
ostring+='%30.15E'%(j)
Opsfile.write(ostring+'\n')
#output which operator is the identity
Opsfile.write('%16i'%(IntOperators['I'])+'\n')
Opsfile.close()
return IntOperators, Opsfilename
[docs] def write_q_complex(self, fileStub, PostProcess=False):
"""
Currently not implemented - Write out a file defining the qoperators for
use by the Fortran routines.
Also passes back a dictionary of operators which are keyed by an increasing
integer instead of the human-readable keys.
**Arguments**
fileStub : string
base for the filename concatenated with 'Ops.dat'
PostProcess: logical, optional
If `PostProcess=True`, the operators are not written to file
and only the mapping between keys and integers is done.
Default to False
"""
forbiddenkeys = ['q_d', 'nqns', 'npns', 'q_numbers', 'q_shifts',
'q_scales', 'transformation']
IntOperators = {}
count = 0
for key in sorted(self.__qoper__.keys()):
if key not in forbiddenkeys:
IntOperators[count]=self.__qoper__[key]
IntOperators[key]=count+1
count=count+1
else:
# do not raise error: forbidden keys are only no matrices
pass
Opsfilename=fileStub+'Ops.dat'
if not PostProcess:
Opsfile=open(Opsfilename,'w')
# number of operators, number of quantum numbers
Opsfile.write('%16i'%(count) + '\n')
for i in range(count):
Operator=IntOperators[i]
# number of index tuple pairs
Opsfile.write('%16i'%(len(Operator['q']))+'\n')
for q in range(len(Operator['q'])):
# write out index tuple pair
tstr=''
for j in Operator['q'][q][0]:
# the j^{th} element of the q^{th} i-tuple
tstr+='%16i'%(j)
for j in Operator['q'][q][1]:
# the j^{th} element of the q^{th} j-tuple
tstr+='%16i'%(j)
Opsfile.write(tstr+'\n')
ostring=''
# write out degeneracy dimensions
Opsfile.write('%16i'%(self.__qoper__['q_d'][Operator['q'][q][0]]) +
'%16i'%(self.__qoper__['q_d'][Operator['q'][q][1]]) +
'\n')
# write out block hami (row major!)
for j in np.reshape(Operator['Blocks'][q],
self.__qoper__['q_d'][Operator['q'][q][0]] *
self.__qoper__['q_d'][Operator['q'][q][1]]):
ostring += '%30.15E'%(np.real(j)) + '%30.15E'%(np.imag(j))
Opsfile.write(ostring + '\n')
# output which operator is the identity
Opsfile.write('%16i'%(IntOperators['I']) + '\n')
Opsfile.close()
return IntOperators, Opsfilename
[docs]def getLocaldim():
"""Return the local dimension"""
global localDimension
return localDimension
[docs]def BuildSpinOperators(spin):
""" Build the spin operators sz, splus, sminus, and identity (I).
Spin is the value of the total spin.
**Arguments**
spin : integer or float
spin value, allowing integer or n * 0.5 with n an integer
**Details**
+-----------------+----------------------+
| Tag operator | Operator |
+=================+======================+
| I | :math:`\hat{I}` |
+-----------------+----------------------+
| sz | :math:`\hat{S^{z}}` |
+-----------------+----------------------+
| splus | :math:`\hat{S^{+}}` |
+-----------------+----------------------+
| sminus | :math:`\hat{S^{-}}` |
+-----------------+----------------------+
"""
if int(2*spin)==2*spin:
spinSize=int(2*spin+1)
else:
raise Exception("2*spin must be an integer in BuildSpinOperators!")
global localDimension
localDimension=spinSize
identity=np.eye(spinSize)
sz=np.zeros((spinSize,spinSize))
splus=np.zeros((spinSize,spinSize))
for i in range(spinSize):
sz[i,i]=spin-i
for i in range(spinSize-1):
splus[i,i+1]=sqrt((spin+(spin-i))*(spin+1-(spin-i)))
sminus=splus.T
Operators={'sz':sz,'splus':splus,'sminus':sminus,'I':identity}
return OperatorList(Operators)
[docs]def BuildBoseOperators(Nmax, Nmin=0, spin=0, nFlavors=1, maxFilling=None):
"""
Build the bose operators b, bdagger, nbtotal, and identity (I).
If spin is nonzero then also build bsz, bsplus, and bsminus.
If nFlavors is greater than one then the bose operators contain
a flavor index Operators['b_flavor'] where flavor=0...nFlavors-1.
If spin is nonzero and the bosons are flavorless then the bose operators
contain a spin index Operators['b_m'] where m=0...2*spin. If both spin and
Flavors are nonzero then a set of spinful+flavorful bose operators are built
with the indexing Operators['b_flavor_m'] where flavor and m are as above.
**Arguments**
Nmax : integer
the maximum total number of bosons allowed on-site.
Nmin : integer, optional
the minimum total number of bosons allowed on-site.
default to 0
spin : integer or float, optional
the spin of the bosons. Multiple of 0.5
default to 0
nFlavors : integer, optional
the number of "flavors" of boson-basically the number
of non-spin internal degrees of freedom.
default to 1
maxFilling : optional
the number of bosons in any one spin or flavor
component allowed on-site. This is equal to Nmax for scalar bosons.
default to ``None``
**Details**
+-----------------+---------------------------+
| Tag operator | Operator |
+=================+===========================+
| I | :math:`\hat{I}` |
+-----------------+---------------------------+
| b | :math:`\hat{b}` |
+-----------------+---------------------------+
| bdagger | :math:`\hat{b}^{\dagger}` |
+-----------------+---------------------------+
| nbtotal | :math:`\hat{n}` |
+-----------------+---------------------------+
"""
global localDimension
if maxFilling is None:
maxFilling=Nmax
if int(2*spin)==2*spin:
spinSize=int(2*spin+1)
else:
raise Exception("2*spin must be an integer in BuildBoseOperators!")
ninternal=spinSize*nFlavors
if ninternal==1:
if maxFilling<Nmax:
raise Exception("maxFilling should always=Nmax for scalar bosons!")
Operators=BuildScalarOperators(Nmax,Nmin)
localDimension=Nmax-Nmin+1
else:
if spinSize==1:
Operators=BuildFlavorfulBoseOperators(Nmax,Nmin,nFlavors,maxFilling)
localDimension=0
for n in range(Nmin,Nmax+1):
localDimension=localDimension+mult(n,nFlavors,maxFilling)
else:
Operators=BuildSpinfulBoseOperators(Nmax,Nmin,spin,nFlavors,maxFilling)
localDimension=0
for n in range(Nmin,Nmax+1):
localDimension=localDimension+mult(n,spinSize*nFlavors,maxFilling)
return OperatorList(Operators)
[docs]def BuildScalarOperators(Nmax, Nmin, phaseOp=False):
"""
Build particle operators with no internal degrees of freedom. Intended
for internal use only.
**Arguments**
Nmax : integer
the maximum total number of bosons allowed on-site.
Nmin : integer, optional
the minimum total number of bosons allowed on-site.
default to 0
phaseOp : bool, optional
Flag for using phases for fermions originating in the Jordan-Wigner
transformation (``True``). If not fermions, use ``False``.
Default to ``False``.
"""
localSize=Nmax-Nmin+1
identity=np.eye(localSize)
b=np.zeros((localSize,localSize))
for i in range(localSize-1):
b[i,i+1]=sqrt((Nmin+i+1))
bdagger=b.T
ntotal=np.dot(bdagger,b)
ntotal[0,0]=Nmin #Correct for the shifted vaccum
if phaseOp:
phaseOp=np.zeros((localSize,localSize))
for i in range(localSize):
phaseOp[i,i]=(-1.0)**(ntotal[i,i])
Operators = {'f' : b, 'fdagger' : bdagger, 'I' : identity,
'nftotal' : ntotal, 'FermiPhase' : phaseOp}
else:
Operators={'b':b,'bdagger':bdagger,'I':identity, 'nbtotal':ntotal}
return Operators
[docs]def steplexVec(lexvec,n):
"""
Step the vector lexvec to the next inferior value in lexicographic order.
**Arguments**
lexvec : tba
n : int
number of particles
"""
#if the last component is n, there are no other components inferior to it
if lexvec[-1]==n:
return True
#find the last k<numFlavors-1, such that n_k!=0
for i in range(1,lexvec.size):
if lexvec[lexvec.size-1-i]!=0:
break
k=lexvec.size-1-i
#n_{1:k-1} is unchanged
#n_{k}=n_{k}-1
lexvec[k]=lexvec[k]-1
#n_{k+1}=N-\sum_{i=1}^{k}n_i
lexvec[k+1]=n-sum(lexvec[0:k+1])
#n_{k+2} onwards zero
lexvec[k+2:]=0
return False
[docs]def mult(n, nFlavors, maxFilling):
"""
Find the multiplicity of a Hilbert space returned as integer.
**Arguments**
n : int
number of particles
nFlavors : int
number of flavors
maxFilling : int
Local Hilbert space has at most maxFilling particles in any
one flavor.
"""
count=0
lexvec=np.zeros(nFlavors,dtype=int)
lexvec[0]=n
while True:
if all(lexvec<=maxFilling):
count=count+1
skip=steplexVec(lexvec,n)
if skip:
break
return count
[docs]def flavorfulHilbertspace(Nmin, Nmax, nFlavors, maxFilling):
"""
Build and index a Hilbert space for multiple flavors of particles.
The function returns the size of the hilbert space and a set of
vectors ``stateList[i,:] = [n_{f1}, n_{f2}, ...]`` giving the
occupations of the various flavors in state i. Intended for internal
use only.
**Arguments**
Nmin : int
Minimum number of particles.
Nmax : int
Maximum number of particles.
nFlavors : int
Number of different types of particles.
maxFilling : int
Local Hilbert space has at most maxFilling particles in any one
flavor component.
"""
#Find multiplicity
localSize=0
for n in range(Nmin,Nmax+1):
localSize=localSize+mult(n,nFlavors,maxFilling)
stateList=np.zeros((localSize,nFlavors))
count=0
for n in range(Nmin,Nmax+1):
#begin with completely superior vector
lexvec=np.zeros(nFlavors,dtype=int)
lexvec[0]=n
#step the vector in lexicographic order
while True:
if all(lexvec<=maxFilling):
stateList[count,:]=lexvec
count=count+1
skip=steplexVec(lexvec,n)
if skip:
break
return localSize,stateList
[docs]def BosedestructionFromList(stateList,i):
"""
Return a bosonic destruction operator from a stateList. Intended for
internal use only.
**Arguments**
stateList : tba
i : tba
"""
ls=stateList.shape[0]
mat=np.zeros((ls,ls))
for j in range(ls):
newstate=deepcopy(stateList[j,:])
newstate[i]=newstate[i]-1
for k in range(ls):
if all(newstate==stateList[k,:]):
mat[k,j]=sqrt(stateList[j,i])
return mat
[docs]def BuildFlavorfulBoseOperators(Nmax, Nmin, nFlavors, maxFilling):
"""
Build Bose operators with flavor degrees of freedom. Intended for internal
use only.
**Arguments**
Nmax : int
Maximum number of particles.
Nmin : int
Minimum number of particles.
nFlavors : int
Number of different types of particles.
maxFilling : int
Local Hilbert space has at most maxFilling particles in any one
flavor component.
"""
#build and index Hilbert space
localSize,stateList=flavorfulHilbertspace(Nmin,Nmax,nFlavors,maxFilling)
Operators={}
identity=np.eye(localSize)
ntotal=np.zeros((localSize,localSize))
for i in range(nFlavors):
Operators['b_'+str(i)]=BosedestructionFromList(stateList,i)
Operators['bdagger_'+str(i)]=Operators['b_'+str(i)].T
Operators['nb_'+str(i)]=np.dot(Operators['bdagger_'+str(i)],
Operators['b_'+str(i)])
Operators['nb_'+str(i)][i,i]=Nmin#correct for shifted vacuum
ntotal=ntotal+Operators['nb_'+str(i)]
Operators['I']=identity
Operators['nbtotal']=ntotal
return Operators
[docs]def BuildSpinfulBoseOperators(Nmax,Nmin,spin,nFlavors,maxFilling):
"""
Build Bose operators with spin and flavor degrees of freedom. Intended
for internal use only.
**Arguments**
Nmax : int
Maximum number of particles.
Nmin : int
Minimum number of particles.
spin : float (half-integer)
Spin for bosons (there is only one spin for all flavors).
nFlavors : int
Number of different types of particles.
maxFilling : int
Local Hilbert space has at most maxFilling particles in any one
flavor component.
"""
SpinOperators=BuildSpinOperators(spin)
spinSize=int(2*spin+1)
idofSize=spinSize*nFlavors
FlavoredOperators=BuildFlavorfulBoseOperators(Nmax, Nmin, idofSize,maxFilling)
Operators={}
bdagger=[]
b=[]
for spinindex in range(spinSize):
for flavorindex in range(nFlavors):
if nFlavors>1:
Operators['b_'+str(flavorindex)+'_'+str(spinindex)]=FlavoredOperators['b_'+str(spinindex*nFlavors+flavorindex)]
b.append(Operators['b_'+str(flavorindex)+'_'+str(spinindex)])
Operators['bdagger_'+str(flavorindex)+'_'+str(spinindex)]=FlavoredOperators['bdagger_'+str(spinindex*nFlavors+flavorindex)]
bdagger.append(Operators['bdagger_'+str(flavorindex)+'_'+str(spinindex)])
Operators['nb_'+str(flavorindex)+'_'+str(spinindex)]=FlavoredOperators['nb_'+str(spinindex*nFlavors+flavorindex)]
else:
Operators['b_'+str(spinindex)]=FlavoredOperators['b_'+str(spinindex*nFlavors+flavorindex)]
b.append(Operators['b_'+str(spinindex)])
Operators['bdagger_'+str(spinindex)]=FlavoredOperators['bdagger_'+str(spinindex*nFlavors+flavorindex)]
bdagger.append(Operators['bdagger_'+str(spinindex)])
Operators['nb_'+str(spinindex)]=FlavoredOperators['nb_'+str(spinindex*nFlavors+flavorindex)]
Operators['nbtotal']=FlavoredOperators['nbtotal']
Operators['I']=FlavoredOperators['I']
localSize=Operators['I'].shape[0]
sz=np.zeros((localSize,localSize))
splus=np.zeros((localSize,localSize))
sminus=np.zeros((localSize,localSize))
for i in range(spinSize):
for k in range(nFlavors):
for j in range(spinSize):
sz=sz+SpinOperators['sz'][i,j]*np.dot(bdagger[i*nFlavors+k],b[j*nFlavors+k])
splus=splus+SpinOperators['splus'][i,j]*np.dot(bdagger[i*nFlavors+k],b[j*nFlavors+k])
sminus=sminus+SpinOperators['sminus'][i,j]*np.dot(bdagger[i*nFlavors+k],b[j*nFlavors+k])
Operators['bsz']=sz
Operators['bsplus']=splus
Operators['bsminus']=sminus
return Operators
[docs]def BuildFermiOperators(Nmax=1, Nmin=0, spin=0, nFlavors=1):
"""
Build the Fermi operators f, fdagger, nftotal, FermiPhase, and identity (I).
Also build the "phased" operators fP=f.FermiPhase and fdaggerP=fdagger.
FermiPhase provided they are linearly independent from f and fdagger.
If spin is nonzero then also build fsz, fsplus, and fsminus.
If nFlavors is greater than one then the fermi operators contain
a flavor index Operators['f_flavor'] where flavor=0...nFlavors-1.
If spin is nonzero and the fermions are flavorless then the fermi operators
contain a spin index Operators['f_spin'] where spin=0...2*spin.
If both spin and nFlavors are nontrivial then a set of spinful+flavorful
fermi operators are built with the indexing Operators['b_flavor_spin'] where
flavor and spin are as above.
**Arguments**
Nmax : integer, optional
the maximum total number of fermions allowed on-site.
default to 1
Nmin : integer, optional
the minimum total number of fermions allowed on-site.
default to 0
spin : integer, float, optional
the spin of the fermions.
default to 0
nFlavors : integer, optional
the number of "flavors" of fermions-basically the number
of non-spin internal degrees of freedom.
default to 1
**Details**
+--------------+-------------------------------------------------------+
| Tag operator | Operator |
+==============+=======================================================+
| I | :math:`\hat{I}` |
+--------------+-------------------------------------------------------+
| nftotal | :math:`\hat{n}` |
+--------------+-------------------------------------------------------+
| FermiPhase | :math:`(-1)^{\hat{n}}` |
+--------------+-------------------------------------------------------+
| f | :math:`\hat{f}` |
+--------------+-------------------------------------------------------+
| fdagger | :math:`\hat{f}^{\dagger}` |
+--------------+-------------------------------------------------------+
| fP | :math:`\hat{f} (-1)^{\hat{n}}` |
| | (only if linearly independent of f) |
+--------------+-------------------------------------------------------+
| fdaggerP | :math:`\hat{f}^{\dagger} (-1)^{\hat{n}}` |
| | (only if linearly independent of fdagger) |
+--------------+-------------------------------------------------------+
"""
global localDimension
if int(2*spin)==2*spin:
spinSize=int(2*spin+1)
else:
raise Exception("2*spin must be an integer in BuildFermiOperators!")
ninternal=spinSize*nFlavors
if ninternal==1:
if Nmax>1:
raise Exception("Spinless fermions cannot have Nmax greater than 1!")
Operators=BuildScalarOperators(Nmax,Nmin,True)
localDimension=Nmax-Nmin+1
else:
if Nmax>ninternal:
print("Warning: Nmax set higher than nFlavors*spinSize in BuildFermiOperators!")
print("Nmax set to nFlavors*spinSize!")
if spinSize==1:
Operators=BuildFlavorfulFermiOperators(Nmax,Nmin,nFlavors)
localDimension=0
for n in range(Nmin,Nmax+1):
localDimension=localDimension+mult(n,nFlavors,1)
else:
Operators=BuildSpinfulFermiOperators(Nmax,Nmin,spin,nFlavors)
localDimension=0
for n in range(Nmin,Nmax+1):
localDimension=localDimension+mult(n,spinSize*nFlavors,1)
if(isinstance(Operators, OperatorList)):
# Already converted
return Operators
else:
# Still convert to OperatorList class
return OperatorList(Operators)
[docs]def FermidestructionFromList(stateList,i):
"""
Return a fermionic destruction operator from a stateList. Intended for
internal use only.
**Arguments**
stateList : tba
i : tba
"""
ls=stateList.shape[0]
mat=np.zeros((ls,ls))
for j in range(ls):
newstate=deepcopy(stateList[j,:])
newstate[i]=newstate[i]-1
for k in range(ls):
if all(newstate==stateList[k,:]):
fphase=0
for l in range(i):
fphase=fphase+stateList[j,l]
mat[k,j]=((-1.0)**fphase)*sqrt(stateList[j,i])
return mat
[docs]def BuildFlavorfulFermiOperators(Nmax, Nmin, nFlavors):
"""
Build Fermi operators with flavor degrees of freedom. Intended for
internal use only.
**Arguments**
Nmax : int
Maximum number of particles.
Nmin : int
Minimum number of particles.
nFlavors : int
Number of different types of particles.
"""
#build and index Hilbert space
localSize,stateList=flavorfulHilbertspace(Nmin,Nmax,nFlavors,maxFilling=1)
b=[]
bdagger=[]
n=[]
identity=np.eye(localSize)
ntotal=np.zeros((localSize,localSize))
Operators={}
for i in range(nFlavors):
Operators['f_'+str(i)]=FermidestructionFromList(stateList,i)
Operators['fdagger_'+str(i)]=Operators['f_'+str(i)].T
Operators['nf_'+str(i)]=np.dot(Operators['fdagger_'+str(i)],Operators['f_'+str(i)])
Operators['nf_'+str(i)][i,i]=Nmin #correct for shifted vacuum
ntotal=ntotal+Operators['nf_'+str(i)]
phaseOp=np.zeros((localSize,localSize))
for i in range(localSize):
phaseOp[i,i]=(-1.0)**(ntotal[i,i])
Operators['I']=identity
Operators['nftotal']=ntotal
Operators['FermiPhase']=phaseOp
Operators = OperatorList(Operators)
for i in range(nFlavors):
Phasedf=np.dot(Operators['f_'+str(i)],phaseOp)
[fp, phasedw] = Operators.check(Phasedf, Add=False)
if fp in Operators and phasedw==1.0:
pass
else:
Operators['f_'+str(i)+'P']=np.dot(Operators['f_'+str(i)],phaseOp)
Phasedf=np.dot(Operators['fdagger_'+str(i)],phaseOp)
[fp, phasedw] = Operators.check(Phasedf, Add=False)
if fp in Operators and phasedw==1.0:
pass
else:
Operators['fdagger_'+str(i)+'P']=np.dot(Operators['fdagger_'+str(i)],phaseOp)
return Operators
[docs]def BuildSpinfulFermiOperators(Nmax, Nmin, spin, nFlavors):
"""
Build Fermi operators with spin and flavor degrees of freedom. Intended
for internal use only.
**Arguments**
Nmax : int
Maximum number of particles.
Nmin : int
Minimum number of particles.
spin : float (half-integer)
Spin for bosons (there is only one spin for all flavors).
nFlavors : int
Number of different types of particles.
"""
SpinOperators=BuildSpinOperators(spin)
spinSize=int(2*spin+1)
idofSize=spinSize*nFlavors
FlavoredOperators=BuildFlavorfulFermiOperators(Nmax, Nmin, idofSize)
Operators = OperatorList({})
bdagger = []
b = []
for spinindex in range(spinSize):
for flavorindex in range(nFlavors):
if nFlavors>1:
Operators['f_'+str(flavorindex)+'_'+str(spinindex)]=FlavoredOperators['f_'+str(spinindex*nFlavors+flavorindex)]
b.append(Operators['f_'+str(flavorindex)+'_'+str(spinindex)])
Operators['fdagger_'+str(flavorindex)+'_'+str(spinindex)]=FlavoredOperators['fdagger_'+str(spinindex*nFlavors+flavorindex)]
bdagger.append(Operators['fdagger_'+str(flavorindex)+'_'+str(spinindex)])
Phasedf=np.dot(Operators['f_'+str(flavorindex)+'_'+str(spinindex)],FlavoredOperators['FermiPhase'])
[fp, phasedw] = Operators.check(Phasedf, Add=False)
if fp in Operators and phasedw==1.0:
pass
else:
Operators['f_'+str(flavorindex)+'_'+str(spinindex)+'P']=Phasedf
Phasedf=np.dot(Operators['fdagger_'+str(flavorindex)+'_'+str(spinindex)],FlavoredOperators['FermiPhase'])
[fp, phasedw] = Operators.check(Phasedf, Add=False)
if fp in Operators and phasedw==1.0:
pass
else:
Operators['fdagger_'+str(flavorindex)+'_'+str(spinindex)+'P']=Phasedf
Operators['nf_'+str(flavorindex)+'_'+str(spinindex)]=FlavoredOperators['nf_'+str(spinindex*nFlavors+flavorindex)]
else:
Operators['f_'+str(spinindex)]=FlavoredOperators['f_'+str(spinindex*nFlavors+flavorindex)]
b.append(Operators['f_'+str(spinindex)])
Operators['fdagger_'+str(spinindex)]=FlavoredOperators['fdagger_'+str(spinindex*nFlavors+flavorindex)]
bdagger.append(Operators['fdagger_'+str(spinindex)])
Operators['nf_'+str(spinindex)]=FlavoredOperators['nf_'+str(spinindex*nFlavors+flavorindex)]
Operators['nftotal']=FlavoredOperators['nftotal']
Operators['FermiPhase']=FlavoredOperators['FermiPhase']
Operators['I']=FlavoredOperators['I']
localSize=Operators['I'].shape[0]
sz=np.zeros((localSize,localSize))
splus=np.zeros((localSize,localSize))
sminus=np.zeros((localSize,localSize))
for i in range(spinSize):
for k in range(nFlavors):
for j in range(spinSize):
sz=sz+SpinOperators['sz'][i,j]*np.dot(bdagger[i*nFlavors+k],b[j*nFlavors+k])
splus=splus+SpinOperators['splus'][i,j]*np.dot(bdagger[i*nFlavors+k],b[j*nFlavors+k])
sminus=sminus+SpinOperators['sminus'][i,j]*np.dot(bdagger[i*nFlavors+k],b[j*nFlavors+k])
Operators['fsz']=sz
Operators['fsplus']=splus
Operators['fsminus']=sminus
return Operators
[docs]def BuildFermiBoseOperators(BNmax, FNmax, BNmin=0, FNmin=0, Bspin=0, Fspin=0,
BnFlavors=1, FnFlavors=1, maxFilling=None):
"""
Build both Bose and Fermi operators in the same Hilbert space.
ntotal is the sum of nbtotal and nftotal.
BNmax : int
Maximum number of bosons.
FNmax : int
Maximum number of fermions.
BNmin : int, optional
Minimum number of bosons.
Default to 0
FNmin : int, optional
Minimum number of fermions.
Default to 0
Bspin : int, optional
Spin of the bosons (eqaul for all flavors).
Default to 0
Fspin : int, optional
Spin of the fermions (equal for all flavors).
Default to 0
BnFlavors : int, optional
Number of different types of bosons.
Default to 0
FnFlavors : int, optional
Number of different types of fermions.
Default to 0
maxFilling : int or ``None``, optional
Local Hilbert space has at most maxFilling particles in any one
boson flavor component.
Default to ``None``.
"""
global localDimension
BoseOperators=BuildBoseOperators(BNmax, BNmin, Bspin, BnFlavors, maxFilling)
FermiOperators=BuildFermiOperators(FNmax, FNmin, Fspin, FnFlavors)
boseSize=BoseOperators['I'].shape[0]
fermiSize=FermiOperators['I'].shape[0]
BoseFermiOperators={}
BoseFermiOperators['I']=np.kron(BoseOperators['I'],FermiOperators['I'])
for key in BoseOperators:
BoseFermiOperators[key]=np.kron(BoseOperators[key],FermiOperators['I'])
for key in FermiOperators:
BoseFermiOperators[key]=np.kron(BoseOperators['I'],FermiOperators[key])
BoseFermiOperators['ntotal']=BoseFermiOperators['nbtotal']+BoseFermiOperators['nftotal']
localDimension=boseSize+fermiSize
return OperatorList(BoseFermiOperators)
[docs]def isprime(x):
"""
Test if an integer is prime using a simple sieve.
**Argument**
x : int
Number to be tested as prime number.
"""
for j in range(2,int(sqrt(x)+1)):
if (x%j)==0:
return False
return True
[docs]def firstnprimes(n):
"""
Return the first n primes.
**Arguments**
int : n
Get n prime numbers.
"""
y=[]
i=0
x=2
while i<n:
if isprime(x):
i+=1
y.append(x)
x+=1
return y
[docs]def SROPHash(q):
"""
Return the square root of primers hash (float).
**Arguments**
q : list, tuple
get a hash of this list or tuple.
**Details**
Given a vector of integers q, compute the square-root of primes hash
:math:`H(q) = \sum_{i} q_i \sqrt{p_i}`, where :math:`p_i` is the
:math:`i^{th}` prime.
"""
myprimes=firstnprimes(len(q))
y=0.0
for i in range(len(q)):
y+=q[i]*sqrt(myprimes[i])
return y