Pathways interface simplification

deerlab/dipolarbackground.py

@@ -1,4 +1,3 @@
-
 # dipolarbackground.py - Multipathway background generator 
 # --------------------------------------------------------
  # This file is a part of DeerLab. License is MIT (see LICENSE.md).
@@ -7,20 +6,25 @@
 import numpy as np
 import types
 
-def dipolarbackground(t,pathinfo,Bmodel,renormalize = True,overtonecoeff = 1):
+def dipolarbackground(t, pathinfo, 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``.
+    pathinfo : 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
@@ -30,8 +34,6 @@ 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
     -----
@@ -48,73 +50,58 @@ 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
+        pathinfo = []
+        pathinfo.append([1-lam]) # unmodulated part, gives offset
+        pathinfo.append([lam, 0]) # main modulation, refocusing at time zero
 
         def Bmodel(t,lam):
             return dl.bg_hom3d(t,conc,lam)
 
         B = dl.dipolarbackground(t,pathinfo,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(pathinfo,list): pathinfo = [pathinfo] 
     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]')
+        pathinfo = [[1-lam], [lam, 0]]
 
-    # 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]
+    pathinfo = [np.atleast_1d(path) for path in pathinfo]
 
+    # Get unmodulated pathways    
+    unmodulated = [pathinfo.pop(i) for i,path in enumerate(pathinfo) 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(pathinfo):
+        if len(path) == 2:
+            # If harmonic is not defined, append default n=1
+            pathinfo[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 pathinfo])
+        Lambda0 /= lamsum
+        for i in range(len(pathinfo)):
+            pathinfo[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 pathinfo:        
+        n,T0,lam = pathway
+        B = B*Bmodel(n*(t-T0),lam)
+        Bnorm = Bnorm*Bmodel(-T0*n,lam)
 
     if renormalize:
         B = B/Bnorm

deerlab/dipolarkernel.py

@@ -24,7 +24,7 @@ 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):
+                  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. 
@@ -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``.
+    pathinfo : 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)``
@@ -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
     --------
@@ -81,8 +84,10 @@ 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]]
-
+        pathinfo = []
+        pathinfo.append([1-lam])      # unmodulated part, gives offset
+        pathinfo.append([lam, 0])     # main modulation, refocusing at time zero
+
         K = dl.dipolarkernel(t,r,pathinfo)
 
 
@@ -97,22 +102,22 @@ 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
+        pathinfo = []
+        pathinfo.append([Lam0])           # unmodulated part, gives offset
+        pathinfo.append([lam1,   0  ])    # main modulation, refocusing at time zero
+        pathinfo.append([lam2,  tau2])    # 2+1 modulation, refocusing at time tau2 
 
         K = dl.dipolarkernel(t,r,pathinfo)
 
 
     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
     """
@@ -123,7 +128,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
@@ -140,51 +144,51 @@ 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(pathinfo,list): pathinfo = [pathinfo] 
     if len(pathinfo) == 1:
         lam = pathinfo[0]
-        pathinfo = np.array([[1-lam, np.NaN], [lam, 0]])
+        pathinfo = [[1-lam], [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.')
+    paths = [np.atleast_1d(path) for path in pathinfo]
 
-    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]')
+    # 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]))
+
+    # 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)

deerlab/ex_models.py

@@ -71,7 +71,7 @@ def ex_4pdeer(param=None):
 
     # Dipolar pathways
     lam = param[0]
-    pathways = np.array([[1-lam, np.NaN], [lam, 0]])
+    pathways = [[1-lam], [lam, 0]]
 
     return pathways
 # ===================================================================
@@ -137,7 +137,7 @@ def ex_ovl4pdeer(param=None):
     # Dipolar pathways
     lam = param[[0,1,2]]
     T02 = param[3]
-    pathways = np.array([[lam[0], np.NaN], [lam[1], 0], [lam[2], T02]])
+    pathways = [[lam[0]], [lam[1], 0], [lam[2], T02]]
 
     return pathways
 # ===================================================================
@@ -203,7 +203,7 @@ def ex_5pdeer(param=None):
     # Dipolar pathways
     lam = param[[0,1,2]]
     T02 = param[3]
-    pathways = np.array([[lam[0], np.NaN], [lam[1], 0], [lam[2], T02]])
+    pathways = [[lam[0]], [lam[1], 0], [lam[2], T02]]
 
     return pathways
 # ===================================================================
@@ -276,7 +276,7 @@ def ex_7pdeer(param=None):
     T02 = param[4]
     T03 = param[5]
 
-    pathways = np.array([[lam[0], np.NaN], [lam[1], 0], [lam[2], T02], [lam[3], T03]])
+    pathways = [[lam[0]], [lam[1], 0], [lam[2], T02], [lam[3], T03]]
 
     return pathways
 # ===================================================================

docsrc/source/basics.rst

@@ -87,10 +87,9 @@ The returned output is
 
 .. code-block:: python
 
-    pathways =[[0.7,  NaN],
-               [0.3,   0 ]]
+    pathways =[[0.7], [0.3,   0 ]]
 
-Each row of this array holds information about one pathway. The first column is modulation amplitude, and the second column is the refocusing point. In the above example, the first row shows a pathway with amplitude 0.7 and no refocusing time, indicating that it represents the unmodulated contribution. The pathway of the second row shows amplitude of 0.3 and refocusing time 0, i.e. this is the primary dipolar pathway.
+Each nested list holds information about one pathway. The first element is modulation amplitude, and the second element is the refocusing point. In the above example, the first list shows a pathway with amplitude 0.7 and no refocusing time, indicating that it represents the unmodulated contribution. The pathway of the second list shows amplitude of 0.3 and refocusing time 0, i.e. this is the primary dipolar pathway.
 
 Kernel matrices
 --------------------------------------------