diff --git a/deerlab/fit.py b/deerlab/fit.py index 967e08f9..b7929b8c 100644 --- a/deerlab/fit.py +++ b/deerlab/fit.py @@ -358,6 +358,8 @@ def fit(model_, y, *constants, par0=None, penalties=None, bootstrap=0, noiselvl= Fitted parameter vector ordered according to the model parameter indices. paramUncert : :ref:`UQResult` Uncertainty quantification of the parameter vector ordered according to the model parameter indices. + paramlist : list + List of the fitted parameter names ordered according to the model parameter indices. model : ndarray Fitted model response. regparam : scalar @@ -505,6 +507,15 @@ def bootstrap_fcn(ysim): FitResult_dict = {key: getattr(fitresults,key) for key in ['y','mask','param','paramUncert','model','cost','plot','residuals','stats','regparam','regparam_stats','__plot_inputs']} _paramlist = model._parameter_list('vector') + param_idx = [[] for _ in _paramlist] + idxprev = 0 + for islinear in [False,True]: + for n,param in enumerate(_paramlist): + if np.all(getattr(model,param).linear == islinear): + N = len(np.atleast_1d(getattr(model,param).idx)) + param_idx[n] = np.arange(idxprev,idxprev + N) + idxprev += N + # Enforce normalization of the linear parameters (if needed) for the final output FitResult_param_,FitResult_paramuq_ = FitResult_param.copy(),FitResult_paramuq.copy() if normalization: @@ -526,7 +537,7 @@ def _scale(x): noiselvl = noiselvl[0] # Generate FitResult object from all the dictionaries - fitresult = FitResult({**FitResult_param_,**FitResult_paramuq_, **FitResult_dict,'penweights':penweights,'noiselvl':noiselvl}) + fitresult = FitResult({**FitResult_param_,**FitResult_paramuq_, **FitResult_dict,'penweights':penweights,'noiselvl':noiselvl,'paramlist':_paramlist, '_param_idx':param_idx}) fitresult._summary = _print_fitresults(fitresult,model) diff --git a/deerlab/fitresult.py b/deerlab/fitresult.py index b36f064c..2dd55ab1 100644 --- a/deerlab/fitresult.py +++ b/deerlab/fitresult.py @@ -87,45 +87,76 @@ def _extarct_params_from_model(self, model): if callable(model): try: modelparam = model._parameter_list('vector') + function_type=False except AttributeError: + function_type=True modelparam = inspect.getfullargspec(model).args if not hasattr(self,'param'): raise ValueError('The fit object does not contain any fitted parameters.') - # Enforce model normalization - normfactor_keys = [] - for key in model._parameter_list(): - param = getattr(model,key) - if np.all(param.linear): - if param.normalization is not None: - normfactor_key = f'{key}_scale' - normfactor_keys.append(normfactor_key) - try: - model.addnonlinear(normfactor_key,lb=-np.inf,ub=np.inf,par0=1,description=f'Normalization factor of {key}') - getattr(model,normfactor_key).freeze(1) - except KeyError: - pass + # # Enforce model normalization + # normfactor_keys = [] + # for key in modelparam: + # param = getattr(model,key) + # if np.all(param.linear): + # if param.normalization is not None: + # normfactor_key = f'{key}_scale' + # normfactor_keys.append(normfactor_key) + # try: + # model.addnonlinear(normfactor_key,lb=-np.inf,ub=np.inf,par0=1,description=f'Normalization factor of {key}') + # getattr(model,normfactor_key).freeze(1) + # except KeyError: + # pass - # Get some basic information on the parameter vector - modelparam = model._parameter_list(order='vector') - param_idx = [[] for _ in model._parameter_list('vector')] - idxprev = 0 - for islinear in [False,True]: - for n,param in enumerate(model._parameter_list('vector')): - if np.all(getattr(model,param).linear == islinear): - N = len(np.atleast_1d(getattr(model,param).idx)) - param_idx[n] = np.arange(idxprev,idxprev + N) - idxprev += N - - params = {key : fitvalue if len(fitvalue)>1 else fitvalue[0] for key,fitvalue in zip(modelparam,[self.param[idx] for idx in param_idx])} + # # Get some basic information on the parameter vector + # modelparam = model._parameter_list(order='vector') + # param_idx = [[] for _ in model._parameter_list('vector')] + # idxprev = 0 + # for islinear in [False,True]: + # for n,param in enumerate(model._parameter_list('vector')): + # if np.all(getattr(model,param).linear == islinear): + # N = len(np.atleast_1d(getattr(model,param).idx)) + # param_idx[n] = np.arange(idxprev,idxprev + N) + # idxprev += N + + fitparams = {key : fitvalue if len(fitvalue)>1 else fitvalue[0] for key, fitvalue in zip(self.paramlist,[self.param[idx] for idx in self._param_idx])} # Check that all parameters are in the fit object for param in modelparam: - if not param in params: + if not param in fitparams: raise KeyError(f'The fit object does not contain the {param} parameter.') + + + params = {param : fitparams[param] for param in modelparam} + params_idx = [self._param_idx[self.paramlist.index(param)] for param in modelparam] + + return modelparam, params, params_idx + + def _extract_params_from_function(self,function): + """ + Extracts the fitted parameters from a callable function. + + Assumes that all parameters are length 1. + + """ + # Get the parameter names from the function definition + modelparam = inspect.getfullargspec(function).args + + fitparam_idx = self._param_idx + + # Get the parameter values from the fit object + fitparams = {param : self.param[i] for i,param in enumerate(self.paramlist)} + params = {param : fitparams[param] for param in modelparam} + params_idx = [fitparam_idx[self.paramlist.index(param)] for param in modelparam] + # fitparams = {key : fitvalue if len(fitvalue)>1 else fitvalue[0] for key, fitvalue in zip(modelparam,[self.param[idx] for idx in fitparam_idx])} + + + + + return modelparam, params, params_idx + - return modelparam, params, param_idx def evaluate(self, model, *constants): # ---------------------------------------------------------------------------- @@ -154,9 +185,14 @@ def evaluate(self, model, *constants): response : array_like Model response at the fitted parameter values. """ - - modelparam, params, _ = self._extarct_params_from_model(model) - parameters = {param: params[param] for param in modelparam} + try: + modelparam = model._parameter_list('vector') + modelparam, fitparams, fitparam_idx = self._extarct_params_from_model(model) + except AttributeError: + modelparam, fitparams, fitparam_idx = self._extract_params_from_function(model) + + + parameters = {param: fitparams[param] for param in modelparam} # Evaluate the input model response = model(*constants,**parameters) @@ -195,11 +231,16 @@ def propagate(self, model, *constants, lb=None, ub=None): Uncertainty quantification of the model's response. """ - modelparam,_, param_idx = self._extarct_params_from_model(model) - # Determine the indices of the subset of parameters the model depends on - subset = [param_idx[np.where(np.asarray(modelparam)==param)[0][0]] for param in modelparam] + try: + modelparam = model._parameter_list('vector') + modelparam, fitparams, fitparam_idx = self._extarct_params_from_model(model) + + except AttributeError: + modelparam, fitparams, fitparam_idx = self._extract_params_from_function(model) + + # Propagate the uncertainty from that subset to the model - modeluq = self.paramUncert.propagate(lambda param: model(*constants,*[param[s] for s in subset]),lb,ub) + modeluq = self.paramUncert.propagate(lambda param: model(*constants,*[param[s] for s in fitparam_idx]),lb,ub) return modeluq diff --git a/examples/advanced/ex_comparing_uncertainties.py b/examples/advanced/ex_comparing_uncertainties.py index 0d0ffd8e..2a34986f 100644 --- a/examples/advanced/ex_comparing_uncertainties.py +++ b/examples/advanced/ex_comparing_uncertainties.py @@ -27,7 +27,7 @@ # Experimental parameters tau1 = 0.3 # First inter-pulse delay, μs tau2 = 4.0 # Second inter-pulse delay, μs -deadtime = 0.1 # Acquisition deadtime, μs +tmin = 0.1 # Start time, μs # Load the experimental data t,Vexp = dl.deerload(path + file) @@ -35,7 +35,8 @@ # Pre-processing Vexp = dl.correctphase(Vexp) # Phase correction Vexp = Vexp/np.max(Vexp) # Rescaling (aesthetic) -t = t + deadtime # Account for deadtime +t = t - t[0] # Account for zerotime +t = t + tmin # Distance vector r = np.arange(2,6,0.05) # nm diff --git a/examples/advanced/ex_dipolarpathways_selection.py b/examples/advanced/ex_dipolarpathways_selection.py index 95b42182..57fb0e18 100644 --- a/examples/advanced/ex_dipolarpathways_selection.py +++ b/examples/advanced/ex_dipolarpathways_selection.py @@ -25,7 +25,7 @@ # Experimental parameters tau1 = 0.5 # First inter-pulse delay, μs tau2 = 3.5 # Second inter-pulse delay, μs -deadtime = 0.1 # Acquisition deadtime, μs +tmin = 0.1 # Start time, μs # Load the experimental data t,Vexp = dl.deerload(path + file) @@ -33,8 +33,8 @@ # Pre-processing Vexp = dl.correctphase(Vexp) # Phase correction Vexp = Vexp/np.max(Vexp) # Rescaling (aesthetic) -t = t + deadtime # Account for deadtime - +t = t - t[0] # Account for zerotime +t = t + tmin # Construct the distance vector r = np.arange(2,5,0.05) diff --git a/examples/advanced/ex_extracting_gauss_constraints.py b/examples/advanced/ex_extracting_gauss_constraints.py index 87d94a5a..759a82ac 100644 --- a/examples/advanced/ex_extracting_gauss_constraints.py +++ b/examples/advanced/ex_extracting_gauss_constraints.py @@ -23,7 +23,7 @@ # Experimental parameters tau1 = 0.3 # First inter-pulse delay, μs tau2 = 4.0 # Second inter-pulse delay, μs -deadtime = 0.1 # Acquisition deadtime, μs +tmin = 0.1 # Start time, μs # Load the experimental data t,Vexp = dl.deerload(path + file) @@ -31,7 +31,8 @@ # Pre-processing Vexp = dl.correctphase(Vexp) # Phase correction Vexp = Vexp/np.max(Vexp) # Rescaling (aesthetic) -t = t + deadtime # Account for deadtime +t = t - t[0] # Account for zerotime +t = t + tmin # Construct the dipolar signal model r = np.arange(2,6,0.02) diff --git a/examples/advanced/ex_forcefield_fit.py b/examples/advanced/ex_forcefield_fit.py index eae7bfb6..c65d4f61 100644 --- a/examples/advanced/ex_forcefield_fit.py +++ b/examples/advanced/ex_forcefield_fit.py @@ -48,7 +48,7 @@ def forcefield_P(c0,c1,c2,c3): # Experimental parameters tau1 = 0.3 # First inter-pulse delay, μs tau2 = 5.0 # Second inter-pulse delay, μs -deadtime = 0.1 # Acquisition deadtime, μs +tmin = 0.1 # Start time, μs # Load the experimental data t,Vexp = dl.deerload(path + file) @@ -56,7 +56,8 @@ def forcefield_P(c0,c1,c2,c3): # Pre-processing Vexp = dl.correctphase(Vexp) # Phase correction Vexp = Vexp/np.max(Vexp) # Rescaling (aesthetic) -t = t + deadtime # Account for deadtime +t = t - t[0] # Account for zerotime +t = t + tmin # Construct the energy and distance distribution models forcefield_energymodel = dl.Model(forcefield_energy) diff --git a/examples/advanced/ex_global_twostates_parametric.py b/examples/advanced/ex_global_twostates_parametric.py index eb0b2726..2f3bf934 100644 --- a/examples/advanced/ex_global_twostates_parametric.py +++ b/examples/advanced/ex_global_twostates_parametric.py @@ -30,7 +30,7 @@ # Experimental parameters tau1 = 0.4 # First inter-pulse delay, μs tau2 = 4.5 # Second inter-pulse delay, μs -deadtime = 0.2 # Acquisition deadtime, μs +tmin = 0.2 # Start time, μs Vmodels, ts, Vexps = [], [], [] for file in files: @@ -39,10 +39,11 @@ t, Vexp = dl.deerload(path + file) # Pre-processing - Vexp = dl.correctphase(Vexp) # Phase correction - Vexp = Vexp / np.max(Vexp) # Rescaling (aesthetic) - t = t + deadtime # Account for deadtime - + Vexp = dl.correctphase(Vexp) # Phase correction + Vexp = Vexp / np.max(Vexp) # Rescaling (aesthetic) + t = t - t[0] # Account for zerotime + t = t + tmin + # Put the datasets into lists ts.append(t) Vexps.append(Vexp) diff --git a/examples/advanced/ex_identifiability_analysis.py b/examples/advanced/ex_identifiability_analysis.py index c874b47f..6ef0fff8 100644 --- a/examples/advanced/ex_identifiability_analysis.py +++ b/examples/advanced/ex_identifiability_analysis.py @@ -23,7 +23,7 @@ # Experimental parameters tau1 = 0.3 # First inter-pulse delay, μs tau2 = 4.0 # Second inter-pulse delay, μs -deadtime = 0.1 # Acquisition deadtime, μs +tmin = 0.1 # Start time, μs # Load the experimental data t,Vexp = dl.deerload(path + file) @@ -31,7 +31,8 @@ # Pre-processing Vexp = dl.correctphase(Vexp) # Phase correction Vexp = Vexp/np.max(Vexp) # Rescaling (aesthetic) -t = t + deadtime # Account for deadtime +t = t - t[0] # Account for zerotime +t = t + tmin # Truncate the signal Vexp_truncated = Vexp[t<=2] diff --git a/examples/advanced/ex_long_threespin_analysis.py b/examples/advanced/ex_long_threespin_analysis.py index 7bc8bf6c..cdbde826 100644 --- a/examples/advanced/ex_long_threespin_analysis.py +++ b/examples/advanced/ex_long_threespin_analysis.py @@ -20,7 +20,7 @@ files = [f'../data/triradical_protein_deer_{dB}dB.DTA' for dB in [0,6,9]] # Experiment information -t0 = 0.280 # Acquisition deadtime, μs +t0 = 0.280 # Start time, μs tau1 = 0.40 # First interpulse delay, μs tau2 = 9.00 # Second interpulse delay, μs @@ -34,7 +34,7 @@ t,Vexp, descriptor = dl.deerload(file,full_output=True) t = t[:-80] Vexp = Vexp[:-80] - # Adjust the start time + # Adjust the Start time t = t - t[0] + t0 # Pre-processing diff --git a/examples/advanced/ex_multigauss_fitting_4pdeer.py b/examples/advanced/ex_multigauss_fitting_4pdeer.py index 073ea087..e2a5cf87 100644 --- a/examples/advanced/ex_multigauss_fitting_4pdeer.py +++ b/examples/advanced/ex_multigauss_fitting_4pdeer.py @@ -24,15 +24,16 @@ # Experimental parameters tau1 = 0.3 # First inter-pulse delay, μs tau2 = 4.0 # Second inter-pulse delay, μs -deadtime = 0.1 # Acquisition deadtime, μ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 + deadtime # Account for deadtime +Vexp = dl.correctphase(Vexp) # Phase correction +Vexp = Vexp / np.max(Vexp) # Rescaling (aesthetic) +t = t - t[0] # Account for zerotime +t = t + tmin # Maximal number of Gaussians in the models Nmax = 5 diff --git a/examples/advanced/ex_profileanalysis.py b/examples/advanced/ex_profileanalysis.py index b0469c4f..82d5fc42 100644 --- a/examples/advanced/ex_profileanalysis.py +++ b/examples/advanced/ex_profileanalysis.py @@ -19,7 +19,7 @@ # Experimental parameters tau1 = 0.3 # First inter-pulse delay, μs tau2 = 4.0 # Second inter-pulse delay, μs -deadtime = 0.1 # Acquisition deadtime, μs +tmin = 0.1 # Start time, μs # Load the experimental data t,Vexp = dl.deerload(path + file) @@ -27,7 +27,8 @@ # Pre-processing Vexp = dl.correctphase(Vexp) # Phase correction Vexp = Vexp/np.max(Vexp) # Rescaling (aesthetic) -t = t + deadtime # Account for deadtime +t = t - t[0] # Account for zerotime +t = t + tmin # Distance vector r = np.arange(2,5,0.05) # nm diff --git a/examples/basic/ex_bootstrapping.py b/examples/basic/ex_bootstrapping.py index ad921fc2..de9a558f 100644 --- a/examples/basic/ex_bootstrapping.py +++ b/examples/basic/ex_bootstrapping.py @@ -30,7 +30,7 @@ # Experimental parameters tau1 = 0.3 # First inter-pulse delay, μs tau2 = 4.0 # Second inter-pulse delay, μs -deadtime = 0.1 # Acquisition deadtime, μs +tmin = 0.1 # Start time, μs # Load the experimental data t,Vexp = dl.deerload(path + file) @@ -38,7 +38,8 @@ # Pre-processing Vexp = dl.correctphase(Vexp) # Phase correction Vexp = Vexp/np.max(Vexp) # Rescaling (aesthetic) -t = t + deadtime # Account for deadtime +t = t - t[0] # Account for zerotime +t = t + tmin # Distance vector r = np.linspace(2,5,100) # nm diff --git a/examples/basic/ex_fitting_4pdeer.py b/examples/basic/ex_fitting_4pdeer.py index fb2791d6..38368f5e 100644 --- a/examples/basic/ex_fitting_4pdeer.py +++ b/examples/basic/ex_fitting_4pdeer.py @@ -21,7 +21,7 @@ # Experimental parameters tau1 = 0.3 # First inter-pulse delay, μs tau2 = 4.0 # Second inter-pulse delay, μs -deadtime = 0.1 # Acquisition deadtime, μs +tmin = 0.1 # Start time, μs # Load the experimental data t,Vexp = dl.deerload(path + file) @@ -29,8 +29,8 @@ # Pre-processing Vexp = dl.correctphase(Vexp) # Phase correction Vexp = Vexp/np.max(Vexp) # Rescaling (aesthetic) -t = t + deadtime # Account for deadtime - +t = t - t[0] # Account for zerotime +t = t + tmin # Distance vector r = np.arange(2.5,5,0.01) # nm diff --git a/examples/basic/ex_fitting_4pdeer_compactness.py b/examples/basic/ex_fitting_4pdeer_compactness.py index 488452c5..b94f782b 100644 --- a/examples/basic/ex_fitting_4pdeer_compactness.py +++ b/examples/basic/ex_fitting_4pdeer_compactness.py @@ -23,7 +23,7 @@ # Experimental parameters tau1 = 0.3 # First inter-pulse delay, μs tau2 = 4.0 # Second inter-pulse delay, μs -deadtime = 0.1 # Acquisition deadtime, μs +tmin = 0.1 # Start time, μs # Load the experimental data t,Vexp = dl.deerload(path + file) @@ -31,7 +31,8 @@ # Pre-processing Vexp = dl.correctphase(Vexp) # Phase correction Vexp = Vexp/np.max(Vexp) # Rescaling (aesthetic) -t = t + deadtime # Account for deadtime +t = t - t[0] # Account for zerotime +t = t + tmin # Distance vector r = np.arange(2,5,0.025) # nm diff --git a/examples/basic/ex_fitting_4pdeer_gauss.py b/examples/basic/ex_fitting_4pdeer_gauss.py index 67da3c42..225827c2 100644 --- a/examples/basic/ex_fitting_4pdeer_gauss.py +++ b/examples/basic/ex_fitting_4pdeer_gauss.py @@ -21,7 +21,7 @@ # Experimental parameters tau1 = 0.3 # First inter-pulse delay, μs tau2 = 4.0 # Second inter-pulse delay, μs -deadtime = 0.1 # Acquisition deadtime, μs +tmin = 0.1 # Start time, μs # Load the experimental data t,Vexp = dl.deerload(path + file) @@ -29,7 +29,8 @@ # Pre-processing Vexp = dl.correctphase(Vexp) # Phase correction Vexp = Vexp/np.max(Vexp) # Rescaling (aesthetic) -t = t + deadtime # Account for deadtime +t = t - t[0] # Account for zerotime +t = t + tmin # Distance vector r = np.arange(1.5,6,0.01) # nm diff --git a/examples/basic/ex_fitting_5pdeer.py b/examples/basic/ex_fitting_5pdeer.py index 54a36e6b..83d7519d 100644 --- a/examples/basic/ex_fitting_5pdeer.py +++ b/examples/basic/ex_fitting_5pdeer.py @@ -29,16 +29,17 @@ file = 'example_5pdeer_1.DTA' # Experimental parameters (reversed 5pDEER) -tau1 = 3.9 # First inter-pulse delay, μs -tau2 = 3.7 # Second inter-pulse delay, μs -tau3 = 0.5 # Third inter-pulse delay, μs -deadtime = 0.3 # Acquisition deadtime, μs +tau1 = 3.9 # First inter-pulse delay, μs +tau2 = 3.7 # Second inter-pulse delay, μs +tau3 = 0.5 # Third inter-pulse delay, μs +tmin = 0.3 # Start time, μs # Load the experimental data t,Vexp = dl.deerload(path + file) Vexp = dl.correctphase(Vexp) # Phase correction -Vexp = Vexp/np.max(Vexp) # Rescaling (aesthetic) -t = t + deadtime # Account for deadtime +Vexp = Vexp/np.max(Vexp) # Rescaling (aesthetic) +t = t - t[0] # Account for zerotime +t = t + tmin # Distance vector r = np.arange(3,5,0.025) # nm diff --git a/examples/basic/ex_fitting_ridme.py b/examples/basic/ex_fitting_ridme.py new file mode 100644 index 00000000..f58ccfa6 --- /dev/null +++ b/examples/basic/ex_fitting_ridme.py @@ -0,0 +1,86 @@ +# %% [markdown] +""" +Basic analysis of a RIDME signal +------------------------------------------------------------------------- + +Fit a simple RIDME signal with a model with a non-parametric +distribution and a stretched exponetial background, using Tikhonov regularization. +""" + +import numpy as np +import matplotlib.pyplot as plt +import deerlab as dl + + +# %% + +# File location +path = '../data/' +file = 'example_ridme_1.DTA' + +# Experimental parameters +tau1 = 0.4 # First inter-pulse delay, μs +tau2 = 4.2 # Second inter-pulse delay, μs +tmin = 0.28 # 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.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) + +# Fit the model to the data +results = dl.fit(Vmodel,Vexp) + +# Print results summary +print(results) + +#%% + +# Extract fitted dipolar signal +Vfit = results.model + +# Extract fitted distance distribution +Pfit = results.P +Pci95 = results.PUncert.ci(95) +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) + +plt.figure(figsize=[6,7]) +violet = '#4550e6' +plt.subplot(211) +# Plot experimental and fitted data +plt.plot(t,Vexp,'.',color='grey',label='Data') +plt.plot(t,Vfit,linewidth=3,color=violet,label='Fit') +plt.plot(t,Bfit,'--',linewidth=3,color=violet,label='Unmodulated contribution') +plt.fill_between(t,Bci[:,0],Bci[:,1],color=violet,alpha=0.3) +plt.legend(frameon=False,loc='best') +plt.xlabel('Time $t$ (μs)') +plt.ylabel('$V(t)$ (arb.u.)') +# Plot the distance distribution +plt.subplot(212) +plt.plot(r,Pfit,color=violet,linewidth=3,label='Fit') +plt.fill_between(r,Pci95[:,0],Pci95[:,1],alpha=0.3,color=violet,label='95%-Conf. Inter.',linewidth=0) +plt.fill_between(r,Pci50[:,0],Pci50[:,1],alpha=0.5,color=violet,label='50%-Conf. Inter.',linewidth=0) +plt.legend(frameon=False,loc='best') +plt.autoscale(enable=True, axis='both', tight=True) +plt.xlabel('Distance $r$ (nm)') +plt.ylabel('$P(r)$ (nm$^{-1}$)') +plt.tight_layout() +plt.show() + +# %% diff --git a/examples/basic/ex_restraints_4pdeer.py b/examples/basic/ex_restraints_4pdeer.py index 7696948c..d58675f3 100644 --- a/examples/basic/ex_restraints_4pdeer.py +++ b/examples/basic/ex_restraints_4pdeer.py @@ -28,7 +28,7 @@ # Experimental parameters tau1 = 0.3 # First inter-pulse delay, μs tau2 = 4.0 # Second inter-pulse delay, μs -deadtime = 0.1 # Acquisition deadtime, μs +tmin = 0.1 # Start time, μs # Load the experimental data t,Vexp = dl.deerload(path + file) @@ -36,7 +36,7 @@ # Pre-processing Vexp = dl.correctphase(Vexp) # Phase correction Vexp = Vexp/np.max(Vexp) # Rescaling (aesthetic) -t = t + deadtime # Account for deadtime +t = t + tmin # Account for zerotime # Distance vector r = np.arange(2,6,0.05) # nm diff --git a/examples/data/example_ridme_1.DSC b/examples/data/example_ridme_1.DSC new file mode 100644 index 00000000..5f5002e4 --- /dev/null +++ b/examples/data/example_ridme_1.DSC @@ -0,0 +1,24 @@ +* Exported from Python using EasySpin, 2023-08-11 12:49:46 +#DESC 1.2 * DESCRIPTOR INFORMATION *********************** +DSRC MAN +BSEQ LIT +IKKF CPLX +IRFMT D +IIFMT D +XTYP IDX +YTYP NODATA +ZTYP NODATA +XPTS 540 +XMIN 280.0 +XWID 4312.000000000004 +* +************************************************************ +* +#SPL 1.2 * STANDARD PARAMETER LAYER +OPER +DATE 11/08/23 +TIME 12:49:46 +CMNT +SAMP +SFOR +MWFQ 3.405e+10 diff --git a/examples/data/example_ridme_1.DTA b/examples/data/example_ridme_1.DTA new file mode 100644 index 00000000..b378fcc8 Binary files /dev/null and b/examples/data/example_ridme_1.DTA differ diff --git a/examples/intermediate/ex_compactness_with_without.py b/examples/intermediate/ex_compactness_with_without.py index a8398878..0b836934 100644 --- a/examples/intermediate/ex_compactness_with_without.py +++ b/examples/intermediate/ex_compactness_with_without.py @@ -21,7 +21,7 @@ # Experimental parameters tau1 = 0.3 # First inter-pulse delay, μs tau2 = 4.0 # Second inter-pulse delay, μs -deadtime = 0.1 # Acquisition deadtime, μs +tmin = 0.1 # Start time, μs # Load the experimental data t,Vexp = dl.deerload(path + file) @@ -29,7 +29,8 @@ # Pre-processing Vexp = dl.correctphase(Vexp) # Phase correction Vexp = Vexp/np.max(Vexp) # Rescaling (aesthetic) -t = t + deadtime # Account for deadtime +t = t - t[0] # Account for zerotime +t = t + tmin # Artifically truncate the signal Vexp = Vexp[t<=3] diff --git a/examples/intermediate/ex_crossing_echoes_masking.py b/examples/intermediate/ex_crossing_echoes_masking.py index db2cba66..dd60009e 100644 --- a/examples/intermediate/ex_crossing_echoes_masking.py +++ b/examples/intermediate/ex_crossing_echoes_masking.py @@ -32,10 +32,10 @@ # Experimental parameters tau1 = 0.5 # First inter-pulse time delay, μs tau2 = 4.5 # Second inter-pulse time delay, μs -t0 = 0.3 # Acquisition deadtime, μs +tmin = 0.3 # Start time, μs -# Adjust for the deadtime -t = t - t[0] + t0 +t = t - t[0] # Account for zerotime +t = t + tmin # Plot the real part of the raw data plt.figure(figsize=[6,4]) diff --git a/examples/intermediate/ex_fitting_4pdeer_pathways.py b/examples/intermediate/ex_fitting_4pdeer_pathways.py index 59d5d993..cef4d058 100644 --- a/examples/intermediate/ex_fitting_4pdeer_pathways.py +++ b/examples/intermediate/ex_fitting_4pdeer_pathways.py @@ -20,7 +20,7 @@ # Experimental parameters tau1 = 0.5 # First inter-pulse delay, μs tau2 = 3.5 # Second inter-pulse delay, μs -deadtime = 0.1 # Acquisition deadtime, μs +tmin = 0.1 # Start time, μs # Load the experimental data t,Vexp = dl.deerload(path + file) @@ -28,8 +28,8 @@ # Pre-processing Vexp = dl.correctphase(Vexp) # Phase correction Vexp = Vexp/np.max(Vexp) # Rescaling (aesthetic) -t = t + deadtime # Account for deadtime - +t = t - t[0] # Account for zerotime +t = t + tmin # Distance vector r = np.arange(2.5,5.5,0.025) # nm diff --git a/examples/intermediate/ex_fitting_5pdeer_pathways.py b/examples/intermediate/ex_fitting_5pdeer_pathways.py index 2d11b555..370a7d93 100644 --- a/examples/intermediate/ex_fitting_5pdeer_pathways.py +++ b/examples/intermediate/ex_fitting_5pdeer_pathways.py @@ -18,17 +18,18 @@ file = 'example_5pdeer_3.DTA' # Experimental parameters (reversed 5pDEER) -tau1 = 3.7 # First inter-pulse delay, μs -tau2 = 3.5 # Second inter-pulse delay, μs -tau3 = 0.3 # Third inter-pulse delay, μs -deadtime = 0.1 # Acquisition deadtime, μs +tau1 = 3.7 # First inter-pulse delay, μs +tau2 = 3.5 # Second inter-pulse delay, μs +tau3 = 0.3 # Third inter-pulse delay, μs +tmin = 0.1 # Start time, μs # Load the experimental data t,Vexp = dl.deerload(path + file) Vexp = dl.correctphase(Vexp) # Phase correction Vexp = Vexp/np.max(Vexp) # Rescaling (aesthetic) -t = t + deadtime # Account for deadtime - +t = t - t[0] # Account for zerotime +t = t + tmin + # Distance vector r = np.arange(2.5,5.5,0.025) # nm diff --git a/examples/intermediate/ex_fitting_sparse_4pdeer.py b/examples/intermediate/ex_fitting_sparse_4pdeer.py index 47b8560b..68adaaa0 100644 --- a/examples/intermediate/ex_fitting_sparse_4pdeer.py +++ b/examples/intermediate/ex_fitting_sparse_4pdeer.py @@ -24,7 +24,7 @@ # Experimental parameters tau1 = 0.400 # First inter-pulse delay, μs tau2 = 8.000 # Second inter-pulse delay, μs -deadtime = 0.482 # Acquisition deadtime, μs +tmin = 0.482 # Start time, μs # Load the experimental data and the grid of recorded timings (this depends on how the data were acquired) _,Vexp = dl.deerload(path + datafile) @@ -34,7 +34,8 @@ # Pre-processing Vexp = dl.correctphase(Vexp) # Phase correction Vexp = Vexp/np.max(Vexp) # Rescaling (aesthetic) -t = t + deadtime # Account for deadtime +t = t - t[0] # Account for zerotime +t = t + tmin # Distance vector r = np.arange(2,10,0.05) # nm diff --git a/examples/intermediate/ex_global_different_deer.py b/examples/intermediate/ex_global_different_deer.py index 7fa0115d..ab7f0ab0 100644 --- a/examples/intermediate/ex_global_different_deer.py +++ b/examples/intermediate/ex_global_different_deer.py @@ -21,23 +21,27 @@ # Experimental parameters (4pDEER) tau1_4p = 0.5 # First inter-pulse delay, μs tau2_4p = 3.5 # Second inter-pulse delay, μs -deadtime_4p = 0.1 # Acquisition deadtime, μs +tmin_4p = 0.1 # Start time, μs # Experimental parameters (reversed 5pDEER) tau1_5p = 2.9 # First inter-pulse delay, μs tau2_5p = 3.3 # Second inter-pulse delay, μs tau3_5p = 0.3 # Third inter-pulse delay, μs -deadtime_5p = 0.1 # Acquisition deadtime, μs +tmin_5p = 0.1 # Start time, μs # Load the experimental data (4pDEER) t4p,V4p = dl.deerload(path + file4p) V4p = dl.correctphase(V4p) # Phase correction V4p = V4p/np.max(V4p) # Rescaling (aesthetic) -t4p = t4p + deadtime_4p # Account for deadtime +t4p = t4p - t4p[0] # Account for zerotime +t4p = t4p + tmin_4p + # Load the experimental data (reversed 5pDEER) t5p,V5p = dl.deerload(path + file5p) V5p = dl.correctphase(V5p) # Phase correction V5p = V5p/np.max(V5p) # Rescaling (aesthetic) -t5p = t5p + deadtime_5p # Account for deadtime +t5p = t5p - t5p[0] # Account for zerotime +t5p = t5p + tmin_5p + # Run fit r = np.arange(2.5,6,0.05) diff --git a/examples/intermediate/ex_global_fitting_4pdeer.py b/examples/intermediate/ex_global_fitting_4pdeer.py index 9cc0cea9..6e6e1c9b 100644 --- a/examples/intermediate/ex_global_fitting_4pdeer.py +++ b/examples/intermediate/ex_global_fitting_4pdeer.py @@ -23,10 +23,10 @@ # Experimental parameters tau1s = [0.3, 0.5] # First inter-pulse delay, μs tau2s = [2.0, 4.0] # Second inter-pulse delay, μs -deadtimes = [0.1, 0.3] # Acquisition deadtime, μs +tmins = [0.1, 0.3] # Start time, μs Vmodels,ts,Vs = [],[],[] -for file, tau1, tau2, deadtime in zip(files, tau1s, tau2s, deadtimes): +for file, tau1, tau2, tmin in zip(files, tau1s, tau2s, tmins): # Load the experimental data t,Vexp = dl.deerload(path + file) @@ -34,7 +34,8 @@ # Pre-processing Vexp = dl.correctphase(Vexp) # Phase correction Vexp = Vexp/np.max(Vexp) # Rescaling (aesthetic) - t = t + deadtime # Account for deadtime + t = t - t[0] # Account for zerotime + t = t + tmin # Distance vector r = np.arange(1.5,7,0.05) # nm diff --git a/examples/intermediate/ex_multipathway_validation.py b/examples/intermediate/ex_multipathway_validation.py index 348e04e4..cf3deaac 100644 --- a/examples/intermediate/ex_multipathway_validation.py +++ b/examples/intermediate/ex_multipathway_validation.py @@ -27,7 +27,7 @@ file = "../data/experimental_mbp_protein_4pdeer.DTA" # Experiment information -t0 = 0.040 +tmin = 0.040 tau1 = 0.4 tau2 = 3.0 @@ -37,7 +37,7 @@ Vexp = Vexp[:-2] Vexp = dl.correctphase(Vexp) Vexp = Vexp/max(Vexp) -t = t- t[0] + t0 +t = t- t[0] + tmin # Define the distance vector r = np.arange(3,4.5,0.05) diff --git a/examples/intermediate/ex_selregparam.py b/examples/intermediate/ex_selregparam.py index e4f99958..417c284d 100644 --- a/examples/intermediate/ex_selregparam.py +++ b/examples/intermediate/ex_selregparam.py @@ -21,7 +21,7 @@ # Experimental parameters tau1 = 0.3 # First inter-pulse delay, μs tau2 = 4.0 # Second inter-pulse delay, μs -deadtime = 0.1 # Acquisition deadtime, μs +tmin = 0.1 # Start time, μs # Load the experimental data t,Vexp = dl.deerload(path + file) @@ -29,7 +29,8 @@ # Pre-processing Vexp = dl.correctphase(Vexp) # Phase correction Vexp = Vexp/np.max(Vexp) # Rescaling (aesthetic) -t = t + deadtime # Account for deadtime +t = t - t[0] # Account for zerotime +t = t + tmin # Distance vector r = np.linspace(1.5,7,50) # nm diff --git a/test/test_model_class.py b/test/test_model_class.py index d8802fba..b8f6b864 100644 --- a/test/test_model_class.py +++ b/test/test_model_class.py @@ -568,6 +568,31 @@ def test_fit_propagate_through_model(mock_data,mock_x,model_type,method): assert np.less_equal(ci_lower,ci_upper).all() # ================================================================ + +@pytest.mark.parametrize('method', ['bootstrap','moment']) +@pytest.mark.parametrize('model_type', ['parametric']) +def test_fit_propagate_through_function(mock_data,mock_x,model_type,method): + + model = _generate_model(model_type, fixed_axis=False) + + if method=='bootstrap': + results = fit(model,mock_data,mock_x, bootstrap=3) + else: + results = fit(model,mock_data,mock_x) + + if model_type=='parametric': + fun = lambda mean1,mean2,std1,std2,amp1,amp2: bigauss(mock_x,mean1,mean2,std1,std2,amp1,amp2) + response_ci = results.propagate(fun, lb=np.zeros_like(mock_x)).ci(95) + + ci_lower = response_ci[:,0] + ci_upper = response_ci[:,1] + + assert np.less_equal(ci_lower,ci_upper).all() + + + +# ================================================================ + def gauss_multiaxis(axis1,axis2,mean,std): return np.exp(-(axis1-mean)**2/std**2/2)