# dipolarmodel.py - DeerLab's dipolar EPR model generator# ---------------------------------------------------------------------------# This file is a part of DeerLab. License is MIT (see LICENSE.md). # Copyright(c) 2019-2022: Luis Fabregas, Stefan Stoll and other contributors.importnumpyasnpfromdeerlab.dipolarkernelimportdipolarkernel,dipolarbackgroundfromdeerlab.regoperatorimportregoperatorfromdeerlab.dd_modelsimportfreedistfromdeerlab.modelimportModel,Penalty,linkfromdeerlabimportbg_hom3dfromdeerlab.utilsimportcholeskycovmatfromdeerlab.constantsimport*frommathimportfactorialfromitertoolsimportpermutationsimportinspectfromscipy.statsimportmultivariate_normal#===============================================================================
[docs]defdipolarmodel(t,r=None,Pmodel=None,Bmodel=bg_hom3d,experiment=None,parametrization='reftimes',npathways=1,spins=2,harmonics=None,excbandwidth=np.inf,orisel=None,g=[ge,ge],gridsize=1000,minamp=1e-3,samespins=True,triangles=None,interp=True,):""" Generate a dipolar EPR signal model. This function generates a model for the dipolar EPR signal arising from the dipolar interactions between the electron spins of the system. The model is defined in terms of the evolution time vector ``t``, the spin-spin distances ``r``, the distance distribution model ``Pmodel``, the intermolecular contribution ``Bmodel``, and an ``experiment`` object containing information about the specific experiment. The number of spins in the system ``spins`` can be specified, as well as the parametrization strategy ``parametrization`` for the dipolar pathway refocusing times. For multi-spin systems (``spins > 2``), the function allows for the inclusion of three-spin interactions by specifying the ``triangles`` list of interspin triangles and enabling the ``samespins`` assumption of spectral permutability. The ``minamp`` parameter can be used to threshold the amplitude of the three-spin interactions. Additional effects such as orientation selection or limited excitation bandwidth can be included via the ``orisel`` and ``excbandwidth`` arguments, respectively. Parameters ---------- t : array_like Vector of dipolar evolution times, in microseconds. r : array_like, optional Vector of spin-spin distances, in nanometers. If not specified, it defaults to the range 1.5nm-8nm with a resolution of 0.05nm. Has no effect on models with ``spins>2``. Pmodel : :ref:`Model`, optional Model for the distance distribution. If not specified, a non-parametric distance distribution defined over ``r`` is used. For models with ``spins>2`` a multivariate Gaussian distance distribution is assumed. Bmodel : :ref:`Model`, optional Model for the intermolecular (background) contribution. If not specified, a background arising from a homogenous 3D distribution of spins is used. experiment : :ref:`ExperimentInfo`, optional Experimental information obtained from experiment models (``ex_``). If specified, the boundaries and start values of the dipolar pathways' refocusing times and amplitudes are refined based on the specific experiment's delays. parametrization : string, optional Parametrization strategy of the dipolar pathway refocusing times. Must be one of the following: * ``'reftimes'`` - Each refocusing time is represented individually as a parameter. * ``'delays'`` - The pulse delays are introduced as parameters from which the refocusing times are computed. Requires ``experiment`` to be specified. * ``'shift'`` - A time shift is introduced as a parameter to represent the variability of the refocusing times from their theoretical values. Requires ``experiment`` to be specified. The default is ``'reftimes'``. npathways : integer scalar, optional Number of dipolar pathways. If not specified, a single dipolar pathway is used. harmonics : list of integers, optional Harmonic prefactors of the dipolar pathways. Must be a list with ``npathways`` harmonics, one for each pathway. spins : integer scalar, optional Number of spins in system. If not specified it defaults to two-spin systems. For multi-spins systems indicated by ``spins>2`` three-spin interactions are accounted for by the model automatically. triangles : list of lists, optional Only required when ``spins>3``. List of ``(Nspins)!/(3!*(Nspins-3)!)`` triangles formed by the different interspin distances. Each triangle is specified by a list of indices of the interspin distances forming the triangle, e.g. for a three-spin system ``triangles=[[1,2,3]]`` where the first, second, and third interspin distances connect to form a triangle. samespins : boolean, optional Only when ``spins>3``. Enables the assumption of spectral permutability, i.e. the assumption that all spins in the system have the same spectral distribution. If enabled, all pairwise pathways of the different dipolar interactions are assumed to have the same amplitude. If disabled, all the individual pathway amplitudes are parametrized separately. Enabled by default. minamp : float scalar, optional Only when ``spins>3``. Threshold amplitude for neglecting three-spin interaction pathways. To enhance performance, all dipolar pathways arising from three-spin interactions with an amplitude ``lam<minamp`` are approximated as having zero amplitude. The default threshold is ``1e-3``. orisel : callable or ``None``, optional Probability distribution of possible orientations of the interspin vector to account for orientation selection. Must be a function taking a value of the angle θ ∈ [0,π/2] between the interspin vector and the external magnetic field and returning the corresponding probability density. If specified as ``None`` (default), a uniform distribution is used. gridsize : scalar, optional Number of points on the grid of powder orientations to be used in the ``'grid'`` kernel calculation method when evaluating models with orientation selection or multi-spin effects. By default set to 1000 points. excbandwidth : scalar, optional Excitation bandwidth of the pulses in MHz to account for limited excitation bandwidth. If not specified, an infinite excitation bandwidth is used. g : scalar, 2-element array, optional Electron g-values of the two spins, ``[g1, g2]``. If a single g is specified, ``[g, g]`` is used. If not specified, g = 2.002319... is used for both spins. interp : boolean, optional Enable dipolar kernel interpolation for computation time reduction. By default enabled. Returns ------- Vmodel : :ref:`Model` Dipolar signal model object. """# If distance axis not specified default to a long range oneifrisNone:r=np.arange(1.5,8,0.05)# Input parsing and validationifnotisinstance(Pmodel,Model)andPmodelisnotNone:raiseTypeError('The argument Pmodel must be a valid Model object')ifnotisinstance(Bmodel,Model)andBmodelisnotNone:raiseTypeError('The argument Bmodel must be a valid Model object or None.')ifnotisinstance(npathways,int)ornpathways<=0:raiseValueError('The number of pathways must be a positive integer.')ifisinstance(harmonics,float)orisinstance(harmonics,int):harmonics=[int(harmonics)]ifnotisinstance(harmonics,list)andharmonicsisnotNone:raiseTypeError('The harmonics must be specified as a list of integer values.')ifexperimentisnotNone:ifnotisinstance(experiment,ExperimentInfo):raiseTypeError('The experiment must be a valid deerlab.ExperimentInfo object.')# Check that the number of requested pathways does not exceed the theoretical limit of the experimentnpathways=experiment.npathwaysmaxpathways=len(experiment.harmonics)ifnpathways>maxpathways:raiseValueError(f'The {experiment.name} experiment can only have up to {maxpathways} dipolar pathways.')#------------------------------------------------------------------------def_importparameter(parameter):""" Private function for importing a parameter's metadata """return{'lb':parameter.lb,'ub':parameter.ub,'par0':parameter.par0,'description':parameter.description,'unit':parameter.unit,'linear':parameter.linear}#------------------------------------------------------------------------# Parse the harmonics of the dipolar pathwaysifharmonicsisNone:ifexperimentisnotNone:harmonics=np.array(experiment.harmonics)[:npathways]else:harmonics=np.ones(npathways)iflen(harmonics)!=npathways:raiseValueError('The number of harmonics must match the number of dipolar pathways.')ifparametrization!='reftimes':ifexperimentisNone:raiseSyntaxError(f"The '{parametrization}' parametrization requires 'experiment' to be specified.")# Determine number of dipolar interactions in the spins systemQ=int(spins*(spins-1)/2)ifspins>2:triangles_required=factorial(spins)/(factorial(3)*factorial(spins-3))ifspins==3andtrianglesisNone:triangles=[[1,2,3]]elifspins>3andtrianglesisNone:raiseSyntaxError(f'For systems with {spins} spins, a list of {triangles_required:n} triangles is required.')elifspins>3andlen(triangles)!=triangles_required:raiseSyntaxError(f'For systems with {spins} spins, a list of {triangles_required:n} triangles is required.')ifnp.min(triangles)!=0:triangles=[np.array(triangle)-1fortriangleintriangles]#------------------------------------------------------------------------defdipolarpathways(*param):""" Generator of dipolar pathways based on a parameter set. Constructs all possible two-spin dipolar pathways based on a parameter set, and then constructs all possible multi-spin dipolar pathways describing three-spin interactions in systems with more than two spins based on the pair-wise dipolar pathways. The true signature of this function is defined dynamically after its definition below. """param=np.atleast_1d(param)ifspins==2:# Construct pair-wise dipolar pathways given by the parameter setifparametrization=='reftimes':λs=param[np.arange(0,len(param),2)]trefs=param[np.arange(1,len(param),2)]elifparametrization=='delays':delays=param[0:len(experiment.delays)]λs=param[len(delays):]trefs=experiment.reftimes(*delays)elifparametrization=='shift':shift=param[0]λs=param[1:]delays=np.array(experiment.delays)+shifttrefs=experiment.reftimes(*delays)# Unmodulated pathwayΛ0=np.maximum(0,1-np.sum(λs))pairpathways=[{'amp':Λ0}]# Modulated pathwaysforλ,tref,δinzip(λs,trefs,harmonics):pairpathways.append({'amp':λ,'reftime':tref,'harmonic':δ})# For two-spin systems these are all the required pathwaysreturnpairpathways# Construct multispin dipolar pathways given by the parameter setλs=[param[np.arange(q,len(param),Q)]forqinrange(Q)]trefs=param[np.arange(Q,len(param),Q)]# Amplitude of the unmodulated pair-wise dipolar pathwayλu=param[-1]# Initialize containerspathways,threespin_pathways=[],[]Λ0=1# Construct all possible multi-spin dipolar pathwaysforidxinrange(npathways):# Construct all the permutations of the two-spin interaction pathways (without repetitions)forperminset(set(permutations([0]+[None]*(Q-1)))):q=int(np.where(np.array(perm)==0)[0])# Compute amplitude of the two-spin interaction pathway λp=λs[q][idx]λ2k=factorial(Q-1)*λp*λu**(Q-1)tref=[trefs[idx]ifn==0elseNoneforninperm]δs=[harmonics[idx]ifn==0else0forninperm]# Add two-spin interaction dipolar pathway to the listpathways.append({'amp':λ2k,'reftime':tuple(tref),'harmonic':tuple(δs)})# Update the unmodulated amplitudeΛ0-=λ2k# Consider now three-spin interactions foridx2innp.arange(idx+1,npathways):# If the pathway has negligible amplitude, ignore it (speed-up) ifλ3k>minamp:# Construct all the permutations of the multi-spin pathway (without repetitions) forperminset(set(permutations([0,1]+[None]*(Q-2)))):# Check if the permutation corresponds to a valid interaction triangleq1,q2=[np.where(perm==n)[0]fornin[0,1]]fortriangleintriangles:ifq1notintriangleandq2notintriangle:continue# If not valid, move on to the next one# Compute amplitude of the three-spin interaction pathway λm=λs[q2][idx]λ3k=2*factorial(Q-2)*λp*λm*λu**(Q-2)tref=[trefs[idx]ifn==0elsetrefs[idx2]ifn==1elseNoneforninperm]δs=[harmonics[idx]ifn==0elseharmonics[idx2]ifn==1else0forninperm]# Add three-spin interaction dipolar pathway to the listthreespin_pathways.append({'amp':λ3k,'reftime':tuple(tref),'harmonic':tuple(δs)})# Update the unmodulated amplitudeΛ0-=λ3k# Add pathway with unmodulated contributionpathways.append({'amp':Λ0})returnpathways,threespin_pathways#------------------------------------------------------------------------# Construct the signature of the dipolarpathways() function dynamicallyifexperimentisnotNone:labels=experiment.labelselse:# Otherwise, just label them in orderlabels=np.arange(1,npathways+1)variables=[]ifparametrization=='delays':pulsedelay_names=inspect.signature(experiment.reftimes).parameters.keys()variables=list(pulsedelay_names)ifparametrization=='shift':variables=['tshift']forninrange(npathways):ifnpathways==1andspins==2:# Use single-pathway notation, with amplitude as modulation depthvariables.append('mod')else:ifspins==2:variables.append(f'lam{labels[n]}')else:forqinrange(Q):variables.append(f'lam{labels[n]}_{q+1}')ifparametrization=='reftimes':ifnpathways==1andspins==2:variables.append(f'reftime')else:variables.append(f'reftime{labels[n]}')# Add unmodulated pairwise pathway amplitude for multi-spin systemsifspins>2:variables.append('lamu')# Create the dipolar pathways model objectPathsModel=Model(dipolarpathways,signature=variables)ifspins>2andsamespins:links={f'lam{labels[n]}':[f'lam{labels[n]}_{q+1}'forqinrange(Q)]forninrange(npathways)}PathsModel=link(PathsModel,**links)#-----------------------------------------------------------------------def_Pmultivar(param):rmeans,cholesky_factors=param[:Q],param[Q:]Σ=choleskycovmat(Q,cholesky_factors)Pmultivar=multivariate_normal(rmeans,cov=Σ)returnPmultivar#-----------------------------------------------------------------------# Determine if distance distributions is non-parametricPnonparametric=PmodelisNoneifPnonparametric:Pmodel=freedist(r)# For multi-spin systems, there is no longer freedom of choice for distance distributionsifspins>2:# Prepare the normal multivariate distance distribution modelvariables=[]variables+=[f'rmean{q+1}'forqinrange(Q)]variables+=[f'chol{q+1}{q+1}'forqinrange(Q)]forqinrange(Q):variables+=[f'chol{j+1}{q+1}'forjinnp.arange(q+1,Q)]Pmodel=Model(_Pmultivar,signature=variables)[getattr(Pmodel,f'rmean{q+1}').set(lb=0,ub=10,par0=3,unit='nm',description=f'Average inter-spin distance #{q+1}')forqinrange(Q)][getattr(Pmodel,f'chol{q+1}{q+1}').set(lb=0,ub=6,par0=0.4,unit='nm',description=f'Cholesky factor ℓ{q+1}{q+1}')forqinrange(Q)]forqinrange(Q):[getattr(Pmodel,f'chol{j+1}{q+1}').set(lb=-1,ub=1,par0=0.0,unit='nm',description=f'Cholesky factor ℓ{j+1}{q+1}')forjinnp.arange(q+1,Q)]Nconstants=len(Pmodel._constantsInfo)# Populate the basic information on the dipolar pathways parametersforninrange(npathways):ifspins==2andnpathways==1:# Special case: use modulation depth notation instead of general pathway amplitudegetattr(PathsModel,f'mod').set(lb=0,ub=1,par0=0.01,description=f'Modulation depth',unit='')else:pairwise=''ifspins>2:pairwise=' pairwise'ifspins==2orsamespins:getattr(PathsModel,f'lam{labels[n]}').set(lb=0,ub=1,par0=0.01,description=f'Amplitude of{pairwise} pathway #{labels[n]}',unit='')else:forqinrange(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='')ifparametrization=='reftimes':ifexperimentisNone:theoretical_reftime=0reftime_variability=20else:theoretical_reftime=experiment.reftimes(*experiment.delays)[n]reftime_variability=3*experiment.pulselengthifspins==2andnpathways==1:# Special case: use reftime notationgetattr(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 notationgetattr(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)ifparametrization=='delays':experimental_delays=experiment.delaysdelay_variability=5*experiment.pulselengthforn,delayinenumerate(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)elifparametrization=='shift':variability=5*experiment.pulselengthgetattr(PathsModel,'tshift').set(par0=0,lb=-variability,ub=variability,description=f'Variability of experimental pulse delays',unit='μs')ifspins>2:getattr(PathsModel,f'lamu').set(lb=0,ub=1,par0=0.5,description='Amplitude of unmodulated pairwise pathway',unit='')# Construct the signature of the dipolar signal model functionsignature=[]parameters,linearparam=[],[]formodelin[PathsModel,Bmodel,Pmodel]:ifmodelisnotNone:forparaminmodel._parameter_list(order='vector'):ifnp.any(getattr(model,param).linear):parameters.append(getattr(model,param))linearparam.append({'name':param,'vec':len(np.atleast_1d(getattr(model,param).idx)),'normalization':getattr(model,param).normalization})elifnot(model==Bmodelandparam=='lam'):signature.append(param)parameters.append(getattr(model,param))# Initialize lists of indices to access subsets of nonlinear parameters Psubset=np.zeros(Pmodel.Nnonlin,dtype=int)PathsSubset=np.zeros(PathsModel.Nnonlin,dtype=int)ifBmodelisNone:Bsubset=[]else:Bsubset=np.zeros(Bmodel.Nnonlin-('lam'inBmodel._parameter_list()),dtype=int)# Determine subset indices based on main function signatureidx=0formodel,subsetinzip([Pmodel,Bmodel,PathsModel],[Psubset,Bsubset,PathsSubset]):ifmodelisnotNone:foridx,paraminenumerate(signature):ifparaminmodel._parameter_list(order='vector'):subset[getattr(model,param).idx]=idxkernelmethod='fresnel'iforiselisNoneandnp.isinf(excbandwidth)else'grid'dt=np.mean(np.diff(t))tinterp=None#------------------------------------------------------------------------defVtwospin_nonlinear_fcn(*nonlin):""" Non-linear part of the dipolar signal function """# Make input arguments as array to access subsets easilynonlin=np.atleast_1d(nonlin)# Construct the basis function of the intermolecular contributionifBmodelisNone:Bfcn=np.ones_like(t)elifhasattr(Bmodel,'lam'):Bfcn=lambdat,lam:Bmodel.nonlinmodel(t,*np.concatenate([nonlin[Bsubset],[lam]]))else:Bfcn=lambdat,_:Bmodel.nonlinmodel(t,*nonlin[Bsubset])# Construct the definition of the dipolar pathwayspathways=PathsModel.nonlinmodel(*nonlin[PathsSubset])# Construct the dipolar kernelKdipolar=dipolarkernel(t,r,pathways=pathways,bg=Bfcn,tinterp=tinterp,excbandwidth=excbandwidth,orisel=orisel,g=g,method=kernelmethod,gridsize=gridsize)# Compute the non-linear part of the distance distributionPnonlin=Pmodel.nonlinmodel(*[r]*Nconstants,*nonlin[Psubset])# Forward calculation of the non-linear part of the dipolar signalVnonlin=Kdipolar@PnonlinreturnVnonlin#------------------------------------------------------------------------#------------------------------------------------------------------------defVmultispin_nonlinear_fcn(*nonlin):""" Non-linear part of the dipolar signal function """# Make input arguments as array to access subsets easilynonlin=np.atleast_1d(nonlin)# Construct the basis function of the intermolecular contributionifBmodelisNone:Bfcn=np.ones_like(t)elifhasattr(Bmodel,'lam'):Bfcn=lambdat,lam:Bmodel.nonlinmodel(t,*np.concatenate([nonlin[Bsubset],[lam]]))else:Bfcn=lambdat,_:Bmodel.nonlinmodel(t,*nonlin[Bsubset])# Sample from the multi-variate distance distributionNsamples=1000000rsamples=Pmodel.nonlinmodel(nonlin[Psubset]).rvs(Nsamples)rsamples=np.maximum(rsamples,1e-16)# Avoid values exactly at zero # Enforce the generalized triangle inequialityfortriangleintriangles:# Loop over all triangle combinationsidx=triangletriangle_condition=np.full_like(rsamples,False)forqinrange(Q):# Evaluate all triangle inequalitiesidx=np.roll(idx,-1,axis=0)triangle_condition[:,q]=np.sum(rsamples[:,idx[:-1]],axis=1)>rsamples[:,idx[-1]]triangle_condition=np.all(triangle_condition,axis=1)# Discard samples that do not satisfy triangle inequalitiesrsamples=rsamples[triangle_condition,:]# Generate all multi-spin dipolar pathwaystwospin_pathways,threespin_pathways=PathsModel.nonlinmodel(*nonlin[PathsSubset])# Start by accounting for the unmodulated contributionΛ0=twospin_pathways.pop(-1)['amp']Vintra=Λ0# Account for all two-spin intramolecular contributionsforqinrange(Q):# Get the marginal univariate distance distributionPunimarg,bins=np.histogram(rsamples[:,q],bins=100)rq=(bins[:-1]+bins[1:])/2rs=[rq]*Qifnp.sum(Punimarg)>0:Punimarg=Punimarg/np.sum(Punimarg)# Get the two-spin pathways modulated by the q-th dipolar interactionpathways_q=[pathwayforpathwayintwospin_pathwaysifpathway['harmonic'][q]!=0]# Construct the corresponding dipolar kernelKtwospin=dipolarkernel(t,rs,pathways=pathways_q,integralop=False,tinterp=tinterp,excbandwidth=excbandwidth,orisel=orisel,g=g,method=kernelmethod,gridsize=gridsize)# Two-spin intramolecular contribution from all pathways modulated by the q-th interactionVintra+=Ktwospin@Punimarg# Account for all three-spin intramolecular contributionsifthreespin_pathways:fortriangleintriangles:# Get the marginal trivariate distance distributionPtrimarg,bins=np.histogramdd(rsamples[:,triangle],bins=10,density=True)r1,r2,r3=[(bins_q[:-1]+bins_q[1:])/2forbins_qinbins]# Get the subset of the most distance combinations with most distribution mass (speed-up)subset=Ptrimarg>np.max(Ptrimarg)/100Ptrimarg=Ptrimarg[subset]ifnp.sum(Ptrimarg)>0:Ptrimarg=Ptrimarg/np.sum(Ptrimarg)r1,r2,r3=[rq[subset]forrqinnp.meshgrid(r1,r2,r3,indexing='ij')]# Get the three-spin pathways modulated by the triangle interactionpathways_triangle=[pathwayforpathwayinthreespin_pathwaysifnp.sum([pathway['harmonic'][q]!=0forqintriangle])==2]# Construct the corresponding dipolar kernelKthreespin=dipolarkernel(t,[r1,r2,r3],pathways=pathways_triangle,integralop=False,excbandwidth=excbandwidth,orisel=orisel,g=g,method='grid',gridsize=gridsize)# Three-spin intramolecular contribution from all pathways modulated by the triangle interactionVintra+=Kthreespin@Ptrimarg# Compute the multi-spin dipolar intermolecular contributionVinter=dipolarbackground(t,twospin_pathways+threespin_pathways,Bfcn)# Construct the full dipolar signal Vnonlin=Vintra*VinterreturnVnonlin#------------------------------------------------------------------------# Create the dipolar model objectifspins==2:DipolarSignal=Model(Vtwospin_nonlinear_fcn,signature=signature)else:DipolarSignal=Model(Vmultispin_nonlinear_fcn,signature=signature)# Add the linear parameters from the subset models forlparaminlinearparam:DipolarSignal.addlinear(lparam['name'],vec=lparam['vec'],normalization=lparam['normalization'])ifPmodelisNone:DipolarSignal.addlinear()# If there are no linear paramters, add a linear scaling parameter ifDipolarSignal.Nlin==0:DipolarSignal.addlinear('scale',lb=0,par0=1,description='Overall echo amplitude/scale')# Import all parameter information from the subset modelsforname,paraminzip(DipolarSignal._parameter_list(order='vector'),parameters):getattr(DipolarSignal,name).set(**_importparameter(param))DipolarSignal.description='Dipolar signal model'# Set prior knowledge on the parameters if experiment is specifiedifexperimentisnotNone:# Specify experiment in model descriptionDipolarSignal.description=f'{experiment.name} dipolar signal model'# Specify start values and boundaries according to experimental timingsreftimes=experiment.reftimes(*experiment.delays)ifinterp:tinterp=np.arange(min(t)-max(reftimes)-3*experiment.pulselength-dt,max(t)+max(reftimes)+3*experiment.pulselength+dt,dt)returnDipolarSignal
[docs]defdipolarpenalty(Pmodel,r,type,selection=None):r""" Construct penalties based on the distance distribution. This function constructs penalties for a given distance distribution that can be used to impose certain desired properties on the distribution. The type of penalty and the function to optimize it are determined by the ``type`` parameter. The ``selection`` parameter allows the user to choose a method for optimizing the penalty weight. Parameters ---------- Pmodel : `Model` object The model of the distance distribution. If ``Pmodel`` is ``None``, a non-parametric distance distribution is assumed. r : array_like Distance axis vector, in nanometers. This vector specifies the distances at which the distance distribution is defined. type : string Type of property to be imposed by the penalty. This parameter determines the type of penalty that will be constructed, as well as the function that will be optimized to determine the penalty weight. The following values are supported for this parameter: - ``'smoothness'`` : Imposes a smoothness constraint on the distance distribution. This penalty is optimized by minimizing the L2 norm of the second derivative of the distribution. - ``'compactness'`` : Imposes a compactness constraint on the distance distribution. This penalty is optimized by minimizing the square root of the product of the distance distribution and the square of the difference between the distances and their mean. selection : string, optional Selection functional for the outer optimization of the penalty weight. Possible values: - ``'aic'`` - Akaike information criterion - ``'bic'`` - Bayesian information criterion - ``'aicc'`` - COrrected Akaike information criterion - ``'icc'`` - Informational complexity criterion If no value is provided, ``'icc'`` is used for ``'compactness'`` and ``'aic'`` for ``'smoothness'``. Returns ------- penalty : ``Penalty`` object Penalty object to be passed to the ``fit`` function. """ifPmodelisNone:Pmodel=freedist(r)Nconstants=len(Pmodel._constantsInfo)# If include compactness penaltyiftype=='compactness':ifselectionisNone:selection='icc'# Define the compactness penalty functiondefcompactness_penalty(*args):P=Pmodel(*[r]*Nconstants,*args)ifnotnp.all(P==0):P=P/np.trapz(P,r)returnnp.sqrt(P*(r-np.trapz(P*r,r))**2*np.mean(np.diff(r)))# Add the penalty to the Pmodelpenalty=Penalty(compactness_penalty,selection,signature=Pmodel._parameter_list(),description='Distance distribution compactness penalty.')penalty.weight.set(lb=1e-6,ub=1e1)# If include smoothness penaltyeliftype=='smoothness':ifselectionisNone:selection='aic'# Define the smoothness penalty functionL=regoperator(r,2)defsmoothness_penalty(*args):returnL@Pmodel(*[r]*Nconstants,*args)# Add the penalty to the Pmodelpenalty=Penalty(smoothness_penalty,selection,signature=Pmodel._parameter_list(),description='Distance distribution smoothness penalty.')penalty.weight.set(lb=1e-9,ub=1e3)else:raiseKeyError(f"The requested {type} is not a valid penalty. Must be 'compactness' or 'smoothness'.")returnpenalty
#===============================================================================#===============================================================================def_checkpathways(pathways,Nmax):# Check that pathways are correctly specifiediflen(pathways)>Nmax:raiseValueError(f"The number of pathways cannot exceed {Nmax}.")ifnp.any(np.array(pathways)<1)ornotnp.all([notp%1forpinpathways])ornp.any(np.array(pathways)>Nmax):raiseValueError(f"The pathways must be specified by integer numbers in the range 1-{Nmax}.")iflen(np.unique(pathways))!=len(pathways):raiseValueError("Some pathways are duplicated.")#===============================================================================#===============================================================================
[docs]classExperimentInfo():r""" Represents information about a dipolar EPR experiment Attributes ---------- name : string Name of the experiment. reftimes : list of floats List of refocusing times, in microseconds. harmonics : list List of harmonics of the refocusing times. pulselength : float Length of the longest microwave pulse in the sequence in microseconds. pathwaylabels : list List of pathway labels. delays : list List of pulse delays. """
[docs]def__init__(self,name,reftimes,harmonics,pulselength,pathwaylabels,delays):""" Populates the attributes of the ExperimentInfo object. Parameters ---------- name : string Name of the experiment. reftimes : list of floats List of refocusing times, in microseconds. harmonics : list List of harmonics of the refocusing times. pulselength : float Length of the longest microwave pulse in the sequence in microseconds. pathwaylabels : list List of pathway labels. delays : list List of pulse delays. """self.npathways=len(harmonics)self.reftimes=reftimesself.harmonics=harmonicsself.pulselength=pulselengthself.labels=pathwaylabelsself.delays=delaysself.name=name
[docs]defex_3pdeer(tau,pathways=[1,2],pulselength=0.016):r""" Generate a 3-pulse DEER dipolar experiment model. The figure below shows the dipolar pathways in 3-pulse DEER. The observer (black) and pump (grey) pulses and their interpulse delays are shown on the top. The middle table summarizes all detectable modulated dipolar pathways `p` along their dipolar phase accumulation factors `\mathbf{s}_p`, harmonics `\delta_p` and refocusing times `t_{\mathrm{ref},p}`. The most commonly encountered pathways are highlighted in color. The bottom panel shows a decomposition of the dipolar signal into the individual intramolecular contributions (shown as colored lines). .. figure:: ../images/ex_3pdeer_pathways.png :width: 350px Source: L. Fábregas-Ibáñez, Advanced data analysis and modeling in dipolar EPR spectroscopy, Doctoral dissertation, 2022 Parameters ---------- tau : float scalar Static interpulse delay `\tau`. pathways : array_like, optional Pathways to include in the model. The pathways are specified as a list of pathways labels `p`. By default, both pathways are included as shown in the table above. pulselength : float scalar, optional Length of the longest microwave pulse in the sequence in microseconds. Used to determine the uncertainty in the boundaries of the pathway refocusing times. Returns ------- experiment : ``ExperimentInfo`` object Experiment object. Can be passed to ``dipolarmodel`` to introduce better constraints into the model. """# Theoretical refocusing timesdefreftimes(tau):tref=[tau,0]# Sort according to pathways orderifpathwaysisnotNone:tref=[tref[pathway-1]forpathwayinpathways]returntref# Pulse delays delays=[tau]# Theoretical dipolar harmonicsharmonics=[1,1]# Pathway labelspathwaylabels=np.arange(1,len(harmonics)+1)# Sort according to pathways orderifpathwaysisnotNone:_checkpathways(pathways,Nmax=len(harmonics))harmonics=[harmonics[pathway-1]forpathwayinpathways]pathwaylabels=[pathwaylabels[pathway-1]forpathwayinpathways]returnExperimentInfo('3-pulse DEER',reftimes,harmonics,pulselength,pathwaylabels,delays)
[docs]defex_4pdeer(tau1,tau2,pathways=[1,2,3,4],pulselength=0.016):r""" Generate a 4-pulse DEER dipolar experiment model. The figure below shows the dipolar pathways in 4-pulse DEER. The observer (black) and pump (grey) pulses and their interpulse delays are shown on the top. The middle table summarizes all detectable modulated dipolar pathways `p` along their dipolar phase accumulation factors `\mathbf{s}_p`, harmonics `\delta_p` and refocusing times `t_{\mathrm{ref},p}`. The most commonly encountered pathways are highlighted in color. The bottom panel shows a decomposition of the dipolar signal into the individual intramolecular contributions (shown as colored lines). The dipolar time axis is defined such that `t=0` right after the second observer pulse. .. figure:: ../images/ex_4pdeer_pathways.png :width: 350px Source: L. Fábregas-Ibáñez, Advanced data analysis and modeling in dipolar EPR spectroscopy, Doctoral dissertation, 2022 Parameters ---------- tau1 : float scalar 1st static interpulse delay `\tau_1`. tau2 : float scalar 2nd static interpulse delay `\tau_2`. pathways : array_like, optional Pathways to include in the model. The pathways are specified as a list of pathways labels `p`. By default, pathways 1-4 are included as shown in the table above. pulselength : float scalar, optional Length of the longest microwave pulse in the sequence in microseconds. Used to determine the uncertainty in the boundaries of the pathway refocusing times. Returns ------- experiment : ``ExperimentInfo`` object Experiment object. Can be passed to ``dipolarmodel`` to introduce better constraints into the model. """# Theoretical refocusing timesdefreftimes(tau1,tau2):tref=[tau1,tau1+tau2,0,tau2]# Sort according to pathways orderifpathwaysisnotNone:tref=[tref[pathway-1]forpathwayinpathways]returntref# Pulse delays delays=[tau1,tau2]# Theoretical dipolar harmonicsharmonics=[1,1,1,1]# Pathway labelspathwaylabels=np.arange(1,len(harmonics)+1)# Sort according to pathways orderifpathwaysisnotNone:_checkpathways(pathways,Nmax=len(harmonics))harmonics=[harmonics[pathway-1]forpathwayinpathways]pathwaylabels=[pathwaylabels[pathway-1]forpathwayinpathways]returnExperimentInfo('4-pulse DEER',reftimes,harmonics,pulselength,pathwaylabels,delays)
[docs]defex_rev5pdeer(tau1,tau2,tau3,pathways=[1,2,3,4,5],pulselength=0.016):r""" Generate a reverse 5-pulse DEER dipolar experiment model. The figure below shows the dipolar pathways in reverse 5-pulse DEER. The observer (black) and pump (grey) pulses and their interpulse delays are shown on the top. The middle table summarizes all detectable modulated dipolar pathways `p` along their dipolar phase accumulation factors `\mathbf{s}_p`, harmonics `\delta_p` and refocusing times `t_{\mathrm{ref},p}`. The most commonly encountered pathways are highlighted in color. The bottom panel shows a decomposition of the dipolar signal into the individual intramolecular contributions (shown as colored lines). .. figure:: ../images/ex_rev5pdeer_pathways.png :width: 350px Source: L. Fábregas-Ibáñez, Advanced data analysis and modeling in dipolar EPR spectroscopy, Doctoral dissertation, 2022 Parameters ---------- tau1 : float scalar 1st static interpulse delay `\tau_1`. tau2 : float scalar 2nd static interpulse delay `\tau_2`. tau3 : float scalar 3rd static interpulse delay `\tau_3`. pathways : array_like, optional Pathways to include in the model. The pathways are specified as a list of pathways labels `p`. By default, pathways 1-5 are included as shown in the table above. pulselength : float scalar, optional Length of the longest microwave pulse in the sequence in microseconds. Used to determine the uncertainty in the boundaries of the pathway refocusing times. Returns ------- experiment : ``ExperimentInfo`` object Experiment object. Can be passed to ``dipolarmodel`` to introduce better constraints into the model. """# Theoretical refocusing timesdefreftimes(tau1,tau2,tau3):tref=[tau3,tau2-tau3,tau1+tau3,tau1+tau2-tau3,tau2,tau1+tau2,0,tau1]# Sort according to pathways orderifpathwaysisnotNone:tref=[tref[pathway-1]forpathwayinpathways]returntref# Pulse delays delays=[tau1,tau2,tau3]# Theoretical dipolar harmonicsharmonics=[1,1,1,1,1,1,1,1]# Pathway labelspathwaylabels=np.arange(1,len(harmonics)+1)# Sort according to pathways orderifpathwaysisnotNone:_checkpathways(pathways,Nmax=len(harmonics))harmonics=[harmonics[pathway-1]forpathwayinpathways]pathwaylabels=[pathwaylabels[pathway-1]forpathwayinpathways]returnExperimentInfo('Reverse 5-pulse DEER',reftimes,harmonics,pulselength,pathwaylabels,delays)
[docs]defex_fwd5pdeer(tau1,tau2,tau3,pathways=[1,2,3,4,5],pulselength=0.016):r""" Generate a forward 5-pulse DEER dipolar experiment model. The figure below shows the dipolar pathways in forward 5-pulse DEER. The observer (black) and pump (grey) pulses and their interpulse delays are shown on the top. The middle table summarizes all detectable modulated dipolar pathways `p` along their dipolar phase accumulation factors `\mathbf{s}_p`, harmonics `\delta_p` and refocusing times `t_{\mathrm{ref},p}`. The most commonly encountered pathways are highlighted in color. The bottom panel shows a decomposition of the dipolar signal into the individual intramolecular contributions (shown as colored lines). .. figure:: ../images/ex_fwd5pdeer_pathways.png :width: 350px Source: L. Fábregas-Ibáñez, Advanced data analysis and modeling in dipolar EPR spectroscopy, Doctoral dissertation, 2022 Parameters ---------- tau1 : float scalar 1st static interpulse delay `\tau_1`. tau2 : float scalar 2nd static interpulse delay `\tau_2`. tau3 : float scalar 3rd static interpulse delay `\tau_3`. pathways : array_like, optional Pathways to include in the model. The pathways are specified as a list of pathways labels `p`. By default, pathways 1-5 are included as shown in the table above. pulselength : float scalar, optional Length of the longest microwave pulse in the sequence in microseconds. Used to determine the uncertainty in the boundaries of the pathway refocusing times. Returns ------- experiment : ``ExperimentInfo`` object Experiment object. Can be passed to ``dipolarmodel`` to introduce better constraints into the model. """# Theoretical refocusing timesdefreftimes(tau1,tau2,tau3):tref=[tau3,tau1-tau3,tau2+tau3,tau1+tau2-tau3,tau1,tau1+tau2,0,tau2]# Sort according to pathways orderifpathwaysisnotNone:tref=[tref[pathway-1]forpathwayinpathways]returntref# Pulse delays delays=[tau1,tau2,tau3]# Theoretical dipolar harmonicsharmonics=[1,1,1,1,1,1,1,1]# Pathway labelspathwaylabels=np.arange(1,len(harmonics)+1)# Sort according to pathways orderifpathwaysisnotNone:_checkpathways(pathways,Nmax=len(harmonics))harmonics=[harmonics[pathway-1]forpathwayinpathways]pathwaylabels=[pathwaylabels[pathway-1]forpathwayinpathways]returnExperimentInfo('Forward 5-pulse DEER',reftimes,harmonics,pulselength,pathwaylabels,delays)
[docs]defex_sifter(tau1,tau2,pathways=[1,2,3],pulselength=0.016):r""" Generate a 4-pulse SIFTER dipolar experiment model. The figure below shows the dipolar pathways in 4-pulse SIFTER. The pulses and the interpulse delays are shown on the top. The middle table summarizes all detectable modulated dipolar pathways `p` along their dipolar phase accumulation factors `\mathbf{s}_p`, harmonics `\delta_p` and refocusing times `t_{\mathrm{ref},p}`. The most commonly encountered pathways are highlighted in color. The bottom panel shows a decomposition of the dipolar signal into the individual intramolecular contributions (shown as colored lines). .. figure:: ../images/ex_sifter_pathways.png :width: 350px Source: L. Fábregas-Ibáñez, Advanced data analysis and modeling in dipolar EPR spectroscopy, Doctoral dissertation, 2022 Parameters ---------- tau1 : float scalar 1st static interpulse delay `\tau_1`. tau2 : float scalar 2nd static interpulse delay `\tau_2`. pathways : array_like, optional Pathways to include in the model. The pathways are specified as a list of pathways labels `p`. By default, all 3 pathways are included as shown in the table above. pulselength : float scalar, optional Length of the longest microwave pulse in the sequence in microseconds. Used to determine the uncertainty in the boundaries of the pathway refocusing times. Returns ------- experiment : ``ExperimentInfo`` object Experiment object. Can be passed to ``dipolarmodel`` to introduce better constraints into the model. """# Theoretical refocusing timesdefreftimes(tau1,tau2):tref=[tau2-tau1,2*tau2,-2*tau1]# Sort according to pathways orderifpathwaysisnotNone:tref=[tref[pathway-1]forpathwayinpathways]returntref# Pulse delays delays=[tau1,tau2]# Theoretical dipolar harmonicsharmonics=[1,1/2,1/2]# Pathway labelspathwaylabels=np.arange(1,len(harmonics)+1)# Sort according to pathways orderifpathwaysisnotNone:_checkpathways(pathways,Nmax=len(harmonics))harmonics=[harmonics[pathway-1]forpathwayinpathways]pathwaylabels=[pathwaylabels[pathway-1]forpathwayinpathways]returnExperimentInfo('4-pulse SIFTER',reftimes,harmonics,pulselength,pathwaylabels,delays)
[docs]defex_ridme(tau1,tau2,pathways=[1,2,3,4],pulselength=0.016):r""" Generate a 5-pulse RIDME dipolar experiment model. The figure below shows the dipolar pathways in 5-pulse RIDME. The pulses and their interpulse delays are shown on the top. The middle table summarizes all detectable modulated dipolar pathways `p` along their dipolar phase accumulation factors `\mathbf{s}_p`, harmonics `\delta_p` and refocusing times `t_{\mathrm{ref},p}`. The most commonly encountered pathways are highlighted in color. The bottom panel shows a decomposition of the dipolar signal into the individual intramolecular contributions (shown as colored lines). The dipolar time axis is defined such that `t=0` right after the second pulse. (Note that the model does not account for any relaxation-induced effects) .. figure:: ../images/ex_ridme_pathways.png :width: 350px Parameters ---------- tau1 : float scalar 1st static interpulse delay `\tau_1`. tau2 : float scalar 2nd static interpulse delay `\tau_2`. pathways : array_like, optional Pathways to include in the model. The pathways are specified as a list of pathways labels `p`. By default, pathways 1-4 are included as shown in the table above. pulselength : float scalar, optional Length of the longest microwave pulse in the sequence, in microseconds. Used to determine the uncertainty in the boundaries of the pathway refocusing times. Returns ------- experiment : ``ExperimentInfo`` object Experiment object. Can be passed to ``dipolarmodel`` to introduce better constraints into the model. """# Theoretical refocusing timesdefreftimes(tau1,tau2):tref=[tau1,tau1+tau2,0,tau2]# Sort according to pathways orderifpathwaysisnotNone:tref=[tref[pathway-1]forpathwayinpathways]returntref# Pulse delays delays=[tau1,tau2]# Theoretical dipolar harmonicsharmonics=[1,1,1,1]# Pathway labelspathwaylabels=np.arange(1,len(harmonics)+1)# Sort according to pathways orderifpathwaysisnotNone:_checkpathways(pathways,Nmax=len(harmonics))harmonics=[harmonics[pathway-1]forpathwayinpathways]pathwaylabels=[pathwaylabels[pathway-1]forpathwayinpathways]returnExperimentInfo('5-pulse RIDME',reftimes,harmonics,pulselength,pathwaylabels,delays)
[docs]defex_dqc(tau1,tau2,tau3,pathways=[1,2,3],pulselength=0.016):r""" Generate a 6-pulse DQC dipolar experiment model. The figure below shows the dipolar pathways in 6-pulse DQC. The pulses and their interpulse delays are shown on the top. The middle table summarizes all detectable modulated dipolar pathways `p` along their dipolar phase accumulation factors `\mathbf{s}_p`, harmonics `\delta_p` and refocusing times `t_{\mathrm{ref},p}`. The most commonly encountered pathways are highlighted in color. The bottom panel shows a decomposition of the dipolar signal into the individual intramolecular contributions (shown as colored lines). .. figure:: ../images/ex_dqc_pathways.png :width: 700px Source: L. Fábregas-Ibáñez, Advanced data analysis and modeling in dipolar EPR spectroscopy, Doctoral dissertation, 2022 Parameters ---------- tau1 : float scalar 1st static interpulse delay `\tau_1`. tau2 : float scalar 2nd static interpulse delay `\tau_2`. tau3 : float scalar 3rd static interpulse delay `\tau_3`. pathways : array_like, optional Pathways to include in the model. The pathways are specified as a list of pathways labels `p`. By default, pathways 1, 2, and 3 are included as shown in the table above. pulselength : float scalar, optional Length of the longest microwave pulse in the sequence in microseconds. Used to determine the uncertainty in the boundaries of the pathway refocusing times. Returns ------- experiment : ``ExperimentInfo`` object Experiment object. Can be passed to ``dipolarmodel`` to introduce better constraints into the model. """# Theoretical refocusing timesdefreftimes(tau1,tau2,tau3):tref=[tau2-tau1,2*tau2,-2*tau1,tau2-tau1-tau3,-2*(tau1-tau3),-2*(tau1+tau3),2*(tau2+tau3),2*(tau1-tau3)]# Sort according to pathways orderifpathwaysisnotNone:tref=[tref[pathway-1]forpathwayinpathways]returntref# Pulse delays delays=[tau1,tau2,tau3]# Theoretical dipolar harmonicsharmonics=[1,1/2,1/2,1,1/2,1/2,1/2,1/2]# Pathway labelspathwaylabels=np.arange(1,len(harmonics)+1)# Sort according to pathways orderifpathwaysisnotNone:_checkpathways(pathways,Nmax=len(harmonics))harmonics=[harmonics[pathway-1]forpathwayinpathways]pathwaylabels=[pathwaylabels[pathway-1]forpathwayinpathways]returnExperimentInfo('6-pulse DQC',reftimes,harmonics,pulselength,pathwaylabels,delays)