#!/usr/bin/python
#
# This file provides an interface from the pyCTQW.MPI python module
# to the libctqwMPI.so FORTRAN library contained in source code src/ctqwMPI.F90
# ----------------------------------------------------------------------------
# pyCTQW - Distributed memory CTQW Fortran library and Python module
# Copyright (C) 2013-2014, Joshua Izaac
#
# pyCTQW is free software: you can redistribute it and/or modify it under the
# terms of version 3 of the GNU General Public License as published by
# the Free Software Foundation.
#
# pyCTQW is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
# more details.
#
# You should have received a copy of the GNU General Public License
# along with pyCTQW. If not, see <http://www.gnu.org/licenses/>.
# ----------------------------------------------------------------------------
"""
Additional ctqw related classes.
"""
import os as _os
import errno as _errno
import sys as _sys
import time as _time
import glob as _glob
from petsc4py import PETSc as _PETSc
from libpyctqw_MPI import ctqwmpi as _ctqwmpi
import io as _io
import plot as _plot
import numpy as _np
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#------------------ Eigensolver object (PETSc matrix input) ------------------
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[docs]class EigSolver(object):
""" Contains methods for setting up and solving for the boundary eigenvalues of a PETSc matrix.
Args:
mat (petsc4py.PETSc.Mat) : input matrix
Keyword Args:
esolver (str): the eigensolver algorithm to use.
* ``'krylovschur'`` *(default)* - Krylov-Schur
* ``'arnoldi'`` - Arnoldi Method
* ``'lanczos'`` - Lanczos Method
* ``'power'`` - Power Iteration/Rayleigh Quotient Iteration
* ``'gd'`` - Generalized Davidson
* ``'jd'`` - Jacobi-Davidson,
* ``'lapack'`` - Uses LAPACK eigenvalue solver routines
* ``'arpack'`` - *only available if SLEPc is\
compiled with ARPACK linking*
workType (str): can be used to set the eigensolver worktype \
(either ``'ncv'`` or ``'mpd'``). The default \
is to let SLEPc decide.
workSize (int): sets the work size **if** :attr:`workType` is set.
tolIn (float): tolerance of the eigenvalue solver \
(*default* ``0.`` (SLEPc decides)).
maxIt (int): maximum number of iterations of the eigenvalue solver \
(*default* ``0`` (SLEPc decides)).
verbose (bool): if ``True``, writes eigensolver information to the console
emax_estimate (float): used to override the calculation \
of the graphs maximum eigenvalue.
emin_estimate (float): used to override the calculation \
of the graphs minimum eigenvalue. For a *finite* graph, \
``emin_estimate=0`` is set by default.
.. caution::
* If supplied, the value of :attr:`emax_estimate` (:math:`\hat{\lambda}_{\max}`) \
**must** satisfy :math:`\hat{\lambda}_{\max}\geq\lambda_{\max}`, \
where :math:`\lambda_{\max}` is the actual maximum eigenvalue of the graph.
* Similarly, the value of :attr:`emin_estimate` (:math:`\hat{\lambda}_{\min}`) \
**must** satisfy :math:`\hat{\lambda}_{\min}\leq\lambda_{\min}`, \
where :math:`\lambda_{\min}` is the actual minimum eigenvalue of the graph.
* The greater the difference value of :math:`|\hat{\lambda}_{\max} \
-\lambda_{\max}|` and :math:`|\hat{\lambda}_{\min} \
-\lambda_{\min}|`, the longer the convergence \
time of the ``chebyshev`` propagator when simulating a CTQW.
Note:
* For more information regarding these properties,refer to \
the `SLEPc documentation\
<http://www.grycap.upv.es/slepc/documentation/manual.htm>`_
"""
def __init__(self,mat,**kwargs):
self.__default = {'esolver' : 'krylovschur',
'workType': 'null',
'workSize': '35',
'tol' : 0.,
'maxIt' : 0,
'verbose' : False,
'emax_estimate' : None,
'emin_estimate' : None,
'Emin_val': None,
'Emax_val': None}
self.__mat = mat
self.__rank = _PETSc.Comm.Get_rank(_PETSc.COMM_WORLD)
for key,default in self.__default.iteritems():
setattr(self, key, kwargs.get(key,default))
[docs] def setEigSolver(self,**kwargs):
"""Set some or all of the eigenvalue solver properties.
Keyword Args:
esolver (str): the eigensolver algorithm to use.
* ``'krylovschur'`` *(default)* - Krylov-Schur
* ``'arnoldi'`` - Arnoldi Method
* ``'lanczos'`` - Lanczos Method
* ``'power'`` - Power Iteration/Rayleigh Quotient Iteration
* ``'gd'`` - Generalized Davidson
* ``'jd'`` - Jacobi-Davidson,
* ``'lapack'`` - Uses LAPACK eigenvalue solver routines
* ``'arpack'`` - *only available if SLEPc is\
compiled with ARPACK linking*
workType (str): can be used to set the eigensolver worktype
(either ``'ncv'`` or ``'mpd'``). The default
is to let SLEPc decide.
workSize (int): sets the work size **if** ``workType`` is set.
tolIn (float): tolerance of the eigenvalue solver
(*default* ``0.`` (SLEPc decides)).
maxIt (int): maximum number of iterations of the eigenvalue solver
(*default* ``0`` (SLEPc decides)).
verbose (bool): if ``True``, writes eigensolver information to the console
emax_estimate (float): used to override the calculation
of the graphs maximum eigenvalue.
.. caution::
* If supplied, the value of :attr:`emax_estimate`:math:`\hat{\lambda}_{\max}` \
**must** satisfy :math:`\hat{\lambda}_{\max}\geq\lambda_{\max}`, \
where :math:`\lambda_{\max}` is the actual maximum eigenvalue of the graph.
* The greater the value of :math:`\hat{\lambda}_{\max} \
-\lambda_{\max}`, the longer the convergence \
time of the ``chebyshev`` propagator when simulating a CTQW.
Note:
* For more information regarding these properties,refer to \
the `SLEPc documentation\
<http://www.grycap.upv.es/slepc/documentation/manual.htm>`_
"""
for key in kwargs:
if hasattr(self, key):
setattr(self, key, kwargs.get(key))
else:
'Property type does not exist!'
[docs] def getEigSolver(self,*args):
""" Get some or all of the GraphISO properties.
For the list of the properties see :func:`setEigSolver`.
"""
for key in args:
if hasattr(self, key):
return getattr(self, key)
else:
'Property type does not exist!'
[docs] def findEmax(self):
""" Returns the maximum Eigenvalue of the matrix, along with associated error.
Example:
>>> Emax, error = EigSolver.findEmax()
.. admonition:: Fortran interface
This function calls the Fortran function :f:func:`min_max_eigs`.
"""
# Calcaulate the eigenvalues
if self.emax_estimate is not None:
self.Emax_val = self.emax_estimate
self.Emax_err = 0.
else:
EmaxStage = _PETSc.Log.Stage('Emax'); EmaxStage.push()
Emax,Emax_error,ierr = _ctqwmpi.min_max_eigs(self.__mat.fortran,self.__rank,'max',
self.esolver,self.workType,self.workSize,self.tol,self.maxIt,self.verbose)
EmaxStage.pop()
if ierr==0:
self.Emax_val = Emax
self.Emax_err = Emax_error
[docs] def findEmin(self):
""" Returns the minimum Eigenvalue of the matrix, along with associated error.
Example:
>>> Emin, error = EigSolver.findEmin()
.. admonition:: Fortran interface
This function calls the Fortran function :f:func:`min_max_eigs`.
"""
# Calcaulate the eigenvalues
if self.emin_estimate is not None:
self.Emin_val = self.emin_estimate
self.Emin_err = 0.
else:
EminStage = _PETSc.Log.Stage('Emin'); EminStage.push()
Emin,Emin_error,ierr = _ctqwmpi.min_max_eigs(self.__mat.fortran,self.__rank,'min',
self.esolver,self.workType,self.workSize,self.tol,self.maxIt,self.verbose)
EminStage.pop()
if ierr==0:
self.Emin_val = Emin
self.Emin_err = Emin_error
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#------------------ Hamiltonian object (grid size input) -----------------------
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[docs]class Hamiltonian(object):
""" Contains methods for initializing, creating and manipulating Hamiltonian matrices.
Args:
N (int) : number of nodes to initialize the Hamiltonian object with.
Returns:
: creates an unallocated PETSc matrix on communicator \
:attr:`petsc4py.PETSc.COMM_WORLD`. This is accessed \
via the attribute :attr:`Hamiltonian.mat`.
:rtype: petsc4py.PETSc.Mat
Returns:
: an associated :class:`EigSolver` object is also created \
at :attr:`Hamiltonian.EigSolver`.
:rtype: :class:`EigSolver`
"""
def __init__(self,N):
self.rank = _PETSc.Comm.Get_rank(_PETSc.COMM_WORLD)
self.N = N
# define matrices
self.mat = _PETSc.Mat()
self.mat.create(_PETSc.COMM_WORLD)
# create eigenvalue solver
self.EigSolver = EigSolver(self.mat)
[docs] def reinitialize(self):
""" Destroy and reinitialize the PETSc matrix :attr:`Hamiltonian.mat`."""
self.destroy()
self.mat = _PETSc.Mat()
self.mat.create(_PETSc.COMM_WORLD)
[docs] def importAdj(self,filename,filetype,d=[0],amp=[0.],layout='spring',delimiter=None):
""" Create a Hamiltonian from an imported adjacency matrix.
Args:
filename (str): path to the file containing the adjacency matrix of the graph
filetype (str): the filetype of the imported adjacency matrix.
* ``'txt'`` - an :math:`N\\times N` dense 2D array in text format.
* ``'bin'`` - an :math:`N\\times N` PETSc binary matrix.
d (array of ints): an array containing *integers* indicating the nodes
where diagonal defects are to be placed (e.g. ``d=[0,1,4]``).
amp (array of floats): an array containing *floats* indicating the diagonal defect
amplitudes corresponding to each element in ``d`` (e.g. ``amp=[0.5,-1,4.2]``).
layout (str): the format to store the position of the nodes.
* ``spring`` *(default)* - spring layout.
* ``circle`` - nodes are arranged in a circle.
* ``spectral`` - nodes are laid out according to the \
spectrum of the graph.
* ``random`` - nodes are arranged in a random pattern.
delimiter (str): this is passed to `numpy.genfromtxt\
<http://docs.scipy.org/doc/numpy/reference/generated/numpy.genfromtxt.html>`_
in the case of strange delimiters in an imported ``txt`` file.
Returns:
:this allocates and builds a Hamiltonian matrix at :attr:`Hamiltonian.mat`.
:rtype: :func:`petsc4py.PETSc.Mat`
Returns:
: the original adjacency matrix is also stored as a PETSc matrix at :attr:`Hamiltonian.Adj`.
:rtype: :func:`petsc4py.PETSc.Mat`
.. important::
* The number of nodes in the imported adjacency matrix \
**must** match the number of nodes the Hamiltonian object \
is initialized with.
* The length of ``amp`` and ``d`` must be identical.
.. note::
* This function calls Python functions located at :mod:`pyCTQW.MPI.io`, \
and is somewhat slower than :func:`importAdjToH`, but with the advantage of \
a more flexible text file importer.
* For the ability to generate 2 or 3 particle Hamiltonians with or without \
interactions, please see :func:`importAdjToH`.
"""
try:
if self.mat.isAssembled():
self.reinitialize()
except: pass
# create the Hamiltonian
Hamiltonian = _PETSc.Log.Stage('Hamiltonian')
Hamiltonian.push()
self.Adj = _io.loadMat(filename,filetype,delimiter=delimiter)
self.mat = _io.adjToH(self.Adj,d=d,amp=amp)
Hamiltonian.pop()
self.nodePos, self.lineX, self.lineY = _plot.getGraphNodes(self.Adj,layout=layout)
[docs] def importAdjToH(self,filename,filetype,d=[0.],amp=[0.],p='1',layout='spring',delimiter=None,interaction=0.,bosonic=False):
""" Create a Hamiltonian from an imported adjacency matrix.
Args:
filename (str): path to the file containing the adjacency matrix of the graph
filetype (str): the filetype of the imported adjacency matrix.
* ``'txt'`` - an :math:`N\\times N` dense 2D array in text format.
* ``'bin'`` - an :math:`N\\times N` PETSc binary matrix.
d (array of ints): an array containing *integers* indicating the nodes
where diagonal defects are to be placed (e.g. ``d=[0,1,4]``).
amp (array of floats): an array containing *floats* indicating the diagonal defect
amplitudes corresponding to each element in ``d`` (e.g. ``amp=[0.5,-1,4.2]``).
p (str) : (``'1'``|``'2'``|``'3'``) - specify whether to create a 1, 2 or 3 \
particle Hamiltonian.
layout (str): the format to store the position of the nodes.
* ``spring`` *(default)* - spring layout.
* ``circle`` - nodes are arranged in a circle.
* ``spectral`` - nodes are laid out according to the \
spectrum of the graph.
* ``random`` - nodes are arranged in a random pattern.
delimiter (str): this is passed to `numpy.genfromtxt \
<http://docs.scipy.org/doc/numpy/reference/generated/numpy.genfromtxt.html>`_
in the case of strange delimiters in an imported ``txt`` file.
interaction (float): the amplitude of interaction between the walkers \
when located on the same vertex.
bosonic (bool): if true, a Bosonic Hamiltonian is generated.
Returns:
: * :attr:`Hamiltonian.mat` *(petsc4py.PETSc.Mat)* - the Hamiltonian matrix.
* :attr:`Hamiltonian.Adj` *(petsc4py.PETSc.Mat)* - the adjacency matrix stored as a PETSc matrix.
* :attr:`Hamiltonian.nodePos` *(array)* - the :math:`(x,y)` coordinates of the graph vertices.
* :attr:`Hamiltonian.lineX` *(array)* - the :math:`x` coordinates of the edges connecting graph vertices.
* :attr:`Hamiltonian.lineY` *(array)* - the :math:`y` coordinates of the edges connecting graph vertices.
.. important::
* The number of nodes in the imported adjacency matrix \
**must** match the number of nodes the Hamiltonian object \
is initialized with.
* The length of ``amp`` and ``d`` must be identical.
.. admonition:: Fortran interface
This function calls the Fortran function :f:func:`importadjtoh` and is \
faster than :func:`importAdjToH`, but with a 'pickier' text file importer.
"""
try:
if self.mat.isAssembled():
self.reinitialize()
except: pass
# create the Hamiltonian
Hamiltonian = _PETSc.Log.Stage('Hamiltonian')
Hamiltonian.push()
_ctqwmpi.importadjtoh(self.mat.fortran,filename,p,d,amp,interaction,bosonic)
Hamiltonian.pop()
# create the adjacency matrix
self.Adj = _io.loadMat(filename,filetype,delimiter=delimiter)
self.nodePos, self.lineX, self.lineY = _plot.getGraphNodes(self.Adj,layout=layout)
[docs] def createLine2P(self,d=[0],amp=[0.],interaction=0.):
""" Create a 2 particle infinite line Hamiltonian.
Args:
d (array of ints): an array containing *integers* indicating the nodes
where diagonal defects are to be placed (e.g. ``d=[0,1,4]``).
amp (array of floats): an array containing *floats* indicating the diagonal defect
amplitudes corresponding to each element in ``d`` (e.g. ``amp=[0.5,-1,4.2]``).
interaction (float): the amplitude of interaction between the walkers \
when located on the same vertex.
Returns:
: :attr:`Hamiltonian.mat`, the Hamiltonian matrix.
:rtype: :func:`petsc4py.PETSc.Mat`
.. important::
* The length of ``amp`` and ``d`` must be identical.
.. admonition:: Fortran interface
This function calls the Fortran function :f:func:`hamiltonian_p2_line`.
"""
try:
if self.mat.isAssembled():
self.reinitialize()
except: pass
self.defectNodes = d
self.defectAmp = amp
# create the Hamiltonian
Hamiltonian = _PETSc.Log.Stage('Hamiltonian')
Hamiltonian.push()
_ctqwmpi.hamiltonian_p2_line(self.mat.fortran,self.defectNodes,self.defectAmp,interaction,self.N)
Hamiltonian.pop()
[docs] def createLine3P(self,d=[0],amp=[0.],interaction=0.):
""" Create a 3 particle infinite line Hamiltonian.
Args:
d (array of ints): an array containing *integers* indicating the nodes
where diagonal defects are to be placed (e.g. ``d=[0,1,4]``).
amp (array of floats): an array containing *floats* indicating the diagonal defect
amplitudes corresponding to each element in ``d`` (e.g. ``amp=[0.5,-1,4.2]``).
interaction (float): the amplitude of interaction between the walkers \
when located on the same vertex.
Returns:
: :attr:`Hamiltonian.mat`, the Hamiltonian matrix.
:rtype: :func:`petsc4py.PETSc.Mat`
.. important::
* The length of ``amp`` and ``d`` must be identical.
.. admonition:: Fortran interface
This function calls the Fortran function :f:func:`hamiltonian_p3_line`.
"""
try:
if self.mat.isAssembled():
self.reinitialize()
except: pass
self.defectNodes = d
self.defectAmp = amp
# create the Hamiltonian
Hamiltonian = _PETSc.Log.Stage('Hamiltonian')
Hamiltonian.push()
_ctqwmpi.hamiltonian_p3_line(self.mat.fortran,self.defectNodes,self.defectAmp,interaction,self.N)
Hamiltonian.pop()
[docs] def createLine(self,d=[0],amp=[0.]):
""" Create a 1 particle infinite line Hamiltonian.
Args:
d (array of ints): an array containing *integers* indicating the nodes
where diagonal defects are to be placed (e.g. ``d=[0,1,4]``).
amp (array of floats): an array containing *floats* indicating the diagonal defect
amplitudes corresponding to each element in ``d`` (e.g. ``amp=[0.5,-1,4.2]``).
interaction (float): the amplitude of interaction between the walkers \
when located on the same vertex.
Returns:
: :attr:`Hamiltonian.mat`, the Hamiltonian matrix.
:rtype: :func:`petsc4py.PETSc.Mat`
.. important::
* The length of ``amp`` and ``d`` must be identical.
.. admonition:: Fortran interface
This function calls the Fortran function :f:func:`hamiltonian_p1_line`.
"""
try:
if self.mat.isAssembled():
self.reinitialize()
except: pass
self.defectNodes = d
self.defectAmp = amp
# create the Hamiltonian
Hamiltonian = _PETSc.Log.Stage('Hamiltonian')
Hamiltonian.push()
_ctqwmpi.hamiltonian_p1_line(self.mat.fortran,self.defectNodes,self.defectAmp,self.N)
Hamiltonian.pop()
[docs] def Emax(self,**kwargs):
""" Returns the maximum Eigenvalue of the matrix.
Example:
>>> Emax = Hamiltonian.Emax()
"""
if self.EigSolver.Emax_val is None:
self.EigSolver.findEmax()
return self.EigSolver.Emax_val
[docs] def Emin(self,**kwargs):
""" Returns the minimum Eigenvalue of the matrix.
Example:
>>> Emin = Hamiltonian.Emin()
"""
if self.EigSolver.Emin_val is None:
self.EigSolver.findEmin()
return self.EigSolver.Emin_val
[docs] def destroy(self):
""" Destroys all associated PETSc matrices."""
self.mat.destroy()
try:
self.Adj.destroy()
except:
pass
[docs]class nodeHandle(object):
""" Creates a handle with the ability to store and updates node probability.
Args:
nodes (array of ints): the nodes to watch (e.g. ``[0,1,4]``).
t (float): the elapsed time at handle initialization.
psi (petsc4py.PETSc.Vec) : vector containing the global probability for particle 1.
psi2 (None or petsc4py.PETSc.Vec) : vector containing the global probability for particle 2.
psi3 (None or petsc4py.PETSc.Vec) : vector containing the global probability for particle 3.
Warning:
Note that the handle attributes/methods are **not** collective;
if running on multiple nodes, only *local* values will be
stored/updated on each node.
"""
def __init__(self,nodes,t,psi,psi2=None,psi3=None):
self.rank = _PETSc.COMM_WORLD.Get_rank()
self.time = [t]
self.all_nodes = nodes
self.local_nodes = []
(Istart,Iend) = psi.getOwnershipRange()
for i in self.all_nodes:
if Istart <= i < Iend:
self.local_nodes.append(i)
self.local_prob = psi.getValues(self.local_nodes)
if psi2 is not None:
self.local_prob2 = psi2.getValues(self.local_nodes)
if psi3 is not None:
self.local_prob3 = psi3.getValues(self.local_nodes)
[docs] def update(self,t,psi,psi2=None,psi3=None):
""" Update the handle with node probability for an additional timestep.
Args:
t (float): the elapsed time at handle update.
psi (petsc4py.PETSc.Vec) : vector containing the global probability for particle 1.
psi2 (None or petsc4py.PETSc.Vec) : vector containing the global probability for particle 2.
psi3 (None or petsc4py.PETSc.Vec) : vector containing the global probability for particle 3.
"""
self.time.append(t)
prob_update = psi.getValues(self.local_nodes)
self.local_prob = _np.vstack([self.local_prob,prob_update])
if psi2 is not None:
prob_update = psi2.getValues(self.local_nodes)
self.local_prob2 = _np.vstack([self.local_prob2,prob_update])
if psi3 is not None:
prob_update = psi3.getValues(self.local_nodes)
self.local_prob3 = _np.vstack([self.local_prob3,prob_update])
[docs] def getLocalNodes(self,t=None,p=1):
""" Get thet stored node probability and timestep array for particle :attr:`p` over all watched nodes.
Keyword Args:
t (None or float): the elapsed time at which to return the node probabilities. \
If ``None``, **all** timesteps will be returned.
p (int) : (``1``,``2``,``3``) - the particle to return probability data for.
Returns:
tuple of arrays : returns arrays corresponding to the times recorded (if ``t=None``) \
and the probabilities over watched nodes.
Example:
>>> timeArray, probXArray, probYArray = nodeHandle.getLocalNodes()
or
>>> probXArray, probYArray = nodeHandle.getLocalNodes(t=4)
"""
if t is not None:
try:
indt = self._time.index(t)
except ValueError:
if self.rank == 0: print '\nERROR: time {} was not handled'.format(t)
return
if p == 1:
return self.local_prob[indt]
elif p == 2:
return self.local_prob[indt], self.local_prob2[indt]
elif p == 3:
return self.local_prob[indt], self.local_prob2[indt], self.local_prob3[indt]
else:
if p == 1:
return _np.array(self.time), self.local_prob.T.tolist()
elif p == 2:
return _np.array(self.time), self.local_prob.T.tolist(), self.local_prob2.T.tolist()
elif p == 3:
return _np.array(self.time), self.local_prob.T.tolist(), self.local_prob2.T.tolist(), self.local_prob3.T.tolist()
[docs] def getLocalNode(self,node,p=1):
""" Get thet stored probability and timestep array for particle :attr:`p` at a specified node.
Args:
node (int) : the node to retrieve data from.
Keyword Args:
p (int) : (``1``,``2``,``3``) - the particle to return probability data for.
Returns:
tuple of arrays : returns arrays corresponding to the times recorded (if ``t=None``) \
and the probabilities over the specified node.
Example:
>>> timeArray, probXArray, probYArray = nodeHandle.getLocalNodes(3)
"""
try:
indN = self.local_nodes.index(node)
except ValueError:
if p == 1:
return _np.array(self.time), []
elif p == 2:
return _np.array(self.time), [], []
elif p == 3:
return _np.array(self.time), [], [], []
if p == 1:
return _np.array(self.time), self.local_prob.T.tolist()[indN]
elif p == 2:
return _np.array(self.time), self.local_prob.T.tolist()[indN], self.local_prob2.T.tolist()[indN]
elif p == 3:
return _np.array(self.time), self.local_prob.T.tolist()[indN], self.local_prob2.T.tolist()[indN], self.local_prob3.T.tolist()[indN]
[docs]class entanglementHandle(object):
""" Creates a handle with the ability to store and updates entanglement.
Args:
t (float): the elapsed time at handle initialization.
psi (petsc4py.PETSc.Vec) : vector containing the statespace.
N (int) : the size of the quantum walker's position space (i.e. number \
of graph vertices).
Keyword Args:
: EigSolver keywords can also be passed;
for more details of the available EigSolver properties,
see :class:`EigSolver`.
.. important::
The Eigsolver *must* be set to lapack, since **all**
the eigenvalues must be found in the calculation of entanglement.
.. admonition:: Fortran interface
This function calls the Fortran function :f:func:`entanglement`.
"""
def __init__(self,t,psi,N,**kwargs):
self.__default = {'esolver' : 'lapack',
'workType': 'null',
'workSize': '35',
'tol' : 0.,
'maxIt' : 0,
'verbose' : False,
'p' : 2}
self._rank = _PETSc.COMM_WORLD.Get_rank()
self.time = [t]
self.N = N
for key,default in self.__default.iteritems():
setattr(self, key, kwargs.get(key,default))
entanglementS = _PETSc.Log.Stage('Entanglement')
entanglementS.push()
if self.p==2:
entInit, ierr = _ctqwmpi.entanglement(psi.fortran,self.N,
self.esolver,self.workType,self.workSize,self.tol,self.maxIt,self.verbose)
entanglementS.pop()
self.entanglement = [entInit]
[docs] def update(self,t,psi):
""" Creates a handle with the ability to store and updates entanglement.
Args:
t (float): the timestep of the handle update.
psi (petsc4py.PETSc.Vec) : vector containing the statespace.
.. admonition:: Fortran interface
This function calls the Fortran function :f:func:`entanglement`.
"""
self.time.append(t)
entanglementS = _PETSc.Log.Stage('Entanglement')
entanglementS.push()
if self.p==2:
entUpdate, ierr = _ctqwmpi.entanglement(psi.fortran,self.N,
self.esolver,self.workType,self.workSize,self.tol,self.maxIt,self.verbose)
entanglementS.pop()
self.entanglement.append(entUpdate)
[docs] def getEntanglement(self,t=None):
""" Returns the stored entanglement.
Keyword Args:
t (None or float): the timestep to retrieve the entanglement. \
If ``None``, the the entanglement at **all** \
recorded timesteps is returned.
Returns:
tuple of arrays : returns arrays corresponding to the times recorded (if ``t=None``) \
and the entanglement over the recorded times.
Example:
>>> timeArray, entArray = entanglementHandle.getEntanglement()
or
>>> entArray = entanglementHandle.getEntanglement(t=4)
"""
if t is not None:
try:
indt = self._time.index(t)
except ValueError:
if self._rank == 0: print '\nERROR: time {} was not handled'.format(t)
return
return self.entanglement[indt]
else:
return _np.array(self.time), _np.array(self.entanglement)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#--------------------------- 1 particle CTQW ---------------------------------
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[docs]class QuantumWalkP1(object):
""" Methods providing a very basic framework of a 1 particle CTQW.
.. caution::
This object constructs some of the *framework* of a 1 particle
quantum walk, but is not feature complete. For a feature complete
construction, see :func:`pyCTQW.MPI.Graph` or :func:`pyCTQW.MPI.Line`
for which this is a base.
"""
def __init__(self,N):
self.rank = _PETSc.Comm.Get_rank(_PETSc.COMM_WORLD)
self.N = N
self.t = 0
self._timestep = False
# define vectors
#: Initial state vector
self.psi0 = _PETSc.Vec()
self.psi0.create(_PETSc.COMM_WORLD)
self.psi0.setSizes(self.N)
self.psi0.setUp()
self.psi = self.psi0.duplicate()
self.prob = self.psi0.duplicate()
# define matrices
self.H = Hamiltonian(self.N)
self.EigSolver = self.H.EigSolver
[docs] def importInitState(self,filename,filetype):
""" Imports the initial state of the quantum walk from a file.
Args:
filename (str): path to the file containing the input state.
filetype (str): the filetype of the imported adjacency matrix.
* ``'txt'`` - an :math:`N` element column vector in text format.
* ``'bin'`` - an :math:`N` element PETSc binary vector.
Returns:
: this creates a PETSc vector containing the initial state, \
accessed via :attr:`Graph.psi0`.
:rtype: :func:`petsc4py.PETSc.Vec`
Warning:
The number of elements in the imported input state vector
**must** match the number of nodes the 1P CTQW object
is initialized with.
"""
self.initState = 'file:'+filename
# create the inital stage
initStateS = _PETSc.Log.Stage('initState')
initStateS.push()
try:
self.psi0 = _io.loadVec(filename,filetype)
self.marginal(self.psi0)
except:
print '\nERROR: incorrect state (is it the correct length?'
_sys.exit()
initStateS.pop()
self.marginal(self.psi0)
[docs] def marginal(self,vec):
# calculate marginal probabilities
Marginal = _PETSc.Log.Stage('Marginal'); Marginal.push()
_ctqwmpi.marginal1(vec.fortran,self.prob.fortran,self.N)
Marginal.pop()
[docs] def watch(self,nodes):
""" Creates a handle that watches node probability during propagation.
Args:
nodes (array of ints): the nodes to watch (e.g. ``[0,1,4]``).
Returns:
: creates a handle that can be
accessed to retrieve node probabilities for various :math:`t`
:rtype: :func:`pyCTQW.MPI.ctqw.nodeHandle`
Example:
>>> walk.watch([0,1,2,3,4])
>>> walk.propagate(5.,method='chebyshev')
>>> timeArray, probArray = walk.handle.getLocalNodes()
Warning:
Note that `walk.handle` attributes are **not** collective;
if running on multiple nodes, only *local* values will be
returned.
"""
self.handle = nodeHandle(nodes,self.t,self.prob)
[docs] def propagate(self,t,method='chebyshev',**kwargs):
""" Propagates the quantum walk for time t.
Args:
t (float): the timestep over which to propagate the state ``walk.psi0`` via the relationship
:math:`\\left|\\psi\\right\\rangle = e^{-iHt}\\left|\\psi0\\right\\rangle`.
method (str): the propagation algorithm to use
(``'chebyshev'`` *(default)* or ``'krylov'``).
Keyword Args:
esolver (str): the eigensolver algorithm to use.
* ``'krylovschur'`` *(default)* - Krylov-Schur
* ``'arnoldi'`` - Arnoldi Method
* ``'lanczos'`` - Lanczos Method
* ``'power'`` - Power Iteration/Rayleigh Quotient Iteration
* ``'gd'`` - Generalized Davidson
* ``'jd'`` - Jacobi-Davidson,
* ``'lapack'`` - Uses LAPACK eigenvalue solver routines
* ``'arpack'`` - *only available if SLEPc is\
compiled with ARPACK linking*
workType (str): can be used to set the eigensolver worktype \
(either ``'ncv'`` or ``'mpd'``). The default \
is to let SLEPc decide.
workSize (int): sets the work size **if** :attr:`workType` is set.
tolIn (float): tolerance of the eigenvalue solver \
(*default* ``0.`` (SLEPc decides)).
maxIt (int): maximum number of iterations of the eigenvalue solver \
(*default* ``0`` (SLEPc decides)).
verbose (bool): if ``True``, writes eigensolver information to the console
emax_estimate (float): used to override the calculation \
of the graphs maximum eigenvalue.
emin_estimate (float): used to override the calculation \
of the graphs minimum eigenvalue. For a *finite* graph, \
``emin_estimate=0`` is set by default.
.. caution::
* If supplied, the value of :attr:`emax_estimate` (:math:`\hat{\lambda}_{\max}`) \
**must** satisfy :math:`\hat{\lambda}_{\max}\geq\lambda_{\max}`, \
where :math:`\lambda_{\max}` is the actual maximum eigenvalue of the graph.
* Similarly, the value of :attr:`emin_estimate` (:math:`\hat{\lambda}_{\min}`) \
**must** satisfy :math:`\hat{\lambda}_{\min}\leq\lambda_{\min}`, \
where :math:`\lambda_{\min}` is the actual minimum eigenvalue of the graph.
* The greater the difference value of :math:`|\hat{\lambda}_{\max} \
-\lambda_{\max}|` and :math:`|\hat{\lambda}_{\min} \
-\lambda_{\min}|`, the longer the convergence \
time of the ``chebyshev`` propagator.
Note:
* The keyword arguments properties only apply if ``propagator='chebyshev'``
* For more information regarding these properties,refer to \
the `SLEPc documentation\
<http://www.grycap.upv.es/slepc/documentation/manual.htm>`_
Returns:
: this creates a PETSc vector containing the
propagated state, accessed via the attribute
:attr:`psi`, as well as a PETSc vector containing
the marginal probabilities, :attr:`prob`.
:rtype: :func:`petsc4py.PETSc.Vec`
Note:
The EigSolver properties only take effect if
`method='chebyshev'`.
.. admonition:: Fortran interface
This function calls the Fortran function :f:func:`expm`. for the Krylov
method, and the function :f:func:`qw_cheby` for the Chebyshev method.
"""
if self._timestep:
self.t = t + self.t
else:
self.t = t
self.EigSolver.setEigSolver(**kwargs)
if method=='krylov':
# SLEPc matrix exponential
expmS = _PETSc.Log.Stage('SLEPc expm'); expmS.push()
_ctqwmpi.qw_krylov(self.H.mat.fortran,t,self.psi0.fortran,self.psi.fortran)
expmS.pop()
elif method=='chebyshev':
# Chebyshev algorithm
chebyS = _PETSc.Log.Stage('Chebyshev'); chebyS.push()
_ctqwmpi.qw_cheby(self.psi0.fortran,self.psi.fortran,t,self.H.mat.fortran,
self.H.Emin(),self.H.Emax(),self.rank,self.N)
chebyS.pop()
self.marginal(self.psi)
try:
self.handle.update(self.t,self.prob)
except:
pass
self._timestep = False
[docs] def plotNodes(self,filename,t=None):
""" Creates a plot of the node probablities over time.
Args:
filename (str): the absolute/relative path to the desired output file.
.. important::
:func:`watch` **must** be called prior to propagation in order for\
the probabilities at the specified nodes to be stored.
Note:
* if multiple nodes are watched, they will **all** be plotted.
* ensure a file extension is present so that filetype is correctly set\
(choose one of png, pdf, ps, eps or svg).
* Timesteps are recorded when propagations occur - to increase the number\
of timesteps, increase the number of propagations.
"""
comm = _PETSc.COMM_WORLD
rank = comm.Get_rank()
node = self.handle.local_nodes
timeArray, probArray = self.handle.getLocalNodes(t=t)
probArray = comm.tompi4py().gather(probArray)
node = comm.tompi4py().gather(node)
if rank == 0:
timeArray = _np.array(timeArray)
nodeArray = _np.array([item for sublist in node for item in sublist])
probArray = _np.array([item for sublist in probArray for item in sublist]).real
_plot.plotNodes(timeArray,nodeArray,probArray,filename)
[docs] def exportState(self,filename,filetype):
""" Exports the final state of the quantum walk to a file.
Args:
filename (str): path to desired output file.
filetype (str): the filetype of the exported adjacency matrix.
* ``'txt'`` - an :math:`N` element column vector in text format.
* ``'bin'`` - an :math:`N` element PETSc binary vector.
"""
_io.exportVec(self.psi,filename,filetype)
[docs] def psiToInit(self):
""" Copies the state :attr:`self.psi` to :attr:`self.psi0`.
Example:
This can be used to iterate the propagation
over discrete timesteps: ::
for i in range(10):
walk.propagate(0.01,method='chebyshev')
walk.psiToInit()
"""
self.psi0 = self.psi
self._timestep = True
[docs] def destroy(self):
""" Destroys the 1 particle quantum walk object, and all \
associated PETSc matrices/vectors.
"""
self.H.destroy()
self.psi.destroy()
self.psi0.destroy()
self.prob.destroy()
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#--------------------------- 2 particle CTQW ---------------------------------
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[docs]class QuantumWalkP2(object):
""" Methods providing a very basic framework of a 2 particle CTQW.
.. caution::
This object constructs some of the *framework* of a 2 particle
quantum walk, but is not feature complete. For a feature complete
construction, see :func:`pyCTQW.MPI.Graph2P` or :func:`pyCTQW.MPI.Line2P`
for which this is a base.
"""
def __init__(self,N):
self.rank = _PETSc.Comm.Get_rank(_PETSc.COMM_WORLD)
self.N = N
self.t = 0
self.timestep = False
self._watchProb = False
self._watchEnt = False
# define vectors
self.psi0 = _PETSc.Vec()
self.psi0.create(_PETSc.COMM_WORLD)
self.psi0.setSizes(self.N**2)
self.psi0.setUp()
self.psi = self.psi0.duplicate()
# create additional vectors
self.psiX = _PETSc.Vec()
self.psiX.create(_PETSc.COMM_WORLD)
self.psiX.setSizes(self.N)
self.psiX.setUp()
self.psiY = self.psiX.duplicate()
# define matrices
self.H = Hamiltonian(self.N)
self.EigSolver = self.H.EigSolver
[docs] def importInitState(self,filename,filetype):
""" Imports the initial state of the quantum walk from a file.
Args:
filename (str): path to the file containing the input state.
filetype (str): the filetype of the imported adjacency matrix.
* ``'txt'`` - an :math:`N\\times N` *2D array* in text format.
* ``'bin'`` - an :math:`N^2` element PETSc binary *vector*.
Returns:
: this creates a PETSc vector containing the initial statespace, \
accessed via :attr:`self.psi0`.
:rtype: :func:`petsc4py.PETSc.Vec`
Warning:
The number of elements in the imported input statespace
**must** equal :math:`N^2` (either as a :math:`N\\times N`
text file or a :math:`N^2` PETSc binary vector),
where :math:`N` is the number of nodes the 2P CTQW object
is initialized with.
"""
self.initState = 'file:'+filename
# create the inital stage
initStateS = _PETSc.Log.Stage('initState')
initStateS.push()
try:
if filetype == 'txt':
self.psi0 = _io.loadMatToVec(filename,filetype)
elif filetype == 'bin':
self.psi0 = _io.loadVec(filename,filetype)
except:
print '\nERROR: incorrect state (is it the correct length?)'
_sys.exit()
initStateS.pop()
self.marginal(self.psi0)
[docs] def marginal(self,vec):
# calculate marginal probabilities
Marginal = _PETSc.Log.Stage('Marginal'); Marginal.push()
_ctqwmpi.marginal2(vec.fortran,self.psiX.fortran,'x',self.N)
_ctqwmpi.marginal2(vec.fortran,self.psiY.fortran,'y',self.N)
Marginal.pop()
[docs] def watch(self,nodes,watchtype='prob',**kwargs):
""" Creates a handle that watches either marginal probability \
or walker entanglement during propagation.
Args:
nodes (array of ints): the nodes to watch marginal probability (e.g. ``[0,1,4]``).
watchtype (str): (`'prob'` , `'entanglement'`).
the type of watch handle to produce.
Keyword Args:
: If ``watchtype='entanglement'``, EigSolver keywords can also be passed;
for more details of the available EigSolver properties,
see :func:`propagate`.
Returns:
: * if ``watchtype='prob'``, creates a handle that can be \
accessed to retrieve marginal node probabilities for various :math:`t`
* if ``watchtype='entanglment'``, creates a handle that can be \
accessed to retrieve entanglement values for various :math:`t`
:rtype: * if ``watchtype='prob'``: :func:`pyCTQW.MPI.ctqw.nodeHandle`
* if ``watchtype='entanglement'``: :func:`pyCTQW.MPI.ctqw.entanglementHandle`
Examples:
To watch the entanglement,
>>> walk.watch(None,watchtype='entanglement')
>>> walk.propagate(5.,method='chebyshev')
>>> timeArray, entArray = walk.entanglementHandle.getEntanglement()
.. note::
* The entanglement measure used in Von Neumann entropy, calculated \
via :math:`S=-\\sum_{i}\\lambda_i\log_{2}\\lambda_i`, where :math:`\lambda_i` \
are the eigenvalues of the reduced density matrix \
:math:`\\rho_2 = \\text{Tr}_1(|\\psi(t)\\rangle\\langle\\psi(t)|)`
* Nodes do not need to be specified, as entanglement is a global measurement.
* As it is a global measurement, there is a large amount of node communication \
which may increase overall program run time.
To watch the probabilities,
>>> walk.watch([0,1,4])
>>> walk.propagate(2.,method='chebyshev')
>>> timeArray, probXArray, probYArray = walk.handle.getLocalNode(4,p=2)
.. warning::
Note that `walk.handle` attributes are **not** collective;
if running on multiple nodes, only *local* values will be
returned.
"""
if watchtype == 'prob':
self._watchProb = True
self.handle = nodeHandle(nodes,self.t,self.psiX,psi2=self.psiY)
elif watchtype == 'entanglement':
self._watchEnt = True
self.entanglementHandle = entanglementHandle(self.t,self.psi0,self.N,**kwargs)
[docs] def propagate(self,t,method='chebyshev',**kwargs):
""" Propagates the quantum walk for time t.
Args:
t (float): the timestep over which to propagate the state ``walk.psi0`` via the relationship
:math:`\\left|\\psi\\right\\rangle = e^{-iHt}\\left|\\psi0\\right\\rangle`.
method (str): the propagation algorithm to use
(``'chebyshev'`` *(default)* or ``'krylov'``).
Keyword Args:
esolver (str): the eigensolver algorithm to use.
* ``'krylovschur'`` *(default)* - Krylov-Schur
* ``'arnoldi'`` - Arnoldi Method
* ``'lanczos'`` - Lanczos Method
* ``'power'`` - Power Iteration/Rayleigh Quotient Iteration
* ``'gd'`` - Generalized Davidson
* ``'jd'`` - Jacobi-Davidson,
* ``'lapack'`` - Uses LAPACK eigenvalue solver routines
* ``'arpack'`` - *only available if SLEPc is\
compiled with ARPACK linking*
workType (str): can be used to set the eigensolver worktype \
(either ``'ncv'`` or ``'mpd'``). The default \
is to let SLEPc decide.
workSize (int): sets the work size **if** :attr:`workType` is set.
tolIn (float): tolerance of the eigenvalue solver \
(*default* ``0.`` (SLEPc decides)).
maxIt (int): maximum number of iterations of the eigenvalue solver \
(*default* ``0`` (SLEPc decides)).
verbose (bool): if ``True``, writes eigensolver information to the console
emax_estimate (float): used to override the calculation \
of the graphs maximum eigenvalue.
emin_estimate (float): used to override the calculation \
of the graphs minimum eigenvalue. For a *finite* graph, \
``emin_estimate=0`` is set by default.
.. caution::
* If supplied, the value of :attr:`emax_estimate` (:math:`\hat{\lambda}_{\max}`) \
**must** satisfy :math:`\hat{\lambda}_{\max}\geq\lambda_{\max}`, \
where :math:`\lambda_{\max}` is the actual maximum eigenvalue of the graph.
* Similarly, the value of :attr:`emin_estimate` (:math:`\hat{\lambda}_{\min}`) \
**must** satisfy :math:`\hat{\lambda}_{\min}\leq\lambda_{\min}`, \
where :math:`\lambda_{\min}` is the actual minimum eigenvalue of the graph.
* The greater the difference value of :math:`|\hat{\lambda}_{\max} \
-\lambda_{\max}|` and :math:`|\hat{\lambda}_{\min} \
-\lambda_{\min}|`, the longer the convergence \
time of the ``chebyshev`` propagator.
Note:
* The keyword arguments properties only apply if ``propagator='chebyshev'``
* For more information regarding these properties,refer to \
the `SLEPc documentation\
<http://www.grycap.upv.es/slepc/documentation/manual.htm>`_
Returns:
: this creates a PETSc vector containing the
propagated state, accessed via the attribute
:attr:`psi`, as well as a PETSc vector containing
the marginal probabilities, :attr:`prob`.
:rtype: :func:`petsc4py.PETSc.Vec`
Note:
The EigSolver properties only take effect if
`method='chebyshev'`.
.. admonition:: Fortran interface
This function calls the Fortran function :f:func:`expm`. for the Krylov
method, and the function :f:func:`qw_cheby` for the Chebyshev method.
"""
if self.timestep:
self.t = t + self.t
else:
self.t = t
self.EigSolver.setEigSolver(**kwargs)
if method=='krylov':
# SLEPc matrix exponential
krylov = _PETSc.Log.Stage('SLEPc krylov'); krylov.push()
_ctqwmpi.qw_krylov(self.H.mat.fortran,t,self.psi0.fortran,self.psi.fortran)
krylov.pop()
elif method=='chebyshev':
# Chebyshev algorithm
chebyS = _PETSc.Log.Stage('Chebyshev'); chebyS.push()
_ctqwmpi.qw_cheby(self.psi0.fortran,self.psi.fortran,t,self.H.mat.fortran,
self.H.Emin(),self.H.Emax(),self.rank,self.N)
chebyS.pop()
self.marginal(self.psi)
if self._watchProb:
self.handle.update(self.t,self.psiX,psi2=self.psiY)
if self._watchEnt:
self.entanglementHandle.update(self.t,self.psi)
self.timestep = False
[docs] def plotEntanglement(self,filename):
""" Creates a plot of the entanglement over time.
Args:
filename (str): the absolute/relative path to the desired output file.
.. important::
:func:`watch` **must** be called prior to propagation in order for\
the entanglement to be stored during propagation.
Note:
* ensure a file extension is present so that filetype is correctly set\
(choose one of png, pdf, ps, eps or svg).
* Timesteps are recorded when propagations occur - to increase the number\
of timesteps, increase the number of propagations prior to plotting.
"""
comm = _PETSc.COMM_WORLD
rank = comm.Get_rank()
if rank == 0:
timeArray, entArray = self.entanglementHandle.getEntanglement()
_plot.plotEntanglement(timeArray,entArray,filename,self.initState,self.defectNodes,self.defectAmp)
[docs] def plotNode(self,filename,node,t=None):
""" Creates a plot of the marginal probablities on a *specified node* over time.
Args:
filename (str): the absolute/relative path to the desired output file.
node (int): the node to plot.
.. important::
:func:`watch` **must** be called prior to propagation in order for\
the probabilities at the specified nodes to be stored.
Note:
* ensure a file extension is present so that filetype is correctly set\
(choose one of png, pdf, ps, eps or svg).
* Timesteps are recorded when propagations occur - to increase the number\
of timesteps, increase the number of propagations.
* See :func:`plotNodes` to plot multiple nodes.
"""
comm = _PETSc.COMM_WORLD
rank = comm.Get_rank()
timeArray, probXArray, probYArray = self.handle.getLocalNode(node,p=2)
probXArray = comm.tompi4py().gather(probXArray)
probYArray = comm.tompi4py().gather(probYArray)
if rank == 0:
timeArray = _np.array(timeArray)
probXArray = _np.array([item for sublist in probXArray for item in sublist]).real
probYArray = _np.array([item for sublist in probYArray for item in sublist]).real
_plot.plotNodes2P(timeArray,node,probXArray,probYArray,filename)
[docs] def plotNodes(self,filename,p=1,t=None):
""" Creates a plot of the node probablities over time.
Args:
filename (str): the absolute/relative path to the desired output file.
p (int): (1|2) - choose whether to plot the marginal probability of particle 1 or 2.
.. important::
:func:`watch` **must** be called prior to propagation in order for\
the probabilities at the specified nodes to be stored.
Note:
* if multiple nodes are watched, they will **all** be plotted.
* if you wish to plot marginal probabilities for both particle \
1 and 2, see :func:`plotNode`.
* ensure a file extension is present so that filetype is correctly set\
(choose one of png, pdf, ps, eps or svg).
* Timesteps are recorded when propagations occur - to increase the number\
of timesteps, increase the number of propagations.
"""
comm = _PETSc.COMM_WORLD
rank = comm.Get_rank()
node = self.handle.local_nodes
timeArray, probXArray, probYArray = self.handle.getLocalNodes(t=t,p=2)
if p == 1:
probArray = probXArray
elif p==2:
probArray = probYArray
else:
print 'p must be either 1 or 2'
return
probArray = comm.tompi4py().gather(probArray)
node = comm.tompi4py().gather(node)
if rank == 0:
timeArray = _np.array(timeArray)
nodeArray = _np.array([item for sublist in node for item in sublist])
probArray = _np.array([item for sublist in probArray for item in sublist]).real
_plot.plotNodes(timeArray,nodeArray,probArray,filename,p=p)
[docs] def exportState(self,filename,filetype):
""" Exports the final state :attr:`psi` of the quantum walk to a file.
Args:
filename (str): path to desired output file.
filetype (str): the filetype of the exported adjacency matrix.
* ``'txt'`` - an :math:`N\\times N` 2D array in text format.
* ``'bin'`` - an :math:`N^2` element PETSc binary vector.
"""
if filetype == 'txt':
_io.exportVecToMat(self.psi,filename,filetype)
elif filetype == 'bin':
_io.exportVec(self.psi,filename,filetype)
[docs] def exportPartialTrace(self,filename,filetype,p=1):
""" Exports the partial trace of the density matrix to a file.
Args:
filename (str): path to desired output file.
filetype (str): the filetype of the exported adjacency matrix.
* ``'txt'`` - an :math:`N\\times N` 2D array in text format.
* ``'bin'`` - an :math:`N\\times N` PETSc binary matrix.
p (int): (1|2) - the particle to trace over.
For example, when ``p=1``, :math:`\\rho_2 = \\text{Tr}_1(\\rho_{12})`,
where :math:`\\rho_{12} = |\\psi(t)\\rangle\\langle\\psi(t)|`, is exported.
.. admonition:: Fortran interface
This function calls the Fortran function :f:func:`partial_trace_mat` for ``bin``
filetypes, and the function :f:func:`partial_trace_array` otherwise.
"""
if p == 1:
if filetype == 'bin':
rho = _PETSc.Mat().create(_PETSc.COMM_WORLD)
_ctqwmpi.partial_trace_mat(self.psi.fortran,rho.fortran,self.N)
_io.exportMat(rho, filename, filetype)
rho.destroy()
else:
rho = _ctqwmpi.partial_trace_array(self.psi.fortran,self.N)
if self.rank == 0: _np.savetxt(filename, rho.real)
if p == 2:
#to do
pass
[docs] def psiToInit(self):
""" Copies the state :attr:`self.psi` to :attr:`self.psi0`.
Example:
This can be used to iterate the propagation
over discrete timesteps: ::
for i in range(10):
walk.propagate(0.01,method='chebyshev')
walk.psiToInit()
"""
self.psi0 = self.psi
self.timestep = True
[docs] def destroy(self):
""" Destroys the 2 particle quantum walk object, and all \
associated PETSc matrices/vectors.
"""
self.H.destroy()
self.psi.destroy()
self.psi0.destroy()
self.psiX.destroy()
self.psiY.destroy()
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#--------------------------- 3 particle CTQW ---------------------------------
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[docs]class QuantumWalkP3(object):
""" Methods providing a very basic framework of a 3 particle CTQW.
.. caution::
This object constructs some of the *framework* of a 3 particle
quantum walk, but is not feature complete. For a feature complete
construction, see :func:`pyCTQW.MPI.Graph3P` or :func:`pyCTQW.MPI.Line3P`
for which this is a base.
"""
def __init__(self,N):
self.rank = _PETSc.Comm.Get_rank(_PETSc.COMM_WORLD)
self.N = N
self.t = 0
self.timestep = False
self._watchProb = False
self._watchEnt = False
# define vectors
self.psi0 = _PETSc.Vec()
self.psi0.create(_PETSc.COMM_WORLD)
self.psi0.setSizes(self.N**3)
self.psi0.setUp()
self.psi = self.psi0.duplicate()
# create additional vectors
self.psiX = _PETSc.Vec()
self.psiX.create(_PETSc.COMM_WORLD)
self.psiX.setSizes(self.N)
self.psiX.setUp()
self.psiY = self.psiX.duplicate()
self.psiZ = self.psiX.duplicate()
# define matrices
self.H = Hamiltonian(self.N)
self.EigSolver = self.H.EigSolver
[docs] def importInitState(self,filename,filetype):
""" Imports the initial state of the quantum walk from a file.
Args:
filename (str): path to the file containing the input state.
filetype (str): the filetype of the imported adjacency matrix.
* ``'txt'`` - TODO: not yet implemented! Please use ``'bin'`` for now.
* ``'bin'`` - an :math:`N^3` element PETSc binary *vector*.
Returns:
: this creates a PETSc vector containing the initial statespace, \
accessed via :attr:`self.psi0`.
:rtype: :func:`petsc4py.PETSc.Vec`
Warning:
The number of elements in the imported input statespace
**must** equal :math:`N^3` where :math:`N` is the number
of nodes the 3P CTQW object is initialized with.
"""
self.initState = 'file:'+filename
# create the inital stage
initStateS = _PETSc.Log.Stage('initState')
initStateS.push()
try:
if filetype == 'txt':
print 'Text file state import not supported for three walkers.'
print 'Please try again with a binary file state import.'
return
elif filetype == 'bin':
self.psi0 = _io.loadVec(filename,filetype)
except:
print '\nERROR: incorrect state (is it the correct length?)'
_sys.exit()
initStateS.pop()
self.marginal(self.psi0)
[docs] def marginal(self,vec):
# calculate marginal probabilities
Marginal = _PETSc.Log.Stage('Marginal'); Marginal.push()
_ctqwmpi.marginal3(vec.fortran,self.psiX.fortran,'x',self.N)
_ctqwmpi.marginal3(vec.fortran,self.psiY.fortran,'y',self.N)
_ctqwmpi.marginal3(vec.fortran,self.psiZ.fortran,'z',self.N)
Marginal.pop()
[docs] def watch(self,nodes,watchtype='prob',**kwargs):
""" Creates a handle that watches either marginal probability \
or walker entanglement during propagation.
Args:
nodes (array of ints): the nodes to watch marginal probability (e.g. ``[0,1,4]``).
watchtype (str): (`'prob'` , `'entanglement'`).
the type of watch handle to produce.
Keyword Args:
: If ``watchtype='entanglement'``, EigSolver keywords can also be passed;
for more details of the available EigSolver properties,
see :func:`propagate`.
Returns:
: * if ``watchtype='prob'``, creates a handle that can be \
accessed to retrieve marginal node probabilities for various :math:`t`
* if ``watchtype='entanglment'``, creates a handle that can be \
accessed to retrieve entanglement values for various :math:`t`
:rtype: * if ``watchtype='prob'``: :func:`pyCTQW.MPI.ctqw.nodeHandle`
* if ``watchtype='entanglement'``: :func:`pyCTQW.MPI.ctqw.entanglementHandle`
Examples:
To watch the entanglement,
>>> walk.watch(None,watchtype='entanglement')
>>> walk.propagate(5.,method='chebyshev')
>>> timeArray, entArray = walk.entanglementHandle.getEntanglement()
.. note::
* The entanglement measure used in Von Neumann entropy, calculated \
via :math:`S=-\\sum_{i}\\lambda_i\log_{2}\\lambda_i`, where :math:`\lambda_i` \
are the eigenvalues of the reduced density matrix \
:math:`\\rho_2 = \\text{Tr}_1(|\\psi(t)\\rangle\\langle\\psi(t)|)`
* Nodes do not need to be specified, as entanglement is a global measurement.
* As it is a global measurement, there is a large amount of node communication \
which may increase overall program run time.
To watch the probabilities,
>>> walk.watch([0,1,4])
>>> walk.propagate(2.,method='chebyshev')
>>> timeArray, probXArray, probYArray, probZArray = walk.handle.getLocalNodes(p=2)
.. warning::
Note that `walk.handle` attributes are **not** collective;
if running on multiple nodes, only *local* values will be
returned.
"""
if watchtype == 'prob':
self._watchProb = True
self.handle = nodeHandle(nodes,self.t,self.psiX,psi2=self.psiY,psi3=self.psiZ)
elif watchtype == 'entanglement':
self._watchEnt = True
self.entanglementHandle = entanglementHandle(self.t,self.psi0,self.N,**kwargs)
[docs] def propagate(self,t,method='expm',**kwargs):
""" Propagates the quantum walk for time t.
Args:
t (float): the timestep over which to propagate the state ``walk.psi0`` via the relationship
:math:`\\left|\\psi\\right\\rangle = e^{-iHt}\\left|\\psi0\\right\\rangle`.
method (str): the propagation algorithm to use
(``'chebyshev'`` *(default)* or ``'krylov'``).
Keyword Args:
esolver (str): the eigensolver algorithm to use.
* ``'krylovschur'`` *(default)* - Krylov-Schur
* ``'arnoldi'`` - Arnoldi Method
* ``'lanczos'`` - Lanczos Method
* ``'power'`` - Power Iteration/Rayleigh Quotient Iteration
* ``'gd'`` - Generalized Davidson
* ``'jd'`` - Jacobi-Davidson,
* ``'lapack'`` - Uses LAPACK eigenvalue solver routines
* ``'arpack'`` - *only available if SLEPc is\
compiled with ARPACK linking*
workType (str): can be used to set the eigensolver worktype \
(either ``'ncv'`` or ``'mpd'``). The default \
is to let SLEPc decide.
workSize (int): sets the work size **if** :attr:`workType` is set.
tolIn (float): tolerance of the eigenvalue solver \
(*default* ``0.`` (SLEPc decides)).
maxIt (int): maximum number of iterations of the eigenvalue solver \
(*default* ``0`` (SLEPc decides)).
verbose (bool): if ``True``, writes eigensolver information to the console
emax_estimate (float): used to override the calculation \
of the graphs maximum eigenvalue.
emin_estimate (float): used to override the calculation \
of the graphs minimum eigenvalue. For a *finite* graph, \
``emin_estimate=0`` is set by default.
.. caution::
* If supplied, the value of :attr:`emax_estimate` (:math:`\hat{\lambda}_{\max}`) \
**must** satisfy :math:`\hat{\lambda}_{\max}\geq\lambda_{\max}`, \
where :math:`\lambda_{\max}` is the actual maximum eigenvalue of the graph.
* Similarly, the value of :attr:`emin_estimate` (:math:`\hat{\lambda}_{\min}`) \
**must** satisfy :math:`\hat{\lambda}_{\min}\leq\lambda_{\min}`, \
where :math:`\lambda_{\min}` is the actual minimum eigenvalue of the graph.
* The greater the difference value of :math:`|\hat{\lambda}_{\max} \
-\lambda_{\max}|` and :math:`|\hat{\lambda}_{\min} \
-\lambda_{\min}|`, the longer the convergence \
time of the ``chebyshev`` propagator.
Note:
* The keyword arguments properties only apply if ``propagator='chebyshev'``
* For more information regarding these properties,refer to \
the `SLEPc documentation\
<http://www.grycap.upv.es/slepc/documentation/manual.htm>`_
Returns:
: this creates a PETSc vector containing the
propagated state, accessed via the attribute
:attr:`psi`, as well as a PETSc vector containing
the marginal probabilities, :attr:`prob`.
:rtype: :func:`petsc4py.PETSc.Vec`
Note:
The EigSolver properties only take effect if
`method='chebyshev'`.
.. admonition:: Fortran interface
This function calls the Fortran function :f:func:`expm`. for the Krylov
method, and the function :f:func:`qw_cheby` for the Chebyshev method.
"""
if self.timestep:
self.t = t + self.t
else:
self.t = t
self.EigSolver.setEigSolver(**kwargs)
if method=='krylov':
# SLEPc matrix exponential
krylov = _PETSc.Log.Stage('SLEPc krylov'); krylov.push()
_ctqwmpi.qw_krylov(self.H.mat.fortran,t,self.psi0.fortran,self.psi.fortran)
krylov.pop()
elif method=='chebyshev':
# Chebyshev algorithm
chebyS = _PETSc.Log.Stage('Chebyshev'); chebyS.push()
_ctqwmpi.qw_cheby(self.psi0.fortran,self.psi.fortran,t,self.H.mat.fortran,
self.H.Emin(),self.H.Emax(),self.rank,self.N)
chebyS.pop()
self.marginal(self.psi)
if self._watchProb:
self.handle.update(self.t,self.psiX,psi2=self.psiY,psi3=self.psiZ)
if self._watchEnt:
self.entanglementHandle.update(self.t,self.psi)
self.timestep = False
[docs] def plotEntanglement(self,filename):
""" Creates a plot of the entanglement over time.
Args:
filename (str): the absolute/relative path to the desired output file.
.. important::
:func:`watch` **must** be called prior to propagation in order for\
the entanglement to be stored during propagation.
Note:
* ensure a file extension is present so that filetype is correctly set\
(choose one of png, pdf, ps, eps or svg).
* Timesteps are recorded when propagations occur - to increase the number\
of timesteps, increase the number of propagations prior to plotting.
"""
comm = _PETSc.COMM_WORLD
rank = comm.Get_rank()
if rank == 0:
timeArray, entArray = self.entanglementHandle.getEntanglement()
_plot.plotEntanglement(timeArray,entArray,filename,self.initState,self.defectNodes,self.defectAmp)
[docs] def plotNode(self,filename,node,t=None):
""" Creates a plot of the marginal probablities on a *specified node* over time.
Args:
filename (str): the absolute/relative path to the desired output file.
node (int): the node to plot.
.. important::
:func:`watch` **must** be called prior to propagation in order for\
the probabilities at the specified nodes to be stored.
Note:
* ensure a file extension is present so that filetype is correctly set\
(choose one of png, pdf, ps, eps or svg).
* Timesteps are recorded when propagations occur - to increase the number\
of timesteps, increase the number of propagations.
* See :func:`plotNodes` to plot multiple nodes.
"""
comm = _PETSc.COMM_WORLD
rank = comm.Get_rank()
timeArray, probXArray, probYArray, probZArray = self.handle.getLocalNode(node,p=3)
probXArray = comm.tompi4py().gather(probXArray)
probYArray = comm.tompi4py().gather(probYArray)
probZArray = comm.tompi4py().gather(probZArray)
if rank == 0:
timeArray = _np.array(timeArray)
probXArray = _np.array([item for sublist in probXArray for item in sublist]).real
probYArray = _np.array([item for sublist in probYArray for item in sublist]).real
probZArray = _np.array([item for sublist in probZArray for item in sublist]).real
_plot.plotNodes3P(timeArray,node,probXArray,probYArray,probZArray,filename)
[docs] def plotNodes(self,filename,p=1,t=None):
""" Creates a plot of the node probablities over time.
Args:
filename (str): the absolute/relative path to the desired output file.
p (int): (1|2|3) - choose whether to plot the marginal probability of particle 1, 2 or 3.
.. important::
:func:`watch` **must** be called prior to propagation in order for\
the probabilities at the specified nodes to be stored.
Note:
* if multiple nodes are watched, they will **all** be plotted.
* if you wish to plot marginal probabilities for all particles, see :func:`plotNode`.
* ensure a file extension is present so that filetype is correctly set\
(choose one of png, pdf, ps, eps or svg).
* Timesteps are recorded when propagations occur - to increase the number\
of timesteps, increase the number of propagations.
"""
comm = _PETSc.COMM_WORLD
rank = comm.Get_rank()
node = self.handle.local_nodes
timeArray, probXArray, probYArray, probZArray = self.handle.getLocalNodes(t=t,p=3)
if p == 1:
probArray = probXArray
elif p==2:
probArray = probYArray
elif p==3:
probArray = probZArray
else:
print 'p must be either 1, 2 or 3'
return
probArray = comm.tompi4py().gather(probArray)
node = comm.tompi4py().gather(node)
if rank == 0:
timeArray = _np.array(timeArray)
nodeArray = _np.array([item for sublist in node for item in sublist])
probArray = _np.array([item for sublist in probArray for item in sublist]).real
_plot.plotNodes(timeArray,nodeArray,probArray,filename,p=p)
[docs] def exportState(self,filename,filetype):
""" Exports the final state :attr:`psi` of the quantum walk to a file.
Args:
filename (str): path to desired output file.
filetype (str): the filetype of the exported adjacency matrix.
* ``'txt'`` - an :math:`N^3` element column vector in text format.
* ``'bin'`` - an :math:`N^3` element PETSc binary vector.
"""
if filetype == 'txt':
_io.exportVec(self.psi,filename,filetype)
elif filetype == 'bin':
_io.exportVec(self.psi,filename,filetype)
[docs] def psiToInit(self):
""" Copies the state :attr:`psi` to :attr:`psi0`.
Example:
This can be used to iterate the propagation
over discrete timesteps: ::
for i in range(10):
walk.propagate(0.01,method='chebyshev')
walk.psiToInit()
"""
self.psi0 = self.psi
self.timestep = True
[docs] def destroy(self):
""" Destroys the 3 particle quantum walk object, and all \
associated PETSc matrices/vectors.
"""
self.H.destroy()
self.psi.destroy()
self.psi0.destroy()
self.psiX.destroy()
self.psiY.destroy()
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#------------------------------- Arbitrary CTQW --------------------------------
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
class Graph(QuantumWalkP1):
""" Performs and analyses 1 particle continuous-time quantum walks on graphs
Args:
N (int) : number of nodes to initialize the walker with. Nodes are \
labeled :math:`j\\in\\{0,N-1\\}`.
Example:
To create a CTQW Graph object for a 10 node graph,
>>> walk = pyCTQW.MPI.Graph(10)
Note:
**if filename AND filetype are provided**, this automatically
creates a PETSc Hamiltonian matrix, neglecting the need to run
:func:`createH`. For details on the other keyword arguments,
see :func:`createH`.
"""
def __init__(self,N,filename=None,filetype=None,d=None,amp=None,layout='spring',delimiter=None):
QuantumWalkP1.__init__(self,N)
self.liveplot = False
self.EigSolver.setEigSolver(emin_estimate=0.)
if (filename and filetype) is not None:
self.createH(filename,filetype,d=d,amp=amp,layout=layout,delimiter=delimiter)
def createH(self,filename,filetype,d=None,amp=None,layout='spring',delimiter=None):
""" Generate the Hamiltonian of the graph.
Args:
filename (str): path to the file containing the adjacency matrix of the graph
filetype (str): the filetype of the imported adjacency matrix.
* ``'txt'`` - an :math:`N\\times N` dense 2D array in text format.
* ``'bin'`` - an :math:`N\\times N` PETSc binary matrix.
d (array of ints): an array containing *integers* indicating the nodes
where diagonal defects are to be placed (e.g. ``d=[0,1,4]``).
amp (array of floats): an array containing *floats* indicating the diagonal defect
amplitudes corresponding to each element in ``d`` (e.g. ``amp=[0.5,-1,4.2]``).
layout (str): the format to store the position of the nodes (only used
when running :func:`plotGraph`).
* ``spring`` *(default)* - spring layout.
* ``circle`` - nodes are arranged in a circle.
* ``spectral`` - nodes are laid out according to the \
spectrum of the graph.
* ``random`` - nodes are arranged in a random pattern.
delimiter (str): this is passed to `numpy.genfromtxt\
<http://docs.scipy.org/doc/numpy/reference/generated/numpy.genfromtxt.html>`_
in the case of strange delimiters in an imported ``txt`` file.
Returns:
:this creates a Hamiltonian object, accessed via the attibute
:attr:`Graph.H`.
:rtype: :func:`pyCTQW.MPI.ctqw.Hamiltonian`
Note:
This needs to be called **only** if the filename and filetype
of the graph were not already called when the Graph object
was initialized.
.. important::
* The number of nodes in the imported adjacency matrix \
**must** match the number of nodes the Graph object \
is initialized with.
* The size of ``amp`` and ``d`` must be identical
>>> amp = [0.5,-1.,4.2]
>>> len(d) == len(amp)
True
"""
if (d and amp) is None:
self.defectNodes = [0]
self.defectAmp = [0.]
else:
self.defectNodes = d
self.defectAmp = amp
self.H.importAdj(filename,filetype,d=self.defectNodes,amp=self.defectAmp,layout=layout,delimiter=delimiter)
def createInitState(self,initState):
""" Generate the initial state of the quantum walk.
Args:
initState (array) : an :math:`n\\times 2` array, containing the initial \
state of the quantum walker in the format ``[[j1,amp1],[j2,amp2],...]``.
Returns:
: this creates a PETSc vector containing the initial state,
accessed via the attribute :attr:`Graph.ps0`.
:rtype: :func:`petsc4py.PETSc.Vec`
Example:
For a CTQW initially located in a
superposition of nodes 1 and 2, e.g.
:math:`\\left|\\psi(0)\\right\\rangle = \\frac{1}{\\sqrt{2}}\\left|1\\right\\rangle \
- \\frac{1}{\\sqrt{2}} \\left|2\\right\\rangle`,
the initial state would be created like so:
>>> import numpy as np
>>> init_state = [[1,1./np.sqrt(2.)], [2,-1./np.sqrt(2.)]]
>>> walk.createInitState(init_state)
"""
self.initState = _np.vstack([_np.array(initState).T[0]-self.N/2+1,
_np.array(initState).T[1]]).T.tolist()
# create the inital stage
initStateS = _PETSc.Log.Stage('initState')
initStateS.push()
_ctqwmpi.p1_init(self.psi0.fortran,self.initState,self.N)
initStateS.pop()
self.marginal(self.psi0)
def plot(self,filename):
""" Creates a plot of probability vs node.
Args:
filename (str): the absolute/relative path to the desired output file.
Note:
* ensure a file extension is present so that filetype is correctly set\
(choose one of png, pdf, ps, eps or svg).
* if :func:`propagate` has not been called, the probability of the \
initial state is plotted (:math:`t=0`). Otherwise, the last propagated \
time (:attr:`Graph.t`) is used.
"""
if _os.path.isabs(filename):
outDir = _os.path.dirname(filename)
else:
outDir = './'+_os.path.dirname(filename)
# create output directory if it doesn't exist
try:
_os.mkdir(outDir)
except OSError as exception:
if exception.errno != _errno.EEXIST:
raise
initstateLabels = []
for i in range(len(self.initState)):
initstateLabels.append([sum(pair) for pair in zip(self.initState[i], [self.N/2-1,0])])
plotStage = _PETSc.Log.Stage('Plotting'); plotStage.push()
_plot.plot(_np.arange(self.N),self.prob,filename,self.t,initstateLabels,
self.defectNodes,self.defectAmp,self.N,self.rank)
plotStage.pop()
def plotGraph(self,output=None,probX=True,size=(12,8),**kwargs):
""" Creates a plot of probability vs node superimposed on a 3D visualisation of the graph vertices.
Args:
output (str) : the absolute/relative path to the desired output file.
probX (bool) : if set to ``True`` *(default)*, the probability is represented
as bars placed on each graph vertex; if ``False``, the graph
is plotted without any probability.
size (tuple) : ``size=(x,y)`` sets the horizontal and vertical size of the output figure.
Keyword Args:
nodesize (float) : size of the vertices in the plot. If left blank, this is
determined automatically.
nodecolor (str) : vertex color (*default* ``'red'``).
For more details on how to specify a color,
see the `matplotlib documentation
<http://matplotlib.org/api/colors_api.html#module-matplotlib.colors>`_.
nodealpha (float) : value between 0 and 1 specifying the vertex opacity (*default* 0.25)
nodetext (bool) : if set ``True``, the vertices are labelled by number.
nodetextcolor (str) : vertex label text color (*default* ``'black'``).
nodetextbg (str) : vertex label background color (*default* ``'None'``).
ntofffset (array of floats): the :math:`(x,y,z)` vertex label offset relative \
to the vertex (*default* ``[0.,0.,-0.15]``).
barscale (float) : scaled height of the probability bars (*default* ``1``).
barcolor (str) : probability bar color (*default* ``'green'``).
baralpha (float) : value between 0 and 1 specifying the opacity (*default* 0.25)
bartext (bool) : if set ``True``, the probability bars are labelled with their value.
bartextcolor (str) : probability label text color (*default* ``'black'``).
bartextbg (str) : probability label background color (*default* ``'None'``).
btoffset (array of floats): the :math:`(x,y,z)` probability label offset relative to the top of
the probability bars (*default* ``[-0.025,-0.025,0.05]``)
Note:
* ensure a file extension is present so that file type is correctly set\
(choose one of png, pdf, ps, eps or svg).
* if :func:`propagate` has not been called, the probability of the \
initial state is plotted (:math:`t=0`). Otherwise, the last propagated \
time (``walk.t``) is used.
"""
rank = _PETSc.COMM_WORLD.Get_rank()
if probX:
# scatter prob to process 0
commX = self.prob.getComm()
rank = _PETSc.COMM_WORLD.getRank()
scatterX, probX0 = _PETSc.Scatter.toZero(self.prob)
scatterX.scatter(self.prob, probX0, False, _PETSc.Scatter.Mode.FORWARD)
if rank==0:
from matplotlib import pyplot as plt
import mpl_toolkits.mplot3d as plt3d
self.fig = plt.figure(figsize=size)
self.ax = plt3d.Axes3D(self.fig)
self.ax.view_init(45, -50)
self.ax.set_axis_off()
if probX:
prob = _np.real(_np.asarray(probX0))
else:
prob = None
_plot.plotGraph(self.ax,self.H.nodePos,self.H.lineX,self.H.lineY,
prob=prob,**kwargs)
self.ax.set_title('$t={}$'.format(self.t))
if type(output) is str:
plt.savefig(output)
else:
plt.show(block=True)
plt.close()
if probX:
# deallocate
commX.barrier()
scatterX.destroy()
probX0.destroy()
def plotLiveGraph(self,dt,size=(12,8),**kwargs):
""" Creates a *live*, updated plot of probability vs node superimposed \
on a 3D visualisation of the graph vertices.
Args:
dt (str) : the amount of time to 'sleep' before resuming the program.
This can be used to 'slow down' propagation, providing
time to view small changes on the graph.
size (tuple) : ``size=(x,y)`` sets the horizontal and vertical size of the output figure.
Keyword Args:
: For available keyword arguments, see :func:`Graph.plotGraph`.
Note:
* Once the live plot is no longer needed, :func:`Graph.clearLiveGraph()` \
should be called in order to destroy the live plot object (this will *not* \
close the plot window).
* if :func:`propagate` has not been called, the probability of the \
initial state is plotted (:math:`t=0`). Otherwise, the last propagated \
time (:attr:`Graph.t`) is used.
Warning:
This feature attempts to uses ``matplotlib`` in a sort of `hackish` way
and is not well tested - your mileage may vary.
"""
rank = _PETSc.COMM_WORLD.Get_rank()
# scatter prob to process 0
commX = self.prob.getComm()
rank = _PETSc.COMM_WORLD.getRank()
scatterX, probX0 = _PETSc.Scatter.toZero(self.prob)
scatterX.scatter(self.prob, probX0, False, _PETSc.Scatter.Mode.FORWARD)
if rank==0:
from matplotlib import pyplot as plt
import mpl_toolkits.mplot3d as plt3d
if not self.liveplot:
plt.ion()
self.figLive = plt.figure(figsize=size)
self.liveplot = True
else:
plt.cla()
self.axLive = plt3d.Axes3D(self.figLive)
self.axLive.view_init(45, -50)
self.axLive.set_axis_off()
prob = _np.real(_np.asarray(probX0))
_plot.plotGraph(self.axLive,self.H.nodePos,self.H.lineX,self.H.lineY,
prob=prob,**kwargs)
self.axLive.set_title('$t={}$'.format(self.t))
plt.draw()#show(block=True)
_time.sleep(dt)
# deallocate
commX.barrier()
scatterX.destroy()
probX0.destroy()
def clearLiveGraph(self):
""" Destroys the live plot object previously created by :func:`Graph.plotLiveGraph()`.
"""
rank = _PETSc.COMM_WORLD.Get_rank()
if rank == 0:
from matplotlib import pyplot as plt
plt.show(block=True)
self.liveplot = False
plt.close()
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#------------------------------- 2P Arbitrary CTQW -----------------------------
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
class Graph2P(QuantumWalkP2):
""" Performs and analyses 2 particle continuous-time quantum walks on graphs
Args:
N (int) : number of nodes to initialize the walker with. Nodes are \
labeled :math:`j\\in\\{0,N-1\\}`.
Example:
To create a 2P CTQW Graph object for a 10 node graph,
>>> walk = pyCTQW.MPI.Graph2P(10)
Note:
**if filename AND filetype are provided**, this automatically
creates a PETSc Hamiltonian matrix, neglecting the need to run
:func:`createH`. For details on the other keyword arguments,
see :func:`createH`.
"""
def __init__(self,N,filename=None,filetype=None,d=None,amp=None,interaction=0.,bosonic=False):
QuantumWalkP2.__init__(self,N)
self.liveplot = False
self.EigSolver.setEigSolver(emin_estimate=0.)
if (filename and filetype) is not None:
self.createH(filename,filetype,d=d,amp=amp,interaction=interaction,bosonic=bosonic)
def createH(self,filename,filetype,d=None,amp=None,layout='spring',delimiter=None,interaction=0.,bosonic=False):
""" Generate the Hamiltonian of the graph.
Args:
filename (str): path to the file containing the adjacency matrix of the graph
filetype (str): the filetype of the imported adjacency matrix.
* ``'txt'`` - an :math:`N\\times N` dense 2D array in text format.
* ``'bin'`` - an :math:`N\\times N` PETSc binary matrix.
d (array of ints): an array containing *integers* indicating the nodes
where diagonal defects are to be placed (e.g. ``d=[0,1,4]``).
amp (array of floats): an array containing *floats* indicating the diagonal defect
amplitudes corresponding to each element in ``d`` (e.g. ``amp=[0.5,-1,4.2]``).
layout (str): the format to store the position of the nodes (only used
when running :func:`plotGraph`).
* ``spring`` *(default)* - spring layout.
* ``circle`` - nodes are arranged in a circle.
* ``spectral`` - nodes are laid out according to the \
spectrum of the graph.
* ``random`` - nodes are arranged in a random pattern.
delimiter (str): this is passed to `numpy.genfromtxt\
<http://docs.scipy.org/doc/numpy/reference/generated/numpy.genfromtxt.html>`_
in the case of strange delimiters in an imported ``txt`` file.
interaction (float): the amplitude of interaction between the two walkers \
when located on the same vertex.
bosonic (bool): if true, a Bosonic Hamiltonian is generated.
Returns:
:this creates a Hamiltonian object, accessed via the attibute
:attr:`Graph2P.H`.
:rtype: :func:`pyCTQW.MPI.ctqw.Hamiltonian`
Note:
This needs to be called **only** if the filename and filetype
of the graph were not already called when the Graph object
was initialized.
.. important::
* The number of nodes in the imported adjacency matrix \
**must** match the number of nodes the :class:`Graph2P` object \
is initialized with.
* The size of ``amp`` and ``d`` must be identical
>>> amp = [0.5,-1.,4.2]
>>> len(d) == len(amp)
True
"""
if (d or amp) is None:
self.defectNodes = [0]
self.defectAmp = [0.]
else:
self.defectNodes = d
self.defectAmp = amp
self.interaction = interaction
self.bosonic = bosonic
self.H.importAdjToH(filename,filetype,
d=self.defectNodes,amp=self.defectAmp,p='2',interaction=self.interaction,bosonic=self.bosonic,
layout=layout,delimiter=delimiter)
# _ctqwmpi.importadjtoh(self.H.mat.fortran,filename,'2',
# d=self.defectNodes,amp=self.defectAmp,layout=layout)
def createInitState(self,initState):
""" Generate the initial state of the quantum walk.
Args:
initState (array) : an :math:`n\\times 3` array, containing the initial \
state of the quantum walker in the format ``[[x1,y1,amp1],[x2,y2,amp2],...]``.
Returns:
: this creates a PETSc vector containing the initial state,
accessed via the attribute :attr:`Graph2P.ps0`.
:rtype: :func:`petsc4py.PETSc.Vec`
Example:
For a CTQW initially located in state
:math:`\\left|\\psi(0)\\right\\rangle = \\frac{1}{\\sqrt{2}}\\left|0\\right\\rangle\\otimes\\left|1\\right\\rangle \
- \\frac{1}{\\sqrt{2}} \\left|1\\right\\rangle\\otimes\\left|0\\right\\rangle`,
the initial state would be created like so:
>>> import numpy as np
>>> init_state = [[0,1,1./np.sqrt(2.)], [1,0,-1./np.sqrt(2.)]]
>>> walk.createInitState(init_state)
.. admonition:: Fortran interface
This function calls the Fortran function :f:func:`p2_init`.
"""
self.initState = _np.vstack([_np.array(initState).T[0]-self.N/2+1,
_np.array(initState).T[1]-self.N/2+1,_np.array(initState).T[2]]).T.tolist()
# create the inital stage
initStateS = _PETSc.Log.Stage('initState')
initStateS.push()
_ctqwmpi.p2_init(self.psi0.fortran,self.initState,self.N)
initStateS.pop()
def plot(self,filename):
""" Creates a plot of probability vs node.
Args:
filename (str): the absolute/relative path to the desired output file.
Note:
* ensure a file extension is present so that filetype is correctly set\
(choose one of png, pdf, ps, eps or svg).
* if :func:`propagate` has not been called, the probability of the \
initial state is plotted (:math:`t=0`). Otherwise, the last propagated \
time (:attr:`Graph2P.t`) is used.
"""
if _os.path.isabs(filename):
outDir = _os.path.dirname(filename)
else:
outDir = './'+_os.path.dirname(filename)
# create output directory if it doesn't exist
try:
_os.mkdir(outDir)
except OSError as exception:
if exception.errno != _errno.EEXIST:
raise
initstateLabels = []
for i in range(len(self.initState)):
initstateLabels.append([sum(pair) for pair in zip(self.initState[i], [self.N/2-1,self.N/2-1,0])])
plotStage = _PETSc.Log.Stage('Plotting'); plotStage.push()
_plot.plot2P(_np.arange(self.N),self.psiX,self.psiY,filename,self.t,initstateLabels,
self.defectNodes,self.defectAmp,self.N,self.rank)
plotStage.pop()
def plotGraph(self,size=(12,8),probX=True,probY=True,output=None,**kwargs):
""" Creates a plot of probability vs node superimposed on a 3D visualisation of the graph vertices.
Args:
output (str) : the absolute/relative path to the desired output file.
probX (bool) : if set to ``True`` *(default)*, the particle 1 marginale probability is represented
as bars placed on each graph vertex.
probY (bool) : if set to ``True`` *(default)*, the particle 2 marginal probability is represented
as bars placed below each graph vertex.
size (tuple) : ``size=(x,y)`` sets the horizontal and vertical size of the output figure.
Keyword Args:
nodesize (float) : size of the vertices in the plot. If left blank, this is
determined automatically.
nodecolor (str) : vertex color (*default* ``'red'``).
For more details on how to specify a color,
see the `matplotlib documentation
<http://matplotlib.org/api/colors_api.html#module-matplotlib.colors>`_.
nodealpha (float) : value between 0 and 1 specifying the vertex opacity (*default* 0.25)
nodetext (bool) : if set ``True``, the vertices are labelled by number.
nodetextcolor (str) : vertex label text color (*default* ``'black'``).
nodetextbg (str) : vertex label background color (*default* ``'None'``).
ntofffset (array of floats): the :math:`(x,y,z)` vertex label offset relative \
to the vertex (*default* ``[0.,0.,-0.15]``).
barscaleP (float) : scaled height of the probability bars (*default* ``1``).
barcolorP (str) : probability bar color (*default* ``'green'``).
baralphaP (float) : value between 0 and 1 specifying the opacity (*default* 0.25)
bartext (bool) : if set ``True``, the probability bars are labelled with their value.
bartextcolorP (str) : probability label text color (*default* ``'black'``).
bartextbgP (str) : probability label background color (*default* ``'None'``).
btoffsetP (array of floats): the :math:`(x,y,z)` probability label offset relative to the top of \
the probability bars (*default* ``[-0.025,-0.025,+-0.05]``)
.. important:
where a keyword argument ends with ``P`` above, substitute ``P=(1|2)`` to
access the probability bar property for particle 1 and 2 respectively.
Note:
* ensure a file extension is present so that file type is correctly set\
(choose one of png, pdf, ps, eps or svg).
* if :func:`propagate` has not been called, the probability of the \
initial state is plotted (:math:`t=0`). Otherwise, the last propagated \
time (``walk.t``) is used.
"""
rank = _PETSc.COMM_WORLD.Get_rank()
if probX:
# scatter prob to process 0
commX = self.psiX.getComm()
rank = _PETSc.COMM_WORLD.getRank()
scatterX, probX0 = _PETSc.Scatter.toZero(self.psiX)
scatterX.scatter(self.psiX, probX0, False, _PETSc.Scatter.Mode.FORWARD)
if probY:
# scatter prob to process 0
commY = self.psiY.getComm()
rank = _PETSc.COMM_WORLD.getRank()
scatterY, probY0 = _PETSc.Scatter.toZero(self.psiY)
scatterX.scatter(self.psiY, probY0, False, _PETSc.Scatter.Mode.FORWARD)
if rank==0:
from matplotlib import pyplot as plt
import mpl_toolkits.mplot3d as plt3d
self.fig = plt.figure(figsize=size)
self.ax = plt3d.Axes3D(self.fig)
self.ax.view_init(45, -50)
self.ax.set_axis_off()
if probX:
prob = _np.real(_np.asarray(probX0))
else:
prob = None
if probY:
prob2 = _np.real(_np.asarray(probY0))
else:
prob2 = None
_plot.plotGraph(self.ax,self.H.nodePos,self.H.lineX,self.H.lineY,
prob=prob,prob2=prob2,**kwargs)
self.ax.set_title('$t={}$'.format(self.t))
if type(output) is str:
plt.savefig(output)
else:
plt.show(block=True)
if probX:
# deallocate
commX.barrier()
scatterX.destroy()
probX0.destroy()
if probY:
# deallocate
commY.barrier()
scatterY.destroy()
probY0.destroy()
def plotLiveGraph(self,dt,size=(12,8),
probX=True,probY=True,**kwargs):
""" Creates a *live*, updated plot of probability vs node superimposed \
on a 3D visualisation of the graph vertices.
Args:
dt (str) : the amount of time to 'sleep' before resuming the program.
This can be used to 'slow down' propagation, providing
time to view small changes on the graph.
size (tuple) : ``size=(x,y)`` sets the horizontal and vertical size of the output figure.
Keyword Args:
: For available keyword arguments, see :func:`Graph2P.plotGraph`.
Note:
* Once the live plot is no longer needed, :func:`Graph.clearLiveGraph()` \
should be called in order to destroy the live plot object (this will *not* \
close the plot window).
* if :func:`propagate` has not been called, the probability of the \
initial state is plotted (:math:`t=0`). Otherwise, the last propagated \
time (:attr:`Graph2P.t`) is used.
Warning:
This feature attempts to uses ``matplotlib`` in a sort of `hackish` way
and is not well tested - your mileage may vary.
"""
rank = _PETSc.COMM_WORLD.Get_rank()
if probX:
# scatter prob to process 0
commX = self.psiX.getComm()
rank = _PETSc.COMM_WORLD.getRank()
scatterX, probX0 = _PETSc.Scatter.toZero(self.psiX)
scatterX.scatter(self.psiX, probX0, False, _PETSc.Scatter.Mode.FORWARD)
if probY:
# scatter prob to process 0
commY = self.psiY.getComm()
rank = _PETSc.COMM_WORLD.getRank()
scatterY, probY0 = _PETSc.Scatter.toZero(self.psiY)
scatterX.scatter(self.psiY, probY0, False, _PETSc.Scatter.Mode.FORWARD)
if rank==0:
from matplotlib import pyplot as plt
import mpl_toolkits.mplot3d as plt3d
if not self.liveplot:
plt.ion()
self.figLive = plt.figure(figsize=size)
self.liveplot = True
else:
plt.cla()
self.axLive = plt3d.Axes3D(self.figLive)
self.axLive.view_init(45, -50)
self.axLive.set_axis_off()
if probX:
prob = _np.real(_np.asarray(probX0))
else:
prob = None
if probY:
prob2 = _np.real(_np.asarray(probY0))
else:
prob2 = None
_plot.plotGraph(self.axLive,self.H.nodePos,self.H.lineX,self.H.lineY,
prob=prob,prob2=prob2,**kwargs)
self.axLive.set_title('$t={}$'.format(self.t))
plt.draw()#show(block=True)
_time.sleep(dt)
if probX:
# deallocate
commX.barrier()
scatterX.destroy()
probX0.destroy()
if probY:
# deallocate
commY.barrier()
scatterY.destroy()
probY0.destroy()
def clearLiveGraph(self):
""" Destroys the live plot object previously created by :func:`Graph2P.plotLiveGraph()`.
"""
rank = _PETSc.COMM_WORLD.Get_rank()
if rank == 0:
from matplotlib import pyplot as plt
plt.show(block=True)
self.liveplot = False
plt.close()
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#------------------------------- 3P Arbitrary CTQW -----------------------------
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
class Graph3P(QuantumWalkP3):
""" Performs and analyses 3 particle continuous-time quantum walks on graphs
Args:
N (int) : number of nodes to initialize the walker with. Nodes are \
labeled :math:`j\\in\\{0,N-1\\}`.
Example:
To create a 3P CTQW Graph object for a 10 node graph,
>>> walk = pyCTQW.MPI.Graph3P(10)
Note:
**if filename AND filetype are provided**, this automatically
creates a PETSc Hamiltonian matrix, neglecting the need to run
:func:`createH`. For details on the other keyword arguments,
see :func:`createH`.
"""
def __init__(self,N,filename=None,filetype=None,d=None,amp=None,interaction=0.):
QuantumWalkP3.__init__(self,N)
self.EigSolver.setEigSolver(emin_estimate=0.)
if (filename and filetype) is not None:
self.createH(filename,filetype,d=d,amp=amp,interaction=interaction)
def createH(self,filename,filetype,d=None,amp=None,layout='spring',delimiter=None,interaction=0.):
""" Generate the Hamiltonian of the graph.
Args:
filename (str): path to the file containing the adjacency matrix of the graph
filetype (str): the filetype of the imported adjacency matrix.
* ``'txt'`` - an :math:`N\\times N` dense 2D array in text format.
* ``'bin'`` - an :math:`N\\times N` PETSc binary matrix.
d (array of ints): an array containing *integers* indicating the nodes
where diagonal defects are to be placed (e.g. ``d=[0,1,4]``).
amp (array of floats): an array containing *floats* indicating the diagonal defect
amplitudes corresponding to each element in ``d`` (e.g. ``amp=[0.5,-1,4.2]``).
layout (str): the format to store the position of the nodes (only used
when running :func:`plotGraph`).
* ``spring`` *(default)* - spring layout.
* ``circle`` - nodes are arranged in a circle.
* ``spectral`` - nodes are laid out according to the \
spectrum of the graph.
* ``random`` - nodes are arranged in a random pattern.
delimiter (str): this is passed to `numpy.genfromtxt\
<http://docs.scipy.org/doc/numpy/reference/generated/numpy.genfromtxt.html>`_
in the case of strange delimiters in an imported ``txt`` file.
interaction (float): the amplitude of interaction between the two walkers \
when located on the same vertex.
Returns:
:this creates a Hamiltonian object, accessed via the attibute
:attr:`Graph3P.H`.
:rtype: :func:`pyCTQW.MPI.ctqw.Hamiltonian`
Note:
This needs to be called **only** if the filename and filetype
of the graph were not already called when the Graph object
was initialized.
.. important::
* The number of nodes in the imported adjacency matrix \
**must** match the number of nodes the :class:`Graph2P` object \
is initialized with.
* The size of ``amp`` and ``d`` must be identical
>>> amp = [0.5,-1.,4.2]
>>> len(d) == len(amp)
True
"""
if (d and amp) is None:
self.defectNodes = [0]
self.defectAmp = [0.]
else:
self.defectNodes = d
self.defectAmp = amp
self.interaction = interaction
self.H.importAdjToH(filename,filetype,
d=self.defectNodes,amp=self.defectAmp,p='3',interaction=self.interaction,
layout=layout,delimiter=delimiter)
# _ctqwmpi.importadjtoh(self.H.mat.fortran,filename,'2',
# d=self.defectNodes,amp=self.defectAmp,layout=layout)
def createInitState(self,initState):
""" Generate the initial state of the quantum walk.
Args:
initState (array) : an :math:`n\\times 4` array, containing the initial \
state of the quantum walker in the format ``[[x1,y1,z1,amp1],[x2,y2,z2,amp2],...]``.
Returns:
: this creates a PETSc vector containing the initial state,
accessed via the attribute :attr:`Graph3P.ps0`.
:rtype: :func:`petsc4py.PETSc.Vec`
Example:
For a CTQW initially located in state
:math:`\\left|\\psi(0)\\right\\rangle = \\frac{1}{\\sqrt{2}}\\left|0\\right\\rangle \
\\otimes\\left|1\\right\\rangle \\otimes\\left|5\\right\\rangle \
- \\frac{i}{\\sqrt{2}} \\left|1\\right\\rangle\\otimes\\left|0\\right\\rangle \
\\otimes\\left|2\\right\\rangle`,
the initial state would be created like so:
>>> import numpy as np
>>> init_state = [[0,1,5,1./np.sqrt(2.)], [1,0,2,-1j./np.sqrt(2.)]]
>>> walk.createInitState(init_state)
.. admonition:: Fortran interface
This function calls the Fortran function :f:func:`p3_init`.
"""
self.initState = _np.vstack( [ _np.array(initState).T[0],
_np.array(initState).T[1],
_np.array(initState).T[2],
_np.array(initState).T[3] ]).T.tolist()
# create the inital stage
initStateS = _PETSc.Log.Stage('initState')
initStateS.push()
_ctqwmpi.p3_init(self.psi0.fortran,self.initState,self.N)
initStateS.pop()
def plot(self,filename):
""" Creates a plot of probability vs node.
Args:
filename (str): the absolute/relative path to the desired output file.
Note:
* ensure a file extension is present so that filetype is correctly set\
(choose one of png, pdf, ps, eps or svg).
* if :func:`propagate` has not been called, the probability of the \
initial state is plotted (:math:`t=0`). Otherwise, the last propagated \
time (:attr:`Graph3P.t`) is used.
"""
if _os.path.isabs(filename):
outDir = _os.path.dirname(filename)
else:
outDir = './'+_os.path.dirname(filename)
# create output directory if it doesn't exist
try:
_os.mkdir(outDir)
except OSError as exception:
if exception.errno != _errno.EEXIST:
raise
plotStage = _PETSc.Log.Stage('Plotting'); plotStage.push()
_plot.plot3P(_np.arange(self.N),self.psiX,self.psiY,self.psiZ,filename,self.t,self.initState,
self.defectNodes,self.defectAmp,self.N,self.rank)
plotStage.pop()
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#--------------------------- 1 particle CTQW on a line -------------------------
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
class Line(QuantumWalkP1):
""" Performs and analyses 1 particle continuous-time \
quantum walks on an infinite line
Args:
N (int) : an **even** number of nodes to initialize the walker with. Nodes are \
labeled :math:`j\\in\\{1-N/2,N/2\\}`.
Example:
To create a CTQW Line object for a 10 node line,
>>> walk = pyCTQW.MPI.Line(10)
"""
def __init__(self,N):
QuantumWalkP1.__init__(self,N)
def createH(self,d=None,amp=None):
""" Generate the Hamiltonian of the graph.
Args:
d (array of ints): an array containing *integers* indicating the nodes
where diagonal defects are to be placed (e.g. ``d=[0,1,4]``).
amp (array of floats): an array containing *floats* indicating the diagonal defect
amplitudes corresponding to each element in ``d`` (e.g. ``amp=[0.5,-1,4.2]``).
Returns:
: this creates a Hamiltonian matrix, accessed via the attribute
:attr:`Line.H`.
:rtype: :func:`pyCTQW.MPI.ctqw.Hamiltonian`
Warning:
The size of ``amp`` and ``d`` must be identical
>>> amp = [0.5,-1.,4.2]
>>> len(d) == len(amp)
True
"""
if (d and amp) is None:
self.defectNodes = [0]
self.defectAmp = [0.]
else:
self.defectNodes = d
self.defectAmp = amp
self.H.createLine(d=self.defectNodes,amp=self.defectAmp)
def createInitState(self,initState):
""" Generate the initial state of the quantum walk.
Args:
initState (array) : an :math:`n\\times 2` array, containing the initial \
state of the quantum walker in the format ``[[j1,amp1],[j2,amp2],...]``.
Returns:
: this creates a PETSc vector containing the initial state,
accessed via the attribute :attr:`Line.psi0`.
:rtype: :func:`petsc4py.PETSc.Vec`
Example:
For a CTQW initially located in a
superposition of nodes -4 and 2, e.g.
:math:`\\left|\\psi(0)\\right\\rangle = \\frac{1}{\\sqrt{2}}\\left|-4\\right\\rangle \
- \\frac{1}{\\sqrt{2}} \\left|2\\right\\rangle`,
the initial state would be created like so:
>>> import numpy as np
>>> init_state = [[-4,1./np.sqrt(2.)], [2,-1./np.sqrt(2.)]]
>>> walk.createInitState(init_state)
.. admonition:: Fortran interface
This function calls the Fortran function :f:func:`p1_init`.
"""
self.initState = initState
# create the inital stage
initStateS = _PETSc.Log.Stage('initState')
initStateS.push()
_ctqwmpi.p1_init(self.psi0.fortran,initState,self.N)
initStateS.pop()
def watch(self,nodes):
""" Creates a handle that watches node probability during propagation.
Args:
nodes (array of ints): the nodes to watch (e.g. ``[0,1,4]``).
Returns:
: creates a handle that can be
accessed to retrieve node probabilities for various :math:`t`
:rtype: :func:`pyCTQW.MPI.ctqw.nodeHandle`
Example:
>>> walk.watch([0,1,2,3,4])
>>> walk.propagate(5.,method='chebyshev')
>>> timeArray, probArray = walk.handle.getLocalNodes()
Warning:
Note that `walk.handle` attributes are **not** collective;
if running on multiple nodes, only *local* values will be
returned.
"""
nodes = [i+self.N/2-1 for i in nodes]
super(Line,self).watch(nodes)
def plot(self,filename):
""" Creates a plot of probability vs node.
Args:
filename (str): the absolute/relative path to the desired output file.
Note:
* ensure a file extension is present so that filetype is correctly set\
(choose one of png, pdf, ps, eps or svg).
* if :func:`propagate` has not been called, the probability of the \
initial state is plotted (:math:`t=0`). Otherwise, the last propagated \
time (``walk.t``) is used.
"""
if _os.path.isabs(filename):
outDir = _os.path.dirname(filename)
else:
outDir = './'+_os.path.dirname(filename)
# create output directory if it doesn't exist
try:
_os.mkdir(outDir)
except OSError as exception:
if exception.errno != _errno.EEXIST:
raise
plotStage = _PETSc.Log.Stage('Plotting'); plotStage.push()
_plot.plot(_np.arange(1-self.N/2,self.N/2+1),self.prob,filename,self.t,self.initState,
self.defectNodes,self.defectAmp,self.N,self.rank)
plotStage.pop()
def plotNodes(self,filename,t=None):
""" Creates a plot of the node probablities over time.
Args:
filename (str): the absolute/relative path to the desired output file.
.. important::
:func:`watch` **must** be called prior to propagation in order for\
the probabilities at the specified nodes to be stored.
Note:
* if multiple nodes are watched, they will **all** be plotted.
* ensure a file extension is present so that filetype is correctly set\
(choose one of png, pdf, ps, eps or svg).
* Timesteps are recorded when propagations occur - to increase the number\
of timesteps, increase the number of propagations.
"""
comm = _PETSc.COMM_WORLD
rank = comm.Get_rank()
node = self.handle.local_nodes
timeArray, probArray = self.handle.getLocalNodes(t=t)
probArray = comm.tompi4py().gather(probArray)
node = comm.tompi4py().gather(node)
if rank == 0:
timeArray = _np.array(timeArray)
nodeArray = _np.array([item-self.N/2+1 for sublist in node for item in sublist])
probArray = _np.array([item for sublist in probArray for item in sublist]).real
_plot.plotNodes(timeArray,nodeArray,probArray,filename)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#--------------------------- 2 particle CTQW on a line -------------------------
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
class Line2P(QuantumWalkP2):
""" Performs and analyses 2 particle continuous-time \
quantum walks on an infinite line
Args:
N (int) : an **even** number of nodes to initialize the walker with. Nodes are \
labeled :math:`j\\in\\{1-N/2,N/2\\}`.
Example:
To create a CTQW Line object for a 10 node line,
>>> walk = pyCTQW.MPI.Line2P(10)
"""
def __init__(self,N,d=None,amp=None,interaction=0.):
QuantumWalkP2.__init__(self,N)
if (d is not None) and (amp is not None):
self.defectNodes = d
self.defectAmp = amp
self.H.createLine2P(d=d,amp=amp,interaction=interaction)
def createH(self,d=None,amp=None,interaction=0.):
""" Generate the Hamiltonian of the graph.
Args:
d (array of ints): an array containing *integers* indicating the nodes
where diagonal defects are to be placed (e.g. ``d=[0,1,4]``).
amp (array of floats): an array containing *floats* indicating the diagonal defect
amplitudes corresponding to each element in ``d`` (e.g. ``amp=[0.5,-1,4.2]``).
interaction (float): the amplitude of interaction between the two walkers \
when located on the same vertex.
Returns:
: this creates a Hamiltonian matrix, accessed via the attribute
:attr:`Line2P.H`.
:rtype: :func:`pyCTQW.MPI.ctqw.Hamiltonian`
Warning:
The size of ``amp`` and ``d`` must be identical
>>> amp = [0.5,-1.,4.2]
>>> len(d) == len(amp)
True
"""
if (d and amp) is None:
self.defectNodes = [0]
self.defectAmp = [0.]
else:
self.defectNodes = d
self.defectAmp = amp
self.interaction = interaction
self.H.createLine2P(d=self.defectNodes,amp=self.defectAmp,interaction=self.interaction)
def createInitState(self,initState):
""" Generate the initial state of the quantum walk.
Args:
initState (array) : an :math:`n\\times 3` array, containing the initial \
state of the quantum walker in the format ``[[x1,y1,amp1],[x2,y2,amp2],...]``.
Returns:
: this creates a PETSc vector containing the initial state,
accessed via the attribute :attr:`Line2P.psi0`.
:rtype: :func:`petsc4py.PETSc.Vec`
Example:
For a CTQW initially located in state
:math:`\\left|\\psi(0)\\right\\rangle = \\frac{1}{\\sqrt{2}}\\left|0\\right\\rangle \
\otimes \\left|0\\right\\rangle \
- \\frac{i}{\\sqrt{2}} \\left|2\\right\\rangle\otimes \\left|2\\right\\rangle`,
the initial state would be created like so:
>>> import numpy as np
>>> init_state = [[0,0,1./np.sqrt(2.)], [2,2,-1j./np.sqrt(2.)]]
>>> walk.createInitState(init_state)
.. admonition:: Fortran interface
This function calls the Fortran function :f:func:`p2_init`.
"""
self.initState = initState
# create the inital stage
initStateS = _PETSc.Log.Stage('initState')
initStateS.push()
_ctqwmpi.p2_init(self.psi0.fortran,initState,self.N)
initStateS.pop()
def watch(self,nodes,watchtype='prob',**kwargs):
""" Creates a handle that watches either marginal probability \
or walker entanglement during propagation.
Args:
nodes (array of ints): the nodes to watch marginal probability (e.g. ``[0,1,4]``).
watchtype (str): (`'prob'` , `'entanglement'`).
the type of watch handle to produce.
Keyword Args:
: If ``watchtype='entanglement'``, EigSolver keywords can also be passed;
for more details of the available EigSolver properties,
see :func:`propagate`.
Returns:
: * if ``watchtype='prob'``, creates a handle that can be \
accessed to retrieve marginal node probabilities for various :math:`t`
* if ``watchtype='entanglment'``, creates a handle that can be \
accessed to retrieve entanglement values for various :math:`t`
:rtype: * if ``watchtype='prob'``: :func:`pyCTQW.MPI.ctqw.nodeHandle`
* if ``watchtype='entanglement'``: :func:`pyCTQW.MPI.ctqw.entanglementHandle`
Examples:
To watch the entanglement,
>>> walk.watch(None,watchtype='entanglement')
>>> walk.propagate(5.,method='chebyshev')
>>> timeArray, entArray = walk.entanglementHandle.getEntanglement()
.. note::
* The entanglement measure used in Von Neumann entropy, calculated \
via :math:`S=-\\sum_{i}\\lambda_i\log_{2}\\lambda_i`, where :math:`\lambda_i` \
are the eigenvalues of the reduced density matrix \
:math:`\\rho_2 = \\text{Tr}_1(|\\psi(t)\\rangle\\langle\\psi(t)|)`
* Nodes do not need to be specified, as entanglement is a global measurement.
* As it is a global measurement, there is a large amount of node communication \
which may increase overall program run time.
To watch the probabilities,
>>> walk.watch([0,1,4])
>>> walk.propagate(2.,method='chebyshev')
>>>timeArray, probXArray, probYArray = self.handle.getLocalNode(0,p=2)
.. warning::
Note that `walk.handle` attributes are **not** collective;
if running on multiple nodes, only *local* values will be
returned.
"""
if watchtype == 'prob':
nodes = [i+self.N/2-1 for i in nodes]
super(Line2P,self).watch(nodes,watchtype=watchtype,**kwargs)
def plot(self,filename):
""" Creates a plot of probability vs node.
Args:
filename (str): the absolute/relative path to the desired output file.
Note:
* ensure a file extension is present so that filetype is correctly set\
(choose one of png, pdf, ps, eps or svg).
* if :func:`propagate` has not been called, the probability of the \
initial state is plotted (:math:`t=0`). Otherwise, the last propagated \
time (:attr:`Line2P.t`) is used.
"""
if _os.path.isabs(filename):
outDir = _os.path.dirname(filename)
else:
outDir = './'+_os.path.dirname(filename)
# create output directory if it doesn't exist
try:
_os.mkdir(outDir)
except OSError as exception:
if exception.errno != _errno.EEXIST:
raise
plotStage = _PETSc.Log.Stage('Plotting'); plotStage.push()
_plot.plot2P(_np.arange(1-self.N/2,self.N/2+1),self.psiX,self.psiY,filename,self.t,self.initState,
self.defectNodes,self.defectAmp,self.N,self.rank)
plotStage.pop()
def plotNode(self,filename,node,t=None):
""" Creates a plot of the marginal probablities on a *specified node* over time.
Args:
filename (str): the absolute/relative path to the desired output file.
node (int): the node to plot.
.. important::
:func:`watch` **must** be called prior to propagation in order for\
the probabilities at the specified nodes to be stored.
Note:
* ensure a file extension is present so that filetype is correctly set\
(choose one of png, pdf, ps, eps or svg).
* Timesteps are recorded when propagations occur - to increase the number\
of timesteps, increase the number of propagations.
* See :func:`plotNodes` to plot multiple nodes.
"""
node = node+self.N/2-1
super(Line2P,self).plotNode(filename,node)
def plotNodes(self,filename,p=1,t=None):
""" Creates a plot of the node probablities over time.
Args:
filename (str): the absolute/relative path to the desired output file.
p (int): (1|2) - choose whether to plot the marginal probability of particle 1 or 2.
.. important::
:func:`watch` **must** be called prior to propagation in order for\
the probabilities at the specified nodes to be stored.
Note:
* if multiple nodes are watched, they will **all** be plotted.
* if you wish to plot marginal probabilities for both particle \
1 and 2, see :func:`plotNode`.
* ensure a file extension is present so that filetype is correctly set\
(choose one of png, pdf, ps, eps or svg).
* Timesteps are recorded when propagations occur - to increase the number\
of timesteps, increase the number of propagations.
"""
comm = _PETSc.COMM_WORLD
rank = comm.Get_rank()
node = self.handle.local_nodes
timeArray, probXArray, probYArray = self.handle.getLocalNodes(t=t,p=2)
if p == 1:
probArray = probXArray
elif p==2:
probArray = probYArray
else:
print 'p must be either 1 or 2'
return
probArray = comm.tompi4py().gather(probArray)
node = comm.tompi4py().gather(node)
if rank == 0:
timeArray = _np.array(timeArray)
nodeArray = _np.array([item-self.N/2+1 for sublist in node for item in sublist])
probArray = _np.array([item for sublist in probArray for item in sublist]).real
_plot.plotNodes(timeArray,nodeArray,probArray,filename,p=p)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#--------------------------- 3 particle CTQW on a line -------------------------
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
class Line3P(QuantumWalkP3):
""" Performs and analyses 3 particle continuous-time \
quantum walks on an infinite line
Args:
N (int) : an **even** number of nodes to initialize the walker with. Nodes are \
labeled :math:`j\\in\\{1-N/2,N/2\\}`.
Example:
To create a CTQW Line object for a 10 node line,
>>> walk = pyCTQW.MPI.Line3P(10)
"""
def __init__(self,N,d=None,amp=None,interaction=0.):
QuantumWalkP3.__init__(self,N)
if (d is not None) and (amp is not None):
self.defectNodes = d
self.defectAmp = amp
self.H.createLine3P(d=d,amp=amp,interaction=interaction)
def createH(self,d=None,amp=None,interaction=0.):
""" Generate the Hamiltonian of the graph.
Args:
d (array of ints): an array containing *integers* indicating the nodes
where diagonal defects are to be placed (e.g. ``d=[0,1,4]``).
amp (array of floats): an array containing *floats* indicating the diagonal defect
amplitudes corresponding to each element in ``d`` (e.g. ``amp=[0.5,-1,4.2]``).
interaction (float): the amplitude of interaction between the two walkers \
when located on the same vertex.
Returns:
: this creates a Hamiltonian matrix, accessed via the attribute
:attr:`Line3P.H`.
:rtype: :func:`pyCTQW.MPI.ctqw.Hamiltonian`
Warning:
The size of ``amp`` and ``d`` must be identical
>>> amp = [0.5,-1.,4.2]
>>> len(d) == len(amp)
True
"""
if (d and amp) is None:
self.defectNodes = [0]
self.defectAmp = [0.]
else:
self.defectNodes = d
self.defectAmp = amp
self.interaction = interaction
self.H.createLine3P(d=self.defectNodes,amp=self.defectAmp,interaction=self.interaction)
def createInitState(self,initState):
""" Generate the initial state of the quantum walk.
Args:
initState (array) : an :math:`n\\times 4` array, containing the initial \
state of the quantum walker in the format ``[[x1,y1,z1,amp1],[x2,y2,z2,amp2],...]``.
Returns:
: this creates a PETSc vector containing the initial state,
accessed via the attribute :attr:`Line3P.ps0`.
:rtype: :func:`petsc4py.PETSc.Vec`
Example:
For a CTQW initially located in state
:math:`\\left|\\psi(0)\\right\\rangle = \\frac{1}{\\sqrt{2}}\\left|-4\\right\\rangle \
\\otimes\\left|1\\right\\rangle \\otimes\\left|5\\right\\rangle \
- \\frac{i}{\\sqrt{2}} \\left|1\\right\\rangle\\otimes\\left|0\\right\\rangle \
\\otimes\\left|2\\right\\rangle`,
the initial state would be created like so:
>>> import numpy as np
>>> init_state = [[-4,1,5,1./np.sqrt(2.)], [1,0,2,-1j./np.sqrt(2.)]]
>>> walk.createInitState(init_state)
.. admonition:: Fortran interface
This function calls the Fortran function :f:func:`p3_init`.
"""
self.initState = _np.vstack( [ _np.array(initState).T[0]+self.N/2-1,
_np.array(initState).T[1]+self.N/2-1,
_np.array(initState).T[2]+self.N/2-1,
_np.array(initState).T[3] ]).T.tolist()
# create the inital stage
initStateS = _PETSc.Log.Stage('initState')
initStateS.push()
_ctqwmpi.p3_init(self.psi0.fortran,self.initState,self.N)
initStateS.pop()
def watch(self,nodes,type='prob'):
""" Creates a handle that watches either marginal probability \
or walker entanglement during propagation.
Args:
nodes (array of ints): the nodes to watch marginal probability (e.g. ``[0,1,4]``).
watchtype (str): (`'prob'` , `'entanglement'`).
the type of watch handle to produce.
Keyword Args:
: If ``watchtype='entanglement'``, EigSolver keywords can also be passed;
for more details of the available EigSolver properties,
see :func:`propagate`.
Returns:
: * if ``watchtype='prob'``, creates a handle that can be \
accessed to retrieve marginal node probabilities for various :math:`t`
* if ``watchtype='entanglment'``, creates a handle that can be \
accessed to retrieve entanglement values for various :math:`t`
:rtype: * if ``watchtype='prob'``: :func:`pyCTQW.MPI.ctqw.nodeHandle`
* if ``watchtype='entanglement'``: :func:`pyCTQW.MPI.ctqw.entanglementHandle`
Examples:
To watch the entanglement,
>>> walk.watch(None,watchtype='entanglement')
>>> walk.propagate(5.,method='chebyshev')
>>> timeArray, entArray = walk.entanglementHandle.getEntanglement()
.. note::
* The entanglement measure used in Von Neumann entropy, calculated \
via :math:`S=-\\sum_{i}\\lambda_i\log_{2}\\lambda_i`, where :math:`\lambda_i` \
are the eigenvalues of the reduced density matrix \
:math:`\\rho_2 = \\text{Tr}_1(|\\psi(t)\\rangle\\langle\\psi(t)|)`
* Nodes do not need to be specified, as entanglement is a global measurement.
* As it is a global measurement, there is a large amount of node communication \
which may increase overall program run time.
To watch the probabilities,
>>> walk.watch([0,1,4])
>>> walk.propagate(2.,method='chebyshev')
>>>timeArray, probXArray, probYArray, probZArray = self.handle.getLocalNode(4,p=2)
.. warning::
Note that `walk.handle` attributes are **not** collective;
if running on multiple nodes, only *local* values will be
returned.
"""
nodes = [i+self.N/2-1 for i in nodes]
super(Line3P,self).watch(nodes,type=type)
def plot(self,filename):
""" Creates a plot of probability vs node.
Args:
filename (str): the absolute/relative path to the desired output file.
Note:
* ensure a file extension is present so that filetype is correctly set\
(choose one of png, pdf, ps, eps or svg).
* if :func:`propagate` has not been called, the probability of the \
initial state is plotted (:math:`t=0`). Otherwise, the last propagated \
time (:attr:`Line3P.t`) is used.
"""
if _os.path.isabs(filename):
outDir = _os.path.dirname(filename)
else:
outDir = './'+_os.path.dirname(filename)
# create output directory if it doesn't exist
try:
_os.mkdir(outDir)
except OSError as exception:
if exception.errno != _errno.EEXIST:
raise
initstateLabels = []
for i in range(len(self.initState)):
initstateLabels.append([sum(pair) for pair in zip(self.initState[i], [1-self.N/2,1-self.N/2,1-self.N/2,0])])
plotStage = _PETSc.Log.Stage('Plotting'); plotStage.push()
_plot.plot3P(_np.arange(1-self.N/2,self.N/2+1),self.psiX,self.psiY,self.psiZ,filename,self.t,initstateLabels,
self.defectNodes,self.defectAmp,self.N,self.rank)
plotStage.pop()
def plotNode(self,filename,node,t=None):
""" Creates a plot of the marginal probablities on a *specified node* over time.
Args:
filename (str): the absolute/relative path to the desired output file.
node (int): the node to plot.
.. important::
:func:`watch` **must** be called prior to propagation in order for\
the probabilities at the specified nodes to be stored.
Note:
* ensure a file extension is present so that filetype is correctly set\
(choose one of png, pdf, ps, eps or svg).
* Timesteps are recorded when propagations occur - to increase the number\
of timesteps, increase the number of propagations.
* See :func:`plotNodes` to plot multiple nodes.
"""
node = node+self.N/2-1
super(Line3P,self).plotNode(filename,node)
def plotNodes(self,filename,p=1,t=None):
""" Creates a plot of the node probablities over time.
Args:
filename (str): the absolute/relative path to the desired output file.
p (int): (1|2|3) - choose whether to plot the marginal probability of particle 1, 2 or 3.
.. important::
:func:`watch` **must** be called prior to propagation in order for\
the probabilities at the specified nodes to be stored.
Note:
* if multiple nodes are watched, they will **all** be plotted.
* if you wish to plot marginal probabilities for all particles, see :func:`plotNode`.
* ensure a file extension is present so that filetype is correctly set\
(choose one of png, pdf, ps, eps or svg).
* Timesteps are recorded when propagations occur - to increase the number\
of timesteps, increase the number of propagations.
"""
comm = _PETSc.COMM_WORLD
rank = comm.Get_rank()
node = self.handle.local_nodes
timeArray, probXArray, probYArray, probZArray = self.handle.getLocalNodes(t=t,p=3)
if p == 1:
probArray = probXArray
elif p==2:
probArray = probYArray
elif p==3:
probArray = probZArray
else:
print 'p must be either 1, 2 or 3'
return
probArray = comm.tompi4py().gather(probArray)
node = comm.tompi4py().gather(node)
if rank == 0:
timeArray = _np.array(timeArray)
nodeArray = _np.array([item-self.N/2+1 for sublist in node for item in sublist])
probArray = _np.array([item for sublist in probArray for item in sublist]).real
_plot.plotNodes(timeArray,nodeArray,probArray,filename,p=p)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#----------------------------------- Graph Isomorphism -------------------------
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
class IncorrectInput(Exception):
pass
class GraphISO(object):
""" A graph isomorphism solver, containing functions for
creating graph certificates and checking isomorphism
of adjacency matrices.
>>> gi = pyCTQW.MPI.GraphISO(p=2,propagator='krylov')
Args:
p (int): number of particles, 2 *(default)* or 3,
to use in constructing the graph certificate.
freqTol (float): the tolerance to use when constructing the
frequency table (*default* ``1.e-2``). See also
:py:func:`GIcert`.
compareTol (float): the tolerance used when comparing two Graph
certificates (*default* ``1.e-10``). See also
:py:func:`isomorphicQ`.
propagator (str): the CTQW propagator algorithm to use
when calculating the graph certificate
(``'chebyshev'`` *(default)* or ``'krylov'``).
bosonic (bool): if ``True``, instructs the solver to use a
Bosonic Hamiltonian (*default* ``False``).
Note:
* For ``freqTol=1.e-2``, probabilities within (100freqTol)% \
are treated as equal.
* Two isomorphic certificates satisfy \
:math:`\max(|cert_1 - cert_2|) < compareTol`.
"""
def __init__(self,**kwargs):
self.__default = {
'p' : 2,
'freqTol' : 5.e-3,
'compareTol' : 1.e-10,
'propagator' : 'chebyshev'
}
self.__eigDefault = {
'esolver' : 'krylovschur',
'emax_estimate' : 0.,
'workType' : 'null',
'workSize' : '35',
'tol' : 0.,
'maxIt' : 0,
'verbose' : False,
'bosonic' : False
}
self.__rank = _PETSc.Comm.Get_rank(_PETSc.COMM_WORLD)
for key, default in self.__default.iteritems():
setattr(self, key, kwargs.get(key,default))
for key, default in self.__eigDefault.iteritems():
setattr(self, key, kwargs.get(key,default))
def setProperties(self,**kwargs):
"""
Set some or all of the GraphISO properties.
For the list of the properties see :class:`GraphISO`.
"""
for key in kwargs:
if self.__default.has_key(key):
setattr(self, key, kwargs.get(key))
else:
print 'Property {} does not exist!'.format(key)
def getProperty(self,*args):
"""
Get some or all of the GraphISO properties.
For the list of the properties see :class:`GraphISO`.
"""
for key in args:
if self.__default.has_key(key):
return getattr(self, key)
else:
print 'Property {} does not exist!'.format(key)
def setEigSolver(self,**kwargs):
"""Set some or all of the eigenvalue solver properties.
Keyword Args:
esolver (str): the eigensolver algorithm to use.
* ``'krylovschur'`` *(default)* - Krylov-Schur
* ``'arnoldi'`` - Arnoldi Method
* ``'lanczos'`` - Lanczos Method
* ``'power'`` - Power Iteration/Rayleigh Quotient Iteration
* ``'gd'`` - Generalized Davidson
* ``'jd'`` - Jacobi-Davidson,
* ``'lapack'`` - Uses LAPACK eigenvalue solver routines
* ``'arpack'`` - *only available if SLEPc is\
compiled with ARPACK linking*
workType (str): can be used to set the eigensolver worktype
(either ``'ncv'`` or ``'mpd'``). The default
is to let SLEPc decide.
workSize (int): sets the work size **if** ``workType`` is set.
tolIn (float): tolerance of the eigenvalue solver
(*default* ``0.`` (SLEPc decides)).
maxIt (int): maximum number of iterations of the eigenvalue solver
(*default* ``0`` (SLEPc decides)).
verbose (bool): if ``True``, writes eigensolver information to the console
emax_estimate (float): used to override the calculation
of the graphs maximum eigenvalue.
.. caution::
* If supplied, the value of :attr:`emax_estimate`:math:`\hat{\lambda}_{\max}` \
**must** satisfy :math:`\hat{\lambda}_{\max}\geq\lambda_{\max}`, \
where :math:`\lambda_{\max}` is the actual maximum eigenvalue of the graph.
* The greater the value of :math:`\hat{\lambda}_{\max} \
-\lambda_{\max}`, the longer the convergence \
time of the ``chebyshev`` propagator.
Note:
* These properties only apply if ``propagator='chebyshev'``
* For more information regarding these properties,refer to \
the `SLEPc documentation\
<http://www.grycap.upv.es/slepc/documentation/manual.htm>`_
"""
for key in kwargs:
if self.__eigDefault.has_key(key):
setattr(self, key, kwargs.get(key))
else:
print 'Eigsolver property {} does not exist!'.format(key)
def getEigSolver(self,*args):
"""
Get some or all of the GraphISO properties.
For the list of the properties see :func:`setEigSolver`.
"""
for key in args:
if self.__eigDefault.has_key(key):
return getattr(self, key)
else:
print 'Eigsolver property {} does not exist!'.format(key)
def GIcert(self,adj):
"""Generate the GI certificate of a graph.
Args:
adj (array or numpy.array): symmetric adjacency matrix in the form
of a dense array of numpy array.
Returns:
graph certificate
:rtype: :func:`numpy.array`
.. admonition:: Fortran interface
This function calls the Fortran function :f:func:`GraphISCert`.
"""
GIcertS = _PETSc.Log.Stage('GICert')
GIcertS.push()
cert, certSize = _ctqwmpi.GraphISCert(adj,self.p,self.freqTol,self.propagator,
self.esolver,self.emax_estimate,self.workType,self.workSize,self.tol,self.maxIt,self.bosonic,self.verbose)
GIcertS.pop()
return _np.array(cert).T[_np.lexsort(_np.array(cert)[:,0:certSize])[::-1]]
def isomorphicQ(self, G1, G2, saveCert1=None, saveCert2=None):
"""Returns the GI certificates in addition to
``True`` if two graphs are isomorphic, ``False`` otherwise.
Args:
G(1|2) (array or numpy.array): \
symmetric adjacency matrices in the form
of a dense array or numpy array, **OR** a
previously calculated two column GI certificate.
Keyword Args:
saveCert(1|2) (str): \
Filename. If provided, the specified GI
certificate is exported to a text file.
Returns:
both graph certificates, boolean value
:rtype: :func:`numpy.array`, *bool*
"""
if _np.all(_np.array(G1).T == _np.array(G1)):
GIcert1 = self.GIcert(G1)
elif _np.array(G1).shape[1] == 2:
GIcert1 = G1
else:
raise IncorrectInput('Argument G1 is not a valid adjacency matrix or GI certificate')
if saveCert1 is not None:
_io.createPath(saveCert1)
_np.savetxt(saveCert1, GIcert1)
if _np.all(_np.array(G2).T == _np.array(G2)):
GIcert2 = self.GIcert(G2)
elif _np.array(G2).shape[1] == 2:
GIcert2 = G2
else:
raise IncorrectInput('Argument G2 is not a valid adjacency matrix or GI certificate')
if saveCert2 is not None:
_io.createPath(saveCert2)
_np.savetxt(saveCert2, GIcert2)
comm = _PETSc.COMM_WORLD
size = comm.getSize()
result = None
if self.__rank == 0:
if GIcert1.shape[0] == GIcert2.shape[0]:
if _np.abs(_np.subtract(GIcert1, GIcert2)).max() < self.compareTol:
result = _np.array([True for i in range(size)])
else:
result = _np.array([False for i in range(size)])
else:
result = _np.array([False for i in range(size)])
result = comm.tompi4py().scatter(result)
return GIcert1, GIcert2, result
def AllIsomorphicQ(self,folder,graphRange=None,info=True,checkSelf=False):
"""Calculates whether each pair of graphs (in a specified set of graphs)
are isomorphic, returning an array :math:`R` with :math:`R_{ij} = 1` if
graphs :math:`i` and :math:`j` are isomorphic, and :math:`R_{ij} = 0` otherwise.
Args:
folder (str): path to a folder containing a collection of adjacency
matrices in dense text file format.
graphRange (array): an array containing graph numbers to test.
By default, all graphs in a folder are tested.
info (bool): if ``True``, information on each :math:`R_{ij}` comparison
is printed to the console.
checkSelf (bool): if ``True``, each graph is also tested against itself.
:rtype: *array of ints*
Note:
The text files must have filenames of the form \*X.txt
where X represents a number (of any number of digits).
These are used to order the graphs.
"""
if not _os.path.isdir(folder):
if self.__rank == 0:
print 'Directory does not exist!'
return
numbers = _glob.re.compile(r'(\d+)')
def numericalSort(value):
parts = numbers.split(value)
parts[1::2] = map(int, parts[1::2])
return parts
filelist = sorted(_glob.glob(folder+'/*[0-9].txt'), key=numericalSort)
adj = []
for graph in filelist:
if graphRange is not None:
if int(_glob.re.findall(r'\d+', graph)[-1]) in graphRange:
adj.append(_np.genfromtxt(graph))
if (self.__rank == 0 and info):
print 'Adding graph ' + graph
else:
adj.append(_np.genfromtxt(graph))
if (self.__rank == 0 and info):
print 'Adding graph ' + graph
NG = len(adj)
comparisonTable = _np.zeros([NG,NG])
cert = [None]*NG
for i in range(NG):
for j in range(i if checkSelf else i+1, NG):
if (self.__rank == 0 and info):
print 'Testing graphs ' + str(i) + ',' + str(j)
if cert[i] is not None and cert[j] is not None:
cert[i], cert[j], checkISO = self.isomorphicQ(cert[i], cert[j])
elif cert[i] is not None:
cert[i], cert[j], checkISO = self.isomorphicQ(cert[i], adj[j])
elif cert[j] is not None:
cert[i], cert[j], checkISO = self.isomorphicQ(adj[i], cert[j])
else:
cert[i], cert[j], checkISO = self.isomorphicQ(adj[i], adj[j])
comparisonTable[i,j] = 1 if checkISO else 0
if (self.__rank == 0 and info):
print '\tIsomorphic' if comparisonTable[i,j] == 1 else '\tNon-isomorphic'
if checkSelf:
return comparisonTable + comparisonTable.T - _np.diag(comparisonTable.diagonal())
else:
return comparisonTable + comparisonTable.T + _np.diag(_np.ones(NG))