From eee415e9075ef57ef614981fe033ea37cf2cdbb0 Mon Sep 17 00:00:00 2001
From: Hugo Karas <hkaras@ethz.ch>
Date: Fri, 10 Jan 2025 16:25:51 +0100
Subject: [PATCH 1/9] Initial push of the new `dipolarbackgroundmodel`function

---
 deerlab/bg_models.py         |   6 +-
 deerlab/dipolarbackground.py | 104 ++++++++++++++++++++++++-
 deerlab/dipolarmodel.py      | 143 +++++++++++++++++++++++------------
 3 files changed, 202 insertions(+), 51 deletions(-)

diff --git a/deerlab/bg_models.py b/deerlab/bg_models.py
index 037f61dcd..107df31a3 100644
--- a/deerlab/bg_models.py
+++ b/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)
+
+
diff --git a/deerlab/dipolarbackground.py b/deerlab/dipolarbackground.py
index 908927cd3..e3bdb9adf 100644
--- a/deerlab/dipolarbackground.py
+++ b/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,105 @@ def basis(t,lam):
                 B *= basis(tdip)
 
     return B
+
+def dipolarbackgroundmodel(experiment,basis=None,parametrization='reftimes'):
+    """
+    
+    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
+    --------
+
+    """
+    from deerlab.dipolarmodel import _populate_dipolar_pathways_parameters
+
+    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
diff --git a/deerlab/dipolarmodel.py b/deerlab/dipolarmodel.py
index 80035e5a6..79fcd3b06 100644
--- a/deerlab/dipolarmodel.py
+++ b/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 = [],[]

From fe4fc17af012bd036a4710059af1ee0331d9e8ca Mon Sep 17 00:00:00 2001
From: Hugo Karas <hkaras@ethz.ch>
Date: Fri, 10 Jan 2025 16:44:33 +0100
Subject: [PATCH 2/9] Documentation update and exampels

---
 deerlab/__init__.py                                  |  2 +-
 deerlab/dipolarbackground.py                         | 12 +++++++++---
 docsrc/source/reference.rst                          |  1 +
 examples/basic/ex_bootstrapping.py                   |  9 +++++----
 examples/basic/ex_fitting_4pdeer.py                  |  9 +++++----
 examples/basic/ex_fitting_4pdeer_compactness.py      |  9 +++++----
 examples/basic/ex_fitting_4pdeer_gauss.py            |  9 +++++----
 examples/basic/ex_fitting_5pdeer.py                  |  6 +++---
 examples/basic/ex_fitting_ridme.py                   | 10 +++++-----
 examples/basic/ex_restraints_4pdeer.py               |  3 ++-
 examples/intermediate/ex_compactness_with_without.py | 10 ++++++----
 examples/intermediate/ex_fitting_4pdeer_pathways.py  |  6 +++---
 examples/intermediate/ex_fitting_5pdeer_pathways.py  |  6 +++---
 examples/intermediate/ex_fitting_dqc_pathways.py     |  6 +++---
 14 files changed, 56 insertions(+), 42 deletions(-)

