Dipolarbackgroundmodel func

deerlab/__init__.py

@@ -5,7 +5,7 @@
 from .deerload import deerload
 from .selregparam import selregparam
 from .dipolarkernel import dipolarkernel
-from .dipolarbackground import dipolarbackground
+from .dipolarbackground import dipolarbackground,dipolarbackgroundmodel
 from .dipolarmodel import dipolarmodel,ExperimentInfo, dipolarpenalty, ex_4pdeer,ex_3pdeer,ex_rev5pdeer,ex_fwd5pdeer,ex_ridme,ex_sifter,ex_dqc
 from .solvers import snlls, fnnls, cvxnnls
 from .regoperator import regoperator

deerlab/bg_models.py

@@ -7,11 +7,11 @@
 import math as m
 from numpy import pi
 import inspect
+from scipy.special import gamma, hyp2f1, sici
+from deerlab.constants import *
 from deerlab.dipolarkernel import dipolarkernel
 from deerlab.utils import formatted_table
 from deerlab.model import Model
-from scipy.special import gamma, hyp2f1, sici
-from deerlab.constants import *
 
 #---------------------------------------------------------------------------------------
 def hyp2f1_repro(a,b,c,z): 
@@ -513,3 +513,5 @@ def _poly3(t,p0,p1,p2,p3):
 bg_poly3.p3.set(description='3rd order weight', lb=-200, ub=200, par0=-1, unit=r'μs\ :sup:`-3`')
 # Add documentation
 bg_poly3.__doc__ = _docstring(bg_poly3,notes)
+
+

deerlab/dipolarbackground.py

@@ -1,7 +1,7 @@
 # dipolarbackground.py - Multipathway background generator
 # --------------------------------------------------------
 # This file is a part of DeerLab. License is MIT (see LICENSE.md).
-# Copyright(c) 2019-2023: Luis Fabregas, Stefan Stoll and other contributors.
+# Copyright(c) 2019-2025: Luis Fabregas, Stefan Stoll, Hugo Karas and other contributors.
 
 import numpy as np
 import inspect
@@ -206,3 +206,111 @@ def basis(t,lam):
                 B *= basis(tdip)
 
     return B
