Scripting with Bonsu¶
Core concepts¶
Bonsu provides an object-oriented scripting interface to each algorithm for ease and interoperability with additional Python modules. To begin using this interface, import the phasing module as follows:
>>> from bonsu import phasing
Once imported, the phasing module provides access to all phase retrieval algorithms. For a complete list see Library Reference – Phasing Algorithms.
Each algorithm will take a number of NumPy arrays as input. Input data is loaded and prepared in wrap around order. It is therefore also necessary to utilise the Wrap Data function. Note that the data must be stored in a “C”-ordered double precision complex Numpy array:
>>> expdata = numpy.array(numpy.load("expdata.npy"), dtype=numpy.cdouble, copy=True, order='C')
>>> numpy.sqrt(expdata, expdata)
>>> expdata[:] = phasing.WrapArray(expdata)
The data mask
in wrap around order is also needed:
>>> mask = numpy.load("mask.npy")
>>> mask[:] = phasing.WrapArray(mask)
An array to contain the reconstructed object is also needed and can be created in a number of ways, with data initialised as needed, i.e. random data or initial starting data from a previous reconstruction.
>>> seqdata = numpy.empty_like(expdata)
>>> seqdata[:] = 1.0 # initialise as 1.0.
To utilise the HIO Mask algorithm, create an instance as follows:
>>> hio = phasing.HIOMask()
Then initialise the remaining parameters for 100 iterations:
>>> hio.SetStartiter(0)
>>> hio.SetNumiter(100)
>>> hio.SetNumthreads(3)
>>> hio.SetBeta(0.9)
>>> hio.SetExpdata(expdata)
>>> hio.SetSupport(support)
>>> hio.SetMask(mask)
>>> hio.SetSeqdata(seqdata)
To begin reconstruction we first prepare the algorithm as follows:
>>> hio.Prepare()
Execution then follows:
>>> hio.Start()
In a threaded or non-blocking use of the algorithm, it is possible the pause, resume and stop the reconstruction using the following methods:
>>> hio.Pause()
>>> hio.Resume()
>>> hio.Stop()
The output of the reconstruction is saved as follows:
>>> numpy.save("output.npy", seqdata)
Each algorithm can use Shrink Wrap
to modify the support during reconstruction. Prefix the algorithm name with SW
to utilise the Shrink Wrap
version of the algorithm. For example, to use the Shrink Wrap
version of the HIO Mask algorithm, import and instantiate the following:
>>> swhio = phasing.SWHIOMask()
Then initialise the remaining parameters as follows:
>>> swhio.SetStartiter(0)
>>> swhio.SetNumiter(100)
>>> swhio.SetNumthreads(3)
>>> swhio.SetExpdata(expdata)
>>> swhio.SetSupport(support)
>>> swhio.SetMask(mask)
>>> swhio.SetSeqdata(seqdata)
>>> swhio.Prepare()
Methods unqiue to the Shrink Wrap
implementation that require initialisation are as follows:
>>> swhio.SetSWCyclelength(20)
>>> swhio.SetSWThreshold(0.2)
>>> swhio.SetSWSigma(3.0)
To begin reconstruction we also first prepare the Shrink Wrap
algorithm as follows:
>>> swhio.SWPrepare()
Execution then follows:
>>> swhio.SWStart()
Custom Algorithms¶
Bonsu provides tools to build custom projection algorithms. To build a custom algorithm, create a script that imports an abstract algorithm class and subclasses the DoIter
and RSCons
methods. DoIter
is called at each iteration of the algorithm. It will apply the real-space and Fourier-space constraints. RSCons
is a method called by DoIter
to apply the real-space constraint.
from bonsu.phasing.abstract import PhaseAbstract
class MyAlgorithm(PhaseAbstract):
"""
My custom HIO algorithm
"""
def DoIter(self):
self.rho_m1[:] = self.seqdata[:]
fftw_stride(self.seqdata,self.seqdata,self.plan,FFTW_TORECIP,1)
self.SetRes()
self.SetAmplitudes()
fftw_stride(self.seqdata,self.seqdata,self.plan,FFTW_TOREAL,1)
n = self.citer_flow[0]
self.residual[n] = self.res/self.sos
amp = numpy.absolute(self.seqdata)
sos1 = numpy.sum(amp*amp)
self.RSCons()
amp = numpy.absolute(self.seqdata)
sos2 = numpy.sum(amp*amp)
norm = sqrt(sos1/sos2)
self.seqdata[:] = norm*self.seqdata[:]
self.citer_flow[0] += 1
def RSCons(self):
realsupport = numpy.real(self.support)
nosupport = realsupport < 0.5
self.seqdata[nosupport] = self.rho_m1[nosupport] - self.seqdata[nosupport] * self.beta
The custom algorithm is then instantiated as before:
myalg = MyAlgorithm()
Initialisation of parameters proceeds as before. To begin reconstruction prepare the algorithm as follows:
myalg.PrepareCustom()
Execution then follows as before:
myalg.Start()
Concurrent Phase Retrieval¶
Bonsu provides tools to reconstruct atomic displacement field information (and subsequently strain field information) via concurrent phase reconstruction using diffraction pattern data that is obtained from multiple Bragg reflections of a single crystal. The method used is described in detail in the following publication: Newton, Phys. Rev. B 102, 014104 (2020) .
The concurrent phase retrieval method can in principle utilise any projection algorithm implemented using the Custom Algorithm interface. The following example illustrates how to use utilise the method using the standard Abstract Algorithm class (HIO). Please note that diffraction pattern data should undergo co-ordinate transformation and interpolation onto a regular grid prior to commencement of the reconstruction process. Please also note that no checking of the data type (numpy.cdouble) and consistency in the array dimensions is performed.
from bonsu.phasing.concurrent import PhaseConcurrent
from bonsu.phasing.abstract import PhaseAbstract
from bonsu.phasing.wrap import WrapArray
# Specify number of diffraction patterns / Q-vectors
NQ = 4
# Instantiate class object
d = PhaseConcurrent()
# Specify class object of phase retrieval algorithm used for reconstruction
# Define your own and use here, if preferred.
d.SetBaseClass(PhaseAbstract)
# Set number of concurrent datasets / Q-vectors
d.SetNQ(NQ)
# Set Q-vectors
d.SetQ(0, [1,0,0])
d.SetQ(1, [1,1,0])
d.SetQ(2, [0,1,0])
d.SetQ(3, [0,0,1])
# Define names of files
expnames = ["diffamp0.npy","diffamp1.npy","diffamp2.npy","diffamp3.npy"]
masknames = ["mask0.npy","mask1.npy","mask2.npy","mask3.npy"]
supportname = "support.npy"
# Use names above to load arrays into list of length NQ
amps = [WrapArray(numpy.load(name)) for name in expnames]
masks = [WrapArray(numpy.load(name)) for name in masknames]
supports = [numpy.load(supportname)]*NQ
seqdatas = [numpy.ones_like(amp) for amp in amps]
# Class methods below are the same as those in the phase retrieval algorithm
# class but will instead take a list of length NQ.
d.SetExpdata(amps)
d.SetSupport(supports)
d.SetMask(masks)
d.SetSeqdata(seqdatas)
d.SetBeta([0.9]*NQ)
# Set iterations
d.SetStartiter(0)
d.SetNumiter(100)
# Number of FFTW threads for each concurrent reconstruction.
d.SetNumthreads(1)
# Prepare
d.Prepare()
# Start
d.Start()
# Save reconstructed phase for each Q-vector
# Filenames are automatically indexed.
d.SaveSequences("output.npy")
# Save displacement field
d.SaveDField("dfield.npy")
Class Reference¶
Abstract Classes¶
- class bonsu.phasing.abstract.PhaseAbstract(parent=None)¶
Bases:
object
Phasing Base Class
- PrepareCustom()¶
Prepare algorithm using FFTW python interface.
- SetSOS()¶
Set Sum of Squares.
- GetSOS()¶
Get Sum of Squares.
- SetStartiter(startiter)¶
Set the starting iteration number.
- SetNumiter(numiter)¶
Set the number of iterations of the algorithm to perform.
- SetNumthreads(nthreads)¶
Set the number of FFTW threads.
- SetSeqdata(seqdata)¶
Set the reconstruction data array.
- GetSeqdata()¶
Get the reconstruction data array.
- SetExpdata(expdata)¶
Set the raw experimental amplitude data array.
- GetExpdata()¶
Get the raw experimental amplitude data array.
- SetSupport(support)¶
Set the support array.
- GetSupport()¶
Get the support array.
- SetMask(mask)¶
Set the mask data array.
- GetMask()¶
Set the mask data array.
- SetBeta(beta)¶
Set beta relaxation parameter.
- GetBeta()¶
Get beta relaxation parameter.
- Pause()¶
Pause the reconstruction process.
- Resume()¶
Resume the reconstruction process.
- Start()¶
Start the reconstruction process.
- Stop()¶
Stop the reconstruction process.
- class bonsu.phasing.abstract.PhaseAbstractPC(parent=None)¶
Bases:
PhaseAbstract
Phasing Partial Coherence Base Class
- SetPSF(psf)¶
Set point-spread function from data array.
- GetPSF()¶
Get point-spread function from data array.
- LorentzFillPSF()¶
Create a point-spread function with Lorentz distribution. This will use the HWHM specified using SetGammaHWHM.
- SetNumiterRLpre(niterrlpre)¶
Set the number of iterations before RL optimisation.
- GetNumiterRLpre()¶
Get the number of iterations before RL optimisation.
- SetNumiterRL(niterrl)¶
Set the number of RL iterations.
- GetNumiterRL()¶
Get the number of RL iterations.
- SetNumiterRLinterval(niterrlinterval)¶
Set number of iterations between each RL optimisation cycle.
- GetNumiterRLinterval()¶
Get number of iterations between each RL optimisation cycle.
- SetGammaHWHM(gammaHWHM)¶
Set the half-width half-maximum of the Lorentz distribution. This is utilised in the LorentzFillPSF method.
- GetGammaHWHM()¶
Get the half-width half-maximum of the Lorentz distribution.
- SetPSFZeroEnd(ze)¶
Voxels upto a distance of [i,j,k] from the perimeter of the PSF array are set to zero. This can improve stability of the algorithm.
- GetPSFZeroEnd()¶
Get zero voxels perimeter size.
- class bonsu.phasing.ShrinkWrap.ShrinkWrap¶
Shrink Wrap algorithm.
- SetSWCyclelength(cycle)¶
Set interval of iterations after which the support is updated.
- GetSWCyclelength()¶
Get interval of iterations after which the support is updated.
- SetSWSigma(sigma)¶
Set standard deviation of the Gaussian smoothing function for the support.
- GetSWSigma()¶
Get standard deviation of the Gaussian smoothing function for the support.
- SetSWThreshold(thresh)¶
Set fractional value below which sequence data is not used when creating the new support.
- GetSWThreshold()¶
Get fractional value below which sequence data is not used when creating the new support.
- SWPrepare()¶
Prepare shrink wrap algorithm.
- SWStart()¶
Start the shrink wrap reconstruction process.
Hybrid Input-Output Classes¶
- class bonsu.phasing.HIO.HIO(parent=None)¶
Bases:
PhaseAbstract
Fienup’s hybrid input-output (HIO) algorithm.
- class bonsu.phasing.HIO.SWHIO¶
Bases:
HIO
,ShrinkWrap
Fienup’s hybrid input-output (HIO) algorithm. Shrink wrapped.
- class bonsu.phasing.HIO.HIOMask(parent=None)¶
Bases:
PhaseAbstract
Fienup’s hybrid input-output (HIO) algorithm with the addition of a Fourier space constraint mask.
- class bonsu.phasing.HIO.SWHIOMask¶
Bases:
HIOMask
,ShrinkWrap
Fienup’s hybrid input-output (HIO) algorithm with the addition of a Fourier space constraint mask. Shrink wrapped.
- class bonsu.phasing.HIO.HIOPlus(parent=None)¶
Bases:
PhaseAbstract
Fienup’s hybrid input-output (HIO) algorithm with non-negativity constraint and with the addition of a Fourier space constraint mask.
- class bonsu.phasing.HIO.SWHIOPlus¶
Bases:
HIOPlus
,ShrinkWrap
Fienup’s hybrid input-output (HIO) algorithm with non-negativity constraint and with the addition of a Fourier space constraint mask. Shrink wrapped.
- class bonsu.phasing.HIO.PCHIO(parent=None)¶
Bases:
PhaseAbstract
Fienup’s hybrid input-output (HIO) algorithm with phase constraint and with the addition of a Fourier space constraint mask.
- SetMaxphase(max)¶
Set phase maximum.
- GetMaxphase()¶
Get phase maximum.
- SetMinphase(min)¶
Set phase minimum.
- GetMinphase()¶
Get phase minimum.
- class bonsu.phasing.HIO.SWPCHIO¶
Bases:
PCHIO
,ShrinkWrap
Fienup’s hybrid input-output (HIO) algorithm with phase constraint and with the addition of a Fourier space constraint mask. Shrink wrapped.
- class bonsu.phasing.HIO.PGCHIO(parent=None)¶
Bases:
PhaseAbstract
Fienup’s hybrid input-output (HIO) algorithm with phase gradient constraint in the Q-vector direction with the addition of a Fourier space constraint mask.
- SetMaxphase(max)¶
Set phase maximum.
- GetMaxphase()¶
Get phase maximum.
- SetMinphase(min)¶
Set phase minimum.
- GetMinphase()¶
Get phase minimum.
- SetQ(q)¶
Set Q-vector tuple
- GetQ()¶
Get Q-vector tuple
- class bonsu.phasing.HIO.SWPGCHIO¶
Bases:
PGCHIO
,ShrinkWrap
Fienup’s hybrid input-output (HIO) algorithm with phase gradient constraint in the Q-vector direction with the addition of a Fourier space constraint mask. Shrink wrapped.
- class bonsu.phasing.HIO.HIOMaskPC(parent=None)¶
Bases:
PhaseAbstractPC
HIO Mask with Partial Coherence Optimisation.
- class bonsu.phasing.HIO.SWHIOMaskPC¶
Bases:
HIOMaskPC
,ShrinkWrap
HIO Mask with Partial Coherence Optimisation. Shrink wrapped.
- Start()¶
Start the reconstruction process.
Error Reduction Classes¶
- class bonsu.phasing.ER.ER(parent=None)¶
Bases:
PhaseAbstract
Error Reduction (ER) algorithm.
- class bonsu.phasing.ER.SWER¶
Bases:
ER
,ShrinkWrap
Error Reduction (ER) algorithm. Shrink wrapped.
- class bonsu.phasing.ER.ERMask(parent=None)¶
Bases:
PhaseAbstract
Error Reduction (ER) algorithm with the addition of a Fourier space constraint mask.
- class bonsu.phasing.ER.SWERMask¶
Bases:
ERMask
,ShrinkWrap
Error Reduction (ER) algorithm with the addition of a Fourier space constraint mask. Shrink wrapped.
- class bonsu.phasing.ER.POER(parent=None)¶
Bases:
PhaseAbstract
Phase Only Error Reduction (POER) algorithm with the addition of a Fourier space constraint mask.
- class bonsu.phasing.ER.ERMaskPC(parent=None)¶
Bases:
PhaseAbstractPC
ER Mask with Partial Coherence Optimisation algorithm.
- class bonsu.phasing.ER.SWERMaskPC¶
Bases:
ERMaskPC
,ShrinkWrap
ER Mask with Partial Coherence Optimisation algorithm. Shrink wrapped.
- Start()¶
Start the reconstruction process.
Hybrid Project-Reflection Classes¶
- class bonsu.phasing.HPR.HPR(parent=None)¶
Bases:
PhaseAbstract
Hybrid Projection Reflection (HPR) algorithm.
- class bonsu.phasing.HPR.SWHPR¶
Bases:
HPR
,ShrinkWrap
Hybrid Projection Reflection (HPR) algorithm. Shrink wrapped.
- class bonsu.phasing.HPR.HPRPC(parent=None)¶
Bases:
PhaseAbstractPC
HPR Mask with Partial Coherence Optimisation algorithm.
- class bonsu.phasing.HPR.SWHPRPC¶
Bases:
HPRPC
,ShrinkWrap
HPR Mask with Partial Coherence Optimisation algorithm. Shrink wrapped.
- Start()¶
Start the reconstruction process.
Relaxed Average Alternating Reflection Classes¶
- class bonsu.phasing.RAAR.RAAR(parent=None)¶
Bases:
PhaseAbstract
Relaxed Average Alternating Reflection (RAAR) algorithm.
- class bonsu.phasing.RAAR.SWRAAR¶
Bases:
RAAR
,ShrinkWrap
Relaxed Average Alternating Reflection (RAAR) algorithm. Shrink wrapped.
- class bonsu.phasing.RAAR.RAARPC(parent=None)¶
Bases:
PhaseAbstractPC
RAAR Mask with Partial Coherence Optimisation algorithm.
- class bonsu.phasing.RAAR.SWRAARPC¶
Bases:
RAARPC
,ShrinkWrap
RAAR Mask with Partial Coherence Optimisation algorithm. Shrink wrapped.
- Start()¶
Start the reconstruction process.
CSHIO Classes¶
- class bonsu.phasing.CSHIO.CSHIO(parent=None)¶
Bases:
PhaseAbstract
Compressed Sensing HIO (CSHIO) algorithm.
- SetPnorm(p)¶
Set p-norm of the Lebesgue space.
- GetPnorm()¶
Get p-norm of the Lebesgue space.
- SetEpsilon(epsilon)¶
Set positive relaxation parameter (epsilon) for the weighted p-norm.
- GetEpsilon()¶
Get positive relaxation parameter (epsilon) for the weighted p-norm.
- SetEpsilonmin(epsilonmin)¶
Set minimum value that epsilon can take.
- GetEpsilonmin()¶
Set minimum value that epsilon can take.
- SetDivisor(d)¶
Set amount by which epsilon is divided when constraint is met.
- GetDivisor()¶
Get amount by which epsilon is divided when constraint is met.
- SetEta(eta)¶
Set parameter in the divisor condition.
- GetEta()¶
Get parameter in the divisor condition.
- class bonsu.phasing.CSHIO.SWCSHIO¶
Bases:
CSHIO
,ShrinkWrap
Compressed Sensing HIO (CSHIO) algorithm. Shrink wrapped.
SO2D Classes¶
- class bonsu.phasing.SO2D.SO2D(parent=None)¶
Bases:
PhaseAbstract
- SetNumsoiter(numsoiter)¶
Set total number of iterations performed when optimising the step length.
- GetNumsoiter()¶
Get total number of iterations performed when optimising the step length.
- SetReweightiter(reweightiter)¶
Set iteration after which reweighting is applied. A negative value implies no reweighting.
- GetReweightiter()¶
Get iteration after which reweighting is applied. A negative value implies no reweighting.
- SetDtaumax(dtaumax)¶
Set maximum amount by which the step can be incremented.
- GetDtaumax()¶
Get maximum amount by which the step can be incremented.
- SetDtaumin(dtaumin)¶
Set minimum amount by which the step can be incremented.
- GetDtaumin()¶
Get minimum amount by which the step can be incremented.
- SetTaumax(taumax)¶
Set maximum size for the total step length.
- GetTaumax()¶
Get maximum size for the total step length.
- SetPsiexitratio(psiexitratio)¶
Set value below |psi|/|psi_0| that will halt step length optimisation.
- SetPsiexiterror(psiexiterror)¶
Set value below (psi^(n+1) - psi^(n))/psi^(n) that will halt step length optimisation.
- GetPsiexiterror()¶
Get value below (psi^(n+1) - psi^(n))/psi^(n) that will halt step length optimisation.
- class bonsu.phasing.SO2D.SWSO2D¶
Bases:
SO2D
,ShrinkWrap