Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Pathways interface simplification #41

Merged
merged 9 commits into from
Oct 20, 2020
Merged

Pathways interface simplification #41

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
a312ab1
dipolarkernel: simplified interface for pathways definition
luisfabib Oct 20, 2020
15deb59
dipolarbackground: simplified interface for pathways definition
luisfabib Oct 20, 2020
2c76a5c
ex_models: adapted to new interface
luisfabib Oct 20, 2020
337706f
docs: adapted and edited docs
luisfabib Oct 20, 2020
da28dbf
dipolarbackground: removed deprecated overtonecoeff option
luisfabib Oct 20, 2020
ac09b31
dipolarbackground & dipolarkernel: fixed normalization of pathway amp…
luisfabib Oct 20, 2020
c5c9dba
dipolarkernel & dipolarbackground: added option to disable pathway am…
luisfabib Oct 20, 2020
0f8017d
dipolarkernel & dipolarbackground: refactored examples in docstrings
luisfabib Oct 20, 2020
a55680a
renamed pathinfo -> pathways throughout
luisfabib Oct 20, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
102 changes: 45 additions & 57 deletions deerlab/dipolarbackground.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,3 @@

# dipolarbackground.py - Multipathway background generator
# --------------------------------------------------------
# This file is a part of DeerLab. License is MIT (see LICENSE.md).
Expand All @@ -7,20 +6,25 @@
import numpy as np
import types

def dipolarbackground(t,pathinfo,Bmodel,renormalize = True,overtonecoeff = 1):
def dipolarbackground(t, pathways, Bmodel, renormalize=True, renormpaths=True):
r""" Constructs background decay functions according to the multi-pathway model.

Parameters
----------
t : array_like
Time axis, in microseconds.
pathinfo : array_like with shape (p,2) or (p,3)
Array of pathway amplitudes (lambda), refocusing points (T0), and harmonics (n)
for multiple (p) dipolar pathways. Each row contains ``[lambda T0 n]`` or ``[lambda T0]``
for one pathway. If n is not given it is assumed to be 1. For a pathway with unmodulated contribution, set the refocusing time to ``numpy.nan``.
pathways : list of lists or scalar
List of pathways. Each pathway is defined as a list of the pathway's amplitude (lambda), refocusing time (T0),
and harmonic (n), i.e. ``[lambda, T0, n]`` or ``[lambda, T0]`` for one pathway. If n is not given it is assumed to be 1.
For a pathway with unmodulated contribution, only the amplitude must be specified, i.e. ``[Lambda0]``.
If a single value is specified, it is interpreted as the 4-pulse DEER pathway amplitude (modulation depth).
Bmodel : callable
Background basis function. A callable function accepting a time-axis array as first input and a pathway amplitude as a second, i.e. ``B = lambda t,lam: bg_model(t,par,lam)``

renormalize : boolean, optional
Re-normalization of the multi-pathway background to ensure the equality ``B(t=0)==1`` is satisfied. Enabled by default.
renormpaths: boolean, optional
Normalization of the pathway amplitudes such that ``Lam0 + lam1 + ... + lamN = 1``. Enabled by default.

Returns
-------
B : ndarray
Expand All @@ -30,12 +34,10 @@ def dipolarbackground(t,pathinfo,Bmodel,renormalize = True,overtonecoeff = 1):
----------------
renormalize : boolean
The multi-pathway background does not necessarily satisfy ``V(0) == 1``. This option enables/disables a re-normalization of the background function to ensure that equality is satisfied, by default it is enabled.
overtonecoeff : array_like
Array containing the overtone coefficients for RIDME experiments.

Notes
-----
Computes the multipathway background ``B`` for the time axis ``t`` corresponding to the dipolar pathways specified in ``pathinfo``. The total background is computed from the basis background function model specified in ``Bmodel``. For a multi-pathway DEER signal (e.g, 4-pulse DEER with 2+1 contribution; 5-pulse DEER with 4-pulse DEER residual signal, and more complicated experiments), ``pathinfo`` is a 2D-array that contains a list of modulation depths (amplitudes), refocusing times (in microseconds), and optional harmonics for all modulated pathway signals
Computes the multipathway background ``B`` for the time axis ``t`` corresponding to the dipolar pathways specified in ``pathways``. The total background is computed from the basis background function model specified in ``Bmodel``. For a multi-pathway DEER signal (e.g, 4-pulse DEER with 2+1 contribution; 5-pulse DEER with 4-pulse DEER residual signal, and more complicated experiments), ``pathways`` is a 2D-array that contains a list of modulation depths (amplitudes), refocusing times (in microseconds), and optional harmonics for all modulated pathway signals

Examples
--------
Expand All @@ -48,73 +50,59 @@ def dipolarbackground(t,pathinfo,Bmodel,renormalize = True,overtonecoeff = 1):
lam = 0.4 # modulation depth main signal
conc = 200 # spin concentration (uM)

pathinfo = [[],[]]
pathinfo[0] = [1-lam, NaN] # unmodulated part, gives offset
pathinfo[1] = [lam, 0] # main modulation, refocusing at time zero
pathways = []
path0 = [1-lam] # unmodulated part, gives offset
path1 = [lam, 0] # main modulation, refocusing at time zero
pathways = [path0, path1]

def Bmodel(t,lam):
return dl.bg_hom3d(t,conc,lam)

B = dl.dipolarbackground(t,pathinfo,Bmodel)


B = dl.dipolarbackground(t,pathways,Bmodel)
"""