diff --git a/deerlab/__init__.py b/deerlab/__init__.py
index 8a8d26e77..7f3b6f580 100644
--- a/deerlab/__init__.py
+++ b/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
diff --git a/deerlab/dipolarbackground.py b/deerlab/dipolarbackground.py
index e3bdb9adf..cb7b583c3 100644
--- a/deerlab/dipolarbackground.py
+++ b/deerlab/dipolarbackground.py
@@ -209,7 +209,9 @@ def basis(t,lam):
 
 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`
@@ -234,10 +236,14 @@ def dipolarbackgroundmodel(experiment,basis=None,parametrization='reftimes'):
 
     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
diff --git a/docsrc/source/reference.rst b/docsrc/source/reference.rst
index afda906d5..e5a6546ff 100644
--- a/docsrc/source/reference.rst
+++ b/docsrc/source/reference.rst
@@ -60,6 +60,7 @@ Reference Index
     deerload
     dipolarkernel 
     dipolarbackground
+    dipolarbackgroundmodel
     fftspec
     distancerange
 
diff --git a/examples/basic/ex_bootstrapping.py b/examples/basic/ex_bootstrapping.py
index 196bd772e..642a03277 100644
--- a/examples/basic/ex_bootstrapping.py
+++ b/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'
diff --git a/examples/basic/ex_fitting_4pdeer.py b/examples/basic/ex_fitting_4pdeer.py
index 38368f5eb..ea29e35bd 100644
--- a/examples/basic/ex_fitting_4pdeer.py
+++ b/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'
diff --git a/examples/basic/ex_fitting_4pdeer_compactness.py b/examples/basic/ex_fitting_4pdeer_compactness.py
index b94f782bb..863eed1f4 100644
--- a/examples/basic/ex_fitting_4pdeer_compactness.py
+++ b/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'
diff --git a/examples/basic/ex_fitting_4pdeer_gauss.py b/examples/basic/ex_fitting_4pdeer_gauss.py
index 225827c22..74f88386a 100644
--- a/examples/basic/ex_fitting_4pdeer_gauss.py
+++ b/examples/basic/ex_fitting_4pdeer_gauss.py
@@ -37,7 +37,8 @@
 
 # Construct the model
 Pmodel= dl.dd_gauss2
-Vmodel = dl.dipolarmodel(t,r,Pmodel, experiment=dl.ex_4pdeer(tau1,tau2, pathways=[1]))
+experimentInfo = dl.ex_4pdeer(tau1,tau2, pathways=[1])
+Vmodel = dl.dipolarmodel(t,r,Pmodel, experiment=experimentInfo)
 
 # Fit the model to the data
 results = dl.fit(Vmodel,Vexp,reg=False)
@@ -59,9 +60,9 @@
 Pci50 = Puncert.ci(50)/scale
 
 # Extract the unmodulated contribution
-Bfcn = lambda mod,conc,reftime: 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 = scale*results.evaluate(Bfcn,t)
+Bci = scale*results.propagate(Bfcn,t).ci(95)
 
 plt.figure(figsize=[6,7])
 violet = '#4550e6'
diff --git a/examples/basic/ex_fitting_5pdeer.py b/examples/basic/ex_fitting_5pdeer.py
index 432d31292..39a056e06 100644
--- a/examples/basic/ex_fitting_5pdeer.py
+++ b/examples/basic/ex_fitting_5pdeer.py
@@ -67,9 +67,9 @@
 Pfit =  Pfit
 
 # Extract the unmodulated contribution
-Bfcn = lambda lam1,lam5,reftime1,reftime5,conc: results.P_scale*(1-lam1-lam5)*dl.bg_hom3d(t-reftime1,conc,lam1)*dl.bg_hom3d(t-reftime5,conc,lam5)
-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'
diff --git a/examples/basic/ex_fitting_ridme.py b/examples/basic/ex_fitting_ridme.py
index f58ccfa6b..56cc75a7f 100644
--- a/examples/basic/ex_fitting_ridme.py
+++ b/examples/basic/ex_fitting_ridme.py
@@ -36,8 +36,8 @@
 r = np.linspace(1.5,6,50) # nm
 
 # Construct the model
-experimentmodel = dl.ex_ridme(tau1,tau2, pathways=[1])
-Vmodel = dl.dipolarmodel(t,r,Bmodel=dl.bg_strexp, experiment =experimentmodel)
+experimentInfo = dl.ex_ridme(tau1,tau2, pathways=[1])
+Vmodel = dl.dipolarmodel(t,r,Bmodel=dl.bg_strexp, experiment =experimentInfo)
 
 # Fit the model to the data
 results = dl.fit(Vmodel,Vexp)
@@ -56,9 +56,9 @@
 Pci50 = results.PUncert.ci(50)
 
 # Extract the unmodulated contribution
-Bfcn = lambda mod,decay,stretch,reftime: results.P_scale*(1-mod)*dl.bg_strexp(t-reftime,decay,stretch)
-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'
diff --git a/examples/basic/ex_restraints_4pdeer.py b/examples/basic/ex_restraints_4pdeer.py
index d58675f37..717ab2e26 100644
--- a/examples/basic/ex_restraints_4pdeer.py
+++ b/examples/basic/ex_restraints_4pdeer.py
@@ -42,7 +42,8 @@
 r = np.arange(2,6,0.05) # nm
 
 # Construct dipolar 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
 fit = dl.fit(Vmodel,Vexp)
diff --git a/examples/intermediate/ex_compactness_with_without.py b/examples/intermediate/ex_compactness_with_without.py
index 0b836934d..39a00d72d 100644
--- a/examples/intermediate/ex_compactness_with_without.py
+++ b/examples/intermediate/ex_compactness_with_without.py
@@ -40,7 +40,8 @@
 r = np.arange(1.5,7,0.05) # 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 with compactness criterion
@@ -69,9 +70,10 @@
     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.subplot(2,2,n+1)
     # Plot experimental and fitted data
diff --git a/examples/intermediate/ex_fitting_4pdeer_pathways.py b/examples/intermediate/ex_fitting_4pdeer_pathways.py
index cef4d058f..b440679f7 100644
--- a/examples/intermediate/ex_fitting_4pdeer_pathways.py
+++ b/examples/intermediate/ex_fitting_4pdeer_pathways.py
@@ -34,8 +34,8 @@
 r = np.arange(2.5,5.5,0.025) # nm
 
 # Construct the model
-experiment = dl.ex_4pdeer(tau1,tau2, pathways=[1,2,3])
-Vmodel = dl.dipolarmodel(t,r,experiment=experiment)
+experimentInfo = dl.ex_4pdeer(tau1,tau2, pathways=[1,2,3])
+Vmodel = dl.dipolarmodel(t,r,experiment=experimentInfo)
 
 # Fit the model to the data
 results = dl.fit(Vmodel,Vexp)
@@ -70,7 +70,7 @@
 lams = [results.lam1, results.lam2, results.lam3]
 reftimes = [results.reftime1, results.reftime2, results.reftime3]
 colors= ['tab:blue',green, red] 
-Vinter = results.P_scale*(1-np.sum(lams))*np.prod([dl.bg_hom3d(t-reftime,results.conc,lam) for lam,reftime in zip(lams,reftimes)],axis=0)
+Vinter = results.P_scale*results.evaluate(dl.dipolarbackgroundmodel(experimentInfo),t)
 for n,(lam,reftime,color) in enumerate(zip(lams,reftimes,colors)):
     Vpath = (1-np.sum(lams) + lam*dl.dipolarkernel(t-reftime,r)@Pfit)*Vinter
     plt.plot(t,Vpath,linewidth=3,label=f'Pathway #{n+1}',color=color)
diff --git a/examples/intermediate/ex_fitting_5pdeer_pathways.py b/examples/intermediate/ex_fitting_5pdeer_pathways.py
index 516e50dce..b1d776ebb 100644
--- a/examples/intermediate/ex_fitting_5pdeer_pathways.py
+++ b/examples/intermediate/ex_fitting_5pdeer_pathways.py
@@ -34,8 +34,8 @@
 r = np.arange(2.5,5.5,0.025) # nm
 
 # Construct the model
-experiment = dl.ex_rev5pdeer(tau1,tau2,tau3, pathways=[1,2,3,5])
-Vmodel = dl.dipolarmodel(t,r,experiment=experiment)
+experimentInfo = dl.ex_rev5pdeer(tau1,tau2,tau3, pathways=[1,2,3,5])
+Vmodel = dl.dipolarmodel(t,r,experiment=experimentInfo)
 
 # Fit the model to the data
 results = dl.fit(Vmodel,Vexp)
@@ -72,7 +72,7 @@
 lams = [results.lam1, results.lam2, results.lam3, results.lam5]
 reftimes = [results.reftime1, results.reftime2, results.reftime3, results.reftime5]
 colors= ['tab:blue','tab:orange', red, green] 
-Vinter = results.P_scale*(1-np.sum(lams))*np.prod([dl.bg_hom3d(t-reftime,results.conc,lam) for lam,reftime in zip(lams,reftimes)],axis=0)
+Vinter = results.P_scale*results.evaluate(dl.dipolarbackgroundmodel(experimentInfo),t)
 for (lam,reftime,color,label) in zip(lams,reftimes,colors,labels):
     Vpath = (1-np.sum(lams) + lam*dl.dipolarkernel(t-reftime,r)@Pfit)*Vinter
     plt.plot(t,Vpath,linewidth=3,label=f'Pathway #{label}',color=color)
diff --git a/examples/intermediate/ex_fitting_dqc_pathways.py b/examples/intermediate/ex_fitting_dqc_pathways.py
index 188ba52b5..56c680d24 100644
--- a/examples/intermediate/ex_fitting_dqc_pathways.py
+++ b/examples/intermediate/ex_fitting_dqc_pathways.py
@@ -37,8 +37,8 @@
 
 # Construct the model
 r = np.arange(2.5,4,0.01) # nm
-experiment = dl.ex_dqc(tau1,tau2,tau3,pathways=[1,2,3])
-Vmodel =  dl.dipolarmodel(t,r,experiment=experiment)
+experimentInfo = dl.ex_dqc(tau1,tau2,tau3,pathways=[1,2,3])
+Vmodel =  dl.dipolarmodel(t,r,experiment=experimentInfo)
 
 # The amplitudes of the second and third pathways must be equal
 Vmodel = dl.link(Vmodel,lam23=['lam2','lam3'])
@@ -76,7 +76,7 @@
 
 # Plot the individual pathway contributions
 plt.subplot(223)
-Vinter = results.P_scale*(1-np.sum(lams))*np.prod([dl.bg_hom3d(tfull-reftime,results.conc,lam) for lam,reftime in zip(lams,reftimes)],axis=0)
+Vinter = results.P_scale*results.evaluate(dl.dipolarbackgroundmodel(experimentInfo),t)
 for n,(lam,reftime,color) in enumerate(zip(lams,reftimes,colors)):
     Vpath = (1-np.sum(lams) + lam*dl.dipolarkernel(tfull-reftime,r)@Pfit)*Vinter
     plt.plot(tfull,Vpath,label=f'Pathway #{n+1}',color=color)

From 553c298b62531847e8f9ae7d3a738ebe7b42dc83 Mon Sep 17 00:00:00 2001
From: Hugo Karas <hkaras@ethz.ch>
Date: Fri, 10 Jan 2025 16:48:08 +0100
Subject: [PATCH 3/9] typo fix

---
 examples/intermediate/ex_fitting_dqc_pathways.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/examples/intermediate/ex_fitting_dqc_pathways.py b/examples/intermediate/ex_fitting_dqc_pathways.py
index 56c680d24..3e349be0f 100644
--- a/examples/intermediate/ex_fitting_dqc_pathways.py
+++ b/examples/intermediate/ex_fitting_dqc_pathways.py
@@ -56,7 +56,7 @@
 
 # Plot the full detectable range
 tfull = np.arange(-2*tau1,2*tau2-4*tau3,0.008)
-Vmodelext =  dl.dipolarmodel(tfull,r,experiment=experiment)
+Vmodelext =  dl.dipolarmodel(tfull,r,experiment=experimentInfo)
 Vmodelext = dl.link(Vmodelext,lam23=['lam2','lam3'])
 
 # Extract results

From aad60009072959d0c5219d6ea7175c184b175768 Mon Sep 17 00:00:00 2001
From: Hugo Karas <hkaras@ethz.ch>
Date: Fri, 10 Jan 2025 16:57:22 +0100
Subject: [PATCH 4/9] Missed Examples

---
 examples/basic/ex_fitting_ridme.py                  |  2 +-
 examples/intermediate/ex_crossing_echoes_masking.py | 13 ++++---------
 examples/intermediate/ex_fitting_sparse_4pdeer.py   | 10 ++++++----
 3 files changed, 11 insertions(+), 14 deletions(-)

diff --git a/examples/basic/ex_fitting_ridme.py b/examples/basic/ex_fitting_ridme.py
index 56cc75a7f..357f1d402 100644
--- a/examples/basic/ex_fitting_ridme.py
+++ b/examples/basic/ex_fitting_ridme.py
@@ -56,7 +56,7 @@
 Pci50 = results.PUncert.ci(50)
 
 # Extract the unmodulated contribution
-Bfcn = dl.dipolarbackgroundmodel(experimentInfo)
+Bfcn = dl.dipolarbackgroundmodel(experimentInfo,dl.bg_strexp)
 Bfit = results.P_scale*results.evaluate(Bfcn,t)
 Bci = results.P_scale*results.propagate(Bfcn,t).ci(95)
 
diff --git a/examples/intermediate/ex_crossing_echoes_masking.py b/examples/intermediate/ex_crossing_echoes_masking.py
index dd60009e2..2defcb860 100644
--- a/examples/intermediate/ex_crossing_echoes_masking.py
+++ b/examples/intermediate/ex_crossing_echoes_masking.py
@@ -110,8 +110,8 @@
 r = np.arange(2,6,0.05) # nm
 
 # Construct the dipolar signal model 
-experiment = dl.ex_4pdeer(tau1,tau2,pathways=[1,2,3])
-Vmodel = dl.dipolarmodel(t,r,experiment=experiment)
+experimentInfo = dl.ex_4pdeer(tau1,tau2,pathways=[1,2,3])
+Vmodel = dl.dipolarmodel(t,r,experiment=experimentInfo)
 
 # Analyze the data while ignoring the crossing echoes
 results = dl.fit(Vmodel,Vexp, mask=mask, noiselvl=noiselevel)
@@ -130,13 +130,8 @@
 Pci50 = results.PUncert.ci(50)
 
 # Extract the unmodulated contribution
-def Vunmodfcn(lam1,lam2,lam3,reftime1,reftime2,reftime3,conc):
-    Lam0 = results.P_scale*(1-lam1-lam2-lam3)
-    Vunmod = Lam0*dl.bg_hom3d(t-reftime1,conc,lam1)
-    Vunmod *= dl.bg_hom3d(t-reftime2,conc,lam2)
-    Vunmod *= dl.bg_hom3d(t-reftime3,conc,lam3)
-    return Vunmod  
-Bfit = results.evaluate(Vunmodfcn)
+Bfcn = dl.dipolarbackgroundmodel(experimentInfo)
+Bfit = results.P_scale*results.evaluate(Bfcn,t)
 
 plt.figure(figsize=[6,7])
 violet = '#4550e6'
diff --git a/examples/intermediate/ex_fitting_sparse_4pdeer.py b/examples/intermediate/ex_fitting_sparse_4pdeer.py
index 68adaaa05..687088b38 100644
--- a/examples/intermediate/ex_fitting_sparse_4pdeer.py
+++ b/examples/intermediate/ex_fitting_sparse_4pdeer.py
@@ -41,7 +41,8 @@
 r = np.arange(2,10,0.05) # nm
 
 # Construct the model with the sparse sampled time vector
-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
@@ -65,9 +66,10 @@
 Pci50 = results.PUncert.ci(50)
 
 # Extract the unmodulated contribution
-Bfcn = lambda mod,conc,reftime: results.P_scale*(1-mod)*dl.bg_hom3d(tuniform-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'

From 6bd18d36fb17e4c76102c79a44b99d922dd48444 Mon Sep 17 00:00:00 2001
From: Hugo Karas <hkaras@ethz.ch>
Date: Wed, 15 Apr 2026 17:02:45 +0200
Subject: [PATCH 5/9] Add multi spin support

---
 deerlab/dipolarbackground.py | 216 +++++++++++++++++++++++++++--------
 1 file changed, 170 insertions(+), 46 deletions(-)

diff --git a/deerlab/dipolarbackground.py b/deerlab/dipolarbackground.py
index cb7b583c3..55123ac67 100644
--- a/deerlab/dipolarbackground.py
+++ b/deerlab/dipolarbackground.py
@@ -207,27 +207,39 @@ def basis(t,lam):
 
     return B
 
-def dipolarbackgroundmodel(experiment,basis=None,parametrization='reftimes'):
+def dipolarbackgroundmodel(experiment=None, basis=None, parametrization='reftimes', spins=2, samespins=True):
     """
-    Construct a dipolar background model for a given dipolar experiment. Currently limited to two-spin systems.
+    Construct a dipolar background model for a given dipolar experiment.
     This model can be used for evaluating and extrapolating the fitted parameters of the dipolar background.
+    Supports both two-spin and multi-spin systems.
 
     Parameters
     ----------
-    experiment : :ref:`ExperimentInfo`
-        The Experimental information obtained from experiment models (ex_).
-    
+    experiment : :ref:`ExperimentInfo` or None, optional
+        The Experimental information obtained from experiment models (ex_). If the experiment is not specified, a single-pathway background model is constructed with a refocusing time of 0.
+
     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'``. 
+        The default is ``'reftimes'``.
+
+    spins : integer, optional
+        Number of spins in the system. If not specified, defaults to two-spin systems. For multi-spin
+        systems (``spins>2``), three-spin interaction contributions are accounted for in the background.
+        The default is ``2``.
+
+    samespins : boolean, optional
+        Only relevant when ``spins>2``. Enables the assumption of spectral permutability, i.e., the
+        assumption that all pairwise interactions have the same amplitude. If enabled, one amplitude
+        parameter ``lam{n}`` is used per pathway. If disabled, individual amplitudes ``lam{n}_{q}``
+        are used for each pathway-interaction combination. Enabled by default.
 
     Returns
     -------
@@ -237,17 +249,36 @@ def dipolarbackgroundmodel(experiment,basis=None,parametrization='reftimes'):
     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)
+
+    To construct a dipolar background model for a three-spin system::
+
+            Bfcn = dl.dipolarbackgroundmodel(experimentInfo, spins=3)
+            Bfit = results.P_scale*results.evaluate(Bfcn, t)
     """
     from deerlab.dipolarmodel import _populate_dipolar_pathways_parameters
     from deerlab.model import Model
+    from itertools import permutations
+    from math import factorial
+
     if basis is None:
         from deerlab.bg_models import bg_hom3d
         basis = bg_hom3d
-    npathways = experiment.npathways
+
+    if experiment is None:
+        npathways = 1
+        harmonics = [1]
+        pulsedelay_names = ['reftime']
+    else:
+        npathways = experiment.npathways
+        harmonics = np.array(experiment.harmonics)[:npathways]
+        pulsedelay_names = list(inspect.signature(experiment.reftimes).parameters.keys())
+
+    Q = spins*(spins-1)//2
+    
 
     signature = []
     basis_signature = basis.signature.copy()
@@ -256,61 +287,154 @@ def dipolarbackgroundmodel(experiment,basis=None,parametrization='reftimes'):
         basis_signature.remove('lam')
     signature.extend(basis_signature)
 
-    if parametrization=='reftimes':
-        if npathways==1:
-            signature.extend(['mod','reftime'])
-        else:
+    if spins == 2:
+        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(pulsedelay_names)
             signature.extend(['lam'+str(i) for i in experiment.labels])
+        elif parametrization=='shift':
+            signature.extend(['tshift'])
+            signature.extend(['lam'+str(i+1) for i in experiment.labels])
+    else:
+        if parametrization=='reftimes':
+            if samespins:
+                signature.extend(['lam'+str(i) for i in experiment.labels])
+            else:
+                signature.extend(['lam'+str(i)+'_'+str(q+1) for i in experiment.labels for q in range(Q)])
             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])
-    
+            signature.append('lamu')
+        elif parametrization=='delays':
+            signature.extend(pulsedelay_names)
+            if samespins:
+                signature.extend(['lam'+str(i) for i in experiment.labels])
+            else:
+                signature.extend(['lam'+str(i)+'_'+str(q+1) for i in experiment.labels for q in range(Q)])
+            signature.append('lamu')
+        elif parametrization=='shift':
+            signature.extend(['tshift'])
+            if samespins:
+                signature.extend(['lam'+str(i) for i in experiment.labels])
+            else:
+                signature.extend(['lam'+str(i)+'_'+str(q+1) for i in experiment.labels for q in range(Q)])
+            signature.append('lamu')
 
     # 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:
+        if spins == 2:
+            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]
+            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
+
+        else:
+            # Multi-spin: parse parameters
+            if parametrization=='reftimes':
+                n_lam = npathways if samespins else npathways*Q
+                λs_flat = pathway_params[:n_lam]
+                trefs = pathway_params[n_lam:n_lam+npathways]
+                lamu = pathway_params[-1]
+            elif parametrization=='delays':
+                delays = pathway_params[0:len(experiment.delays)]
+                remaining = pathway_params[len(experiment.delays):]
+                trefs = experiment.reftimes(*delays)
+                λs_flat = remaining[:-1]
+                lamu = remaining[-1]
+            elif parametrization=='shift':
+                shift = pathway_params[0]
+                remaining = pathway_params[1:]
+                delays = np.array(experiment.delays) + shift
+                trefs = experiment.reftimes(*delays)
+                λs_flat = remaining[:-1]
+                lamu = remaining[-1]
+
+            # Build λs as [Q x npathways] array: λs_matrix[q, idx]
+            if samespins:
+                λs_matrix = np.tile(λs_flat, (Q, 1))
+            else:
+                λs_matrix = np.array(λs_flat).reshape(npathways, Q).T
+
+            # Wrap basis into a callable for dipolarbackground
             if lam_in_basis:
-                basis_constants[-1] = λs[i]
-            prod *= basis(t-trefs[i],*basis_constants[1:])
-        return scale*prod
+                def basis_fcn(tdip, lam):
+                    bc = basis_constants.copy()
+                    bc.append(lam)
+                    return basis(tdip, *bc[1:])
+            else:
+                def basis_fcn(tdip):
+                    return basis(tdip, *basis_constants[1:])
+
+            # Construct all multi-spin dipolar pathways
+            paths, threespin_paths = [], []
+            Λ0 = 1.0
+
+            for idx in range(npathways):
+                tref_idx = trefs[idx]
+                δ_idx = harmonics[idx]
+
+                for perm in set(permutations([0]+[None]*(Q-1))):
+                    q = int(np.where(np.array(perm)==0)[0][0])
+                    λp = λs_matrix[q, idx]
+                    λ2k = factorial(Q-1)*λp*lamu**(Q-1)
+                    tref_path = tuple([tref_idx if n==0 else None for n in perm])
+                    δs_path = tuple([δ_idx if n==0 else 0 for n in perm])
+                    paths.append({'amp': λ2k, 'reftime': tref_path, 'harmonic': δs_path})
+                    Λ0 -= λ2k
+
+                    for idx2 in range(idx+1, npathways):
+                        tref_idx2 = trefs[idx2]
+                        δ_idx2 = harmonics[idx2]
+                        for perm2 in set(permutations([0,1]+[None]*(Q-2))):
+                            perm2_arr = np.array(perm2)
+                            q2 = int(np.where(perm2_arr==1)[0][0])
+                            λm = λs_matrix[q2, idx]
+                            λ3k = 2*factorial(Q-2)*λp*λm*lamu**(Q-2)
+                            tref_path2 = tuple([tref_idx if n==0 else tref_idx2 if n==1 else None for n in perm2])
+                            δs_path2 = tuple([δ_idx if n==0 else δ_idx2 if n==1 else 0 for n in perm2])
+                            threespin_paths.append({'amp': λ3k, 'reftime': tref_path2, 'harmonic': δs_path2})
+                            Λ0 -= λ3k
+
+            paths.append({'amp': Λ0})
+            return dipolarbackground(t, paths+threespin_paths, basis_fcn)
 
     # 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)
-    
+    bckgrnd_model = _populate_dipolar_pathways_parameters(bckgrnd_model, npathways, spins=spins, samespins=samespins,
+                                                          experiment=experiment, parametrization=parametrization,
+                                                          Q=Q if spins>2 else None)
+
     # 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

From 0074e3af0312c23e09c39406ac36abc5252e14c8 Mon Sep 17 00:00:00 2001
From: Hugo Karas <hkaras@ethz.ch>
Date: Wed, 15 Apr 2026 17:03:40 +0200
Subject: [PATCH 6/9] `dl.dipolarmodel` now exports with the `Bmodel`,
 `Pmodel`, `r` and `t`

---
 deerlab/dipolarmodel.py | 20 ++++++++++++++++++--
 1 file changed, 18 insertions(+), 2 deletions(-)

diff --git a/deerlab/dipolarmodel.py b/deerlab/dipolarmodel.py
index 179de16ed..5d2b91042 100644
--- a/deerlab/dipolarmodel.py
+++ b/deerlab/dipolarmodel.py
@@ -5,7 +5,7 @@
 
 import numpy as np
 from deerlab.dipolarkernel import dipolarkernel
-from deerlab.dipolarbackground import dipolarbackground
+from deerlab.dipolarbackground import dipolarbackground, dipolarbackgroundmodel
 from deerlab.regoperator import regoperator
 from deerlab.dd_models import freedist
 from deerlab.model import Model,Penalty, link
@@ -212,7 +212,13 @@ def dipolarmodel(t, r=None, Pmodel=None, Bmodel=bg_hom3d, experiment=None, param
     Returns
     -------
     Vmodel : :ref:`Model`
-        Dipolar signal model object.
+        Dipolar signal model object. With these additional attributes:
+            * ``t`` - The time axis of the model.
+            * ``r`` - The distance axis of the model, only for two-spin systems.
+            * ``Bmodel`` - The multi-pathway background model object used in the construction of the dipolar signal model.
+            * ``Pmodel`` - The distance distribution model object used in the construction of the dipolar signal model.
+            
+
     """
 
     # If distance axis not specified default to a long range one
@@ -629,6 +635,16 @@ def Vmultispin_nonlinear_fcn(*nonlin):
                     max(t)+max(reftimes) + 3*experiment.pulselength + dt,
                     dt)
 
+    # 
+    if Bmodel is not None:
+        DipolarSignal.Bmodel = dipolarbackgroundmodel(experiment=experiment, basis=Bmodel, parametrization=parametrization, spins=spins, samespins=samespins)
+    
+    if Pmodel is not None:
+        DipolarSignal.Pmodel = Pmodel
+    DipolarSignal.t = t
+    if spins==2:
+        DipolarSignal.r = r
+
     return DipolarSignal
 #===============================================================================
 

From a4cc84cb7cd8c779d5e4ae0d2ec230d1a2cb3ffa Mon Sep 17 00:00:00 2001
From: Hugo Karas <hkaras@ethz.ch>
Date: Wed, 15 Apr 2026 17:04:09 +0200
Subject: [PATCH 7/9] Fit returns `bg` and `bgUncert`

---
 deerlab/fit.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/deerlab/fit.py b/deerlab/fit.py
index 71d262310..c1c27da30 100644
--- a/deerlab/fit.py
+++ b/deerlab/fit.py
@@ -548,9 +548,15 @@ def _scale(x):
     if len(noiselvl)==1: 
         noiselvl = noiselvl[0]
     
+
     # Generate FitResult object from all the dictionaries
     fitresult = FitResult({**FitResult_param_,**FitResult_paramuq_, **FitResult_dict,'penweights':penweights,'noiselvl':noiselvl,'paramlist':_paramlist, '_param_idx':param_idx}) 
 
+    # Add background if it is part of the model
+    if hasattr(model,'Bmodel'):
+        fitresult['bg'] = fitresult.evaluate(model.Bmodel,model.t)
+        fitresult['bgUncert'] = fitresult.propagate(model.Bmodel,model.t, samples=bootstrap)
+
     fitresult._summary = _print_fitresults(fitresult,model)
 
     return fitresult

From dca05749b9e1ba37f32f2be7bb3d6d745f629153 Mon Sep 17 00:00:00 2001
From: Hugo Karas <hkaras@ethz.ch>
Date: Wed, 15 Apr 2026 17:04:23 +0200
Subject: [PATCH 8/9] Added tests and an example

---
 examples/intermediate/ex_background_model.py |  75 ++++++++
 test/test_dipolarbackground.py               | 177 ++++++++++++++++++-
 test/test_dipolarmodel.py                    |  22 ++-
 3 files changed, 269 insertions(+), 5 deletions(-)
 create mode 100644 examples/intermediate/ex_background_model.py

diff --git a/examples/intermediate/ex_background_model.py b/examples/intermediate/ex_background_model.py
new file mode 100644
index 000000000..9d47cb60e
--- /dev/null
+++ b/examples/intermediate/ex_background_model.py
@@ -0,0 +1,75 @@
+# %% [markdown]
+""" 
+Generating a background model for a 4-pulse DEER experiment
+-------------------------------------------------------------------------
+
+Since DeerLab v1.2, the background model can be generated for a given experiment using the :func:`deerlab.dipolarbackgroundmodel` function, as is normally automatically calculated in the fit function as `FitResult.bg` and  `FitResult.bgUncert`.
+
+However, it can still be useful to generate the background model separately, and so this example shows how to do that for a 4-pulse DEER experiment.
+
+First we will fit a simple 4-pulse DEER dataset with multiple-pathways. 
+""" 
+
+
+import numpy as np
+import matplotlib.pyplot as plt
+import deerlab as dl
+# %%
+# File location
+path = '../data/'
+file = 'example_4pdeer_1.DTA'
+
+# Experimental parameters
+tau1 = 0.3      # First inter-pulse delay, μs
+tau2 = 4.0      # Second inter-pulse delay, μs
+tmin = 0.1      # Start time, μs
+
+# Load the experimental data
+t,Vexp = dl.deerload(path + file)
+
+# Pre-processing
+Vexp = dl.correctphase(Vexp) # Phase correction
+Vexp = Vexp/np.max(Vexp)     # Rescaling (aesthetic)
+t = t - t[0]                  # Account for zerotime
+t = t + tmin    
+# Distance vector
+r = np.arange(2.5,5,0.01) # nm
+
+# Construct the model
+experimentInfo = dl.ex_4pdeer(tau1,tau2, pathways=[1,2,3])
+Vmodel = dl.dipolarmodel(t,r, experiment = experimentInfo, Bmodel=dl.bg_hom3d)
+
+# Fit the model to the data
+results = dl.fit(Vmodel,Vexp)
+
+# Print results summary
+print(results)
+
+
+# %%
+"""
+Now using DeerLab > v1.2, the background model can be directly plotted from the fit results.
+"""
+plt.figure(figsize=[6,7])
+violet = '#4550e6'
+plt.plot(t,Vexp,'.',color='grey',label='Data')
+plt.plot(t,results.model,linewidth=3,color=violet,label='Fit')
+plt.plot(t,results.bg,'--',linewidth=3,color=violet,label='Unmodulated contribution')
+
+
+# %%
+"""
+Alternatively we could generate the background model separately using the `dl.dipolarbackgroundmodel` function, which takes the experimental information as input.
+"""
+
+bg_model = dl.dipolarbackgroundmodel(experimentInfo, basis = dl.bg_hom3d)
+bg = results.evaluate(bg_model,t)
+bgUncert = results.propagate(bg_model,t)
+
+"""
+In DeerLab <1.2, the multi-pathway background model needs to be generated manually from a function.
+"""
+
+bf_func = lambda lam1,lam2,lam3,reftime1,reftime2,reftime3,conc: results.P_scale*(1-lam1-lam2-lam3)*dl.bg_hom3d(t-reftime1,conc,lam1)*dl.bg_hom3d(t-reftime2,conc,lam2)*dl.bg_hom3d(t-reftime3,conc,lam3)
+bf = results.evaluate(bf_func, t)
+bfUncert = results.propagate(bf_func, t)
diff --git a/test/test_dipolarbackground.py b/test/test_dipolarbackground.py
index f25f7c6e7..73a2420d4 100644
--- a/test/test_dipolarbackground.py
+++ b/test/test_dipolarbackground.py
@@ -1,9 +1,12 @@
 
 import numpy as np
-from deerlab.bg_models import _hom3d, bg_exp, _homfractal
-from deerlab.dipolarbackground import dipolarbackground
+from itertools import permutations
+from math import factorial
+from deerlab.bg_models import _hom3d, bg_exp, _homfractal, bg_hom3d
+from deerlab.dipolarbackground import dipolarbackground, dipolarbackgroundmodel
+from deerlab.dipolarmodel import ex_3pdeer, ex_4pdeer
 from deerlab.utils import assert_docstring
-import pytest 
+import pytest
 
 # Fixtures 
 # -----------------------------------------------------------------------------------
@@ -132,3 +135,171 @@ def test_docstring():
     "Check that the docstring includes all variables and keywords."
     assert_docstring(dipolarbackground)
 # ======================================================================
+
+
+# ==================================================================================
+# Tests for dipolarbackgroundmodel
+# ==================================================================================
+
+tau = 0.5
+tau1, tau2 = 0.5, 3.0
+
+# ==================================================================================
+def test_dipolarbackgroundmodel_twospin_single_pathway():
+    "Test dipolarbackgroundmodel for a two-spin system with a single pathway"
+    mod = 0.4
+    tref = tau
+    Bref = (1-mod) * _hom3d(t-tref, conc, mod)
+    experiment = ex_3pdeer(tau, pathways=[1])
+    Bmodel = dipolarbackgroundmodel(experiment, bg_hom3d)
+    Bval = Bmodel(t, conc, mod, tref)
+    assert np.allclose(Bval, Bref)
+# ==================================================================================
+
+# ==================================================================================
+def test_dipolarbackgroundmodel_twospin_single_pathway_no_experiment():
+    "Test dipolarbackgroundmodel for a two-spin system with a single pathway with no experiment"
+    mod = 0.4
+    tref = tau
+    Bref = (1-mod) * _hom3d(t-tref, conc, mod)
+    Bmodel = dipolarbackgroundmodel(None, bg_hom3d)
+    Bval = Bmodel(t, conc, mod, tref)
+    assert np.allclose(Bval, Bref)
+# ==================================================================================
+
+
+# ==================================================================================
+def test_dipolarbackgroundmodel_twospin_multi_pathway():
+    "Test dipolarbackgroundmodel for a two-spin system with multiple pathways"
+    lam1, lam2 = 0.3, 0.1
+    tref1, tref2 = tau1, tau1+tau2
+    Bref = (1-lam1-lam2) * _hom3d(t-tref1, conc, lam1) * _hom3d(t-tref2, conc, lam2)
+    experiment = ex_4pdeer(tau1, tau2, pathways=[1,2])
+    Bmodel = dipolarbackgroundmodel(experiment, bg_hom3d)
+    Bval = Bmodel(t, conc, lam1, lam2, tref1, tref2)
+    assert np.allclose(Bval, Bref)
+# ==================================================================================
+
+# ==================================================================================
+def test_dipolarbackgroundmodel_twospin_delays():
+    "Test dipolarbackgroundmodel for a two-spin system with delays parametrization"
+    lam1, lam2 = 0.3, 0.1
+    tref1, tref2 = tau1, tau1+tau2
+    Bref = (1-lam1-lam2) * _hom3d(t-tref1, conc, lam1) * _hom3d(t-tref2, conc, lam2)
+    experiment = ex_4pdeer(tau1, tau2, pathways=[1,2])
+    Bmodel = dipolarbackgroundmodel(experiment, bg_hom3d, parametrization='delays')
+    Bval = Bmodel(t, conc, tau1, tau2, lam1, lam2)
+    assert np.allclose(Bval, Bref)
+# ==================================================================================
+
+# ==================================================================================
+def test_dipolarbackgroundmodel_threespin_single_pathway():
+    "Test dipolarbackgroundmodel for a three-spin system with a single pathway"
+    lam, lamu = 0.4, 0.8
+    Q = 3
+    # Construct reference pathways manually
+    Λ0 = 1.0
+    paths_ref = []
+    for perm in set(permutations([0]+[None]*(Q-1))):
+        λ2k = factorial(Q-1)*lam*lamu**(Q-1)
+        tref_path = tuple([tau if n==0 else None for n in perm])
+        δs_path = tuple([1 if n==0 else 0 for n in perm])
+        paths_ref.append({'amp': λ2k, 'reftime': tref_path, 'harmonic': δs_path})
+        Λ0 -= λ2k
+    paths_ref.append({'amp': Λ0})
+    Bref = dipolarbackground(t, paths_ref, lambda td, l: _hom3d(td, conc, l))
+    experiment = ex_3pdeer(tau, pathways=[1])
+    Bmodel = dipolarbackgroundmodel(experiment, bg_hom3d, spins=3)
+    Bval = Bmodel(t, conc, lam, tau, lamu)
+    assert np.allclose(Bval, Bref)
+# ==================================================================================
+
+# ==================================================================================
+def test_dipolarbackgroundmodel_threespin_multi_pathway():
+    "Test dipolarbackgroundmodel for a three-spin system with multiple pathways"
+    lam1, lam2, lamu = 0.3, 0.1, 0.7
+    Q = 3
+    harmonics_exp = [1, 1]
+    trefs_exp = [tau1, tau1+tau2]
+
+    # Construct reference pathways manually
+    Λ0 = 1.0
+    paths_ref, threespin_ref = [], []
+    for idx, (tref_idx, δ_idx) in enumerate(zip(trefs_exp, harmonics_exp)):
+        lam_idx = [lam1, lam2][idx]
+        for perm in set(permutations([0]+[None]*(Q-1))):
+            q = int(np.where(np.array(perm)==0)[0][0])
+            λp = lam_idx
+            λ2k = factorial(Q-1)*λp*lamu**(Q-1)
+            paths_ref.append({'amp': λ2k,
+                               'reftime': tuple([tref_idx if n==0 else None for n in perm]),
+                               'harmonic': tuple([δ_idx if n==0 else 0 for n in perm])})
+            Λ0 -= λ2k
+            for idx2 in range(idx+1, 2):
+                tref_idx2 = trefs_exp[idx2]
+                δ_idx2 = harmonics_exp[idx2]
+                for perm2 in set(permutations([0,1]+[None]*(Q-2))):
+                    perm2_arr = np.array(perm2)
+                    q2 = int(np.where(perm2_arr==1)[0][0])
+                    λm = lam_idx
+                    λ3k = 2*factorial(Q-2)*λp*λm*lamu**(Q-2)
+                    threespin_ref.append({'amp': λ3k,
+                                          'reftime': tuple([tref_idx if n==0 else tref_idx2 if n==1 else None for n in perm2]),
+                                          'harmonic': tuple([δ_idx if n==0 else δ_idx2 if n==1 else 0 for n in perm2])})
+                    Λ0 -= λ3k
+    paths_ref.append({'amp': Λ0})
+    Bref = dipolarbackground(t, paths_ref+threespin_ref, lambda td, l: _hom3d(td, conc, l))
+
+    experiment = ex_4pdeer(tau1, tau2, pathways=[1,2])
+    Bmodel = dipolarbackgroundmodel(experiment, bg_hom3d, spins=3)
+    Bval = Bmodel(t, conc, lam1, lam2, tau1, tau1+tau2, lamu)
+    assert np.allclose(Bval, Bref)
+# ==================================================================================
+
+# ==================================================================================
+def test_dipolarbackgroundmodel_threespin_nosamespins():
+    "Test dipolarbackgroundmodel for a three-spin system with samespins=False"
+    lam_q = [0.3, 0.4, 0.2]   # different amplitudes per interaction
+    lamu = 0.8
+    Q = 3
+    # Construct reference pathways with per-interaction amplitudes
+    Λ0 = 1.0
+    paths_ref = []
+    for perm in set(permutations([0]+[None]*(Q-1))):
+        q = int(np.where(np.array(perm)==0)[0][0])
+        λ2k = factorial(Q-1)*lam_q[q]*lamu**(Q-1)
+        tref_path = tuple([tau if n==0 else None for n in perm])
+        δs_path = tuple([1 if n==0 else 0 for n in perm])
+        paths_ref.append({'amp': λ2k, 'reftime': tref_path, 'harmonic': δs_path})
+        Λ0 -= λ2k
+    paths_ref.append({'amp': Λ0})
+    Bref = dipolarbackground(t, paths_ref, lambda td, l: _hom3d(td, conc, l))
+
+    experiment = ex_3pdeer(tau, pathways=[1])
+    Bmodel = dipolarbackgroundmodel(experiment, bg_hom3d, spins=3, samespins=False)
+    Bval = Bmodel(t, conc, *lam_q, tau, lamu)
+    assert np.allclose(Bval, Bref)
+# ==================================================================================
+
+# ==================================================================================
+def test_dipolarbackgroundmodel_threespin_delays():
+    "Test dipolarbackgroundmodel for a three-spin system with delays parametrization"
+    lam, lamu = 0.4, 0.8
+    Q = 3
+    Λ0 = 1.0
+    paths_ref = []
+    for perm in set(permutations([0]+[None]*(Q-1))):
+        λ2k = factorial(Q-1)*lam*lamu**(Q-1)
+        tref_path = tuple([tau if n==0 else None for n in perm])
+        δs_path = tuple([1 if n==0 else 0 for n in perm])
+        paths_ref.append({'amp': λ2k, 'reftime': tref_path, 'harmonic': δs_path})
+        Λ0 -= λ2k
+    paths_ref.append({'amp': Λ0})
+    Bref = dipolarbackground(t, paths_ref, lambda td, l: _hom3d(td, conc, l))
+
+    experiment = ex_3pdeer(tau, pathways=[1])
+    Bmodel = dipolarbackgroundmodel(experiment, bg_hom3d, spins=3, parametrization='delays')
+    # signature: t(const), conc, delay1, lam1, lamu
+    Bval = Bmodel(t, conc, tau, lam, lamu)
+    assert np.allclose(Bval, Bref)
+# ==================================================================================
diff --git a/test/test_dipolarmodel.py b/test/test_dipolarmodel.py
index f7b987e03..f0ebd8f39 100644
--- a/test/test_dipolarmodel.py
+++ b/test/test_dipolarmodel.py
@@ -4,7 +4,7 @@
 from deerlab.model import Model
 from deerlab.fit import fit
 from deerlab.dipolarmodel import ExperimentInfo,dipolarpenalty, dipolarmodel, ex_4pdeer, ex_3pdeer,ex_fwd5pdeer, ex_rev5pdeer, ex_sifter, ex_ridme, ex_dqc
-from deerlab import dd_gauss,dd_gauss2,bg_hom3d,bg_exp
+from deerlab import dd_gauss,dd_gauss2,bg_hom3d,bg_exp,UQResult
 from deerlab.utils import assert_docstring
 from pytest import fixture 
 
@@ -153,6 +153,7 @@ def test_phenomenological_Bmodel(V1path_phenoB):
     Vsim = Vmodel(mod=0.3,reftime=0.0,mean=3,std=0.2,decay=0.1,scale=1e5)
 
     assert np.allclose(Vsim,V1path_phenoB)
+   
 # ======================================================================
 
 # ======================================================================
@@ -164,6 +165,7 @@ def test_no_Bmodel(V1path_noB):
     Vsim = Vmodel(mod=0.3,reftime=0.0,mean=3,std=0.2,scale=1e5)
 
     assert np.allclose(Vsim,V1path_noB)
+    assert not hasattr(Vmodel,'Bmodel')
 # ======================================================================
 
 # ======================================================================
@@ -177,6 +179,19 @@ def test_model_1pathways(V1path):
     assert np.allclose(Vsim,V1path)
 # ======================================================================
 
+# ======================================================================
+def test_model_1pathways_background(V1path): 
+    "Check that the model with one dipolar pathway is correct"
+
+    Vmodel = dipolarmodel(t,r,dd_gauss,bg_hom3d,npathways=1)
+    
+    Vsim = Vmodel(mod=0.3,reftime=0.0,conc=50,mean=3,std=0.2,scale=1e5)
+
+    assert np.allclose(Vsim,V1path)
+    assert hasattr(Vmodel,'Bmodel') and Vmodel.Bmodel is not None
+    assert isinstance(Vmodel.Bmodel,Model)
+# ======================================================================
+
 # ======================================================================
 def test_model_2pathways(V2path): 
     "Check that the model with two dipolar pathways is correct"
@@ -212,6 +227,8 @@ def test_fit_1pathways(V1path):
     result = fit(Vmodel,V1path,ftol=1e-4)
 
     assert np.allclose(result.model,V1path)
+    assert hasattr(result,'bg') and result.bg is not None and isinstance(result.bg,np.ndarray)
+    assert hasattr(result,'bgUncert') and result.bgUncert is not None and isinstance(result.bgUncert,UQResult)
 # ======================================================================
 
 # ======================================================================
@@ -449,7 +466,8 @@ def test_ex_fwd5pdeer_fit(Vfwd5pdeer,Pr):
     result = fit(Vmodel,Vfwd5pdeer,ftol=1e-4)
 
     assert np.allclose(Vfwd5pdeer,result.model,rtol=1e-3) and ovl(result.P/1e5,Pr)>0.975
-# ======================================================================
+# ========
+# ==============================================================
 
 # ======================================================================
 def test_ex_fwd5pdeer_pulselength(): 

From 89eb9902dec8449fe4ae5c9f476939cef47648fa Mon Sep 17 00:00:00 2001
From: Hugo Karas <hkaras@ethz.ch>
Date: Wed, 15 Apr 2026 17:44:54 +0200
Subject: [PATCH 9/9] Bug fixes

---
 deerlab/dipolarbackground.py | 36 ++++++++++++++++++++++--------------
 deerlab/dipolarmodel.py      |  2 +-
 2 files changed, 23 insertions(+), 15 deletions(-)

diff --git a/deerlab/dipolarbackground.py b/deerlab/dipolarbackground.py
index 55123ac67..070ca8014 100644
--- a/deerlab/dipolarbackground.py
+++ b/deerlab/dipolarbackground.py
@@ -207,7 +207,7 @@ def basis(t,lam):
 
     return B
 
-def dipolarbackgroundmodel(experiment=None, basis=None, parametrization='reftimes', spins=2, samespins=True):
+def dipolarbackgroundmodel(experiment=None, basis=None, parametrization='reftimes', spins=2, samespins=True, **kwargs):
     """
     Construct a dipolar background model for a given dipolar experiment.
     This model can be used for evaluating and extrapolating the fitted parameters of the dipolar background.
@@ -241,6 +241,9 @@ def dipolarbackgroundmodel(experiment=None, basis=None, parametrization='reftime
         parameter ``lam{n}`` is used per pathway. If disabled, individual amplitudes ``lam{n}_{q}``
         are used for each pathway-interaction combination. Enabled by default.
 
+    **kwargs:
+        Additional keyword arguments. If ``experiment`` is not specified, the number of pathways can be specified with the keyword argument ``npathways`` (default is 1).
+
     Returns
     -------
 
@@ -269,13 +272,18 @@ def dipolarbackgroundmodel(experiment=None, basis=None, parametrization='reftime
         basis = bg_hom3d
 
     if experiment is None:
-        npathways = 1
+        npathways = kwargs.get('npathways', 1)
         harmonics = [1]
-        pulsedelay_names = ['reftime']
+        if npathways>1:
+            pulsedelay_names = ['reftime'+str(i+1) for i in range(npathways)]
+        else:
+            pulsedelay_names = ['reftime']
+        exp_labels = list(np.arange(1,npathways+1,dtype=int))
     else:
         npathways = experiment.npathways
         harmonics = np.array(experiment.harmonics)[:npathways]
         pulsedelay_names = list(inspect.signature(experiment.reftimes).parameters.keys())
+        exp_labels = experiment.labels
 
     Q = spins*(spins-1)//2
     
@@ -292,35 +300,35 @@ def dipolarbackgroundmodel(experiment=None, basis=None, parametrization='reftime
             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])
+                signature.extend(['lam'+str(i) for i in exp_labels])
+                signature.extend(['reftime'+str(i) for i in exp_labels])
         elif parametrization=='delays':
             signature.extend(pulsedelay_names)
-            signature.extend(['lam'+str(i) for i in experiment.labels])
+            signature.extend(['lam'+str(i) for i in exp_labels])
         elif parametrization=='shift':
             signature.extend(['tshift'])
-            signature.extend(['lam'+str(i+1) for i in experiment.labels])
+            signature.extend(['lam'+str(i) for i in exp_labels])
     else:
         if parametrization=='reftimes':
             if samespins:
-                signature.extend(['lam'+str(i) for i in experiment.labels])
+                signature.extend(['lam'+str(i) for i in exp_labels])
             else:
-                signature.extend(['lam'+str(i)+'_'+str(q+1) for i in experiment.labels for q in range(Q)])
-            signature.extend(['reftime'+str(i) for i in experiment.labels])
+                signature.extend(['lam'+str(i)+'_'+str(q+1) for i in exp_labels for q in range(Q)])
+            signature.extend(['reftime'+str(i) for i in exp_labels])
             signature.append('lamu')
         elif parametrization=='delays':
             signature.extend(pulsedelay_names)
             if samespins:
-                signature.extend(['lam'+str(i) for i in experiment.labels])
+                signature.extend(['lam'+str(i) for i in exp_labels])
             else:
-                signature.extend(['lam'+str(i)+'_'+str(q+1) for i in experiment.labels for q in range(Q)])
+                signature.extend(['lam'+str(i)+'_'+str(q+1) for i in exp_labels for q in range(Q)])
             signature.append('lamu')
         elif parametrization=='shift':
             signature.extend(['tshift'])
             if samespins:
-                signature.extend(['lam'+str(i) for i in experiment.labels])
+                signature.extend(['lam'+str(i) for i in exp_labels])
             else:
-                signature.extend(['lam'+str(i)+'_'+str(q+1) for i in experiment.labels for q in range(Q)])
+                signature.extend(['lam'+str(i)+'_'+str(q+1) for i in exp_labels for q in range(Q)])
             signature.append('lamu')
 
     # Construct the dipolar background model
diff --git a/deerlab/dipolarmodel.py b/deerlab/dipolarmodel.py
index 5d2b91042..1551d1b59 100644
--- a/deerlab/dipolarmodel.py
+++ b/deerlab/dipolarmodel.py
@@ -637,7 +637,7 @@ def Vmultispin_nonlinear_fcn(*nonlin):
 
     # 
     if Bmodel is not None:
-        DipolarSignal.Bmodel = dipolarbackgroundmodel(experiment=experiment, basis=Bmodel, parametrization=parametrization, spins=spins, samespins=samespins)
+        DipolarSignal.Bmodel = dipolarbackgroundmodel(experiment=experiment, basis=Bmodel, parametrization=parametrization, spins=spins, samespins=samespins, npathways=npathways)
     
     if Pmodel is not None:
         DipolarSignal.Pmodel = Pmodel