# dipolarmodel.py - DeerLab's dipolar EPR model generator
# ---------------------------------------------------------------------------
# This file is a part of DeerLab. License is MIT (see LICENSE.md).
# Copyright(c) 2019-2022: Luis Fabregas, Stefan Stoll and other contributors.
import numpy as np
from deerlab.dipolarkernel import dipolarkernel, dipolarbackground
from deerlab.regoperator import regoperator
from deerlab.dd_models import freedist
from deerlab.model import Model,Penalty, link
from deerlab import bg_hom3d
from deerlab.utils import choleskycovmat
from deerlab.constants import *
from math import factorial
from itertools import permutations
import inspect
from scipy.stats import multivariate_normal
#===============================================================================
[docs]
def dipolarmodel(t, r=None, Pmodel=None, Bmodel=bg_hom3d, experiment=None, parametrization='reftimes', npathways=1, spins=2, harmonics=None,
excbandwidth=np.inf, orisel=None, g=[ge,ge], gridsize=1000, minamp=1e-3, samespins=True, triangles=None, interp=True,):
"""
Generate a dipolar EPR signal model.
This function generates a model for the dipolar EPR signal arising from the dipolar interactions
between the electron spins of the system. The model is defined in terms of the evolution time vector
``t``, the spin-spin distances ``r``, the distance distribution model ``Pmodel``, the intermolecular contribution
``Bmodel``, and an ``experiment`` object containing information about the specific experiment.
The number of spins in the system ``spins`` can be specified, as well as the parametrization strategy
``parametrization`` for the dipolar pathway refocusing times.
For multi-spin systems (``spins > 2``), the function allows for the inclusion of three-spin interactions by
specifying the ``triangles`` list of interspin triangles and enabling the ``samespins`` assumption of spectral
permutability. The ``minamp`` parameter can be used to threshold the amplitude of the three-spin interactions.
Additional effects such as orientation selection or limited excitation bandwidth can be included via the ``orisel``
and ``excbandwidth`` arguments, respectively.
Parameters
----------
t : array_like
Vector of dipolar evolution times, in microseconds.
r : array_like, optional
Vector of spin-spin distances, in nanometers. If not specified, it defaults to the range
1.5nm-8nm with a resolution of 0.05nm. Has no effect on models with ``spins>2``.
Pmodel : :ref:`Model`, optional
Model for the distance distribution. If not specified, a non-parametric
distance distribution defined over ``r`` is used. For models with ``spins>2`` a multivariate
Gaussian distance distribution is assumed.
Bmodel : :ref:`Model`, optional
Model for the intermolecular (background) contribution. If not specified,
a background arising from a homogenous 3D distribution of spins is used.
experiment : :ref:`ExperimentInfo`, optional
Experimental information obtained from experiment models (``ex_``). If specified, the
boundaries and start values of the dipolar pathways' refocusing times and amplitudes
are refined based on the specific experiment's delays.
parametrization : string, optional
Parametrization strategy of the dipolar pathway refocusing times. Must be one of the following:
* ``'reftimes'`` - Each refocusing time is represented individually as a parameter.
* ``'delays'`` - The pulse delays are introduced as parameters from which the refocusing times are computed. Requires ``experiment`` to be specified.
* ``'shift'`` - A time shift is introduced as a parameter to represent the variability of the refocusing times from their theoretical values. Requires ``experiment`` to be specified.
The default is ``'reftimes'``.
npathways : integer scalar, optional
Number of dipolar pathways. If not specified, a single dipolar pathway is used.
harmonics : list of integers, optional
Harmonic prefactors of the dipolar pathways. Must be a list with ``npathways`` harmonics, one
for each pathway.
spins : integer scalar, optional
Number of spins in system. If not specified it defaults to two-spin systems. For multi-spins
systems indicated by ``spins>2`` three-spin interactions are accounted for by the model automatically.
triangles : list of lists, optional
Only required when ``spins>3``. List of ``(Nspins)!/(3!*(Nspins-3)!)`` triangles formed by the different interspin distances.
Each triangle is specified by a list of indices of the interspin
distances forming the triangle, e.g. for a three-spin system ``triangles=[[1,2,3]]`` where the first,
second, and third interspin distances connect to form a triangle.
samespins : boolean, optional
Only when ``spins>3``. Enables the assumption of spectral permutability, i.e. the assumption that all spins in the system
have the same spectral distribution. If enabled, all pairwise pathways of the different dipolar interactions are assumed to
have the same amplitude. If disabled, all the individual pathway amplitudes are parametrized separately. Enabled by default.
minamp : float scalar, optional
Only when ``spins>3``. Threshold amplitude for neglecting three-spin interaction pathways. To enhance performance, all dipolar
pathways arising from three-spin interactions with an amplitude ``lam<minamp`` are approximated as having zero amplitude.
The default threshold is ``1e-3``.
orisel : callable or ``None``, optional
Probability distribution of possible orientations of the interspin vector to account
for orientation selection. Must be a function taking a value of the angle θ ∈ [0,π/2]
between the interspin vector and the external magnetic field and returning the corresponding
probability density. If specified as ``None`` (default), a uniform distribution is used.
gridsize : scalar, optional
Number of points on the grid of powder orientations to be used in the ``'grid'`` kernel
calculation method when evaluating models with orientation selection or multi-spin effects.
By default set to 1000 points.
excbandwidth : scalar, optional
Excitation bandwidth of the pulses in MHz to account for limited excitation bandwidth.
If not specified, an infinite excitation bandwidth is used.
g : scalar, 2-element array, optional
Electron g-values of the two spins, ``[g1, g2]``. If a single g is specified, ``[g, g]`` is used.
If not specified, g = 2.002319... is used for both spins.
interp : boolean, optional
Enable dipolar kernel interpolation for computation time reduction. By default enabled.
Returns
-------
Vmodel : :ref:`Model`
Dipolar signal model object.
"""
# If distance axis not specified default to a long range one
if r is None:
r = np.arange(1.5,8,0.05)
# Input parsing and validation
if not isinstance(Pmodel,Model) and Pmodel is not None:
raise TypeError('The argument Pmodel must be a valid Model object')
if not isinstance(Bmodel,Model) and Bmodel is not None:
raise TypeError('The argument Bmodel must be a valid Model object or None.')
if not isinstance(npathways,int) or npathways<=0:
raise ValueError('The number of pathways must be a positive integer.')
if isinstance(harmonics,float) or isinstance(harmonics,int):
harmonics = [int(harmonics)]
if not isinstance(harmonics,list) and harmonics is not None:
raise TypeError('The harmonics must be specified as a list of integer values.')
if experiment is not None:
if not isinstance(experiment,ExperimentInfo):
raise TypeError('The experiment must be a valid deerlab.ExperimentInfo object.')
# Check that the number of requested pathways does not exceed the theoretical limit of the experiment
npathways = experiment.npathways
maxpathways = len(experiment.harmonics)
if npathways>maxpathways:
raise ValueError(f'The {experiment.name} experiment can only have up to {maxpathways} dipolar pathways.')
#------------------------------------------------------------------------
def _importparameter(parameter):
""" Private function for importing a parameter's metadata """
return {
'lb' : parameter.lb,
'ub' : parameter.ub,
'par0' : parameter.par0,
'description' : parameter.description,
'unit' : parameter.unit,
'linear' : parameter.linear
}
#------------------------------------------------------------------------
# Parse the harmonics of the dipolar pathways
if harmonics is None:
if experiment is not None:
harmonics = np.array(experiment.harmonics)[:npathways]
else:
harmonics = np.ones(npathways)
if len(harmonics)!=npathways:
raise ValueError('The number of harmonics must match the number of dipolar pathways.')
if parametrization!='reftimes':
if experiment is None:
raise SyntaxError(f"The '{parametrization}' parametrization requires 'experiment' to be specified.")
# Determine number of dipolar interactions in the spins system
Q = int(spins*(spins-1)/2)
if spins>2:
triangles_required = factorial(spins)/(factorial(3)*factorial(spins-3))
if spins==3 and triangles is None:
triangles = [[1,2,3]]
elif spins>3 and triangles is None:
raise SyntaxError(f'For systems with {spins} spins, a list of {triangles_required:n} triangles is required.')
elif spins>3 and len(triangles)!=triangles_required:
raise SyntaxError(f'For systems with {spins} spins, a list of {triangles_required:n} triangles is required.')
if np.min(triangles)!=0:
triangles = [np.array(triangle)-1 for triangle in triangles]
#------------------------------------------------------------------------
def dipolarpathways(*param):
"""
Generator of dipolar pathways based on a parameter set.
Constructs all possible two-spin dipolar pathways based on a parameter set,
and then constructs all possible multi-spin dipolar pathways describing three-spin
interactions in systems with more than two spins based on the pair-wise dipolar pathways.
The true signature of this function is defined dynamically after its definition below.
"""
param = np.atleast_1d(param)
if spins==2:
# Construct pair-wise dipolar pathways given by the parameter set
if parametrization=='reftimes':
λs = param[np.arange(0,len(param),2)]
trefs = param[np.arange(1,len(param),2)]
elif parametrization=='delays':
delays = param[0:len(experiment.delays)]
λs = param[len(delays):]
trefs = experiment.reftimes(*delays)
elif parametrization=='shift':
shift = param[0]
λs = param[1:]
delays = np.array(experiment.delays) + shift
trefs = experiment.reftimes(*delays)
# Unmodulated pathway
Λ0 = np.maximum(0,1 - np.sum(λs))
pairpathways = [{'amp':Λ0}]
# Modulated pathways
for λ,tref,δ in zip(λs,trefs,harmonics):
pairpathways.append({'amp': λ, 'reftime': tref, 'harmonic': δ})
# For two-spin systems these are all the required pathways
return pairpathways
# Construct multispin dipolar pathways given by the parameter set
λs = [param[np.arange(q,len(param),Q)] for q in range(Q)]
trefs = param[np.arange(Q,len(param),Q)]
# Amplitude of the unmodulated pair-wise dipolar pathway
λu = param[-1]
# Initialize containers
pathways,threespin_pathways = [],[]
Λ0 = 1
# Construct all possible multi-spin dipolar pathways
for idx in range(npathways):
# Construct all the permutations of the two-spin interaction pathways (without repetitions)
for perm in set(set(permutations([0]+[None]*(Q-1)))):
q = int(np.where(np.array(perm)==0)[0])
# Compute amplitude of the two-spin interaction pathway
λp = λs[q][idx]
λ2k = factorial(Q-1)*λp*λu**(Q-1)
tref = [trefs[idx] if n==0 else None for n in perm]
δs = [harmonics[idx] if n==0 else 0 for n in perm]
# Add two-spin interaction dipolar pathway to the list
pathways.append({'amp':λ2k, 'reftime': tuple(tref), 'harmonic': tuple(δs)})
# Update the unmodulated amplitude
Λ0 -= λ2k
# Consider now three-spin interactions
for idx2 in np.arange(idx + 1,npathways):
# If the pathway has negligible amplitude, ignore it (speed-up)
if λ3k>minamp:
# Construct all the permutations of the multi-spin pathway (without repetitions)
for perm in set(set(permutations([0,1]+[None]*(Q-2)))):
# Check if the permutation corresponds to a valid interaction triangle
q1,q2 = [np.where(perm==n)[0] for n in [0,1]]
for triangle in triangles:
if q1 not in triangle and q2 not in triangle:
continue # If not valid, move on to the next one
# Compute amplitude of the three-spin interaction pathway
λm = λs[q2][idx]
λ3k = 2*factorial(Q-2)*λp*λm*λu**(Q-2)
tref = [trefs[idx] if n==0 else trefs[idx2] if n==1 else None for n in perm]
δs = [harmonics[idx] if n==0 else harmonics[idx2] if n==1 else 0 for n in perm]
# Add three-spin interaction dipolar pathway to the list
threespin_pathways.append({'amp':λ3k, 'reftime': tuple(tref), 'harmonic': tuple(δs)})
# Update the unmodulated amplitude
Λ0 -= λ3k
# Add pathway with unmodulated contribution
pathways.append({'amp':Λ0})
return pathways,threespin_pathways
#------------------------------------------------------------------------
# Construct the signature of the dipolarpathways() function dynamically
if experiment is not None:
labels = experiment.labels
else:
# Otherwise, just label them in order
labels = np.arange(1,npathways+1)
variables = []
if parametrization=='delays':
pulsedelay_names = inspect.signature(experiment.reftimes).parameters.keys()
variables = list(pulsedelay_names)
if parametrization=='shift':
variables = ['tshift']
for n in range(npathways):
if npathways == 1 and spins==2:
# Use single-pathway notation, with amplitude as modulation depth
variables.append('mod')
else:
if spins==2:
variables.append(f'lam{labels[n]}')
else:
for q in range(Q):
variables.append(f'lam{labels[n]}_{q+1}')
if parametrization=='reftimes':
if npathways == 1 and spins==2:
variables.append(f'reftime')
else:
variables.append(f'reftime{labels[n]}')
# Add unmodulated pairwise pathway amplitude for multi-spin systems
if spins>2:
variables.append('lamu')
# Create the dipolar pathways model object
PathsModel = Model(dipolarpathways,signature=variables)
if spins>2 and samespins:
links = {f'lam{labels[n]}': [f'lam{labels[n]}_{q+1}' for q in range(Q)] for n in range(npathways)}
PathsModel = link(PathsModel, **links)
#-----------------------------------------------------------------------
def _Pmultivar(param):
rmeans,cholesky_factors = param[:Q], param[Q:]
Σ = choleskycovmat(Q,cholesky_factors)
Pmultivar = multivariate_normal(rmeans, cov=Σ)
return Pmultivar
#-----------------------------------------------------------------------
# Determine if distance distributions is non-parametric
Pnonparametric = Pmodel is None
if Pnonparametric:
Pmodel = freedist(r)
# For multi-spin systems, there is no longer freedom of choice for distance distributions
if spins>2:
# Prepare the normal multivariate distance distribution model
variables = []
variables += [f'rmean{q+1}' for q in range(Q)]
variables += [f'chol{q+1}{q+1}' for q in range(Q)]
for q in range(Q):
variables += [f'chol{j+1}{q+1}' for j in np.arange(q+1,Q)]
Pmodel = Model(_Pmultivar,signature=variables)
[getattr(Pmodel,f'rmean{q+1}').set(lb=0,ub=10,par0=3,unit='nm',description=f'Average inter-spin distance #{q+1}') for q in range(Q)]
[getattr(Pmodel,f'chol{q+1}{q+1}').set(lb=0,ub=6,par0=0.4,unit='nm',description=f'Cholesky factor ℓ{q+1}{q+1}') for q in range(Q)]
for q in range(Q):
[getattr(Pmodel,f'chol{j+1}{q+1}').set(lb=-1,ub=1,par0=0.0,unit='nm',description=f'Cholesky factor ℓ{j+1}{q+1}') for j in np.arange(q+1,Q)]
Nconstants = len(Pmodel._constantsInfo)
# Populate the basic information on the dipolar pathways parameters
for n in range(npathways):
if spins==2 and npathways==1:
# Special case: use modulation depth notation instead of general pathway amplitude
getattr(PathsModel,f'mod').set(lb=0,ub=1,par0=0.01,description=f'Modulation depth',unit='')
else:
pairwise = ''
if spins>2:
pairwise = ' pairwise'
if spins==2 or samespins:
getattr(PathsModel,f'lam{labels[n]}').set(lb=0,ub=1,par0=0.01,description=f'Amplitude of{pairwise} pathway #{labels[n]}',unit='')
else:
for q in range(Q):
getattr(PathsModel,f'lam{labels[n]}_{q+1}').set(lb=0,ub=1,par0=0.01,description=f'Amplitude of{pairwise} pathway #{labels[n]} on interaction #{q+1}',unit='')
if parametrization=='reftimes':
if experiment is None:
theoretical_reftime = 0
reftime_variability = 20
else:
theoretical_reftime = experiment.reftimes(*experiment.delays)[n]
reftime_variability = 3*experiment.pulselength
if spins==2 and npathways==1:
# Special case: use reftime notation
getattr(PathsModel,f'reftime').set(description=f'Refocusing time',unit='μs',
par0 = theoretical_reftime,
lb = theoretical_reftime - reftime_variability,
ub = theoretical_reftime + reftime_variability)
else:
# Special case: use reftime notation
getattr(PathsModel,f'reftime{labels[n]}').set(description=f'Refocusing time of{pairwise} pathway #{labels[n]}',unit='μs',
par0 = theoretical_reftime,
lb = theoretical_reftime - reftime_variability,
ub = theoretical_reftime + reftime_variability)
if parametrization=='delays':
experimental_delays = experiment.delays
delay_variability = 5*experiment.pulselength
for n,delay in enumerate(pulsedelay_names):
getattr(PathsModel,delay).set(description=f'Pulse delay {delay} of the experiment',unit='μs',
par0 = experimental_delays[n],
lb = experimental_delays[n] - delay_variability,
ub = experimental_delays[n] + delay_variability)
elif parametrization=='shift':
variability = 5*experiment.pulselength
getattr(PathsModel,'tshift').set(par0=0,lb=-variability,ub=variability,description=f'Variability of experimental pulse delays',unit='μs')
if spins>2:
getattr(PathsModel,f'lamu').set(lb=0,ub=1,par0=0.5,description='Amplitude of unmodulated pairwise pathway',unit='')
# Construct the signature of the dipolar signal model function
signature = []
parameters,linearparam = [],[]
for model in [PathsModel,Bmodel,Pmodel]:
if model is not None:
for param in model._parameter_list(order='vector'):
if np.any(getattr(model,param).linear):
parameters.append(getattr(model,param))
linearparam.append({'name':param,'vec':len(np.atleast_1d(getattr(model,param).idx)),'normalization':getattr(model,param).normalization})
elif not (model==Bmodel and param=='lam'):
signature.append(param)
parameters.append(getattr(model,param))
# Initialize lists of indices to access subsets of nonlinear parameters
Psubset = np.zeros(Pmodel.Nnonlin,dtype=int)
PathsSubset = np.zeros(PathsModel.Nnonlin,dtype=int)
if Bmodel is None:
Bsubset = []
else:
Bsubset = np.zeros(Bmodel.Nnonlin-('lam' in Bmodel._parameter_list()),dtype=int)
# Determine subset indices based on main function signature
idx = 0
for model,subset in zip([Pmodel,Bmodel,PathsModel],[Psubset,Bsubset,PathsSubset]):
if model is not None:
for idx,param in enumerate(signature):
if param in model._parameter_list(order='vector'):
subset[getattr(model,param).idx] = idx
kernelmethod = 'fresnel' if orisel is None and np.isinf(excbandwidth) else 'grid'
dt = np.mean(np.diff(t))
tinterp = None
#------------------------------------------------------------------------
def Vtwospin_nonlinear_fcn(*nonlin):
""" Non-linear part of the dipolar signal function """
# Make input arguments as array to access subsets easily
nonlin = np.atleast_1d(nonlin)
# Construct the basis function of the intermolecular contribution
if Bmodel is None:
Bfcn = np.ones_like(t)
elif hasattr(Bmodel,'lam'):
Bfcn = lambda t,lam: Bmodel.nonlinmodel(t,*np.concatenate([nonlin[Bsubset],[lam]]))
else:
Bfcn = lambda t,_: Bmodel.nonlinmodel(t,*nonlin[Bsubset])
# Construct the definition of the dipolar pathways
pathways = PathsModel.nonlinmodel(*nonlin[PathsSubset])
# Construct the dipolar kernel
Kdipolar = dipolarkernel(t, r, pathways=pathways, bg=Bfcn, tinterp=tinterp,
excbandwidth=excbandwidth, orisel=orisel,
g=g, method=kernelmethod, gridsize=gridsize)
# Compute the non-linear part of the distance distribution
Pnonlin = Pmodel.nonlinmodel(*[r]*Nconstants,*nonlin[Psubset])
# Forward calculation of the non-linear part of the dipolar signal
Vnonlin = Kdipolar@Pnonlin
return Vnonlin
#------------------------------------------------------------------------
#------------------------------------------------------------------------
def Vmultispin_nonlinear_fcn(*nonlin):
""" Non-linear part of the dipolar signal function """
# Make input arguments as array to access subsets easily
nonlin = np.atleast_1d(nonlin)
# Construct the basis function of the intermolecular contribution
if Bmodel is None:
Bfcn = np.ones_like(t)
elif hasattr(Bmodel,'lam'):
Bfcn = lambda t,lam: Bmodel.nonlinmodel(t,*np.concatenate([nonlin[Bsubset],[lam]]))
else:
Bfcn = lambda t,_: Bmodel.nonlinmodel(t,*nonlin[Bsubset])
# Sample from the multi-variate distance distribution
Nsamples = 1000000
rsamples = Pmodel.nonlinmodel(nonlin[Psubset]).rvs(Nsamples)
rsamples = np.maximum(rsamples,1e-16) # Avoid values exactly at zero
# Enforce the generalized triangle inequiality
for triangle in triangles: # Loop over all triangle combinations
idx = triangle
triangle_condition = np.full_like(rsamples,False)
for q in range(Q): # Evaluate all triangle inequalities
idx = np.roll(idx,-1,axis=0)
triangle_condition[:,q] = np.sum(rsamples[:,idx[:-1]],axis=1) > rsamples[:,idx[-1]]
triangle_condition = np.all(triangle_condition,axis=1)
# Discard samples that do not satisfy triangle inequalities
rsamples = rsamples[triangle_condition,:]
# Generate all multi-spin dipolar pathways
twospin_pathways, threespin_pathways = PathsModel.nonlinmodel(*nonlin[PathsSubset])
# Start by accounting for the unmodulated contribution
Λ0 = twospin_pathways.pop(-1)['amp']
Vintra = Λ0
# Account for all two-spin intramolecular contributions
for q in range(Q):
# Get the marginal univariate distance distribution
Punimarg,bins = np.histogram(rsamples[:,q], bins=100)
rq = (bins[:-1] + bins[1:])/2
rs = [rq]*Q
if np.sum(Punimarg)>0:
Punimarg = Punimarg/np.sum(Punimarg)
# Get the two-spin pathways modulated by the q-th dipolar interaction
pathways_q = [pathway for pathway in twospin_pathways if pathway['harmonic'][q]!=0]
# Construct the corresponding dipolar kernel
Ktwospin = dipolarkernel(t,rs, pathways=pathways_q, integralop=False, tinterp=tinterp,
excbandwidth=excbandwidth, orisel=orisel,
g=g, method=kernelmethod, gridsize=gridsize)
# Two-spin intramolecular contribution from all pathways modulated by the q-th interaction
Vintra += Ktwospin@Punimarg
# Account for all three-spin intramolecular contributions
if threespin_pathways:
for triangle in triangles:
# Get the marginal trivariate distance distribution
Ptrimarg,bins = np.histogramdd(rsamples[:,triangle],bins=10, density=True)
r1,r2,r3 = [ (bins_q[:-1] + bins_q[1:])/2 for bins_q in bins ]
# Get the subset of the most distance combinations with most distribution mass (speed-up)
subset = Ptrimarg > np.max(Ptrimarg)/100
Ptrimarg = Ptrimarg[subset]
if np.sum(Ptrimarg)>0:
Ptrimarg = Ptrimarg/np.sum(Ptrimarg)
r1,r2,r3 = [rq[subset] for rq in np.meshgrid(r1,r2,r3, indexing='ij')]
# Get the three-spin pathways modulated by the triangle interaction
pathways_triangle = [pathway for pathway in threespin_pathways if np.sum([pathway['harmonic'][q]!=0 for q in triangle])==2]
# Construct the corresponding dipolar kernel
Kthreespin = dipolarkernel(t,[r1,r2,r3], pathways=pathways_triangle, integralop=False,
excbandwidth=excbandwidth, orisel=orisel,
g=g, method='grid', gridsize=gridsize)
# Three-spin intramolecular contribution from all pathways modulated by the triangle interaction
Vintra += Kthreespin@Ptrimarg
# Compute the multi-spin dipolar intermolecular contribution
Vinter = dipolarbackground(t, twospin_pathways+threespin_pathways, Bfcn)
# Construct the full dipolar signal
Vnonlin = Vintra*Vinter
return Vnonlin
#------------------------------------------------------------------------
# Create the dipolar model object
if spins==2:
DipolarSignal = Model(Vtwospin_nonlinear_fcn,signature=signature)
else:
DipolarSignal = Model(Vmultispin_nonlinear_fcn,signature=signature)
# Add the linear parameters from the subset models
for lparam in linearparam:
DipolarSignal.addlinear(lparam['name'],vec=lparam['vec'],normalization=lparam['normalization'])
if Pmodel is None:
DipolarSignal.addlinear()
# If there are no linear paramters, add a linear scaling parameter
if DipolarSignal.Nlin==0:
DipolarSignal.addlinear('scale',lb=0,par0=1,description='Overall echo amplitude/scale')
# Import all parameter information from the subset models
for name,param in zip(DipolarSignal._parameter_list(order='vector'),parameters):
getattr(DipolarSignal,name).set(**_importparameter(param))
DipolarSignal.description = 'Dipolar signal model'
# Set prior knowledge on the parameters if experiment is specified
if experiment is not None:
# Specify experiment in model description
DipolarSignal.description = f'{experiment.name} dipolar signal model'
# Specify start values and boundaries according to experimental timings
reftimes = experiment.reftimes(*experiment.delays)
if interp:
tinterp = np.arange(
min(t)-max(reftimes) - 3*experiment.pulselength - dt,
max(t)+max(reftimes) + 3*experiment.pulselength + dt,
dt)
return DipolarSignal
#===============================================================================
#===============================================================================
[docs]
def dipolarpenalty(Pmodel, r, type, selection=None):
r"""
Construct penalties based on the distance distribution.
This function constructs penalties for a given distance distribution that can be used to impose
certain desired properties on the distribution. The type of penalty and the function to optimize
it are determined by the ``type`` parameter. The ``selection`` parameter allows the user to choose
a method for optimizing the penalty weight.
Parameters
----------
Pmodel : `Model` object
The model of the distance distribution. If ``Pmodel`` is ``None``, a non-parametric distance distribution is assumed.
r : array_like
Distance axis vector, in nanometers. This vector specifies the distances at which the distance distribution is defined.
type : string
Type of property to be imposed by the penalty. This parameter determines the type of penalty
that will be constructed, as well as the function that will be optimized to determine the
penalty weight. The following values are supported for this parameter:
- ``'smoothness'`` : Imposes a smoothness constraint on the distance distribution. This penalty
is optimized by minimizing the L2 norm of the second derivative of the distribution.
- ``'compactness'`` : Imposes a compactness constraint on the distance distribution. This penalty
is optimized by minimizing the square root of the product of the distance distribution and
the square of the difference between the distances and their mean.
selection : string, optional
Selection functional for the outer optimization of the penalty weight. Possible values:
- ``'aic'`` - Akaike information criterion
- ``'bic'`` - Bayesian information criterion
- ``'aicc'`` - COrrected Akaike information criterion
- ``'icc'`` - Informational complexity criterion
If no value is provided, ``'icc'`` is used for ``'compactness'`` and ``'aic'`` for ``'smoothness'``.
Returns
-------
penalty : ``Penalty`` object
Penalty object to be passed to the ``fit`` function.
"""
if Pmodel is None:
Pmodel = freedist(r)
Nconstants = len(Pmodel._constantsInfo)
# If include compactness penalty
if type=='compactness':
if selection is None:
selection = 'icc'
# Define the compactness penalty function
def compactness_penalty(*args):
P = Pmodel(*[r]*Nconstants,*args)
if not np.all(P==0):
P = P/np.trapz(P,r)
return np.sqrt(P*(r - np.trapz(P*r,r))**2*np.mean(np.diff(r)))
# Add the penalty to the Pmodel
penalty = Penalty(compactness_penalty,selection,
signature = Pmodel._parameter_list(),
description = 'Distance distribution compactness penalty.')
penalty.weight.set(lb=1e-6, ub=1e1)
# If include smoothness penalty
elif type=='smoothness':
if selection is None:
selection = 'aic'
# Define the smoothness penalty function
L = regoperator(r,2)
def smoothness_penalty(*args):
return L@Pmodel(*[r]*Nconstants,*args)
# Add the penalty to the Pmodel
penalty = Penalty(smoothness_penalty,selection,
signature = Pmodel._parameter_list(),
description = 'Distance distribution smoothness penalty.')
penalty.weight.set(lb=1e-9, ub=1e3)
else:
raise KeyError(f"The requested {type} is not a valid penalty. Must be 'compactness' or 'smoothness'.")
return penalty
#===============================================================================
#===============================================================================
def _checkpathways(pathways,Nmax):
# Check that pathways are correctly specified
if len(pathways)>Nmax:
raise ValueError(f"The number of pathways cannot exceed {Nmax}.")
if np.any(np.array(pathways)<1) or not np.all([not p%1 for p in pathways]) or np.any(np.array(pathways)>Nmax):
raise ValueError(f"The pathways must be specified by integer numbers in the range 1-{Nmax}.")
if len(np.unique(pathways))!=len(pathways):
raise ValueError("Some pathways are duplicated.")
#===============================================================================
#===============================================================================
[docs]
class ExperimentInfo():
r"""
Represents information about a dipolar EPR experiment
Attributes
----------
name : string
Name of the experiment.
reftimes : list of floats
List of refocusing times, in microseconds.
harmonics : list
List of harmonics of the refocusing times.
pulselength : float
Length of the longest microwave pulse in the sequence in microseconds.
pathwaylabels : list
List of pathway labels.
delays : list
List of pulse delays.
"""
[docs]
def __init__(self,name,reftimes,harmonics,pulselength,pathwaylabels,delays):
""" Populates the attributes of the ExperimentInfo object.
Parameters
----------
name : string
Name of the experiment.
reftimes : list of floats
List of refocusing times, in microseconds.
harmonics : list
List of harmonics of the refocusing times.
pulselength : float
Length of the longest microwave pulse in the sequence in microseconds.
pathwaylabels : list
List of pathway labels.
delays : list
List of pulse delays.
"""
self.npathways = len(harmonics)
self.reftimes = reftimes
self.harmonics = harmonics
self.pulselength = pulselength
self.labels = pathwaylabels
self.delays = delays
self.name = name
#===============================================================================
#===============================================================================
[docs]
def ex_3pdeer(tau, pathways=[1,2], pulselength=0.016):
r"""
Generate a 3-pulse DEER dipolar experiment model.
The figure below shows the dipolar pathways in 3-pulse DEER. The observer (black) and pump (grey) pulses and their interpulse delays are shown on the top.
The middle table summarizes all detectable modulated dipolar pathways `p` along their dipolar phase accumulation factors `\mathbf{s}_p`,
harmonics `\delta_p` and refocusing times `t_{\mathrm{ref},p}`. The most commonly encountered pathways are highlighted in color.
The bottom panel shows a decomposition of the dipolar signal into the individual intramolecular contributions (shown as colored lines).
.. figure:: ../images/ex_3pdeer_pathways.png
:width: 350px
Source: L. Fábregas-Ibáñez, Advanced data analysis and modeling in dipolar EPR spectroscopy, Doctoral dissertation, 2022
Parameters
----------
tau : float scalar
Static interpulse delay `\tau`.
pathways : array_like, optional
Pathways to include in the model. The pathways are specified as a list of pathways labels `p`.
By default, both pathways are included as shown in the table above.
pulselength : float scalar, optional
Length of the longest microwave pulse in the sequence in microseconds. Used to determine the uncertainty in the
boundaries of the pathway refocusing times.
Returns
-------
experiment : ``ExperimentInfo`` object
Experiment object. Can be passed to ``dipolarmodel`` to introduce better
constraints into the model.
"""
# Theoretical refocusing times
def reftimes(tau):
tref = [ tau, 0]
# Sort according to pathways order
if pathways is not None:
tref = [tref[pathway-1] for pathway in pathways]
return tref
# Pulse delays
delays = [tau]
# Theoretical dipolar harmonics
harmonics = [ 1, 1]
# Pathway labels
pathwaylabels = np.arange(1,len(harmonics)+1)
# Sort according to pathways order
if pathways is not None:
_checkpathways(pathways,Nmax=len(harmonics))
harmonics = [harmonics[pathway-1] for pathway in pathways]
pathwaylabels = [pathwaylabels[pathway-1] for pathway in pathways]
return ExperimentInfo('3-pulse DEER', reftimes, harmonics, pulselength, pathwaylabels,delays)
#===============================================================================
#===============================================================================
[docs]
def ex_4pdeer(tau1, tau2, pathways=[1,2,3,4], pulselength=0.016):
r"""
Generate a 4-pulse DEER dipolar experiment model.
The figure below shows the dipolar pathways in 4-pulse DEER. The observer (black) and pump (grey) pulses and their interpulse delays are shown on the top.
The middle table summarizes all detectable modulated dipolar pathways `p` along their dipolar phase accumulation factors `\mathbf{s}_p`,
harmonics `\delta_p` and refocusing times `t_{\mathrm{ref},p}`. The most commonly encountered pathways are highlighted in color.
The bottom panel shows a decomposition of the dipolar signal into the individual intramolecular contributions (shown as colored lines).
The dipolar time axis is defined such that `t=0` right after the second observer pulse.
.. figure:: ../images/ex_4pdeer_pathways.png
:width: 350px
Source: L. Fábregas-Ibáñez, Advanced data analysis and modeling in dipolar EPR spectroscopy, Doctoral dissertation, 2022
Parameters
----------
tau1 : float scalar
1st static interpulse delay `\tau_1`.
tau2 : float scalar
2nd static interpulse delay `\tau_2`.
pathways : array_like, optional
Pathways to include in the model. The pathways are specified as a list of pathways labels `p`.
By default, pathways 1-4 are included as shown in the table above.
pulselength : float scalar, optional
Length of the longest microwave pulse in the sequence in microseconds. Used to determine the uncertainty in the
boundaries of the pathway refocusing times.
Returns
-------
experiment : ``ExperimentInfo`` object
Experiment object. Can be passed to ``dipolarmodel`` to introduce better constraints into the model.
"""
# Theoretical refocusing times
def reftimes(tau1,tau2):
tref = [tau1,
tau1+tau2,
0,
tau2]
# Sort according to pathways order
if pathways is not None:
tref = [tref[pathway-1] for pathway in pathways]
return tref
# Pulse delays
delays = [tau1,tau2]
# Theoretical dipolar harmonics
harmonics = [1, 1, 1, 1]
# Pathway labels
pathwaylabels = np.arange(1,len(harmonics)+1)
# Sort according to pathways order
if pathways is not None:
_checkpathways(pathways,Nmax=len(harmonics))
harmonics = [harmonics[pathway-1] for pathway in pathways]
pathwaylabels = [pathwaylabels[pathway-1] for pathway in pathways]
return ExperimentInfo('4-pulse DEER', reftimes, harmonics, pulselength, pathwaylabels, delays)
#===============================================================================
#===============================================================================
[docs]
def ex_rev5pdeer(tau1, tau2, tau3, pathways=[1,2,3,4,5], pulselength=0.016):
r"""
Generate a reverse 5-pulse DEER dipolar experiment model.
The figure below shows the dipolar pathways in reverse 5-pulse DEER. The observer (black) and pump (grey) pulses and their interpulse delays are shown on the top.
The middle table summarizes all detectable modulated dipolar pathways `p` along their dipolar phase accumulation factors `\mathbf{s}_p`,
harmonics `\delta_p` and refocusing times `t_{\mathrm{ref},p}`. The most commonly encountered pathways are highlighted in color.
The bottom panel shows a decomposition of the dipolar signal into the individual intramolecular contributions (shown as colored lines).
.. figure:: ../images/ex_rev5pdeer_pathways.png
:width: 350px
Source: L. Fábregas-Ibáñez, Advanced data analysis and modeling in dipolar EPR spectroscopy, Doctoral dissertation, 2022
Parameters
----------
tau1 : float scalar
1st static interpulse delay `\tau_1`.
tau2 : float scalar
2nd static interpulse delay `\tau_2`.
tau3 : float scalar
3rd static interpulse delay `\tau_3`.
pathways : array_like, optional
Pathways to include in the model. The pathways are specified as a list of pathways labels `p`.
By default, pathways 1-5 are included as shown in the table above.
pulselength : float scalar, optional
Length of the longest microwave pulse in the sequence in microseconds. Used to determine the uncertainty in the
boundaries of the pathway refocusing times.
Returns
-------
experiment : ``ExperimentInfo`` object
Experiment object. Can be passed to ``dipolarmodel`` to introduce better
constraints into the model.
"""
# Theoretical refocusing times
def reftimes(tau1,tau2,tau3):
tref = [ tau3,
tau2-tau3,
tau1+tau3,
tau1+tau2-tau3,
tau2,
tau1+tau2,
0,
tau1]
# Sort according to pathways order
if pathways is not None:
tref = [tref[pathway-1] for pathway in pathways]
return tref
# Pulse delays
delays = [tau1,tau2,tau3]
# Theoretical dipolar harmonics
harmonics = [ 1, 1, 1, 1, 1, 1, 1, 1]
# Pathway labels
pathwaylabels = np.arange(1,len(harmonics)+1)
# Sort according to pathways order
if pathways is not None:
_checkpathways(pathways,Nmax=len(harmonics))
harmonics = [harmonics[pathway-1] for pathway in pathways]
pathwaylabels = [pathwaylabels[pathway-1] for pathway in pathways]
return ExperimentInfo('Reverse 5-pulse DEER', reftimes, harmonics, pulselength, pathwaylabels, delays)
#===============================================================================
#===============================================================================
[docs]
def ex_fwd5pdeer(tau1, tau2, tau3, pathways=[1,2,3,4,5], pulselength=0.016):
r"""
Generate a forward 5-pulse DEER dipolar experiment model.
The figure below shows the dipolar pathways in forward 5-pulse DEER. The observer (black) and pump (grey) pulses and their interpulse delays are shown on the top.
The middle table summarizes all detectable modulated dipolar pathways `p` along their dipolar phase accumulation factors `\mathbf{s}_p`,
harmonics `\delta_p` and refocusing times `t_{\mathrm{ref},p}`. The most commonly encountered pathways are highlighted in color.
The bottom panel shows a decomposition of the dipolar signal into the individual intramolecular contributions (shown as colored lines).
.. figure:: ../images/ex_fwd5pdeer_pathways.png
:width: 350px
Source: L. Fábregas-Ibáñez, Advanced data analysis and modeling in dipolar EPR spectroscopy, Doctoral dissertation, 2022
Parameters
----------
tau1 : float scalar
1st static interpulse delay `\tau_1`.
tau2 : float scalar
2nd static interpulse delay `\tau_2`.
tau3 : float scalar
3rd static interpulse delay `\tau_3`.
pathways : array_like, optional
Pathways to include in the model. The pathways are specified as a list of pathways labels `p`.
By default, pathways 1-5 are included as shown in the table above.
pulselength : float scalar, optional
Length of the longest microwave pulse in the sequence in microseconds. Used to determine the uncertainty in the
boundaries of the pathway refocusing times.
Returns
-------
experiment : ``ExperimentInfo`` object
Experiment object. Can be passed to ``dipolarmodel`` to introduce better constraints into the model.
"""
# Theoretical refocusing times
def reftimes(tau1,tau2,tau3):
tref = [ tau3,
tau1-tau3,
tau2+tau3,
tau1+tau2-tau3,
tau1,
tau1+tau2,
0,
tau2]
# Sort according to pathways order
if pathways is not None:
tref = [tref[pathway-1] for pathway in pathways]
return tref
# Pulse delays
delays = [tau1,tau2,tau3]
# Theoretical dipolar harmonics
harmonics = [1, 1, 1, 1, 1, 1, 1, 1]
# Pathway labels
pathwaylabels = np.arange(1,len(harmonics)+1)
# Sort according to pathways order
if pathways is not None:
_checkpathways(pathways,Nmax=len(harmonics))
harmonics = [harmonics[pathway-1] for pathway in pathways]
pathwaylabels = [pathwaylabels[pathway-1] for pathway in pathways]
return ExperimentInfo('Forward 5-pulse DEER', reftimes, harmonics, pulselength, pathwaylabels, delays)
#===============================================================================
#===============================================================================
[docs]
def ex_sifter(tau1, tau2, pathways=[1,2,3], pulselength=0.016):
r"""
Generate a 4-pulse SIFTER dipolar experiment model.
The figure below shows the dipolar pathways in 4-pulse SIFTER. The pulses and the interpulse delays are shown on the top.
The middle table summarizes all detectable modulated dipolar pathways `p` along their dipolar phase accumulation factors `\mathbf{s}_p`,
harmonics `\delta_p` and refocusing times `t_{\mathrm{ref},p}`. The most commonly encountered pathways are highlighted in color.
The bottom panel shows a decomposition of the dipolar signal into the individual intramolecular contributions (shown as colored lines).
.. figure:: ../images/ex_sifter_pathways.png
:width: 350px
Source: L. Fábregas-Ibáñez, Advanced data analysis and modeling in dipolar EPR spectroscopy, Doctoral dissertation, 2022
Parameters
----------
tau1 : float scalar
1st static interpulse delay `\tau_1`.
tau2 : float scalar
2nd static interpulse delay `\tau_2`.
pathways : array_like, optional
Pathways to include in the model. The pathways are specified as a list of pathways labels `p`.
By default, all 3 pathways are included as shown in the table above.
pulselength : float scalar, optional
Length of the longest microwave pulse in the sequence in microseconds. Used to determine the uncertainty in the
boundaries of the pathway refocusing times.
Returns
-------
experiment : ``ExperimentInfo`` object
Experiment object. Can be passed to ``dipolarmodel`` to introduce better constraints into the model.
"""
# Theoretical refocusing times
def reftimes(tau1,tau2):
tref = [ tau2-tau1,
2*tau2,
-2*tau1]
# Sort according to pathways order
if pathways is not None:
tref = [tref[pathway-1] for pathway in pathways]
return tref
# Pulse delays
delays = [tau1,tau2]
# Theoretical dipolar harmonics
harmonics = [ 1, 1/2, 1/2]
# Pathway labels
pathwaylabels = np.arange(1,len(harmonics)+1)
# Sort according to pathways order
if pathways is not None:
_checkpathways(pathways,Nmax=len(harmonics))
harmonics = [harmonics[pathway-1] for pathway in pathways]
pathwaylabels = [pathwaylabels[pathway-1] for pathway in pathways]
return ExperimentInfo('4-pulse SIFTER', reftimes, harmonics, pulselength, pathwaylabels, delays)
#===============================================================================
#===============================================================================
[docs]
def ex_ridme(tau1, tau2, pathways=[1,2,3,4], pulselength=0.016):
r"""
Generate a 5-pulse RIDME dipolar experiment model.
The figure below shows the dipolar pathways in 5-pulse RIDME. The pulses and their interpulse delays are shown on the top.
The middle table summarizes all detectable modulated dipolar pathways `p` along their dipolar phase accumulation factors `\mathbf{s}_p`,
harmonics `\delta_p` and refocusing times `t_{\mathrm{ref},p}`. The most commonly encountered pathways are highlighted in color.
The bottom panel shows a decomposition of the dipolar signal into the individual intramolecular contributions (shown as colored lines).
The dipolar time axis is defined such that `t=0` right after the second pulse.
(Note that the model does not account for any relaxation-induced effects)
.. figure:: ../images/ex_ridme_pathways.png
:width: 350px
Parameters
----------
tau1 : float scalar
1st static interpulse delay `\tau_1`.
tau2 : float scalar
2nd static interpulse delay `\tau_2`.
pathways : array_like, optional
Pathways to include in the model. The pathways are specified as a list of pathways labels `p`.
By default, pathways 1-4 are included as shown in the table above.
pulselength : float scalar, optional
Length of the longest microwave pulse in the sequence, in microseconds. Used to determine
the uncertainty in the boundaries of the pathway refocusing times.
Returns
-------
experiment : ``ExperimentInfo`` object
Experiment object. Can be passed to ``dipolarmodel`` to introduce better
constraints into the model.
"""
# Theoretical refocusing times
def reftimes(tau1, tau2):
tref = [tau1,
tau1+tau2,
0,
tau2]
# Sort according to pathways order
if pathways is not None:
tref = [tref[pathway-1] for pathway in pathways]
return tref
# Pulse delays
delays = [tau1, tau2]
# Theoretical dipolar harmonics
harmonics = [1, 1, 1, 1]
# Pathway labels
pathwaylabels = np.arange(1, len(harmonics)+1)
# Sort according to pathways order
if pathways is not None:
_checkpathways(pathways, Nmax=len(harmonics))
harmonics = [harmonics[pathway-1] for pathway in pathways]
pathwaylabels = [pathwaylabels[pathway-1] for pathway in pathways]
return ExperimentInfo('5-pulse RIDME', reftimes, harmonics, pulselength, pathwaylabels, delays)
#===============================================================================
#===============================================================================
[docs]
def ex_dqc(tau1, tau2, tau3, pathways=[1,2,3], pulselength=0.016):
r"""
Generate a 6-pulse DQC dipolar experiment model.
The figure below shows the dipolar pathways in 6-pulse DQC. The pulses and their interpulse delays are shown on the top.
The middle table summarizes all detectable modulated dipolar pathways `p` along their dipolar phase accumulation factors `\mathbf{s}_p`,
harmonics `\delta_p` and refocusing times `t_{\mathrm{ref},p}`. The most commonly encountered pathways are highlighted in color.
The bottom panel shows a decomposition of the dipolar signal into the individual intramolecular contributions (shown as colored lines).
.. figure:: ../images/ex_dqc_pathways.png
:width: 700px
Source: L. Fábregas-Ibáñez, Advanced data analysis and modeling in dipolar EPR spectroscopy, Doctoral dissertation, 2022
Parameters
----------
tau1 : float scalar
1st static interpulse delay `\tau_1`.
tau2 : float scalar
2nd static interpulse delay `\tau_2`.
tau3 : float scalar
3rd static interpulse delay `\tau_3`.
pathways : array_like, optional
Pathways to include in the model. The pathways are specified as a list of pathways labels `p`.
By default, pathways 1, 2, and 3 are included as shown in the table above.
pulselength : float scalar, optional
Length of the longest microwave pulse in the sequence in microseconds. Used to determine the uncertainty in the
boundaries of the pathway refocusing times.
Returns
-------
experiment : ``ExperimentInfo`` object
Experiment object. Can be passed to ``dipolarmodel`` to introduce better constraints into the model.
"""
# Theoretical refocusing times
def reftimes(tau1,tau2,tau3):
tref = [ tau2-tau1,
2*tau2,
-2*tau1,
tau2-tau1-tau3,
-2*(tau1-tau3),
-2*(tau1+tau3),
2*(tau2+tau3),
2*(tau1-tau3)]
# Sort according to pathways order
if pathways is not None:
tref = [tref[pathway-1] for pathway in pathways]
return tref
# Pulse delays
delays = [tau1,tau2,tau3]
# Theoretical dipolar harmonics
harmonics = [ 1,
1/2,
1/2,
1,
1/2,
1/2,
1/2,
1/2]
# Pathway labels
pathwaylabels = np.arange(1,len(harmonics)+1)
# Sort according to pathways order
if pathways is not None:
_checkpathways(pathways,Nmax=len(harmonics))
harmonics = [harmonics[pathway-1] for pathway in pathways]
pathwaylabels = [pathwaylabels[pathway-1] for pathway in pathways]
return ExperimentInfo('6-pulse DQC', reftimes, harmonics, pulselength, pathwaylabels, delays)
#===============================================================================