# Ensure that all inputs are numpy arrays
t = np.atleast_1d(t)
pathinfo = np.atleast_1d(pathinfo)
overtonecoeff = np.atleast_1d(overtonecoeff)


if type(Bmodel) is not types.LambdaType:
raise TypeError('For a model with multiple modulated pathways, B must be a function handle of the type: @(t,lambda) bg_model(t,par,lambda)')

if not np.isreal(pathinfo).all:
raise TypeError('lambda/pathinfo must be a numeric array.')
if not isinstance(pathways,list): pathways = [pathways]
if len(pathways) == 1:
lam = pathways[0]
pathways = [[1-lam], [lam, 0]]

if len(pathinfo) == 1:
lam = pathinfo[0]
pathinfo = np.array([[1-lam, np.NaN], [lam, 0]])

if not np.any(np.shape(pathinfo)[1] != np.array([2, 3])):
raise TypeError('pathinfo must be a numeric array with two or three columns.')

if np.any(np.isnan(pathinfo[:, 0])):
raise ValueError('In pathinfo, NaN can only appear in the second column (refocusing time) e.g. path[1,:] = [Lam0 NaN]')

# Normalize the pathway amplitudes to unity
pathinfo[:, 0] = pathinfo[:, 0]/sum(pathinfo[:, 0])
lam = pathinfo[:, 0]
T0 = pathinfo[:, 1]
if np.shape(pathinfo)[1] == 2:
n = np.ones(np.shape(T0))
else:
n = pathinfo[:, 2]
pathways = [np.atleast_1d(path) for path in pathways]

# Get unmodulated pathways
unmodulated = [pathways.pop(i) for i,path in enumerate(pathways) if len(path)==1]
# Combine all unmodulated contributions
Lambda0 = sum(np.concatenate([path for path in unmodulated]))

# Check structure of pathways
for i,path in enumerate(pathways):
if len(path) == 2:
# If harmonic is not defined, append default n=1
pathways[i] = np.append(path,1)
elif len(path) != 3:
# Otherwise paths are not correctly defined
raise KeyError('The pathway #{} must be a list of two or three elements [lam, T0] or [lam, T0, n]'.format(i))

# Combine all unmodulated components, and eliminate from list
unmodulated = np.where(np.isnan(T0))
lam = np.delete(lam, unmodulated)
T0 = np.delete(T0, unmodulated)
n = np.delete(n, unmodulated)

# Fold overtones into pathway list
nCoeffs = len(overtonecoeff)
lam_,T0_,n_ = (np.empty(0) for _ in range(3))
for i in range(nCoeffs):
lam_ = np.concatenate((lam_, lam*overtonecoeff[i]))
T0_ = np.concatenate((T0_,T0))
n_ = np.concatenate((n_,n*(i+1)))
lam,T0,n = (lam_,T0_,n_)
nModPathways = len(lam)
# Normalize the pathway amplitudes to unity
if renormpaths:
lamsum = Lambda0 + sum([path[0] for path in pathways])
Lambda0 /= lamsum
for i in range(len(pathways)):
pathways[i][0] /= lamsum