+
+def dipolarbackgroundmodel(experiment,basis=None,parametrization='reftimes'):
+    """
+    Construct a dipolar background model for a given dipolar experiment. Currently limited to two-spin systems.
+    This model can be used for evaluating and extrapolating the fitted parameters of the dipolar background.
+
+    Parameters
+    ----------
+    experiment : :ref:`ExperimentInfo`
+        The Experimental information obtained from experiment models (ex_).
+
+    basis : :ref:`Model`, optional
+        The basis model for the intermolecular (background) contribution. If not specified, a background arising from a homogenous 3D distribution of spins is used.
+
+    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'``. 
+
+    Returns
+    -------
+
+    BModel : :ref:`Model`
+
+    Examples
+    --------
+    To construct a dipolar background model from the fit results using the specified :ref:`ExperimentInfo`::
+
+            Bfcn = dl.dipolarbackgroundmodel(experimentInfo)
+            Bfit = results.P_scale*results.evaluate(Bfcn,t)
+            Bci = results.P_scale*results.propagate(Bfcn,t).ci(95)
+    """
+    from deerlab.dipolarmodel import _populate_dipolar_pathways_parameters
+    from deerlab.model import Model
+    if basis is None:
+        from deerlab.bg_models import bg_hom3d
+        basis = bg_hom3d
+    npathways = experiment.npathways
+
+    signature = []
+    basis_signature = basis.signature.copy()
+    lam_in_basis = 'lam' in basis_signature
+    if lam_in_basis:
+        basis_signature.remove('lam')
+    signature.extend(basis_signature)
+
+    if parametrization=='reftimes':
+        if npathways==1:
+            signature.extend(['mod','reftime'])
+        else:
+            signature.extend(['lam'+str(i) for i in experiment.labels])
+            signature.extend(['reftime'+str(i) for i in experiment.labels])
+
+    elif parametrization=='delays':
+        signature.extend(['delay'+str(i+1) for i in range(len(experiment.delays))])
+        signature.extend(['lam'+str(i) for i in experiment.labels])
+
+    elif parametrization=='shift':
+        signature.extend(['shift'])
+        signature.extend(['lam'+str(i+1) for i in experiment.labels])
+
+
+    # Construct the dipolar background model
+    def bckgrnd_fun(*param:tuple):
+        basis_constants = list(param[0:len(basis_signature)])
+        t = basis_constants[0] # t is always the first basis parameter
+        pathway_params = np.array(param[len(basis_signature):])
+        if parametrization=='reftimes':
+            λs = pathway_params[np.arange(0, len(pathway_params)//2, 1,)]
+            trefs = pathway_params[np.arange(len(pathway_params)//2,len(pathway_params),1)]
+        elif parametrization=='delays':
+            delays = pathway_params[0:len(experiment.delays)]
+            λs = pathway_params[len(delays):]
+            trefs = experiment.reftimes(*delays)
+        elif parametrization=='shift':
+            shift = pathway_params[0]
+            λs = pathway_params[1:]
+            delays = np.array(experiment.delays) + shift
+            trefs = experiment.reftimes(*delays)
+
+        if lam_in_basis:
+                basis_constants.append(λs[0])
+
+        prod = np.ones_like(t)
+        scale = 1
+
+        for i in range(npathways):
+            scale -= λs[i]
+            if lam_in_basis:
+                basis_constants[-1] = λs[i]
+            prod *= basis(t-trefs[i],*basis_constants[1:])
+        return scale*prod
+
+    # Define the dipolar background model
+    bckgrnd_model = Model(bckgrnd_fun, signature=signature,constants=['t'])
+    bckgrnd_model = _populate_dipolar_pathways_parameters(bckgrnd_model,npathways,experiment=experiment,parametrization=parametrization)
+
+    # Copy the basis parameters to the background model
+    for param in basis_signature:
+        if param =='t':
+            continue
+        getattr(bckgrnd_model,param).setas(getattr(basis,param))
+
+    return bckgrnd_model

deerlab/dipolarmodel.py

@@ -4,7 +4,8 @@
 # Copyright(c) 2019-2022: Luis Fabregas, Stefan Stoll and other contributors.
 
 import numpy as np
-from deerlab.dipolarkernel import dipolarkernel, dipolarbackground
+from deerlab.dipolarkernel import dipolarkernel
+from deerlab.dipolarbackground import dipolarbackground
 from deerlab.regoperator import regoperator
 from deerlab.dd_models import freedist
 from deerlab.model import Model,Penalty, link
@@ -16,7 +17,97 @@
 import inspect
 from scipy.stats import multivariate_normal
 
+
 #===============================================================================
+def _populate_dipolar_pathways_parameters(model,npathways,spins=2,samespins=True,experiment=None,parametrization='reftimes',Q=None):
+    """
+    A private function to populate the dipolar pathways parameters in the model.
+
+    Parameters
+    ----------
+    model : :ref:`Model`
+        Model object to be populated with the dipolar pathways parameters.
+    npathways : int
+        Number of dipolar pathways.
+    spins : int
+        Number of spins in the system.
+    samespins : bool    
+        Enables the assumption of spectral permutability.
+    experiment : :ref:`ExperimentInfo`
+        Experimental information obtained from experiment models (``ex_``).
+    parametrization : str
+        Parametrization strategy of the dipolar pathway refocusing times.
+    Q : int
+        Number of interactions in the spins system.
+
+    Returns
+    -------
+    model : :ref:`Model`
+        Model object populated with the dipolar pathways parameters.
+    """
+
+    # Populate the basic information on the dipolar pathways parameters
+    if experiment is not None:
+        labels = experiment.labels
+        pulsedelay_names = inspect.signature(experiment.reftimes).parameters.keys()
+    else:
+        # Otherwise, just label them in order
+        labels = np.arange(1,npathways+1)
+
+    if spins>2:
+        if Q is None:
+            raise ValueError('If spins>2, the number of interactions Q must be specified.')
+
+
+    for n in range(npathways):
+        if spins==2 and npathways==1:
+            # Special case: use modulation depth notation instead of general pathway amplitude
+            getattr(model,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(model,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(model,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(model,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(model,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(model,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(model,'tshift').set(par0=0,lb=-variability,ub=variability,description=f'Variability of experimental pulse delays',unit='μs')
+    if spins>2: 
+        getattr(model,f'lamu').set(lb=0,ub=1,par0=0.5,description='Amplitude of unmodulated pairwise pathway',unit='')
+
+    return model
+
+
 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,):
     """
@@ -342,53 +433,9 @@ def _Pmultivar(param):
             [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='')
-
+    PathsModel = _populate_dipolar_pathways_parameters(
+        PathsModel,npathways,spins=spins,samespins=samespins,experiment=experiment,parametrization=parametrization,Q=Q)
+
     # Construct the signature of the dipolar signal model function
     signature = []
     parameters,linearparam = [],[]

docsrc/source/reference.rst

@@ -60,6 +60,7 @@ Reference Index
     deerload
     dipolarkernel 
     dipolarbackground
+    dipolarbackgroundmodel
     fftspec
     distancerange
 

examples/basic/ex_bootstrapping.py

@@ -45,7 +45,8 @@
 r = np.linspace(2,5,100) # nm
 
 # Construct the model
-Vmodel = dl.dipolarmodel(t,r, experiment = dl.ex_4pdeer(tau1,tau2, pathways=[1]))
+experimentInfo = dl.ex_4pdeer(tau1,tau2, pathways=[1])
+Vmodel = dl.dipolarmodel(t,r, experiment = experimentInfo)
 
 # Fit the model to the data
 results = dl.fit(Vmodel,Vexp,bootstrap=20)
@@ -67,9 +68,9 @@
 Pci50 = results.PUncert.ci(50)
 
 # Extract the unmodulated contribution
-Bfcn = lambda mod,conc,reftime: results.P_scale*(1-mod)*dl.bg_hom3d(t-reftime,conc,mod)
-Bfit = results.evaluate(Bfcn)
-Bci = results.propagate(Bfcn).ci(95)
+Bfcn = dl.dipolarbackgroundmodel(experimentInfo)
+Bfit = results.P_scale*results.evaluate(Bfcn,t)
+Bci = results.P_scale*results.propagate(Bfcn,t).ci(95)
 
 plt.figure(figsize=[6,7])
 violet = '#4550e6'

examples/basic/ex_fitting_4pdeer.py

@@ -35,7 +35,8 @@
 r = np.arange(2.5,5,0.01) # nm
 
 # Construct the model
-Vmodel = dl.dipolarmodel(t,r, experiment = dl.ex_4pdeer(tau1,tau2, pathways=[1]))
+experimentInfo = dl.ex_4pdeer(tau1,tau2, pathways=[1])
+Vmodel = dl.dipolarmodel(t,r, experiment = experimentInfo)
 
 # Fit the model to the data
 results = dl.fit(Vmodel,Vexp)
@@ -54,9 +55,9 @@
 Pci50 = results.PUncert.ci(50)
 
 # Extract the unmodulated contribution
-Bfcn = lambda mod,conc,reftime: results.P_scale*(1-mod)*dl.bg_hom3d(t-reftime,conc,mod)
-Bfit = results.evaluate(Bfcn)
-Bci = results.propagate(Bfcn).ci(95)
+Bfcn = dl.dipolarbackgroundmodel(experimentInfo)
+Bfit = results.P_scale*results.evaluate(Bfcn,t)
+Bci = results.P_scale*results.propagate(Bfcn,t).ci(95)
 
 plt.figure(figsize=[6,7])
 violet = '#4550e6'

examples/basic/ex_fitting_4pdeer_compactness.py

@@ -38,7 +38,8 @@
 r = np.arange(2,5,0.025) # nm
 
 # Construct the model
-Vmodel = dl.dipolarmodel(t, r, experiment = dl.ex_4pdeer(tau1,tau2, pathways=[1]))
+experimentInfo = dl.ex_4pdeer(tau1,tau2, pathways=[1])
+Vmodel = dl.dipolarmodel(t, r, experiment = experimentInfo)
 compactness = dl.dipolarpenalty(Pmodel=None, r=r, type='compactness')
 
 # Fit the model to the data
@@ -58,9 +59,9 @@
 Pci50 = results.PUncert.ci(50)
 
 # Extract the unmodulated contribution
-Bfcn = lambda mod,conc,reftime: results.P_scale*(1-mod)*dl.bg_hom3d(t-reftime,conc,mod)
-Bfit = results.evaluate(Bfcn)
-Bci = results.propagate(Bfcn).ci(95)
+Bfcn = dl.dipolarbackgroundmodel(experimentInfo)
+Bfit = results.P_scale*results.evaluate(Bfcn,t)
+Bci = results.P_scale*results.propagate(Bfcn,t).ci(95)
 
 plt.figure(figsize=[6,7])
 violet = '#4550e6'