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.

GetPsiexitratio()

Get 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.

SetPsiresetratio(psiresetratio)

Set value above |psi|/|psi_0| that will reset the step length to that of the initial value.

GetPsiresetratio()

Get value above |psi|/|psi_0| that will reset the step length to that of the initial value.

class bonsu.phasing.SO2D.SWSO2D

Bases: SO2D, ShrinkWrap