# Construction of multi-pathway background function
#-------------------------------------------------------------------------------
Bnorm = 1
B = 1
for pathway in range(nModPathways):
B = B*Bmodel(n[pathway]*(t-T0[pathway]),lam[pathway])
Bnorm = Bnorm*Bmodel(-T0[pathway]*n[pathway],lam[pathway])
for pathway in pathways:
n,T0,lam = pathway
B = B*Bmodel(n*(t-T0),lam)
Bnorm = Bnorm*Bmodel(-T0*n,lam)

if renormalize:
B = B/Bnorm
Expand Down
109 changes: 56 additions & 53 deletions deerlab/dipolarkernel.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,8 +23,8 @@
def w0(g):
return (mu0/2)*muB**2*g[0]*g[1]/h*1e21 # Hz m^3 -> MHz nm^3 -> Mrad s^-1 nm^3

def dipolarkernel(t,r,pathinfo = 1, B = 1, method = 'fresnel', excbandwidth = inf, g = [ge, ge],
integralop = True, nKnots = 5001, renormalize = True, clearcache = False):
def dipolarkernel(t, r, pathways = 1, B = 1, method = 'fresnel', excbandwidth = inf, g = [ge, ge],
integralop = True, nKnots = 5001, renormalize = True, renormpaths=True, clearcache = False):
#===================================================================================================
r"""Compute the dipolar kernel operator which enables the linear transformation from
distance-domain to time-domain data.
Expand All @@ -35,10 +35,11 @@ def dipolarkernel(t,r,pathinfo = 1, B = 1, method = 'fresnel', excbandwidth = in
Dipolar time axis, in microseconds.
r : array_like
Distance axis, in nanometers.
pathinfo : array_like with shape (p,2) or (p,3)
Array of pathway amplitudes (lambda), refocusing points (T0), and harmonics (n)
for multiple (p) dipolar pathways. Each row contains ``[lambda T0 n]`` or ``[lambda T0]``
for one pathway. If n is not given it is assumed to be 1. For a pathway with unmodulated contribution, set the refocusing time to ``numpy.nan``.
pathways : list of lists or scalar
List of pathways. Each pathway is defined as a list of the pathway's amplitude (lambda), refocusing time (T0),
and harmonic (n), i.e. ``[lambda, T0, n]`` or ``[lambda, T0]`` for one pathway. If n is not given it is assumed to be 1.
For a pathway with unmodulated contribution, only the amplitude must be specified, i.e. ``[Lambda0]``.
If a single value is specified, it is interpreted as the 4-pulse DEER pathway amplitude (modulation depth).
B : callable or array_like
For a single-pathway model, the numerical background decay can be passed as an array.
For multiple pathways, a callable function must be passed, accepting a time-axis array as first input and a pathway amplitude as a second, i.e. ``B = lambda t,lam: bg_model(t,par,lam)``
Expand All @@ -60,9 +61,11 @@ def dipolarkernel(t,r,pathinfo = 1, B = 1, method = 'fresnel', excbandwidth = in
nKnots : scalar, optional
Number of knots for the grid of powder orientations to be used in the 'grid' kernel calculation method.
renormalize : boolean, optional
Re-normalization of multi-pathway kernels to ensure the equality ``K(t=0,r)==1`` is satisfied. Default is True.
Re-normalization of multi-pathway kernels to ensure the equality ``K(t=0,r)==1`` is satisfied. Enabled by default.
renormpaths: boolean, optional
Normalization of the pathway amplitudes such that ``Lam0 + lam1 + ... + lamN = 1``. Enabled by default.
clearcache : boolean, optional
Clear the cached dipolar kernels at the beginning of the function. Default is False.
Clear the cached dipolar kernels at the beginning of the function. Disabled by default.

Returns
--------
Expand All @@ -71,7 +74,7 @@ def dipolarkernel(t,r,pathinfo = 1, B = 1, method = 'fresnel', excbandwidth = in

Notes
-----
For a multi-pathway DEER [1]_ signal (e.g, 4-pulse DEER with 2+1 contribution 5-pulse DEER with 4-pulse DEER residual signal, and more complicated experiments), ``pathinfo`` contains a list of modulation depths (amplitudes) and refocusing times (in microseconds).
For a multi-pathway DEER [1]_ signal (e.g, 4-pulse DEER with 2+1 contribution 5-pulse DEER with 4-pulse DEER residual signal, and more complicated experiments), ``pathways`` contains a list of modulation depths (amplitudes) and refocusing times (in microseconds).
The background function specified as ''B'' is used as basis function, and the actual multipathway background included into the kernel is compued using :ref:`dipolarbackground`. The background in included in the dipolar kernel definition [2]_.
Optionally, the harmonic (1 = fundamental, 2 = first overtone, etc.) can be given as a third value in each row. This can be useful for modeling RIDME signals [3]_. If not given, the harmonic is 1 for all pathways.

Expand All @@ -81,9 +84,9 @@ def dipolarkernel(t,r,pathinfo = 1, B = 1, method = 'fresnel', excbandwidth = in
To specify the standard model for 4-pulse DEER with an unmodulated offset and a single dipolar pathway that refocuses at time 0, use::

lam = 0.4 # modulation depth main signal
pathinfo = [[1-lam, nan], [lam, 0]]

K = dl.dipolarkernel(t,r,pathinfo)
pathways = [[1-lam], [lam, 0]]
luisfabib marked this conversation as resolved.
Show resolved Hide resolved
K = dl.dipolarkernel(t,r,pathways)


A shorthand input syntax equivalent to this input::
Expand All @@ -97,22 +100,23 @@ def dipolarkernel(t,r,pathinfo = 1, B = 1, method = 'fresnel', excbandwidth = in
Lam0 = 0.5 # unmodulated part
lam = 0.4 # modulation depth main signal
lam21 = 0.1 # modulation depth 2+1 contribution
tau2 = 4 # refocusing time of 2+1 contribution
pathinfo = [[],[],[]]
pathinfo[0] = [Lam0, nan] # unmodulated part, gives offset
pathinfo[1] = [lama, 0 ] # main modulation, refocusing at time zero
pathinfo[2] = [lam21, tau2] # 2+1 modulation, refocusing at time tau2
tau2 = 4 # refocusing time (us) of 2+1 contribution

path0 = Lam0 # unmodulated pathway
path1 = [lam1, 0] # main dipolar pathway, refocusing at time zero
path1 = [lam2, tau2] # 2+1 dipolar pathway, refocusing at time tau2
pathways = [path0, path1, path2]

luisfabib marked this conversation as resolved.
Show resolved Hide resolved
K = dl.dipolarkernel(t,r,pathinfo)
K = dl.dipolarkernel(t,r,pathways)


References
----------
.. [1] L. Fábregas Ibáñez, G. Jeschke, and S. Stoll: "DeerLab: A comprehensive toolbox for analyzing dipolar EPR spectroscopy data", Magn. Reson. Discuss., 2020
.. [1] L. Fábregas Ibáñez, G. Jeschke, and S. Stoll. DeerLab: A comprehensive toolbox for analyzing dipolar EPR spectroscopy data, Magn. Reson., 1, 209–224, 2020

.. [2] Fábregas Ibáñez, L. and Jeschke, G.: Optimal background treatment in dipolar spectroscopy, Physical Chemistry Chemical Physics, 22, 1855–1868, 2020.
.. [2] L. Fábregas Ibáñez, and G. Jeschke: Optimal background treatment in dipolar spectroscopy, Physical Chemistry Chemical Physics, 22, 1855–1868, 2020.

.. [3] Keller, K., Mertens, V., Qi, M., Nalepa, A. I., Godt, A., Savitsky, A., Jeschke, G., and Yulikov, M.: Computing distance distributions from dipolar evolution data with overtones: RIDME spectroscopy with Gd(III)-based spin labels, Physical Chemistry Chemical Physics, 19
.. [3] K. Keller, V. Mertens, M. Qi, A. I. Nalepa, A. Godt, A. Savitsky, G. Jeschke, and M. Yulikov: Computing distance distributions from dipolar evolution data with overtones: RIDME spectroscopy with Gd(III)-based spin labels, Physical Chemistry Chemical Physics, 19

.. [4] J.E. Banham, C.M. Baker, S. Ceola, I.J. Day, G.H. Grant, E.J.J. Groenen, C.T. Rodgers, G. Jeschke, C.R. Timmel, Distance measurements in the borderline region of applicability of CW EPR and DEER: A model study on a homologous series of spin-labelled peptides, Journal of Magnetic Resonance, 191, 2, 2008, 202-218
"""
Expand All @@ -123,7 +127,6 @@ def dipolarkernel(t,r,pathinfo = 1, B = 1, method = 'fresnel', excbandwidth = in
# Ensure that all inputs are numpy arrays
r = np.atleast_1d(r)
t = np.atleast_1d(t)
pathinfo = np.atleast_1d(pathinfo)

if type(B) is types.LambdaType:
ismodelB = True
Expand All @@ -140,54 +143,54 @@ def dipolarkernel(t,r,pathinfo = 1, B = 1, method = 'fresnel', excbandwidth = in
if np.any(r<=0):
raise ValueError("All elements in r must be nonnegative and nonzero.")

if not np.isreal(pathinfo).all:
raise TypeError('lambda/pathinfo must be a numeric array.')
if not isinstance(pathways,list): pathways = [pathways]
if len(pathways) == 1:
lam = pathways[0]
pathways = [[1-lam], [lam, 0]]

if len(pathinfo) == 1:
lam = pathinfo[0]
pathinfo = np.array([[1-lam, np.NaN], [lam, 0]])
paths = [np.atleast_1d(path) for path in pathways]

if not np.any(np.shape(pathinfo)[1] != np.array([2, 3])):
raise TypeError('pathinfo must be a numeric array with two or three columns.')
# Get unmodulated pathways
unmodulated = [paths.pop(i) for i,path in enumerate(paths) if len(path)==1]
# Combine all unmodulated contributions
Lambda0 = sum(np.concatenate([path for path in unmodulated]))

if np.any(np.isnan(pathinfo[:,0])):
raise ValueError('In pathinfo, NaN can only appear in the second column (refocusing time) e.g. path[1,:] = [Lam0 NaN]')
# Check structure of pathways
for i,path in enumerate(paths):
if len(path) == 2:
# If harmonic is not defined, append default n=1
paths[i] = np.append(path,1)
elif len(path) != 3:
# Otherwise paths are not correctly defined
raise KeyError('The pathway #{} must be a list of two or three elements [lam, T0] or [lam, T0, n]'.format(i))

# Normalize the pathway amplitudes to unity
pathinfo[:,0] = pathinfo[:,0]/sum(pathinfo[:,0])
lam = pathinfo[:,0]
T0 = pathinfo[:,1]
if np.shape(pathinfo)[1] == 2:
n = np.ones(np.shape(T0))
else:
n = pathinfo[:,2]

# Combine all unmodulated components into Lambda0, and eliminate from list
unmodulated = np.where(np.isnan(T0))
Lambda0 = np.sum(lam[unmodulated])
lam = np.delete(lam, unmodulated)
T0 = np.delete(T0, unmodulated)
n = np.delete(n, unmodulated)
nModPathways = len(lam)
if renormpaths:
lamsum = Lambda0 + sum([path[0] for path in paths])
Lambda0 /= lamsum
for i in range(len(paths)):
paths[i][0] /= lamsum

# Define kernel matrix auxiliary function
kernelmatrix = lambda t: calckernelmatrix(t,r,method,excbandwidth,nKnots,g)

# Build dipolar kernel matrix, summing over all pathways
K = Lambda0
for pathway in range(nModPathways):
K = K + lam[pathway]*kernelmatrix(n[pathway]*(t-T0[pathway]))
for pathway in paths:
lam,T0,n = pathway
K = K + lam*kernelmatrix(n*(t-T0))

# Renormalize if requested
if renormalize:
Knorm = Lambda0
for pathway in range(nModPathways):
Knorm = Knorm + lam[pathway]*kernelmatrix(-T0[pathway]*n[pathway])
for pathway in paths:
lam,T0,n = pathway
Knorm = Knorm + lam*kernelmatrix(-T0*n)
K = K/Knorm


# Multiply by background
if ismodelB:
B = dipolarbackground(t,pathinfo,B)
B = dipolarbackground(t,pathways,B)
K = K*B[:,np.newaxis]

# Include delta-r factor for integration
Expand Down
Loading