diff --git a/.github/ISSUE_TEMPLATE/---bug-report.md b/.github/ISSUE_TEMPLATE/---bug-report.md deleted file mode 100644 index f4f998475..000000000 --- a/.github/ISSUE_TEMPLATE/---bug-report.md +++ /dev/null @@ -1,34 +0,0 @@ ---- -name: "\U0001F41B Bug report" -about: Something is no working and needs a fix. -title: "[BUG] A short description of the bug" -labels: bug -assignees: '' - ---- - -### Describe the bug -A clear and concise description of what the bug is. - -### Environment -**DeerLab Version:** X.X.X -**Python Version:** X.X.X -**OS:** Windows, Max, Linux? - -### How to reproduce it? -Steps to reproduce the behavior: -1. Go to '...' -2. Click on '....' -3. Scroll down to '....' -4. See error - -### Expected behavior -A clear and concise description of what you expected to happen. - -### Code Snippet -```` -Copy-paste your code here -```` - -### Files -If there are any files required to run your code please upload them. diff --git a/.github/workflows/deploy_ghpages.yml b/.github/workflows/deploy_ghpages.yml index 57cfa2bea..e9dbd05f2 100644 --- a/.github/workflows/deploy_ghpages.yml +++ b/.github/workflows/deploy_ghpages.yml @@ -7,9 +7,6 @@ on: - master paths: - 'docsrc/**' - schedule: - # Run once a week on Sunday at 12:00 PM - - cron: '0 12 * * 0' jobs: diff --git a/CHANGELOG.md b/CHANGELOG.md index d65727874..a76fba134 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,4 +1,43 @@ +------------------------------- + +Version 0.11.0 - September 2020 +----------------------------- + +#### Overall changes + +* All Gauss models (``dd_gauss``,etc.) now use the standard deviation ``sigma`` instead of the FWHM as the width parameter for consistency with other method such as Rice distributions. (#19) + +* All hard-wired random seeds have been removed. + +* A new method ``plot()`` has been added to the ``FitResult`` class returned by all fit functions. This will create a basic plot of the fit results. (#7) + +#### Specific changes +* ``snlls``: + - Renamed option ``penalty`` as ``reg`` and improved its interface (#13). + - The regularization parameter of the optimal solution is returned now (#20). + +* ``whitegaussnoise``: Added a ``seed`` option to select static noise realizations. + +* ``correctzerotime``: + - Fixed bug when zero-time is at start/end of array (#24). + - Function no longer rescales the experimental data passed on to the function. + +* ``fitsignal``: + - The regularization parameter of the optimal solution is returned now (#20). + - Bug fixed when fitting dipolar evolution functions (no background and no experiment models) with a parametric distance distribution. + +* ``fitmultimodel``: Start points are now spread over constrained parameter space grid instead of being randomble initiated (#22). + +* ``deerload``: + - Now returns the time axis in microseconds instead of nanoseconds (#21). + - The bug appearing when loading certain BES3T files has been fixed (#14). + +* ``fitregmodel``: Now returns the fitted dipolar signal in the ``FitResult`` output + +* ``correctscale``: The parameter fit ranges have been adjusted. + + ------------------------------- Version 0.10.0 - August 2020 @@ -6,31 +45,31 @@ Version 0.10.0 - August 2020 As of this version, DeerLab is written in Python in contrast to older versions based on MATLAB. -Deprecated functions - The following functions have been deprecated due to limited usability or due to functionality overlap with other DeerLab functions: ``aptkernel``, ``backgroundstart``, ``fitbackground``, ``paramodel``, and ``time2freq``. +#### Deprecated functions +The following functions have been deprecated due to limited usability or due to functionality overlap with other DeerLab functions: ``aptkernel``, ``backgroundstart``, ``fitbackground``, ``paramodel``, and ``time2freq``. -Overall changes - * All fit functions now return a single ``FitResult`` output which will contain all results. +#### Overall changes +* All fit functions now return a single ``FitResult`` output which will contain all results. - * All functions are now compatible with non-uniformly increasing distance axes. +* All functions are now compatible with non-uniformly increasing distance axes. - * All fit functions are completely agnostic with respect of the abolute values of the signal amplitude. This is automatically fitted by all function and return as part of the results. +* All fit functions are completely agnostic with respect of the abolute values of the signal amplitude. This is automatically fitted by all function and return as part of the results. - * Uncertainty quantification for all fit functions is returned as a ``UncertQuant`` object from which confidence intervals, parameter distributions, etc. can be generated generalizing the uncertainty interface for all DeerLab. Uncertainty can now be propagated to arbitrary functions. +* Uncertainty quantification for all fit functions is returned as a ``UncertQuant`` object from which confidence intervals, parameter distributions, etc. can be generated generalizing the uncertainty interface for all DeerLab. Uncertainty can now be propagated to arbitrary functions. -Specific changes - * ``fitparamodel``: the functionality has been streamlined. Function now fits arbitrary parametric models using non-linear leas-squares without consideration of whether it is a time-domain or distance-domain model. The models no longer need to take two inputs (axis+parameters) and now only tk the parameters as input. +#### Specific changes +* ``fitparamodel``: the functionality has been streamlined. Function now fits arbitrary parametric models using non-linear leas-squares without consideration of whether it is a time-domain or distance-domain model. The models no longer need to take two inputs (axis+parameters) and now only tk the parameters as input. - * ``fitregmodel``: goodness-of-fit estimators are now computed using the proper estimation the degrees of freedom. +* ``fitregmodel``: goodness-of-fit estimators are now computed using the proper estimation the degrees of freedom. - * ``fitmultimodel``: added internal measures to avoid situations where one or several components are suppressed by fitting zero-amplitudes making the method more stable. +* ``fitmultimodel``: added internal measures to avoid situations where one or several components are suppressed by fitting zero-amplitudes making the method more stable. - * ``uqst``: the uncertainty distributions of the parameters are now returned as properly normalized probability density functions. +* ``uqst``: the uncertainty distributions of the parameters are now returned as properly normalized probability density functions. - * ``fftspec``: frequency axis construction has been corrected. +* ``fftspec``: frequency axis construction has been corrected. - * ``regoperator``: now calculates the numerically exact finite-difference matrix using Fornberg's method. +* ``regoperator``: now calculates the numerically exact finite-difference matrix using Fornberg's method. - * ``correctphase``: now can handle 2D-datasets. +* ``correctphase``: now can handle 2D-datasets. ------------------------------- \ No newline at end of file diff --git a/VERSION b/VERSION index f78dc3652..5a9560e05 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v0.10.0 \ No newline at end of file +v0.11.0-dev diff --git a/deerlab/bg_models.py b/deerlab/bg_models.py index a3c45279f..c0434f7a5 100644 --- a/deerlab/bg_models.py +++ b/deerlab/bg_models.py @@ -1,3 +1,8 @@ +# bg_models.py - Background parametric models +# --------------------------------------------------------------------------- +# This file is a part of DeerLab. License is MIT (see LICENSE.md). +# Copyright(c) 2019-2020: Luis Fabregas, Stefan Stoll and other contributors. + import numpy as np import math as m import scipy as scp @@ -46,7 +51,7 @@ def bg_exp(*args): ------------------------------------------------- Parameter Units Lower Upper Start ------------------------------------------------- - Decay Rate us-1 0 200 0.35 + Decay Rate μs⁻¹ 0 200 0.35 ------------------------------------------------- @@ -76,7 +81,7 @@ def bg_exp(*args): if not args: info = dict( Parameters = ['Decay Rate'], - Units = ['us-1'], + Units = ['μs⁻¹'], Start = np.asarray([0.35]), Lower = np.asarray([0]), Upper = np.asarray([200]) @@ -115,7 +120,7 @@ def bg_hom3d(*args): ---------------------------------------------------------------------- Parameter Units Lower Upper Start ---------------------------------------------------------------------- - Concentration of pumped spins uM 0.01 5000 50 + Concentration of pumped spins μM 0.01 5000 50 ---------------------------------------------------------------------- Parameters @@ -144,7 +149,7 @@ def bg_hom3d(*args): if not args: info = dict( Parameters = ['Concentration of pumped spins'], - Units = ['uM'], + Units = ['μM'], Start = np.asarray([50]), Lower = np.asarray([0.01]), Upper = np.asarray([5000]) @@ -193,7 +198,7 @@ def bg_hom3dex(*args): ----------------------------------------------------------------------------- Parameter Units Lower Upper Start ----------------------------------------------------------------------------- - Fractal Concentration of pumped spins umol/dm^d 0.01 5000 50 + Fractal Concentration of pumped spins μmol/dmᵈ 0.01 5000 50 Exclusion distance nm 0.10 20 1 ----------------------------------------------------------------------------- @@ -223,7 +228,7 @@ def bg_hom3dex(*args): if not args: info = dict( Parameters = ['Fractal Concentration of pumped spins','Fractal dimensionality'], - Units = ['umol/dm^d',''], + Units = ['μmol/dmᵈ',''], Start = np.asarray([50, 1]), Lower = np.asarray([0.01, 0.01]), Upper = np.asarray([5000, 20]) @@ -292,7 +297,7 @@ def bg_homfractal(*args): ----------------------------------------------------------------------------- Parameter Units Lower Upper Start ----------------------------------------------------------------------------- - Fractal Concentration of pumped spins umol/dm^d 0.01 5000 50 + Fractal Concentration of pumped spins μmol/dmᵈ 0.01 5000 50 Fractal dimensionality 0 6 3 ----------------------------------------------------------------------------- @@ -322,7 +327,7 @@ def bg_homfractal(*args): if not args: info = dict( Parameters = ['Fractal Concentration of pumped spins','Fractal dimensionality'], - Units = ['umol/dm^d',''], + Units = ['μmol/dmᵈ',''], Start = np.asarray([50, 3]), Lower = np.asarray([0.01, 0+np.finfo(float).eps]), Upper = np.asarray([5000, 6-np.finfo(float).eps]) @@ -383,7 +388,7 @@ def bg_strexp(*args): ---------------------------------------------------- Parameter Units Lower Upper Start ---------------------------------------------------- - Decay Rate us-1 0 200 0.25 + Decay Rate μs⁻¹ 0 200 0.25 Stretch factor 0 6 1 ---------------------------------------------------- @@ -413,7 +418,7 @@ def bg_strexp(*args): if not args: info = dict( Parameters = ['Decay Rate','Stretch factor'], - Units = ['us-1',''], + Units = ['μs⁻¹',''], Start = np.asarray([0.25, 1]), Lower = np.asarray([0, 0]), Upper = np.asarray([200, 6]) @@ -422,7 +427,7 @@ def bg_strexp(*args): t,param,lam = _parsargs(args,npar=2) # Unpack model paramters - kappa = param[0] # decay rate, us-1 + kappa = param[0] # decay rate, µs^-1 d = param[1] # fractal dimension B = np.exp(-lam*kappa*abs(t)**d) @@ -452,9 +457,9 @@ def bg_prodstrexp(*args): ----------------------------------------------------------------- Parameter Units Lower Upper Start ----------------------------------------------------------------- - Decay Rate of 1st component us-1 0 200 0.25 + Decay Rate of 1st component μs⁻¹ 0 200 0.25 Stretch factor 1st component 0 6 1 - Decay Rate of 2nd component us-1 0 200 0.25 + Decay Rate of 2nd component μs⁻¹ 0 200 0.25 Stretch factor 2nd component 0 6 1 ----------------------------------------------------------------- @@ -484,7 +489,7 @@ def bg_prodstrexp(*args): if not args: info = dict( Parameters = ['Decay Rate of 1st component','Stretch factor of 1st component','Decay Rate of 2nd component','Stretch factor of 2nd component'], - Units = ['us-1','','us-1',''], + Units = ['µs^-1','','µs^-1',''], Start = np.asarray([0.25, 1, 0.25, 1]), Lower = np.asarray([ 0, 0, 0, 0]), Upper = np.asarray([200, 6, 200, 6]) @@ -527,10 +532,10 @@ def bg_sumstrexp(*args): ----------------------------------------------------------------- Parameter Units Lower Upper Start ----------------------------------------------------------------- - Decay Rate of 1st component us-1 0 200 0.25 + Decay Rate of 1st component μs⁻¹ 0 200 0.25 Stretch factor 1st component 0 6 1 Amplitude of 1st component 0 1 0.50 - Decay Rate of 2nd component us-1 0 200 0.25 + Decay Rate of 2nd component μs⁻¹ 0 200 0.25 Stretch factor 2nd component 0 6 1 ----------------------------------------------------------------- @@ -560,7 +565,7 @@ def bg_sumstrexp(*args): if not args: info = dict( Parameters = ['Decay Rate of 1st component','Stretch factor of 1st component','Amplitude of 1st component','Decay Rate of 2nd component','Stretch factor of 2nd component'], - Units = ['us-1','','','us-1',''], + Units = ['μs⁻¹','','','μs⁻¹',''], Start = np.asarray([0.25, 1, 0.5, 0.25, 1]), Lower = np.asarray([ 0, 0, 0, 0, 0]), Upper = np.asarray([200, 6, 1, 200, 6]) @@ -606,7 +611,7 @@ def bg_poly1(*args): Parameter Units Lower Upper Start ---------------------------------------------------------- Intercept 0 200 1 - 1st-order coefficient us^-1 -200 200 -1 + 1st-order coefficient μs⁻¹ -200 200 -1 ---------------------------------------------------------- Parameters @@ -635,7 +640,7 @@ def bg_poly1(*args): if not args: info = dict( Parameters = ['Intercept','1st-order coefficient'], - Units = ['','us-1'], + Units = ['','μs⁻¹'], Start = np.asarray([ 1, -1 ]), Lower = np.asarray([ 0, -200]), Upper = np.asarray([200, 200]) @@ -678,8 +683,8 @@ def bg_poly2(*args): Parameter Units Lower Upper Start ---------------------------------------------------------- Intercept 0 200 1 - 1st-order coefficient us^-1 -200 200 -1 - 2nd-order coefficient us^-2 -200 200 -1 + 1st-order coefficient μs⁻¹ -200 200 -1 + 2nd-order coefficient μs⁻² -200 200 -1 ---------------------------------------------------------- Parameters @@ -708,7 +713,7 @@ def bg_poly2(*args): if not args: info = dict( Parameters = ['Intercept','1st-order coefficient','2nd-order coefficient'], - Units = ['','us-1','us-2'], + Units = ['','μs⁻¹','μs⁻²'], Start = np.asarray([ 1, -1 , -1]), Lower = np.asarray([ 0, -200, -200]), Upper = np.asarray([200, 200, 200]) @@ -747,9 +752,9 @@ def bg_poly3(*args): Parameter Units Lower Upper Start ---------------------------------------------------------- Intercept 0 200 1 - 1st-order coefficient us^-1 -200 200 -1 - 2nd-order coefficient us^-2 -200 200 -1 - 3rd-order coefficient us^-3 -200 200 -1 + 1st-order coefficient μs⁻¹ -200 200 -1 + 2nd-order coefficient μs⁻² -200 200 -1 + 3rd-order coefficient μs⁻³ -200 200 -1 ---------------------------------------------------------- Parameters @@ -778,7 +783,7 @@ def bg_poly3(*args): if not args: info = dict( Parameters = ['Intercept','1st-order coefficient','2nd-order coefficient','3rd-order coefficient'], - Units = ['','us-1','us-2','us-3'], + Units = ['','μs⁻¹','μs⁻²','μs⁻³'], Start = np.asarray([ 1, -1 , -1, -1 ]), Lower = np.asarray([ 0, -200, -200, -200]), Upper = np.asarray([200, 200, 200, 200]) diff --git a/deerlab/bootan.py b/deerlab/bootan.py index 9bb2b29f6..736c51bd8 100644 --- a/deerlab/bootan.py +++ b/deerlab/bootan.py @@ -15,7 +15,7 @@ def bootan(fcn,Vexp,Vfit, samples=1000, resampling='gaussian', verbose = False): fcn : callable Function to be analyzed. Must be a callable function accepting a signal array as input and returning a tuple with all variables to be analyzed. All variables must be numerical arrays (no strings or booleans) and must preserve shape between calls. - V : array_like or list of array_like + Vexp : array_like or list of array_like Experimental dataset(s). Vfit : array or list of array_like Fit of the dataset(s). @@ -97,9 +97,6 @@ def myfcn(V): Vresample = [0]*nSignals for iSample in range(nSamples-1): - # Change the random seed each iteration, but keeping reproducibility between runs - np.random.seed(iSample) - # Inform of progress if requested if verbose: print('Bootstrapping: #{}/#{} samples finished'.format(iSample+1,nSamples), end='\r', flush=True) @@ -110,11 +107,11 @@ def myfcn(V): if iSample>0: #Determine requested re-sampling method - if resampling is 'gaussian': + if resampling == 'gaussian': # Resample from a Gaussian distribution with variance estimated from the residuals Vresample[i] = Vfit[i] + np.random.normal(0, sigma[i], len(Vfit[i])) - elif resampling is 'residual': + elif resampling == 'residual': # Resample from the residual directly Vresample[i] = Vfit[i] + residuals[i][np.random.permutation(len(Vfit[i]))] else: diff --git a/deerlab/classes.py b/deerlab/classes.py index 87509f1f5..cfb16d4d3 100644 --- a/deerlab/classes.py +++ b/deerlab/classes.py @@ -1,3 +1,12 @@ +# This file is a part of DeerLab. License is MIT (see LICENSE.md). +# Copyright(c) 2019-2020: Luis Fabregas, Stefan Stoll and other contributors. + +import numpy as np +import numdifftools as nd +from scipy.stats import norm +from scipy.signal import fftconvolve +import copy + class FitResult(dict): r""" Represents the results of a fit. @@ -21,6 +30,13 @@ class FitResult(dict): * ``stats['aicc']`` - Corrected Akaike information criterion * ``stats['bic']`` - Bayesian information criterion + Methods + ------- + plot() + Display the fit results on a Matplotlib window. The script returns a + `matplotlib.axes `_ object. + All graphical parameters can be adjusted from this object. + Notes ----- There may be additional attributes not listed above depending of the @@ -49,48 +65,13 @@ def __repr__(self): def __dir__(self): return list(self.keys()) - - - - - - - - - - - - - - - - - - - - - - - - - - -# This file is a part of DeerLab. License is MIT (see LICENSE.md). -# Copyright(c) 2019-2020: Luis Fabregas, Stefan Stoll and other contributors. - -import numpy as np -import numdifftools as nd -from scipy.stats import norm -from scipy.signal import fftconvolve -import copy - class UncertQuant: - r""" Represens the uncertainty quantification of fit results. + r""" Represents the uncertainty quantification of fit results. Attributes ---------- type : string - Uncertainty quanification approach: + Uncertainty quantification approach: * 'covariance' - Covariance-based uncertainty analysis * 'bootstrap' - Bootstrapped uncertainty analysis @@ -101,9 +82,9 @@ class UncertQuant: Median values of the uncertainty distribution of the parameters. std : ndarray Standard deviations of the uncertainty distribution of the parameters. - covmat: ndarray + covmat : ndarray Covariance matrix - nparam: int scalar + nparam : int scalar Number of parameters in the analysis. Methods @@ -125,7 +106,7 @@ def __init__(self,uqtype,data,covmat=[],lb=[],ub=[]): elif uqtype == 'bootstrap': # Scheme 2: UncertQuant('bootstrap',samples) samples = data - self.__samples = samples + self.samples = samples nParam = np.shape(samples)[1] else: @@ -194,7 +175,7 @@ def pardist(self,n): if self.type == 'bootstrap': # Get bw using silverman's rule (1D only) - samplen = self.__samples[:, n] + samplen = self.samples[:, n] sigma = np.std(samplen, ddof=1) bw = sigma * (len(samplen) * 3 / 4.0) ** (-1 / 5) @@ -292,8 +273,8 @@ def ci(self,coverage): elif self.type=='bootstrap': # Compute bootstrap-based confidence intervals # Clip possible artifacts from the percentile estimation - x[:,0] = np.minimum(self.percentile(p*100), np.amax(self.__samples)) - x[:,1] = np.maximum(self.percentile((1-p)*100), np.amin(self.__samples)) + x[:,0] = np.minimum(self.percentile(p*100), np.amax(self.samples)) + x[:,1] = np.maximum(self.percentile((1-p)*100), np.amin(self.samples)) return x diff --git a/deerlab/correctphase.py b/deerlab/correctphase.py index c6dc438ab..a86948d29 100644 --- a/deerlab/correctphase.py +++ b/deerlab/correctphase.py @@ -82,7 +82,7 @@ def correctphase(V, phase = [], imagoffset=False, full_output=False): # Fit only imaginary offset for i in range(Ntraces): par0 = 0 - fun = lambda offset: _imaginarynorm([Phase, offset],V[:,i]) + fun = lambda offset: _imaginarynorm([phase, offset],V[:,i]) ImagOffset[i] = opt.fmin(fun,par0) ImagOffset = ImagOffset*1j diff --git a/deerlab/correctscale.py b/deerlab/correctscale.py index facd62c36..669775e93 100644 --- a/deerlab/correctscale.py +++ b/deerlab/correctscale.py @@ -55,16 +55,17 @@ def correctscale(V,t,tmax=[],model='deer'): V_ = V_[idx] t_ = t[ idx] - if model is 'deer': + if model == 'deer': # DEER model with Gaussian distibution if len(V_)<5: raise KeyError('Number of points in fit range cannot be smaller than number of model parameters. Increase tmax.') r = np.linspace(0.5,20,300) fitmodel = lambda p: p[0]*(dipolarkernel(t_,r,p[1],bg_exp(t_,p[4]))@dd_gauss(r,p[[2,3]])) + # parameters: amplitude, mod depth,Gauss position, Gauss std, conc par0 = [1, 0.5, 3, 0.3, 0.2] - lb = [1e-3, 0.5, 0, 0, 0] + lb = [1e-3, 0.01, 0, 0.01, 0] ub = [10, 1, 20, 5, 100] - elif model is 'gauss': + elif model == 'gauss': # Gaussian function if len(V_)<3: raise KeyError('Number of points in fit range cannot be smaller than number of model parameters. Increase tmax.') diff --git a/deerlab/correctzerotime.py b/deerlab/correctzerotime.py index 41b0348dc..fed41508f 100644 --- a/deerlab/correctzerotime.py +++ b/deerlab/correctzerotime.py @@ -32,17 +32,18 @@ def correctzerotime(V,t): if len(t)!=len(V): raise TypeError('The dipolar signal (V) and time axis (t) must have the size.') - # Rescale to avoid overflow during integration later - V /= max(V) + # Generate finely-grained interpolated signal and time axis resolution = 10 t_ = np.linspace(min(t), max(t), (len(t)-1)*resolution+1) V_ = scp.interpolate.interp1d(t, np.real(V), kind='cubic')(t_) + # Rescale to avoid overflow during integration later + V_ /= max(V_) + # Determine location of maximum idxmax = np.argmax(V_) - - if idxmax!=1 and idxmax!=len(V): + if idxmax!=0 and idxmax!=len(V_)-1: # If maximum is not the first or last point, then do moment analysis # (i.e. minimize the magntitude of the first-moment integral) diff --git a/deerlab/dd_models.py b/deerlab/dd_models.py index 954346a9c..45fdc1cb2 100644 --- a/deerlab/dd_models.py +++ b/deerlab/dd_models.py @@ -1,8 +1,11 @@ +# dd_models.py - Distance distribtion parametric models +# --------------------------------------------------------------------------- +# This file is a part of DeerLab. License is MIT (see LICENSE.md). +# Copyright(c) 2019-2020: Luis Fabregas, Stefan Stoll and other contributors. + import math as m import numpy as np import scipy.special as spc -__docformat__ = 'reStructuredText' - def _parsargs(args,npar): #================================================================= @@ -21,14 +24,13 @@ def _normalize(r,P): return P #================================================================= -def _multigaussfun(r,r0,fwhm,a): +def _multigaussfun(r,r0,sig,a): #================================================================= "Compute a distribution with multiple Gaussians" n = len(r0) P = np.zeros_like(r) for k in range(n): - sig = fwhm[k]/2/m.sqrt(2*m.log(2)) - P += a[k]*m.sqrt(1/(2*m.pi))*1/sig*np.exp(-0.5*((r-r0[k])/sig)**2) + P += a[k]*m.sqrt(1/(2*m.pi))*1/sig[k]*np.exp(-0.5*((r-r0[k])/sig[k])**2) P = _normalize(r,P) return P #================================================================= @@ -88,28 +90,28 @@ def dd_gauss(*args): Model parameters: ------------------- - ------------------------------------------------ - Parameter Units Lower Upper Start - ------------------------------------------------ - Center nm 1 20 3.5 - FWHM nm 0.2 5 0.5 - ------------------------------------------------ + ------------------------------------------------------------- + Parameter Units Lower Upper Start + ------------------------------------------------------------- + Mean nm 1 20 3.5 + Standard deviation nm 0.05 2.5 0.2 + ------------------------------------------------------------- """ if not args: info = dict( - Parameters = ('Center','FWHM'), + Parameters = ('Mean','Standard deviation'), Units = ('nm','nm'), - Start = np.asarray([3.5, 0.5]), - Lower = np.asarray([1, 0.2]), - Upper = np.asarray([20, 5]) + Start = np.asarray([3.5, 0.2]), + Lower = np.asarray([1, 0.05]), + Upper = np.asarray([20, 2.5]) ) return info r,p = _parsargs(args,npar=2) r0 = [p[0]] - fwhm = [p[1]] + sigma = [p[1]] a = [1.0] - P = _multigaussfun(r,r0,fwhm,a) + P = _multigaussfun(r,r0,sigma,a) return P #================================================================= @@ -149,34 +151,34 @@ def dd_gauss2(*args): P : ndarray Distance distribution. - --------------------------------------------------------------- - Parameter Units Lower Upper Start - --------------------------------------------------------------- - Center of 1st Gaussian nm 1 20 2.5 - FWHM of 1st Gaussian nm 0.2 5 0.5 - Amplitude of 1sd Gaussian 0 1 0.5 - Center of 2nd Gaussian nm 1 20 3.5 - FWHM of 2nd Gaussian nm 0.2 5 0.5 - Amplitude of 2nd Gaussian 0 1 0.5 - -------------------------------------------------------------- + ------------------------------------------------------------------------- + Parameter Units Lower Upper Start + ------------------------------------------------------------------------- + Mean of 1st Gaussian nm 1 20 2.5 + Standard deviation of 1st Gaussian nm 0.05 2.5 0.2 + Amplitude of 1st Gaussian 0 1 0.5 + Mean of 2nd Gaussian nm 1 20 3.5 + Standard deviation of 2nd Gaussian nm 0.05 2.5 0.2 + Amplitude of 2nd Gaussian 0 1 0.5 + ------------------------------------------------------------------------- """ if not args: info = dict( - Parameters = ('Center of 1st Gaussian', 'FWHM of 1st Gaussian', 'Amplitude of 1st Gaussian', - 'Center of 2nd Gaussian', 'FWHM of 2nd Gaussian', 'Amplitude of 2nd Gaussian'), + Parameters = ('Mean of 1st Gaussian', 'Standard deviation of 1st Gaussian', 'Amplitude of 1st Gaussian', + 'Mean of 2nd Gaussian', 'Standard deviation of 2nd Gaussian', 'Amplitude of 2nd Gaussian'), Units = ('nm','nm','','nm','nm',''), - Start = np.asarray([2.5, 0.5, 0.5, 3.5, 0.5, 0.5]), - Lower = np.asarray([1, 0.2, 0, 1, 0.2, 0]), - Upper = np.asarray([20, 5, 1, 20, 5, 1]) + Start = np.asarray([2.5, 0.2, 0.5, 3.5, 0.2, 0.5]), + Lower = np.asarray([1, 0.05, 0, 1, 0.05, 0]), + Upper = np.asarray([20, 2.5, 1, 20, 2.5, 1]) ) return info r,p = _parsargs(args,npar=6) r0 = [p[0], p[3]] - fwhm = [p[1], p[4]] + sigma = [p[1], p[4]] a = [p[2], p[5]] - P = _multigaussfun(r,r0,fwhm,a) + P = _multigaussfun(r,r0,sigma,a) return P #================================================================= @@ -220,38 +222,38 @@ def dd_gauss3(*args): Model parameters: ------------------- - --------------------------------------------------------------- - Parameter Units Lower Upper Start - --------------------------------------------------------------- - Center of 1st Gaussian nm 1 20 2.5 - FWHM of 1st Gaussian nm 0.2 5 0.5 - Amplitude of 1sd Gaussian 0 1 0.3 - Center of 2nd Gaussian nm 1 20 3.5 - FWHM of 2nd Gaussian nm 0.2 5 0.5 - Amplitude of 2nd Gaussian 0 1 0.3 - Center of 3rd Gaussian nm 1 20 5.0 - FWHM of 3rd Gaussian nm 0.2 5 0.5 - Amplitude of 3rd Gaussian 0 1 0.3 - -------------------------------------------------------------- + ------------------------------------------------------------------------- + Parameter Units Lower Upper Start + ------------------------------------------------------------------------- + Mean of 1st Gaussian nm 1 20 2.5 + Standard deviation of 1st Gaussian nm 0.05 2.5 0.2 + Amplitude of 1sd Gaussian 0 1 0.3 + Mean of 2nd Gaussian nm 1 20 3.5 + Standard deviation of 2nd Gaussian nm 0.05 2.5 0.2 + Amplitude of 2nd Gaussian 0 1 0.3 + Mean of 3rd Gaussian nm 1 20 5.0 + Standard deviation of 3rd Gaussian nm 0.05 2.5 0.2 + Amplitude of 3rd Gaussian 0 1 0.3 + ------------------------------------------------------------------------- """ if not args: info = dict( - Parameters = ('Center of 1st Gaussian', 'FWHM of 1st Gaussian', 'Amplitude of 1st Gaussian', - 'Center of 2nd Gaussian', 'FWHM of 2nd Gaussian', 'Amplitude of 2nd Gaussian', - 'Center of 3rd Gaussian', 'FWHM of 3rd Gaussian', 'Amplitude of 3rd Gaussian'), + Parameters = ('Mean of 1st Gaussian', 'Standard deviation of 1st Gaussian', 'Amplitude of 1st Gaussian', + 'Mean of 2nd Gaussian', 'Standard deviation of 2nd Gaussian', 'Amplitude of 2nd Gaussian', + 'Mean of 3rd Gaussian', 'Standard deviation of 3rd Gaussian', 'Amplitude of 3rd Gaussian'), Units = ('nm','nm','','nm','nm','','nm','nm',''), - Start = np.asarray([2.5, 0.5, 0.3, 3.5, 0.5, 0.3, 5, 0.5, 0.3]), - Lower = np.asarray([1, 0.2, 0, 1, 0.2, 0, 1, 0.2, 0]), - Upper = np.asarray([20, 5, 1, 20, 5, 1, 20, 5, 1]) + Start = np.asarray([2.5, 0.2, 0.3, 3.5, 0.2, 0.3, 5, 0.2, 0.3]), + Lower = np.asarray([1, 0.05, 0, 1, 0.05, 0, 1, 0.05, 0]), + Upper = np.asarray([20, 2.5, 1, 20, 2.5, 1, 20, 2.5, 1]) ) return info r,p = _parsargs(args,npar=9) r0 = [p[0], p[3], p[6]] - fwhm = [p[1], p[4], p[7]] + sigma = [p[1], p[4], p[7]] a = [p[2], p[5], p[8]] - P = _multigaussfun(r,r0,fwhm,a) + P = _multigaussfun(r,r0,sigma,a) return P #================================================================= @@ -296,14 +298,14 @@ def dd_gengauss(*args): ------------------------------------------- Parameter Units Lower Upper Start ------------------------------------------- - Center nm 1 20 3.5 - FWHM nm 0.2 5 0.5 - Kurtosis 0.25 15 5 + Mean nm 1 20 3.5 + Spread nm 0.05 2.5 0.2 + Kurtosis 0.25 15 5 -------------------------------------------- """ if not args: info = dict( - Parameters = ('Center','FWHM','Kurtosis'), + Parameters = ('Mean','Spread','Kurtosis'), Units = ('nm','nm',''), Start = np.asarray([3.5, 0.5,0.5]), Lower = np.asarray([1, 0.2, 0.25]), @@ -314,9 +316,8 @@ def dd_gengauss(*args): # Compute the model distance distribution r0 = p[0] - width = p[1] + sigma = p[1] beta = p[2] - sigma = width/(2*np.sqrt(2*np.log(2))) x = abs(r-r0)/sigma P = beta/(2*sigma*spc.gamma(1/beta))*np.exp(-x**beta) P = _normalize(r,P) @@ -366,14 +367,14 @@ def dd_skewgauss(*args): ------------------------------------------- Parameter Units Lower Upper Start ------------------------------------------- - Location nm 1 20 3.5 - FWHM nm 0.2 5 0.5 + Center nm 1 20 3.5 + Spread nm 0.05 2.5 0.2 Skewness -25 25 5 -------------------------------------------- """ if not args: info = dict( - Parameters = ('Center','FWHM','Kurtosis'), + Parameters = ('Center','Spread','Kurtosis'), Units = ('nm','nm',''), Start = np.asarray([3.5, 0.5, 5]), Lower = np.asarray([1, 0.2, -25]), @@ -384,9 +385,8 @@ def dd_skewgauss(*args): # Compute the model distance distribution r0 = p[0] - width = p[1] + sigma = p[1] alpha = p[2] - sigma = width/(2*np.sqrt(2*np.log(2))) x = abs(r-r0)/sigma/np.sqrt(2) P = 1/np.sqrt(2*np.pi)*np.exp(-x**2)*(1 + spc.erf(alpha*x)) P = _normalize(r,P) @@ -436,13 +436,13 @@ def dd_rice(*args): ------------------------------------------------ Parameter Units Lower Upper Start ------------------------------------------------ - Center nm 1 10 3.5 - Width nm 0.1 5 0.7 + Location nm 1 10 3.5 + Spread nm 0.1 5 0.7 ------------------------------------------------ """ if not args: info = dict( - Parameters = ('Center','Width'), + Parameters = ('Location','Spread'), Units = ('nm','nm'), Start = np.asarray([3.5, 0.7]), Lower = np.asarray([1, 0.1]), @@ -500,19 +500,19 @@ def dd_rice2(*args): --------------------------------------------------------------- Parameter Units Lower Upper Start --------------------------------------------------------------- - Center of 1st Gaussian nm 1 10 2.5 - Width of 1st Gaussian nm 0.1 5 0.7 - Amplitude of 1sd Gaussian 0 1 0.5 - Center of 2nd Gaussian nm 1 10 4.0 - Width of 2nd Gaussian nm 0.1 5 0.7 - Amplitude of 2nd Gaussian 0 1 0.5 + Location of 1st Rician nm 1 10 2.5 + Spread of 1st Rician nm 0.1 5 0.7 + Amplitude of 1sd Rician 0 1 0.5 + Location of 2nd Rician nm 1 10 4.0 + Spread of 2nd Rician nm 0.1 5 0.7 + Amplitude of 2nd Rician 0 1 0.5 -------------------------------------------------------------- """ if not args: info = dict( - Parameters = ('Center of 1st Gaussian', 'Width of 1st Gaussian', 'Amplitude of 1st Gaussian', - 'Center of 2nd Gaussian', 'Width of 2nd Gaussian', 'Amplitude of 2nd Gaussian'), + Parameters = ('Location of 1st Rician', 'Spread of 1st Rician', 'Amplitude of 1st Rician', + 'Location of 2nd Rician', 'Spread of 2nd Rician', 'Amplitude of 2nd Rician'), Units = ('nm','nm','','nm','nm',''), Start = np.asarray([2.5, 0.7, 0.5, 4.0, 0.7, 0.5]), Lower = np.asarray([1, 0.1, 0, 1, 0.1, 0]), @@ -570,23 +570,23 @@ def dd_rice3(*args): --------------------------------------------------------------- Parameter Units Lower Upper Start --------------------------------------------------------------- - Center of 1st Gaussian nm 1 10 2.5 - Width of 1st Gaussian nm 0.1 5 0.7 - Amplitude of 1sd Gaussian 0 1 0.3 - Center of 2nd Gaussian nm 1 10 3.5 - Width of 2nd Gaussian nm 0.1 5 0.7 - Amplitude of 2nd Gaussian 0 1 0.3 - Center of 3rd Gaussian nm 1 10 5.0 - Width of 3rd Gaussian nm 0.1 5 0.7 - Amplitude of 3rd Gaussian 0 1 0.3 + Location of 1st Rician nm 1 10 2.5 + Spread of 1st Rician nm 0.1 5 0.7 + Amplitude of 1sd Rician 0 1 0.3 + Location of 2nd Rician nm 1 10 3.5 + Spread of 2nd Rician nm 0.1 5 0.7 + Amplitude of 2nd Rician 0 1 0.3 + Location of 3rd Rician nm 1 10 5.0 + Spread of 3rd Rician nm 0.1 5 0.7 + Amplitude of 3rd Rician 0 1 0.3 -------------------------------------------------------------- """ if not args: info = dict( - Parameters = ('Center of 1st Gaussian', 'Width of 1st Gaussian', 'Amplitude of 1st Gaussian', - 'Center of 2nd Gaussian', 'Width of 2nd Gaussian', 'Amplitude of 2nd Gaussian', - 'Center of 3rd Gaussian', 'Width of 3rd Gaussian', 'Amplitude of 3rd Gaussian'), + Parameters = ('Location of 1st Rician', 'Spread of 1st Rician', 'Amplitude of 1st Rician', + 'Location of 2nd Rician', 'Spread of 2nd Rician', 'Amplitude of 2nd Rician', + 'Location of 3rd Rician', 'Spread of 3rd Rician', 'Amplitude of 3rd Rician'), Units = ('nm','nm','','nm','nm','','nm','nm',''), Start = np.asarray([2.5, 0.7, 0.3, 3.5, 0.7, 0.3, 5, 0.7, 0.3]), Lower = np.asarray([1, 0.1, 0, 1, 0.1, 0, 1, 0.1, 0]), diff --git a/deerlab/deerload.py b/deerlab/deerload.py index 17562c212..4afe503b0 100644 --- a/deerlab/deerload.py +++ b/deerlab/deerload.py @@ -5,7 +5,7 @@ import matplotlib.pyplot as plt #------------------------------------------------------------------------------- -def deerload(fullbasename,Scaling=None,plot=False,*args,**kwargs): +def deerload(fullbasename,Scaling=None,plot=False,full_output=False,*args,**kwargs): r""" Load file in BES3T format (Bruker EPR Standard for Spectrum Storage and Transfer) @@ -22,17 +22,14 @@ def deerload(fullbasename,Scaling=None,plot=False,*args,**kwargs): Returns ------- t : ndarray - Time axis. + Time axis in microseconds. V : ndarray Experimental signal. pars : dict - Parameter file entries. + Parameter file entries, returned if ``full_output`` is ``True``. Notes ----- - Most commercial spectrometers save their data in nanoseconds. Since the - required time unit in DeerLab is microseconds, it is important to check - the values of ``t`` returned by ``deerload`` and convert to microseconds if necessary. Code based on BES3T version 1.2 (Xepr >=2.1). """ filename = fullbasename[:-4] @@ -109,7 +106,7 @@ def deerload(fullbasename,Scaling=None,plot=False,*args,**kwargs): if parDESC["IIFMT"] != parDESC["IRFMT"]: raise ValueError("IRFMT and IIFMT in DSC file must be identical.") - # Preallocation of theabscissa + # Preallocation of the abscissa maxlen = max(nx,ny,nz) abscissa = np.full((maxlen,3),np.nan) # Construct abscissa vectors @@ -258,6 +255,14 @@ def deerload(fullbasename,Scaling=None,plot=False,*args,**kwargs): else: warn('Cannot scale by temperature, since STMP in the DSC file is missing.') + # ns -> us converesion + abscissa /= 1e3 + + # Ensue proper numpy formatting + abscissa,data = np.atleast_1d(abscissa,data) + abscissa = np.squeeze(abscissa) + data = np.squeeze(data) + if plot: plt.plot(abscissa/1e3,np.real(data),abscissa/1e3,np.imag(data)) plt.xlabel("time (μs)") @@ -265,8 +270,11 @@ def deerload(fullbasename,Scaling=None,plot=False,*args,**kwargs): plt.grid() plt.show() - return abscissa, data, parameters - + if full_output: + return abscissa, data, parameters + else: + return abscissa, data + def read_description_file(DSCFileName): """ diff --git a/deerlab/ex_models.py b/deerlab/ex_models.py index 6189b85dd..df8262d98 100644 --- a/deerlab/ex_models.py +++ b/deerlab/ex_models.py @@ -1,10 +1,11 @@ - - -import numpy as np -from deerlab.utils import isempty +# ex_models.py - Experiment parametric models +# --------------------------------------------------------------------------- # This file is a part of DeerLab. License is MIT (see LICENSE.md). # Copyright(c) 2019-2020: Luis Fabregas, Stefan Stoll and other contributors. +import numpy as np + + def _parsargs(param,npar): #================================================================= param = np.atleast_1d(param) @@ -14,7 +15,7 @@ def _parsargs(param,npar): #================================================================= -def ex_4pdeer(param=[]): +def ex_4pdeer(param=None): # =================================================================== r""" Single-pathway 4-pulse DEER experiment model @@ -57,7 +58,7 @@ def ex_4pdeer(param=[]): Modulation depth 0 1 0.3 ------------------------------------------------------- """ - if isempty(param): + if param is None: info = dict( Parameters = ['Modulation depth'], Units = [''], @@ -76,7 +77,7 @@ def ex_4pdeer(param=[]): # =================================================================== -def ex_ovl4pdeer(param=[]): +def ex_ovl4pdeer(param=None): # =================================================================== r""" 4-pulse DEER with band overlap experiment model @@ -118,13 +119,14 @@ def ex_ovl4pdeer(param=[]): Amplitude of unmodulated components 0 1 0.3 Amplitude of 1st modulated pathway 0 1 0.3 Amplitude of 2nd modulated pathway 0 1 0.3 - Refocusing time of 2nd modulated pathway us 0 20 5 + Refocusing time of 2nd modulated pathway μs 0 20 5 --------------------------------------------------------------------------- """ - if isempty(param): + if param is None: info = dict( - Parameters = ['Amplitude of unmodulated components','Amplitude of 1st modulated pathway','Amplitude of 2nd modulated pathway','Refocusing time of 2nd modulated pathway'], - Units = ['','','','us'], + Parameters = ['Amplitude of unmodulated components','Amplitude of 1st modulated pathway', + 'Amplitude of 2nd modulated pathway','Refocusing time of 2nd modulated pathway'], + Units = ['','','','μs'], Start = np.asarray([0.3, 0.3, 0.3,5]), Lower = np.asarray([0, 0, 0, 0]), Upper = np.asarray([1, 1, 1, 20]) @@ -142,7 +144,7 @@ def ex_ovl4pdeer(param=[]): -def ex_5pdeer(param=[]): +def ex_5pdeer(param=None): # =================================================================== r""" 5-pulse DEER experiment model @@ -184,13 +186,13 @@ def ex_5pdeer(param=[]): Amplitude of unmodulated components 0 1 0.3 Amplitude of 1st modulated pathway 0 1 0.3 Amplitude of 2nd modulated pathway 0 1 0.3 - Refocusing time of 2nd modulated pathway us 0 20 5 + Refocusing time of 2nd modulated pathway μs 0 20 5 --------------------------------------------------------------------------- """ - if isempty(param): + if param is None: info = dict( Parameters = ['Amplitude of unmodulated components','Amplitude of 1st modulated pathway','Amplitude of 2nd modulated pathway','Refocusing time of 2nd modulated pathway'], - Units = ['','','','us'], + Units = ['','','','μs'], Start = np.asarray([0.3, 0.3, 0.3,5]), Lower = np.asarray([0, 0, 0, 0]), Upper = np.asarray([1, 1, 1, 20]) @@ -208,7 +210,7 @@ def ex_5pdeer(param=[]): -def ex_7pdeer(param=[]): +def ex_7pdeer(param=None): # =================================================================== r""" 7-pulse DEER experiment model @@ -252,14 +254,16 @@ def ex_7pdeer(param=[]): Amplitude of 1st modulated pathway 0 1 0.5 Amplitude of 2nd modulated pathway 0 1 0.3 Amplitude of 3rd modulated pathway 0 1 0.2 - Refocusing time of 2nd modulated pathway us 0 20 1.5 - Refocusing time of 3rd modulated pathway us 0 20 3.5 + Refocusing time of 2nd modulated pathway μs 0 20 1.5 + Refocusing time of 3rd modulated pathway μs 0 20 3.5 --------------------------------------------------------------------------- """ - if isempty(param): + if param is None: info = dict( - Parameters = ['Amplitude of unmodulated components','Amplitude of 1st modulated pathway','Amplitude of 2nd modulated pathway','Refocusing time of 2nd modulated pathway'], - Units = ['','','','us'], + Parameters = ['Amplitude of unmodulated components','Amplitude of 1st modulated pathway', + 'Amplitude of 2nd modulated pathway','Amplitude of 3rd modulated pathway', + 'Refocusing time of 2nd modulated pathway','Refocusing time of 3rd modulated pathway'], + Units = ['','','','','μs','μs'], Start = np.asarray([0.3, 0.5, 0.3, 0.2, 1.5, 3.5]), Lower = np.asarray([0, 0, 0, 0, 0, 0]), Upper = np.asarray([1, 1, 1, 1, 20, 20]) diff --git a/deerlab/fftspec.py b/deerlab/fftspec.py index 3a94d9197..629e6d109 100644 --- a/deerlab/fftspec.py +++ b/deerlab/fftspec.py @@ -35,7 +35,7 @@ def fftspec(V,t,mode='abs',zerofilling='auto',apodization=True): Use of a Hamming apodization window, the default is True. """ - if zerofilling is 'auto': + if zerofilling == 'auto': zerofilling = 2*len(V) #If requested apply Hamming apodization window @@ -48,11 +48,11 @@ def fftspec(V,t,mode='abs',zerofilling='auto',apodization=True): spec = fftshift(fft(V,zerofilling)) #Get the requested component/type of spectrum - if mode is 'abs': + if mode == 'abs': spec = np.abs(spec) - elif mode is 'real': + elif mode == 'real': spec = spec.real - elif mode is 'imag': + elif mode == 'imag': spec = spec.imag else: raise KeyError("Invalid spectrum mode. Must be 'abs', 'real', or 'imag'. ") diff --git a/deerlab/fitmultimodel.py b/deerlab/fitmultimodel.py index 53787a833..9db726de8 100644 --- a/deerlab/fitmultimodel.py +++ b/deerlab/fitmultimodel.py @@ -6,6 +6,7 @@ import copy import numpy as np import numdifftools as nd +import matplotlib.pyplot as plt from types import FunctionType import deerlab as dl from deerlab.utils import hccm, goodness_of_fit @@ -74,6 +75,11 @@ def fitmultimodel(V, Kmodel, r, model, maxModels, method='aic', lb=None, ub=None Selection functional values (as specified as ``method``) for the all fitted multi-component models. scale : float int or list of float int Amplitude scale(s) of the dipolar signal(s). + plot : callable + Function to display the results. It will display the + fitted signals, the distance distribution with confidence intervals, + and the values of the selection functional. If requested, the function + returns the `matplotlib.axes` object as output. stats : dict Goodness of fit statistical estimators: @@ -180,7 +186,7 @@ def K4pdeer(Kpar): if len(lbK) is not nKparam or len(ubK) is not nKparam: raise ValueError('The upper/lower bounds of the kernel parameters must be ',nKparam,'-element arrays') - areCenterDistances = [str in ['Center','Location'] for str in paramnames] + areCenterDistances = [str in ['Mean','Location'] for str in paramnames] if any(areCenterDistances): # If the center of the basis function is a parameter limit it # to the distance axis range (stabilizes parameter search) @@ -234,7 +240,7 @@ def Pmodel(nlinpar,linpar): def logestimators(V,Vfit,plin,pnonlin,functionals): #=============================================================================== """ - Log-Likelihood Estimators + Log-Likelihood Estimators --------------------------- Computes the estimated likelihood of a multi-component model being the optimal choice. @@ -260,6 +266,24 @@ def logestimators(V,Vfit,plin,pnonlin,functionals): return functionals #=============================================================================== + def spread_within_box(lb,ub,Nmodels): + #=============================================================================== + """ + Start values within boundary box + --------------------------------- + Computes the start values of the parameters of the multiple basis functions + spread equally within the box constraints. + """ + par = [] + for lbi,ubi in zip(lb,ub): + par.append(np.linspace(lbi,ubi,Nmodels+2)[1:-1]) + par0 = [] + for i in range(Nmodels): + for pari in par: + par0.append(pari[i]) + return par0 + #=============================================================================== + # Pre-allocate containers fits,Vfit,Pfit,plin_,pnonlin_,nlin_ub_,nlin_lb_,lin_ub_,lin_lb_ = ([] for _ in range(9)) logest = [] @@ -272,6 +296,11 @@ def logestimators(V,Vfit,plin,pnonlin,functionals): # =========================================== Knonlin = lambda par: nonlinmodel(par,Nmodels) + # Start values of non-linear parameters + par0_K = (ubK - lbK)/2 # start in the middle + par0_P = spread_within_box(lb,ub,Nmodels) # start spread within boundaries + par0 = np.concatenate((par0_K, par0_P),axis=None) + # Box constraints for the model parameters (non-linear parameters) nlin_lb = np.tile(lb,(1,Nmodels)) nlin_ub = np.tile(ub,(1,Nmodels)) @@ -280,10 +309,6 @@ def logestimators(V,Vfit,plin,pnonlin,functionals): nlin_lb = np.concatenate((lbK, nlin_lb),axis=None) nlin_ub = np.concatenate((ubK, nlin_ub),axis=None) - # Start values of non-linear parameters - np.random.seed(1) - par0 = np.random.uniform(size=(len(nlin_lb)),low=nlin_lb, high=nlin_ub) - # Box constraints for the components amplitudes (linear parameters) lin_lb = np.ones(Nmodels) lin_ub = np.full(Nmodels,np.inf) @@ -291,7 +316,7 @@ def logestimators(V,Vfit,plin,pnonlin,functionals): # Separable non-linear least-squares (SNLLS) fit scale = 1e2 fit = dl.snlls(V*scale,Knonlin,par0,nlin_lb,nlin_ub,lin_lb,lin_ub, - weights=weights, penalty=False, uqanalysis=False, lin_tol=[], lin_maxiter=[],**kwargs) + weights=weights, reg=False, uqanalysis=False, lin_tol=[], lin_maxiter=[],**kwargs) pnonlin = fit.nonlin plin = fit.lin @@ -384,7 +409,54 @@ def logestimators(V,Vfit,plin,pnonlin,functionals): if len(scales)==1: scales = scales[0] + # Results display function + plotfcn = lambda: _plot(Vsubsets,V,Vfit,r,Pfit,Puq,fcnals,maxModels,method) + return FitResult(P=Pfit, Pparam=fitparam_P, Kparam=fitparam_K, amps=fitparam_amp, Puncert=Puq, - paramUncert=paramuq, selfun=fcnals, Nopt=Nopt, Pn=Peval, scale=scales, + paramUncert=paramuq, selfun=fcnals, Nopt=Nopt, Pn=Peval, scale=scales, plot=plotfcn, stats=stats, cost=fit.cost, residuals=fit.residuals, success=fit.success) - # ========================================================================= \ No newline at end of file +# ========================================================================= + + +def _plot(Vsubsets,V,Vfit,r,Pfit,Puq,fcnals,maxModels,method): +# ========================================================================= + nSignals = len(Vsubsets) + _,axs = plt.subplots(nSignals+1,figsize=[7,3+3*nSignals]) + for i in range(nSignals): + subset = Vsubsets[i] + # Plot the experimental signal and fit + axs[i].plot(V[subset],'.',color='grey',alpha=0.5) + axs[i].plot(Vfit[subset],'tab:blue') + axs[i].grid(alpha=0.3) + axs[i].set_xlabel('Array Elements') + axs[i].set_ylabel('V[{}]'.format(i)) + axs[i].legend(('Data','Fit')) + + # Confidence intervals of the fitted distance distribution + Pci95 = Puq.ci(95) # 95#-confidence interval + Pci50 = Puq.ci(50) # 50#-confidence interval + + ax = plt.subplot(nSignals+1,2,2*(nSignals+1)-1) + ax.plot(r,Pfit,color='tab:blue',linewidth=1.5) + ax.fill_between(r,Pci50[:,0],Pci50[:,1],color='tab:blue',linestyle='None',alpha=0.45) + ax.fill_between(r,Pci95[:,0],Pci95[:,1],color='tab:blue',linestyle='None',alpha=0.25) + ax.grid(alpha=0.3) + ax.legend(['truth','optimal fit','95%-CI']) + ax.set_xlabel('Distance [nm]') + ax.set_ylabel('P [nm⁻¹]') + axs = np.append(axs,ax) + + # Compute the Akaike weights + dfcnals = fcnals - min(fcnals) + ax = plt.subplot(nSignals+1,2,2*(nSignals+1)) + ax.bar(np.arange(maxModels)+1,dfcnals + abs(min(dfcnals)),facecolor='tab:blue',alpha=0.6) + ax.grid(alpha=0.3) + ax.set_ylabel('Δ{}'.format(method.upper())) + ax.set_xlabel('Number of components') + axs = np.append(axs,ax) + + plt.tight_layout() + plt.show() + return axs +# ========================================================================= + diff --git a/deerlab/fitparamodel.py b/deerlab/fitparamodel.py index 555a13b36..29afa3c7c 100644 --- a/deerlab/fitparamodel.py +++ b/deerlab/fitparamodel.py @@ -7,6 +7,7 @@ import numdifftools as nd from deerlab.utils import isempty, multistarts, hccm, parse_multidatasets, goodness_of_fit from deerlab.classes import UncertQuant, FitResult +import matplotlib.pyplot as plt from scipy.optimize import least_squares import warnings @@ -37,6 +38,10 @@ def fitparamodel(V, model, par0=[],lb=[],ub=[], weights = 1, Covariance-based uncertainty quantification of the fitted parameters. scale : float int or list of float int Amplitude scale(s) of the dipolar signal(s). + plot : callable + Function to display the results. It will + display the fitted signals. If requested, the function returns + the `matplotlib.axes` object as output. stats : dict Goodness of fit statistical estimators: @@ -193,13 +198,13 @@ def lsqresiduals(p): atUpper = abs(parfit-ub)0: print() info = dd_model() - ci = paruq['dd'].ci(95) + if uqanalysis: + ci = paruq['dd'].ci(95) + else: + ci = np.full((len(parfit['dd']),2),np.nan) for j in range(len(parfit['dd'])): c = parfit['dd'][j] print(pstr.format('ddparam',j,c,ci[j,0],ci[j,1],info['Parameters'][j],info['Units'][j])) @@ -571,25 +588,25 @@ def _display_results(): if includeBackground[i]: if len(parfit['bg'])>0: info = bg_model[i]() + if uqanalysis: ci = paruq['bg'][i].ci(95) - for j in range(len(parfit['bg'][i])): + else: + ci = np.full((len(parfit['bg'][i]),2),np.nan) + for j in range(len(parfit['bg'][i])): c = parfit['bg'][i][j] print(pstr.format('bgparam',j,c,ci[j,0],ci[j,1],info['Parameters'][j],info['Units'][j])) if includeExperiment[i]: if len(parfit['ex'])>0: info = ex_model[i]() + if uqanalysis: ci = paruq['ex'][i].ci(95) + else: + ci = np.full((len(parfit['ex'][i]),2),np.nan) for j in range(len(parfit['ex'][i])): c = parfit['ex'][i][j] print(pstr.format('exparam',j,c,ci[j,0],ci[j,1],info['Parameters'][j],info['Units'][j])) print('----------------------------------------------------------------------------') - # ========================================================================= - - # Plotting - # -------- - if display: - _display_results() # Return numeric arrays and not lists if there is only one signal if nSignals==1: @@ -607,7 +624,7 @@ def _display_results(): return FitResult(V=Vfit, P=Pfit, B=Bfit, exparam=parfit['ex'], bgparam=parfit['bg'], ddparam=parfit['dd'], Vuncert = modfituq['Vfit'], Puncert = modfituq['Pfit'], Buncert = modfituq['Bfit'], exparamUncert = paruq['ex'], bgparamUncert = paruq['bg'], - ddparamUncert = paruq['dd'], scale=scales, stats=stats, cost=fit.cost, + ddparamUncert = paruq['dd'], regparam = alphaopt, plot=_display_results, scale=scales, stats=stats, cost=fit.cost, residuals=fit.residuals, success=fit.success) diff --git a/deerlab/nnls.py b/deerlab/nnls.py index 2c36af3da..0e09b33cd 100644 --- a/deerlab/nnls.py +++ b/deerlab/nnls.py @@ -121,6 +121,14 @@ def cvxnnls(AtA, Atb, tol=[], maxiter=[]): #===================================================================================== """ NNLS problem solved via CVXOPT + + + References: + ----------- + [1] + Rein, Lewe, Andrade, Kacprzak and Weber. J. Magn. Reson., 295 (2018) 17–26. + Global analysis of complex PELDOR time traces + https://doi.org/10.1016/j.jmr.2018.07.015 """ N = np.shape(AtA)[1] @@ -259,4 +267,4 @@ def nnlsbpp(AtA,AtB,x0=None): nInfeasible = sum(xFnegative | yGnegative) return x -#===================================================================================== \ No newline at end of file +#===================================================================================== diff --git a/deerlab/noiselevel.py b/deerlab/noiselevel.py index fd00a395c..7b7c6c082 100644 --- a/deerlab/noiselevel.py +++ b/deerlab/noiselevel.py @@ -92,18 +92,18 @@ def noiselevel(V,*args): # Estimation of the noise level # ----------------------------- - if estimationMethod is '2D': + if estimationMethod == '2D': # Estimate standard deviations for all time point, and average over scans if shape(V)[1] < 10: raise Warning('Only a few scans are given. Noise standard deviation estimate will be inaccurate.') sigma = std(V,1) sigma = mean(sigma) - elif estimationMethod is 'filtering': + elif estimationMethod == 'filtering': # Filter the noise in the signal - if filterType is 'movmean': + if filterType == 'movmean': Vfilt = movmean(V,3) - elif filterType is 'savgol': + elif filterType == 'savgol': with warnings.catch_warnings(): warnings.simplefilter('ignore') Vfilt = savgol_filter(V,11,3) @@ -112,13 +112,13 @@ def noiselevel(V,*args): # And estimate the noiselevel from the resulting residual sigma = std(V - Vfilt) - elif estimationMethod is 'complex': + elif estimationMethod == 'complex': # Optimize the phase of the signal _,Vim,_,_ = correctphase(V,full_output=True) # And estimate the noiselevel from the imaginary part sigma = std(Vim) - elif estimationMethod is 'reference': + elif estimationMethod == 'reference': sigma = std(V - Vref) return sigma diff --git a/deerlab/snlls.py b/deerlab/snlls.py index 0d18c727c..514f1be63 100644 --- a/deerlab/snlls.py +++ b/deerlab/snlls.py @@ -1,8 +1,13 @@ +# snlls.py - Separable non-linear least-squares solver +# --------------------------------------------------------------- +# This file is a part of DeerLab. License is MIT (see LICENSE.md). +# Copyright(c) 2019-2020: Luis Fabregas, Stefan Stoll and other contributors. import copy import numpy as np import numdifftools as nd from scipy.optimize import least_squares, lsq_linear +import matplotlib.pyplot as plt from numpy.linalg import solve # Import DeerLab depencies @@ -11,7 +16,7 @@ from deerlab.nnls import cvxnnls, fnnls, nnlsbpp from deerlab.classes import UncertQuant, FitResult -def snlls(y, Amodel, par0, lb=None, ub=None, lbl=None, ubl=None, nnlsSolver='cvx', penalty=None, weights=1, +def snlls(y, Amodel, par0, lb=None, ub=None, lbl=None, ubl=None, nnlsSolver='cvx', reg='auto', weights=1, regtype='tikhonov', regparam='aic', multistart=1, regorder=2, alphareopt=1e-3, nonlin_tol=1e-9, nonlin_maxiter=1e8, lin_tol=1e-15, lin_maxiter=1e4, huberparam=1.35, uqanalysis=True): @@ -50,6 +55,11 @@ def snlls(y, Amodel, par0, lb=None, ub=None, lbl=None, ubl=None, nnlsSolver='cvx * ``paramuq.ci(n)`` - n%-CI of the full parameter set * ``paramuq.ci(n,'lin')`` - n%-CI of the linear parameter set * ``paramuq.ci(n,'nonlin')`` - n%-CI of the non-linear parameter set + regparam : scalar + Regularization parameter value used for the regularization of the linear parameters. + plot : callable + Function to display the results. It will display the fitted data. + If requested, the function returns the `matplotlib.axes` object as output. stats : dict Goodness of fit statistical estimators @@ -68,10 +78,14 @@ def snlls(y, Amodel, par0, lb=None, ub=None, lbl=None, ubl=None, nnlsSolver='cvx Other parameters ---------------- - penalty : boolean - Forces the use of a regularization penalty on the solution of the linear problem. - If not specified it is determined automatically based con the condition number of - the non-linear model ``Amodel``. + reg : boolean or string + Determines the use of regularization on the solution of the linear problem. + + * ``'auto'`` - Automatic decision based con the condition number of the non-linear model ``Amodel``. + * ``True`` - Forces regularization regardless of the condition number + * ``False`` - Disables regularization regardless of the condition number + The default is ``'auto'``. + regType : string Regularization penalty type: @@ -171,20 +185,21 @@ def snlls(y, Amodel, par0, lb=None, ub=None, lbl=None, ubl=None, nnlsSolver='cvx par0 = np.atleast_1d(par0) # Parse multiple datsets and non-linear operators into a single concatenated vector/matrix - y, Amodel, weights, subsets, prescales = dl.utils.parse_multidatasets(y, Amodel, weights, precondition=True) + y, Amodel, weights, subsets, _ = dl.utils.parse_multidatasets(y, Amodel, weights, precondition=True) # Get info on the problem parameters and non-linear operator A0 = Amodel(par0) Nnonlin = len(par0) Nlin = np.shape(A0)[1] linfit = np.zeros(Nlin) + alpha = 0 # Determine whether to use regularization penalty illConditioned = np.linalg.cond(A0) > 10 - if illConditioned and penalty is None: - includePenalty = True + if reg == 'auto': + includePenalty = illConditioned else: - includePenalty = penalty + includePenalty = reg # Checks for bounds constraints # ---------------------------------------------------------- @@ -247,11 +262,11 @@ def snlls(y, Amodel, par0, lb=None, ub=None, lbl=None, ubl=None, nnlsSolver='cvx elif linearConstrained and nonNegativeOnly: # Non-negative linear LSQ - if nnlsSolver is 'fnnls': + if nnlsSolver == 'fnnls': linSolver = lambda AtA, Aty: fnnls(AtA, Aty, tol=lin_tol, maxiter=lin_maxiter) - elif nnlsSolver is 'nnlsbpp': + elif nnlsSolver == 'nnlsbpp': linSolver = lambda AtA, Aty: nnlsbpp(AtA, Aty, np.linalg.solve(AtA, Aty)) - elif nnlsSolver is 'cvx': + elif nnlsSolver == 'cvx': linSolver = lambda AtA, Aty: cvxnnls(AtA, Aty, tol=lin_tol, maxiter=lin_maxiter) parseResult = lambda result: result # ---------------------------------------------------------- @@ -270,7 +285,7 @@ def ResidualsFcn(p): non-linear least-squares solver. """ - nonlocal par_prev, check, regparam_prev, linfit + nonlocal par_prev, check, regparam_prev, linfit, alpha # Non-linear model evaluation A = Amodel(p) @@ -397,7 +412,10 @@ def ci(coverage,ptype='full'): if len(stats) == 1: stats = stats[0] - return FitResult(nonlin=nonlinfit, lin=linfit, uncertainty=paramuq, + # Display function + plotfcn = lambda: _plot(subsets,y,yfit) + + return FitResult(nonlin=nonlinfit, lin=linfit, uncertainty=paramuq, regparam=alpha, plot=plotfcn, stats=stats, cost=fvals, residuals=sol.fun, success=sol.success) # =========================================================================================== @@ -414,13 +432,13 @@ def _augment(res, J, regtype, alpha, L, x, eta, Nnonlin): """ eps = np.finfo(float).eps # Compute the regularization penalty augmentation for the residual and the Jacobian - if regtype is 'tikhonov': + if regtype == 'tikhonov': resreg = L@x Jreg = L - elif regtype is 'tv': + elif regtype == 'tv': resreg =((L@x)**2 + eps)**(1/4) Jreg = 2/4*((( ( (L@x)**2 + eps)**(-3/4) )*(L@x))[:, np.newaxis]*L) - elif regtype is 'huber': + elif regtype == 'huber': resreg = np.sqrt(np.sqrt((L@x/eta)**2 + 1) - 1) Jreg = 0.5/(eta**2)*((((np.sqrt((L@x/eta)**2 + 1) - 1 + eps)**(-1/2)*(((L@x/eta)**2 + 1+ eps)**(-1/2)))*(L@x))[:, np.newaxis]*L) @@ -436,3 +454,23 @@ def _augment(res, J, regtype, alpha, L, x, eta, Nnonlin): return res, J # =========================================================================================== + + +def _plot(subsets,y,yfit): +# =========================================================================================== + nSignals = len(subsets) + _,axs = plt.subplots(nSignals+1,figsize=[7,3*nSignals]) + for i in range(nSignals): + subset = subsets[i] + # Plot the experimental signal and fit + axs[i].plot(y[subset],'.',color='grey',alpha=0.5) + axs[i].plot(yfit[subset],'tab:blue') + axs[i].grid(alpha=0.3) + axs[i].set_xlabel('Array elements') + axs[i].set_ylabel('Data #{}'.format(i)) + axs[i].legend(('Data','Fit')) + plt.tight_layout() + plt.show() + return axs +# =========================================================================================== + diff --git a/deerlab/whitegaussnoise.py b/deerlab/whitegaussnoise.py index b8366bd29..68c9914af 100644 --- a/deerlab/whitegaussnoise.py +++ b/deerlab/whitegaussnoise.py @@ -1,37 +1,47 @@ import numpy as np -def whitegaussnoise(t,level,rescale=False): - r"""White Gaussian noise generator - +def whitegaussnoise(t,level=1,rescale=False,seed=None): + r""" + Generates a vector of white Gaussian (normal) noise + Parameters ---------- t : array_like Time axis. level : float scalar - Noise level (standard deviation of noise vector) + Noise level, i.e. standard deviation of underlying Gaussian distribution. rescale : boolean - Enable/disable rescaling of noise vector to enforce the exact noise level. - + If true, rescales the noise vector such that its standard deviation is exactly equal + to ``level``. If false (default), the standard deviation of the noise vector can deviate + slightly from ``level``, particularly for short vectors. + seed : None or integer scalar + If None (default), do not seed the random number generator. If an integer scalar is + given (e.g. ``seed=137``), seed the random number generator with this number. + Returns ------- noise : ndarray Noise vector. - + Notes ----- - The noise vector is generated from a vector of pseudo-random numbers generated with numpy. - Each call of ``whitegaussnoise`` will return a different realization of the noise vector. To set the output - to a fixed noise realization, the random number generator must be fixed by calling ``numpy.random.seed(k)`` - before ``whitegaussnoise``, where ``k`` is some integer number. - + The noise vector is generated using pseudo-random numbers generated with NumPy. + Without seeding the random number generator, subsequent calls of ``whitegaussnoise`` return + different realizations of the noise vector. To obtain a reproducible noise realization, seed + the random number generator by using the ``seed`` kewyword arguement, or call ``numpy.random.seed(k)`` + with some integer number ``k`` before calling ``whitegaussnoise``. """ - t = np.atleast_1d(t) - N = len(t) - + + # Seed RNG if wanted + if seed is not None: + np.random.seed(seed) + + # Draw from standard normal distribution + N = len(np.atleast_1d(t)) noise = np.random.normal(0,level,N) - + + # Rescale to sample std if wanted if rescale: - noise = level*noise/np.std(noise) - + noise = level/np.std(noise)*noise + return noise - diff --git a/docsrc/source/basics.rst b/docsrc/source/basics.rst index 7a8a09d5c..fd06c5bc7 100644 --- a/docsrc/source/basics.rst +++ b/docsrc/source/basics.rst @@ -24,7 +24,7 @@ DeerLab distinguishes between **non-parametric** and **parametric** distance dis .. code-block:: python r = np.linspace(1.5,6,201) # distance, in nanometers - P = dd_gauss(r,[3 0.5]) # Gaussian distribution, center 3 nm, fwhm 0.5 nm + P = dd_gauss(r,[3 0.2]) # Gaussian distribution, center 3 nm, fwhm 0.5 nm plt.plot(r,P) The first line generates a distance axis ``r`` as a row vector. The second line generates the distance distribution ``P`` as a vector defined over ``r``. The function ``dd_gauss`` is a distance distribution model function provided by DeerLab. It takes as inputs the distance axis and a two-element parameter array (center and full width at half maximum) and returns the calculated distance distribution consisting of one Gaussian, truncated to the range of ``r`` and normalized. diff --git a/docsrc/source/firstscript.rst b/docsrc/source/firstscript.rst index 6d72c2ff2..756110034 100644 --- a/docsrc/source/firstscript.rst +++ b/docsrc/source/firstscript.rst @@ -7,9 +7,9 @@ We start by importing the required function/libraries: .. code-block:: python - from numpy import np - from matplotlib.pyplot import plt - from deerlab import * + import numpy as np + import matplotlib.pyplot as plt + import deerlab as dl We continue by generating a distance distribution consisting of a single Gaussian: @@ -17,7 +17,7 @@ We continue by generating a distance distribution consisting of a single Gaussia N = 201 # number of points r = np.linspace(1.5,7,N) # distance range, in nanometers - P = dd_gauss(r,[3.5, 0.4]) # single-Gaussian distance distribution + P = dl.dd_gauss(r,[3.5, 0.15]) # single-Gaussian distance distribution Next, we calculate the background decay function due to a homogeneus 3D distribution of spins: @@ -26,42 +26,33 @@ Next, we calculate the background decay function due to a homogeneus 3D distribu t = np.linspace(0,3,N) # time axis, in microseconds conc = 100 # spin concentration, in micromolar lam = 0.4 # modulation depth - B = bg_hom3d(t,conc,lam) # homogeneous 3D background decay function + B = dl.bg_hom3d(t,conc,lam) # homogeneous 3D background decay function Next, we combine the distance distribution and the background into a full 4-pulse DEER signal and add some noise: .. code-block:: python - K = dipolarkernel(t,r,lam,B) # DEER kernel - V = K@P # DEER signal - sig = 0.01 # noise level - Vexp = V + whitegaussnoise(V,sig) # add noise + K = dl.dipolarkernel(t,r,lam,B) # DEER kernel + V = K@P # DEER signal + sig = 0.01 # noise level + Vexp = V + dl.whitegaussnoise(V,sig) # add noise We can look at the result: .. code-block:: python - plt.plot(t,V,t,Vexp,t,(1-lam)*B) # plotting + plt.plot(t,V,t,Vexp) # plotting plt.show() Now that we have a noisy DEER trace, we fit it (in a single step) with a non-parametric distance distribution and a homogeneous 3D background. .. code-block:: python - fit = fitsignal(Vexp,t,r,'P',bg_hom3d,ex_4pdeer) # fitting + fit = dl.fitsignal(Vexp,t,r,'P',dl.bg_hom3d,dl.ex_4pdeer) # fitting Finally, we plot the results .. code-block:: python - Pfit = fit.P # fitted distance distribution - Bfit = fit.B # fitted background - Vfit = fit.V # fitted dipolar signal - lamfit = fit.exparam # fitted modulation depth - - # plotting - plt.subplot(2,1,1) - plt.plot(t,Vexp,t,Vfit,t,(1-lamfit)*Bfit) # plot fitted model and background - plt.subplot(2,1,2) - plt.plot(r,P,r,Pfit) # plot fitted distribution - plt.show() + # Plotting + fit.plot() diff --git a/docsrc/source/models/bg_exp.rst b/docsrc/source/models/bg_exp.rst index e8c1fe223..69eafd2e6 100644 --- a/docsrc/source/models/bg_exp.rst +++ b/docsrc/source/models/bg_exp.rst @@ -14,11 +14,11 @@ Model B(t) = \exp\left(-\lambda\kappa \vert t \vert\right) -============== =============== ========= ============= ============= ============================== - Variable Symbol Default Lower bound Upper bound Description -============== =============== ========= ============= ============= ============================== -``param(1)`` :math:`\kappa` 0.35 0 200 Decay rate -============== =============== ========= ============= ============= ============================== +============== =============== ============= ============= ============= ================================ + Variable Symbol Start Value Lower bound Upper bound Description +============== =============== ============= ============= ============= ================================ +``param[0]`` :math:`\kappa` 0.35 0 200 Decay rate (μs\ :sup:`-1`) +============== =============== ============= ============= ============= ================================ Although the ``bg_exp`` model has the same functional form as ``bg_hom3d``, it is distinct since its parameter is a decay rate constant and not a spin concentration like for ``bg_hom3d``. diff --git a/docsrc/source/models/bg_hom3d.rst b/docsrc/source/models/bg_hom3d.rst index 73c45db3a..7cc9e73c1 100644 --- a/docsrc/source/models/bg_hom3d.rst +++ b/docsrc/source/models/bg_hom3d.rst @@ -25,8 +25,9 @@ where `c_p` is the pumped-spin concentration (entered in spins/m\ :sup:`3` into .. math:: D = \frac{\mu_0}{4\pi}\frac{(g_\mathrm{e}\mu_\mathrm{B})^2}{\hbar} -============= ============= ========= ============= ============= ============================================= - Variable Symbol Default Lower bound Upper bound Description -============= ============= ========= ============= ============= ============================================= -``param(1)`` :math:`c_p` 50 0.01 5000 Pumped spin concentration (μM) -============= ============= ========= ============= ============= ============================================= +============== =============== ============= ============= ============= ================================= + Variable Symbol Start Value Lower bound Upper bound Description +============== =============== ============= ============= ============= ================================= +``param[0]`` :math:`c_p` 50 0.01 5000 Pumped spin concentration (μM) +============== =============== ============= ============= ============= ================================= + diff --git a/docsrc/source/models/bg_hom3dex.rst b/docsrc/source/models/bg_hom3dex.rst index 44c19f02f..27b9614fb 100644 --- a/docsrc/source/models/bg_hom3dex.rst +++ b/docsrc/source/models/bg_hom3dex.rst @@ -27,9 +27,9 @@ where :math:`c` is the spin concentration (entered in spins/m\ :sup:`3` into thi The function :math:`\alpha(R)` of the exclusion distance :math:`R` captures the excluded-volume effect. It is a smooth function, but doesn't have an analytical representation. For details, see `Kattnig et al, J.Phys.Chem.B 2013, 117, 16542 `_. -============= =================== ========= ============= ============= ================================================ - Variable Symbol Default Lower bound Upper bound Description -============= =================== ========= ============= ============= ================================================ -``param(1)`` :math:`c` 50 0.01 1000 spin concentration (μM) -``param(2)`` :math:`R` 1 0.1 20 exclusion distance (nm) -============= =================== ========= ============= ============= ================================================ +============== =============== ============= ============= ============= ================================= + Variable Symbol Start Value Lower bound Upper bound Description +============== =============== ============= ============= ============= ================================= +``param[0]`` :math:`c` 50 0.01 1000 Spin concentration (μM) +``param[1]`` :math:`R` 1 0.1 20 Exclusion distance (nm) +============== =============== ============= ============= ============= ================================= \ No newline at end of file diff --git a/docsrc/source/models/bg_homfractal.rst b/docsrc/source/models/bg_homfractal.rst index f5e134695..a7476990c 100644 --- a/docsrc/source/models/bg_homfractal.rst +++ b/docsrc/source/models/bg_homfractal.rst @@ -14,9 +14,9 @@ Model This implements the background due to a homogeneous distribution of spins in a d-dimensional space, with d-dimensional spin concentration ``c_d``. -============= ============= ========= ============= ============= ========================================================== - Variable Symbol Default Lower bound Upper bound Description -============= ============= ========= ============= ============= ========================================================== -``param(1)`` :math:`c_d` 50 0.01 5000 Pumped spin fractal concentration (μmol/dm\ :sup:`d`) -``param(2)`` :math:`d` 3 0 6 Fractal dimension -============= ============= ========= ============= ============= ========================================================== +============= ============= ============= ============= ============= ========================================================== + Variable Symbol Start Value Lower bound Upper bound Description +============= ============= ============= ============= ============= ========================================================== +``param[0]`` :math:`c_d` 50 0.01 5000 Pumped spin fractal concentration (μmol/dm\ :sup:`d`) +``param[1]`` :math:`d` 3 0 6 Fractal dimension +============= ============= ============= ============= ============= ========================================================== diff --git a/docsrc/source/models/bg_poly1.rst b/docsrc/source/models/bg_poly1.rst index 8eb155f02..0a17af713 100644 --- a/docsrc/source/models/bg_poly1.rst +++ b/docsrc/source/models/bg_poly1.rst @@ -12,9 +12,9 @@ Model :math:`B(t) = p_0 + p_1(\lambda t)` -============= ============= ========= ============= ============= ============================== - Variable Symbol Default Lower bound Upper bound Description -============= ============= ========= ============= ============= ============================== -``param(1)`` :math:`p_0` 1 0 200 Intercept -``param(2)`` :math:`p_1` -1 -200 200 1st order weight -============= ============= ========= ============= ============= ============================== +============== =============== ============= ============= ============= ==================================== + Variable Symbol Start Value Lower bound Upper bound Description +============== =============== ============= ============= ============= ==================================== +``param[0]`` :math:`p_0` 1 0 200 Intercept +``param[1]`` :math:`p_1` -1 -200 200 1st order weight (μs\ :sup:`-1`) +============== =============== ============= ============= ============= ==================================== \ No newline at end of file diff --git a/docsrc/source/models/bg_poly2.rst b/docsrc/source/models/bg_poly2.rst index 207b34069..dc59b4ba4 100644 --- a/docsrc/source/models/bg_poly2.rst +++ b/docsrc/source/models/bg_poly2.rst @@ -13,10 +13,10 @@ Model :math:`B(t) = p_0 + p_1(\lambda t) + p_2(\lambda t)^2` -============= ============= ========= ============= ============= ============================== - Variable Symbol Default Lower bound Upper bound Description -============= ============= ========= ============= ============= ============================== -``param(1)`` :math:`p_0` 1 0 200 Intercept -``param(2)`` :math:`p_1` -1 -200 200 1st order weight -``param(3)`` :math:`p_2` -1 -200 200 2nd order weight -============= ============= ========= ============= ============= ============================== +============== =============== ============= ============= ============= =================================== + Variable Symbol Start Value Lower bound Upper bound Description +============== =============== ============= ============= ============= =================================== +``param[0]`` :math:`p_0` 1 0 200 Intercept +``param[1]`` :math:`p_1` -1 -200 200 1st order weight (μs\ :sup:`-1`) +``param[2]`` :math:`p_2` -1 -200 200 2nd order weight (μs\ :sup:`-2`) +============== =============== ============= ============= ============= =================================== \ No newline at end of file diff --git a/docsrc/source/models/bg_poly3.rst b/docsrc/source/models/bg_poly3.rst index a9bf604c4..e2e3321d7 100644 --- a/docsrc/source/models/bg_poly3.rst +++ b/docsrc/source/models/bg_poly3.rst @@ -13,11 +13,11 @@ Model :math:`B(t) = p_0 + p_1(\lambda t) + p_2(\lambda t)^2 + p_3(\lambda t)^3` -============= ============= ========= ============= ============= ============================== - Variable Symbol Default Lower bound Upper bound Description -============= ============= ========= ============= ============= ============================== -``param(1)`` :math:`p_0` 1 0 200 Intercept -``param(2)`` :math:`p_1` -1 -200 200 1st order weight -``param(3)`` :math:`p_2` -1 -200 200 2nd order weight -``param(3)`` :math:`p_3` -1 -200 200 3rd order weight -============= ============= ========= ============= ============= ============================== +============== =============== ============= ============= ============= =================================== + Variable Symbol Start Value Lower bound Upper bound Description +============== =============== ============= ============= ============= =================================== +``param[0]`` :math:`p_0` 1 0 200 Intercept +``param[1]`` :math:`p_1` -1 -200 200 1st order weight (μs\ :sup:`-1`) +``param[2]`` :math:`p_2` -1 -200 200 2nd order weight (μs\ :sup:`-2`) +``param[3]`` :math:`p_3` -1 -200 200 3rd order weight (μs\ :sup:`-3`) +============== =============== ============= ============= ============= =================================== \ No newline at end of file diff --git a/docsrc/source/models/bg_prodstrexp.rst b/docsrc/source/models/bg_prodstrexp.rst index 7b3a4ad15..5e99a942e 100644 --- a/docsrc/source/models/bg_prodstrexp.rst +++ b/docsrc/source/models/bg_prodstrexp.rst @@ -13,11 +13,11 @@ Model :math:`B(t) = \exp\left(-\lambda\kappa_1 \vert t \vert^{d_1}\right) \exp\left(-\lambda\kappa_2 \vert t\vert^{d_2}\right)` -============= ================== ========= ============= ============= ============================== - Variable Symbol Default Lower bound Upper bound Description -============= ================== ========= ============= ============= ============================== -``param(1)`` :math:`\kappa_1` 3.5 0 200 1st strexp decay rate -``param(2)`` :math:`d_1` 1 0 6 1st strexp stretch factor -``param(3)`` :math:`\kappa_2` 3.5 0 200 2nd strexp decay rate -``param(4)`` :math:`d_2` 1 0 6 2nd strexp stretch factor -============= ================== ========= ============= ============= ============================== +============== ================= ============= ============= ============= ================================= + Variable Symbol Start Value Lower bound Upper bound Description +============== ================= ============= ============= ============= ================================= +``param[0]`` :math:`\kappa_1` 3.5 0 200 1st strexp decay rate +``param[1]`` :math:`d_1` 1 0 6 1st strexp stretch factor +``param[2]`` :math:`\kappa_2` 3.5 0 200 2nd strexp decay rate +``param[3]`` :math:`d_2` 1 0 6 2nd strexp stretch factor +============== ================= ============= ============= ============= ================================= diff --git a/docsrc/source/models/bg_strexp.rst b/docsrc/source/models/bg_strexp.rst index 1dc80bad1..78ba1768c 100644 --- a/docsrc/source/models/bg_strexp.rst +++ b/docsrc/source/models/bg_strexp.rst @@ -15,11 +15,11 @@ Model B(t) = \exp\left(-\lambda \kappa \vert t\vert^{d}\right) -============= ================= ========= ============= ============= ======================== - Variable Symbol Default Lower bound Upper bound Description -============= ================= ========= ============= ============= ======================== -``param(1)`` :math:`\kappa` 3.5 0 200 decay rate (1/us) -``param(2)`` :math:`d` 1 0 6 stretch factor -============= ================= ========= ============= ============= ======================== +============== ================= ============= ============= ============= ================================= + Variable Symbol Start Value Lower bound Upper bound Description +============== ================= ============= ============= ============= ================================= +``param[0]`` :math:`\kappa` 3.5 0 200 Decay rate (μs\ :sup:`-d`) +``param[1]`` :math:`d` 1 0 6 Stretch factor +============== ================= ============= ============= ============= ================================= Although the ``bg_strexp`` model has the same functional form as ``bg_homfractal``, it is distinct since its first parameter is a decay rate constant and not a spin concentration like for ``bg_homfractal``. diff --git a/docsrc/source/models/bg_sumstrexp.rst b/docsrc/source/models/bg_sumstrexp.rst index b077de8c3..959754153 100644 --- a/docsrc/source/models/bg_sumstrexp.rst +++ b/docsrc/source/models/bg_sumstrexp.rst @@ -12,12 +12,12 @@ Model :math:`B(t) = A_1\exp \left(-\lambda\kappa_1 \vert t \vert^{d_1}\right) + (1-A_1)\exp\left(-\lambda\kappa_2 \vert t \vert^{d_2}\right)` -============= ================= ========= ============= ============= ============================== - Variable Symbol Default Lower bound Upper bound Description -============= ================= ========= ============= ============= ============================== -``param(1)`` :math:`\kappa_1` 3.5 0 200 1st strexp decay rate -``param(2)`` :math:`d_1` 1 0 6 1st strexp stretch factor -``param(3)`` :math:`\kappa_2` 3.5 0 200 2nd strexp decay rate -``param(4)`` :math:`d_2` 1 0 6 2nd strexp stretch factor -``param(5)`` :math:`A_1` 0.5 0 1 Relative amplitude -============= ================= ========= ============= ============= ============================== +============== ================= ============= ============= ============= ======================================== + Variable Symbol Start Value Lower bound Upper bound Description +============== ================= ============= ============= ============= ======================================== +``param[0]`` :math:`\kappa_1` 3.5 0 200 1st strexp decay rate (μs\ :sup:`-d`) +``param[1]`` :math:`d_1` 1 0 6 1st strexp stretch factor +``param[2]`` :math:`\kappa_2` 3.5 0 200 2nd strexp decay rate (μs\ :sup:`-d`) +``param[3]`` :math:`d_2` 1 0 6 2nd strexp stretch factor +``param[4]`` :math:`A_1` 0.5 0 1 Relative amplitude +============== ================= ============= ============= ============= ======================================== \ No newline at end of file diff --git a/docsrc/source/models/dd_circle.rst b/docsrc/source/models/dd_circle.rst index 903a86b2e..cd994ee59 100644 --- a/docsrc/source/models/dd_circle.rst +++ b/docsrc/source/models/dd_circle.rst @@ -15,15 +15,15 @@ Model This provides a `semi-circle distribution `_, defined by :math:`P(r) = 2\pi\sqrt{(r-r_0)^2/R^2+1}` for :math:`r_0-R\le r\le r_0+R` and zero otherwise. -============== ======================== ========= ============= ============= ======================== - Variable Symbol Default Lower bound Upper bound Description -============== ======================== ========= ============= ============= ======================== -``param(1)`` :math:`r_0` 3.0 0.1 20 center -``param(2)`` :math:`R` 0.5 0.1 5 radius -============== ======================== ========= ============= ============= ======================== +============== ================= ============= ============= ============= ================================= + Variable Symbol Start Value Lower bound Upper bound Description +============== ================= ============= ============= ============= ================================= +``param[0]`` :math:`r_0` 3.0 0.1 20 Center (nm) +``param[1]`` :math:`R` 0.5 0.1 5 Radius (nm) +============== ================= ============= ============= ============= ================================= -Example using default parameters: +Example using the start values: .. image:: ../images/model_dd_circle.png :width: 650px diff --git a/docsrc/source/models/dd_cos.rst b/docsrc/source/models/dd_cos.rst index f2c3e9f4b..5e81cd7ff 100644 --- a/docsrc/source/models/dd_cos.rst +++ b/docsrc/source/models/dd_cos.rst @@ -16,15 +16,15 @@ Model This provides a `raised-cosine distribution `_, defined by :math:`P(r) = \frac{1}{2w}\cos\left(\frac{r-r_0}{w}\pi\right)` for :math:`r_0-w \le r \le r_0+w`, and zero otherwise. -============== ======================== ========= ============= ============= ======================== - Variable Symbol Default Lower bound Upper bound Description -============== ======================== ========= ============= ============= ======================== -``param(1)`` :math:`r_0` 3.0 0.1 20 center, in nm -``param(2)`` :math:`w` 0.5 0.1 5 fwhm, in nm -============== ======================== ========= ============= ============= ======================== +============== ================= ============= ============= ============= ================================= + Variable Symbol Start Value Lower bound Upper bound Description +============== ================= ============= ============= ============= ================================= +``param[0]`` :math:`r_0` 3.0 0.1 20 center, in nm +``param[1]`` :math:`w` 0.5 0.1 5 fwhm, in nm +============== ================= ============= ============= ============= ================================= -Example using default parameters: +Example using the start values: .. image:: ../images/model_dd_cos.png :width: 650px diff --git a/docsrc/source/models/dd_gauss.rst b/docsrc/source/models/dd_gauss.rst index 7674e263d..5442d211f 100644 --- a/docsrc/source/models/dd_gauss.rst +++ b/docsrc/source/models/dd_gauss.rst @@ -13,17 +13,15 @@ Model :math:`P(r) = \sqrt{\frac{2}{\pi}}\frac{1}{\sigma}\exp\left(-\frac{(r-\left)^2}{\sigma^2}\right)` -with :math:`\sigma = \mathrm{FWHM}/\sqrt{2ln(2)}` +============== ======================== ============= ============= ============= =========================== + Variable Symbol Start value Lower bound Upper bound Description +============== ======================== ============= ============= ============= =========================== +``param[0]`` :math:`\left` 3.5 1.0 20 Mean (nm) +``param[1]`` :math:`\sigma` 0.2 0.05 2.5 Standard deviation (nm) +============== ======================== ============= ============= ============= =========================== -============== ======================== ========= ============= ============= ======================== - Variable Symbol Default Lower bound Upper bound Description -============== ======================== ========= ============= ============= ======================== -``param(1)`` :math:`r0` 3.5 1.0 20 center distance (nm) -``param(2)`` :math:`\mathrm{FWHM}` 0.5 0.2 5 FWHM (nm) -============== ======================== ========= ============= ============= ======================== - -Example using default parameters: +Example using the start values: .. image:: ../images/model_dd_gauss.png :width: 650px diff --git a/docsrc/source/models/dd_gauss2.rst b/docsrc/source/models/dd_gauss2.rst index 29f23fef5..3ed447a11 100644 --- a/docsrc/source/models/dd_gauss2.rst +++ b/docsrc/source/models/dd_gauss2.rst @@ -13,21 +13,19 @@ Model :math:`P(r) = a_1\sqrt{\frac{2}{\pi}}\frac{1}{\sigma_1}\exp\left(-\frac{(r-\left)^2}{\sigma_1^2}\right) + a_2\sqrt{\frac{2}{\pi}}\frac{1}{\sigma_2}\exp\left(-\frac{(r-\left)^2}{\sigma_2^2}\right)` -with :math:`\sigma_i = \mathrm{FWHM}_i/\sqrt{2ln(2)}` - -============== ======================== ========= ======== ========= =================================== - Variable Symbol Default Lower Upper Description -============== ======================== ========= ======== ========= =================================== -``param(1)`` :math:`\left` 2.5 1.0 20 1st Gaussian center distance -``param(2)`` :math:`\mathrm{FWHM}_1` 0.5 0.2 5 1st Gaussian FWHM -``param(3)`` :math:`a_1` 0.5 0 1 1st Gaussian amplitude -``param(4)`` :math:`\left` 3.5 1.0 20 2nd Gaussian center distance -``param(5)`` :math:`\mathrm{FWHM}_2` 0.5 0.2 5 2nd Gaussian FWHM -``param(6)`` :math:`a_2` 0.5 0 1 2nd Gaussian amplitude -============== ======================== ========= ======== ========= =================================== - -Example using default parameters: +============== ========================= ============= ============= ============= ====================================== + Variable Symbol Start Value Lower bound Upper bound Description +============== ========================= ============= ============= ============= ====================================== +``param[0]`` :math:`\left` 2.5 1.0 20 1st Gaussian mean distance (nm) +``param[1]`` :math:`\sigma_1` 0.2 0.05 2.5 1st Gaussian standard deviation (nm) +``param[2]`` :math:`a_1` 0.5 0 1 1st Gaussian amplitude +``param[3]`` :math:`\left` 3.5 1.0 20 2nd Gaussian mean distance (nm) +``param[4]`` :math:`\sigma_2` 0.2 0.05 2.5 2nd Gaussian standard deviation (nm) +``param[5]`` :math:`a_2` 0.5 0 1 2nd Gaussian amplitude +============== ========================= ============= ============= ============= ====================================== + +Example using the start values: .. image:: ../images/model_dd_gauss2.png :width: 650px diff --git a/docsrc/source/models/dd_gauss3.rst b/docsrc/source/models/dd_gauss3.rst index f9dc1b3f0..9f2056388 100644 --- a/docsrc/source/models/dd_gauss3.rst +++ b/docsrc/source/models/dd_gauss3.rst @@ -13,24 +13,22 @@ Model :math:`P(r) = a_1\sqrt{\frac{2}{\pi}}\frac{1}{\sigma_1}\exp\left(-\frac{(r-\left)^2}{\sigma_1^2}\right) + a_2\sqrt{\frac{2}{\pi}}\frac{1}{\sigma_2}\exp\left(-\frac{(r-\left)^2}{\sigma_2^2}\right) + a_3\sqrt{\frac{2}{\pi}}\frac{1}{\sigma_3}\exp\left(-\frac{(r-\left)^2}{\sigma_3^2}\right)` -with :math:`\sigma_i = \mathrm{FWHM}_i/\sqrt{2ln(2)}` - -================ ======================== ========= ======== ========= =================================== - Variable Symbol Default Lower Upper Description -================ ======================== ========= ======== ========= =================================== -``param(1)`` :math:`\left` 2.5 1.0 20 1st Gaussian center distance -``param(2)`` :math:`\mathrm{FWHM}_1` 0.5 0.2 5 1st Gaussian FWHM -``param(3)`` :math:`a_1` 0.3 0 1 1st Gaussian amplitude -``param(4)`` :math:`\left` 3.5 1.0 20 2nd Gaussian center distance -``param(5)`` :math:`\mathrm{FWHM}_2` 0.5 0.2 5 2nd Gaussian FWHM -``param(6)`` :math:`a_2` 0.3 0 1 2nd Gaussian amplitude -``param(7)`` :math:`\left` 5.0 1.0 20 3rd Gaussian center distance -``param(8)`` :math:`\mathrm{FWHM}_3` 0.5 0.2 5 3rd Gaussian FWHM -``param(9)`` :math:`a_3` 0.3 0 1 3rd Gaussian amplitude -================ ======================== ========= ======== ========= =================================== - - -Example using default parameters: +============== ========================== ============= ============= ============= ======================================= + Variable Symbol Start Value Lower bound Upper bound Description +============== ========================== ============= ============= ============= ======================================= +``param[0]`` :math:`\left` 2.5 1.0 20 1st Gaussian mean distance (nm) +``param[1]`` :math:`\sigma_1` 0.2 0.05 2.5 1st Gaussian standard deviation (nm) +``param[2]`` :math:`a_1` 0.3 0 1 1st Gaussian amplitude +``param[3]`` :math:`\left` 3.5 1.0 20 2nd Gaussian mean distance (nm) +``param[4]`` :math:`\sigma_2` 0.2 0.05 2.5 2nd Gaussian standard deviation (nm) +``param[5]`` :math:`a_2` 0.3 0 1 2nd Gaussian amplitude +``param[6]`` :math:`\left` 5.0 1.0 20 3rd Gaussian mean distance (nm) +``param[7]`` :math:`\sigma_3` 0.2 0.05 2.5 3rd Gaussian standard deviation (nm) +``param[8]`` :math:`a_3` 0.3 0 1 3rd Gaussian amplitude +============== ========================== ============= ============= ============= ======================================= + + +Example using the start values: .. image:: ../images/model_dd_gauss3.png :width: 650px diff --git a/docsrc/source/models/dd_gengauss.rst b/docsrc/source/models/dd_gengauss.rst index dff3b471c..c0bac1a72 100644 --- a/docsrc/source/models/dd_gengauss.rst +++ b/docsrc/source/models/dd_gengauss.rst @@ -17,18 +17,16 @@ Model :math:`P(r) = \frac{\beta}{2\sigma\Gamma(1/\beta)}\exp\left(-\left(\frac{(r-\left)}{\sigma}\right)^\beta \right)` -with :math:`\sigma = w/(2\sqrt{2ln(2)})` +============== ========================== ============= ============= ============= ======================================= + Variable Symbol Start Value Lower bound Upper bound Description +============== ========================== ============= ============= ============= ======================================= +``param[0]`` :math:`r_0` 3.5 1.0 20 Mean (nm) +``param[1]`` :math:`\sigma` 0.2 0.05 2.5 Spread (nm) +``param[2]`` :math:`\beta` 5.0 0.25 15 kurtosis +============== ========================== ============= ============= ============= ======================================= -============== ======================== ========= ============= ============= ======================== - Variable Symbol Default Lower bound Upper bound Description -============== ======================== ========= ============= ============= ======================== -``param(1)`` :math:`r_0` 3.5 1.0 20 center distance (nm) -``param(2)`` :math:`w` 0.5 0.2 5 FWHM (nm) -``param(2)`` :math:`\beta` 5.0 0.25 15 kurtosis -============== ======================== ========= ============= ============= ======================== - -Example using default parameters: +Example using start values: .. image:: ../images/model_dd_gengauss.png :width: 650px diff --git a/docsrc/source/models/dd_randcoil.rst b/docsrc/source/models/dd_randcoil.rst index 36651ad1d..a3e310f29 100644 --- a/docsrc/source/models/dd_randcoil.rst +++ b/docsrc/source/models/dd_randcoil.rst @@ -18,13 +18,13 @@ Model where :math:`\nu_0 = 3/(12\pi r_0 N \nu)^{3/2}` -============== =========== ======== ======== ======== ================================== - Variable Symbol Default Lower Upper Description -============== =========== ======== ======== ======== ================================== -``param(1)`` :math:`N` 50 2 1000 Number of residues -``param(2)`` :math:`R_0` 0.20 0.10 0.40 Segment length -``param(3)`` :math:`\nu` 0.60 0.33 1.00 Scaling exponent -============== =========== ======== ======== ======== ================================== +============== ============= ============= ============= ============= ======================================= + Variable Symbol Start Value Lower bound Upper bound Description +============== ============= ============= ============= ============= ======================================= +``param[0]`` :math:`N` 50 2 1000 Number of residues +``param[1]`` :math:`R_0` 0.20 0.10 0.40 Segment length (nm) +``param[2]`` :math:`\nu` 0.60 0.33 1.00 Scaling exponent +============== ============= ============= ============= ============= ======================================= Example using default parameters: diff --git a/docsrc/source/models/dd_rice.rst b/docsrc/source/models/dd_rice.rst index e6b57f489..ef99308e9 100644 --- a/docsrc/source/models/dd_rice.rst +++ b/docsrc/source/models/dd_rice.rst @@ -17,15 +17,15 @@ Model where :math:`n=3` and :math:`I_{n/2-1}(x)` is the modified Bessel function of the first kind with order :math:`n/2-1`. This is a three-dimensional non-central chi distribution, the 3D generalization of the 2D Rice distribution. -============== ======================== ========= ============= ============= ======================== - Variable Symbol Default Lower bound Upper bound Description -============== ======================== ========= ============= ============= ======================== -``param(1)`` :math:`\nu` 3.5 1.0 10 Center -``param(2)`` :math:`\sigma` 0.7 0.1 5 Width -============== ======================== ========= ============= ============= ======================== +============== ======================== ============= ============= ============= ======================================= + Variable Symbol Start Value Lower bound Upper bound Description +============== ======================== ============= ============= ============= ======================================= +``param[0]`` :math:`\nu` 3.5 1.0 10 Location (nm) +``param[1]`` :math:`\sigma` 0.7 0.1 5 Spread (nm) +============== ======================== ============= ============= ============= ======================================= -Example using default parameters: +Example using start values: .. image:: ../images/model_dd_rice.png :width: 650px diff --git a/docsrc/source/models/dd_rice2.rst b/docsrc/source/models/dd_rice2.rst index dc93a09aa..665c37567 100644 --- a/docsrc/source/models/dd_rice2.rst +++ b/docsrc/source/models/dd_rice2.rst @@ -17,19 +17,19 @@ Model where :math:`n=3` and :math:`I_{n/2-1}(x)` is the modified Bessel function of the first kind with order :math:`n/2-1`. This is a three-dimensional non-central chi distribution, the 3D generalization of the 2D Rice distribution. -============== ======================== ========= ======== ======== =============================== - Variable Symbol Default Lower Upper Description -============== ======================== ========= ======== ======== =============================== -``param(1)`` :math:`\nu_1` 2.5 1.0 10 1st Rician center distance -``param(2)`` :math:`\sigma_1` 0.7 0.1 5 1st Rician width -``param(3)`` :math:`a_1` 0.5 0 1 1st Rician amplitude -``param(4)`` :math:`\nu_2` 4.0 1.0 10 2nd Rician center distance -``param(5)`` :math:`\sigma_2` 0.7 0.1 5 2nd Rician width -``param(6)`` :math:`a_2` 0.5 0 1 2nd Rician amplitude -============== ======================== ========= ======== ======== =============================== - - -Example using default parameters: +============== ======================== ============= ============= ============= ======================================= + Variable Symbol Start Value Lower bound Upper bound Description +============== ======================== ============= ============= ============= ======================================= +``param[0]`` :math:`\nu_1` 2.5 1.0 10 1st Rician location (nm) +``param[1]`` :math:`\sigma_1` 0.7 0.1 5 1st Rician spread (nm) +``param[2]`` :math:`a_1` 0.5 0 1 1st Rician amplitude +``param[3]`` :math:`\nu_2` 4.0 1.0 10 2nd Rician location (nm) +``param[4]`` :math:`\sigma_2` 0.7 0.1 5 2nd Rician spread (nm) +``param[5]`` :math:`a_2` 0.5 0 1 2nd Rician amplitude +============== ======================== ============= ============= ============= ======================================= + + +Example using start values: .. image:: ../images/model_dd_rice2.png :width: 650px diff --git a/docsrc/source/models/dd_rice3.rst b/docsrc/source/models/dd_rice3.rst index 2e633c7b5..903162622 100644 --- a/docsrc/source/models/dd_rice3.rst +++ b/docsrc/source/models/dd_rice3.rst @@ -20,22 +20,22 @@ where :math:`n=3` and :math:`I_{n/2-1}(x)` is the modified Bessel function of th This is a three-dimensional non-central chi distribution, the 3D generalization of the 2D Rice distribution. -============== ======================== ========= ======== ======== =============================== - Variable Symbol Default Lower Upper Description -============== ======================== ========= ======== ======== =============================== -``param(1)`` :math:`\nu_1` 2.5 1.0 10 1st Rician center distance -``param(2)`` :math:`\sigma_1` 0.7 0.1 5 1st Rician width -``param(3)`` :math:`a_1` 0.3 0 1 1st Rician amplitude -``param(4)`` :math:`\nu_2` 4.0 1.0 10 2nd Rician center distance -``param(5)`` :math:`\sigma_2` 0.7 0.1 5 2nd Rician width -``param(6)`` :math:`a_2` 0.3 0 1 2nd Rician amplitude -``param(7)`` :math:`\nu_3` 5.0 1.0 10 3rd Rician center distance -``param(8)`` :math:`\sigma_3` 0.7 0.1 5 3rd Rician width -``param(9)`` :math:`a_3` 0.3 0 1 3rd Rician amplitude -============== ======================== ========= ======== ======== =============================== - - -Example using default parameters: +============== ======================== ============= ============= ============= ======================================= + Variable Symbol Start Value Lower bound Upper bound Description +============== ======================== ============= ============= ============= ======================================= +``param[0]`` :math:`\nu_1` 2.5 1.0 10 1st Rician location (nm) +``param[1]`` :math:`\sigma_1` 0.7 0.1 5 1st Rician spread (nm) +``param[2]`` :math:`a_1` 0.3 0 1 1st Rician amplitude +``param[3]`` :math:`\nu_2` 4.0 1.0 10 2nd Rician location (nm) +``param[4]`` :math:`\sigma_2` 0.7 0.1 5 2nd Rician spread (nm) +``param[5]`` :math:`a_2` 0.3 0 1 2nd Rician amplitude +``param[6]`` :math:`\nu_3` 5.0 1.0 10 3rd Rician location (nm) +``param[7]`` :math:`\sigma_3` 0.7 0.1 5 3rd Rician spread (nm) +``param[8]`` :math:`a_3` 0.3 0 1 3rd Rician amplitude +============== ======================== ============= ============= ============= ======================================= + + +Example using start values: .. image:: ../images/model_dd_rice3.png :width: 650px diff --git a/docsrc/source/models/dd_shell.rst b/docsrc/source/models/dd_shell.rst index 25f07994a..7b094899f 100644 --- a/docsrc/source/models/dd_shell.rst +++ b/docsrc/source/models/dd_shell.rst @@ -29,15 +29,15 @@ and :math:`R_2 = R + w` -================ ============== ========= ======== ========= =================================== - Variable Symbol Default Lower Upper Description -================ ============== ========= ======== ========= =================================== -``param(1)`` :math:`R` 1.5 0.1 20 Inner shell radius -``param(2)`` :math:`w` 0.5 0.1 20 Shell thickness -================ ============== ========= ======== ========= =================================== +============== ============== ============= ============= ============= ======================================= + Variable Symbol Start Value Lower bound Upper bound Description +============== ============== ============= ============= ============= ======================================= +``param[0]`` :math:`R` 1.5 0.1 20 Inner shell radius (nm) +``param[1]`` :math:`w` 0.5 0.1 20 Shell thickness (nm) +============== ============== ============= ============= ============= ======================================= -Example using default parameters: +Example using start values: .. image:: ../images/model_dd_shell.png :width: 650px diff --git a/docsrc/source/models/dd_shellshell.rst b/docsrc/source/models/dd_shellshell.rst index 0dbeea92e..2a08ce747 100644 --- a/docsrc/source/models/dd_shellshell.rst +++ b/docsrc/source/models/dd_shellshell.rst @@ -31,16 +31,16 @@ and -================ ============== ========= ======== ========= =================================== - Variable Symbol Default Lower Upper Description -================ ============== ========= ======== ========= =================================== -``param(1)`` :math:`R` 1.5 0.1 20 Inner shell radius -``param(2)`` :math:`w_1` 0.5 0.1 20 1st Shell thickness -``param(3)`` :math:`w_2` 0.5 0.1 20 2nd Shell thickness -================ ============== ========= ======== ========= =================================== +============== ============== ============= ============= ============= ======================================= + Variable Symbol Start Value Lower bound Upper bound Description +============== ============== ============= ============= ============= ======================================= +``param[0]`` :math:`R` 1.5 0.1 20 Inner shell radius (nm) +``param[1]`` :math:`w_1` 0.5 0.1 20 1st Shell thickness (nm) +``param[2]`` :math:`w_2` 0.5 0.1 20 2nd Shell thickness (nm) +============== ============== ============= ============= ============= ======================================= -Example using default parameters: +Example using start values: .. image:: ../images/model_dd_shellshell.png :width: 650px diff --git a/docsrc/source/models/dd_shellsphere.rst b/docsrc/source/models/dd_shellsphere.rst index f16ceeb13..8b45ad3d5 100644 --- a/docsrc/source/models/dd_shellsphere.rst +++ b/docsrc/source/models/dd_shellsphere.rst @@ -23,15 +23,15 @@ with :math:`R_2 = R + w_1` -================ ============== ========= ======== ========= =================================== - Variable Symbol Default Lower Upper Description -================ ============== ========= ======== ========= =================================== -``param(1)`` :math:`R` 1.5 0.1 20 Sphere radius -``param(2)`` :math:`w` 0.5 0.1 20 Shell thickness -================ ============== ========= ======== ========= =================================== +============== ============== ============= ============= ============= ========================= + Variable Symbol Start Value Lower bound Upper bound Description +============== ============== ============= ============= ============= ========================= +``param[0]`` :math:`R` 1.5 0.1 20 Sphere radius (nm) +``param[1]`` :math:`w` 0.5 0.1 20 Shell thickness (nm) +============== ============== ============= ============= ============= ========================= -Example using default parameters: +Example using start values: .. image:: ../images/model_dd_shellsphere.png :width: 650px diff --git a/docsrc/source/models/dd_shellvoidshell.rst b/docsrc/source/models/dd_shellvoidshell.rst index ece3f9e3c..3d2d45a2f 100644 --- a/docsrc/source/models/dd_shellvoidshell.rst +++ b/docsrc/source/models/dd_shellvoidshell.rst @@ -31,17 +31,17 @@ and :math:`R_4 = R + w_1 + d + w_2` -================ ============== ========= ======== ========= =================================== - Variable Symbol Default Lower Upper Description -================ ============== ========= ======== ========= =================================== -``param(1)`` :math:`R` 0.75 0.1 20 Sphere radius -``param(2)`` :math:`w_1` 1.00 0.1 20 1st Shell thickness -``param(3)`` :math:`w_2` 1.00 0.1 20 2nd Shell thickness -``param(4)`` :math:`d` 0.50 0.1 20 Shell-Shell separation -================ ============== ========= ======== ========= =================================== +============== ============== ============= ============= ============= ============================= + Variable Symbol Start Value Lower bound Upper bound Description +============== ============== ============= ============= ============= ============================= +``param[0]`` :math:`R` 0.75 0.1 20 Sphere radius (nm) +``param[1]`` :math:`w_1` 1.00 0.1 20 1st Shell thickness (nm) +``param[2]`` :math:`w_2` 1.00 0.1 20 2nd Shell thickness (nm) +``param[3]`` :math:`d` 0.50 0.1 20 Shell-Shell separation (nm) +============== ============== ============= ============= ============= ============================= -Example using default parameters: +Example using start values: .. image:: ../images/model_dd_shellvoidshell.png :width: 650px diff --git a/docsrc/source/models/dd_shellvoidsphere.rst b/docsrc/source/models/dd_shellvoidsphere.rst index 3e5653c11..b38a1cca5 100644 --- a/docsrc/source/models/dd_shellvoidsphere.rst +++ b/docsrc/source/models/dd_shellvoidsphere.rst @@ -29,16 +29,16 @@ and :math:`R_3 = R + d + w` -================ ============== ========= ======== ========= =================================== - Variable Symbol Default Lower Upper Description -================ ============== ========= ======== ========= =================================== -``param(1)`` :math:`R` 1.5 0.1 20 Sphere radius -``param(2)`` :math:`w` 1.0 0.1 20 Shell thickness -``param(2)`` :math:`d` 0.5 0.1 20 Shell-Sphere separation -================ ============== ========= ======== ========= =================================== +============== ============== ============= ============= ============= ============================== + Variable Symbol Start Value Lower bound Upper bound Description +============== ============== ============= ============= ============= ============================== +``param[0]`` :math:`R` 1.5 0.1 20 Sphere radius (nm) +``param[1]`` :math:`w` 1.0 0.1 20 Shell thickness (nm) +``param[2]`` :math:`d` 0.5 0.1 20 Shell-Sphere separation (nm) +============== ============== ============= ============= ============= ============================== -Example using default parameters: +Example using start values: .. image:: ../images/model_dd_shellvoidsphere.png :width: 650px diff --git a/docsrc/source/models/dd_skewgauss.rst b/docsrc/source/models/dd_skewgauss.rst index 706c2271b..d629b1c3a 100644 --- a/docsrc/source/models/dd_skewgauss.rst +++ b/docsrc/source/models/dd_skewgauss.rst @@ -17,18 +17,16 @@ Model :math:`P(r) = \sqrt{\frac{2}{\pi}}\frac{1}{\sigma}\exp\left(-\frac{(r-\left)^2}{\sqrt(2)\sigma^2}\right)\frac{1}{2}\left(1 + erf\left(\frac{(r-\left)}{\sqrt{2}\sigma}\right) \right)` -with :math:`\sigma = w/(2\sqrt{2ln(2)})` +============== ============== ============= ============= ============= ========================= + Variable Symbol Start Value Lower bound Upper bound Description +============== ============== ============= ============= ============= ========================= +``param[0]`` :math:`r_0` 3.5 1.0 20 Location (nm) +``param[1]`` :math:`\sigma` 0.5 0.2 5 Spread (nm) +``param[2]`` :math:`\alpha` 5.0 -15 15 Skewness +============== ============== ============= ============= ============= ========================= -============== ======================== ========= ============= ============= ======================== - Variable Symbol Default Lower bound Upper bound Description -============== ======================== ========= ============= ============= ======================== -``param(1)`` :math:`r_0` 3.5 1.0 20 Center distance (nm) -``param(2)`` :math:`w` 0.5 0.2 5 FWHM (nm) -``param(2)`` :math:`\alpha` 5.0 -15 15 Skewness -============== ======================== ========= ============= ============= ======================== - -Example using default parameters: +Example using start values: .. image:: ../images/model_dd_skewgauss.png :width: 650px diff --git a/docsrc/source/models/dd_sphere.rst b/docsrc/source/models/dd_sphere.rst index f0c5f015e..5ee91928e 100644 --- a/docsrc/source/models/dd_sphere.rst +++ b/docsrc/source/models/dd_sphere.rst @@ -17,14 +17,14 @@ Model :math:`P(r) = \begin{cases} \frac{3r^5}{16R^6} - \frac{9r^3}{4R^4} + \frac{3r^2}{R^3} \quad \text{for} \quad 0 \leq r < 2R \\ 0 \quad \text{for} \quad \text{otherwise} \end{cases}` -================ ============== ========= ======== ========= =================================== - Variable Symbol Default Lower Upper Description -================ ============== ========= ======== ========= =================================== -``param(1)`` :math:`R` 2.5 0.1 20 Sphere radius -================ ============== ========= ======== ========= =================================== +============== ============== ============= ============= ============= ========================= + Variable Symbol Start Value Lower bound Upper bound Description +============== ============== ============= ============= ============= ========================= +``param[0]`` :math:`R` 2.5 0.1 20 Sphere radius (nm) +============== ============== ============= ============= ============= ========================= -Example using default parameters: +Example using start values: .. image:: ../images/model_dd_sphere.png :width: 650px diff --git a/docsrc/source/models/dd_spherepoint.rst b/docsrc/source/models/dd_spherepoint.rst index 132ecfc4a..069efd646 100644 --- a/docsrc/source/models/dd_spherepoint.rst +++ b/docsrc/source/models/dd_spherepoint.rst @@ -17,15 +17,15 @@ Model :math:`P(r) = \begin{cases} \frac{3r(R^2-(d-r)^2)}{4dR^3} \quad \text{for} \quad d-R \leq r < d+R \\ 0 \quad \text{for} \quad \text{otherwise} \end{cases}` -================ ============== ========= ======== ========= =================================== - Variable Symbol Default Lower Upper Description -================ ============== ========= ======== ========= =================================== -``param(1)`` :math:`R` 1.5 0.1 20 Sphere radius -``param(2)`` :math:`d` 3.5 0.1 20 Distance to point -================ ============== ========= ======== ========= =================================== +============== ============== ============= ============= ============= ========================= + Variable Symbol Start Value Lower bound Upper bound Description +============== ============== ============= ============= ============= ========================= +``param[0]`` :math:`R` 1.5 0.1 20 Sphere radius (nm) +``param[1]`` :math:`d` 3.5 0.1 20 Distance to point (nm) +============== ============== ============= ============= ============= ========================= -Example using default parameters: +Example using start values: .. image:: ../images/model_dd_spherepoint.png :width: 650px diff --git a/docsrc/source/models/dd_spheresurf.rst b/docsrc/source/models/dd_spheresurf.rst index 844ca02c4..04ffa56fb 100644 --- a/docsrc/source/models/dd_spheresurf.rst +++ b/docsrc/source/models/dd_spheresurf.rst @@ -18,14 +18,14 @@ Model :math:`P(r) = \begin{cases} \frac{r}{2R^2} \quad \text{for} \quad 0 \leq r < 2R \\ 0 \quad \text{for} \quad \text{otherwise} \end{cases}` -================ ============== ========= ======== ========= =================================== - Variable Symbol Default Lower Upper Description -================ ============== ========= ======== ========= =================================== -``param(1)`` :math:`R` 2.5 0.1 20 Sphere radius -================ ============== ========= ======== ========= =================================== +============== ============== ============= ============= ============= ========================= + Variable Symbol Start Value Lower bound Upper bound Description +============== ============== ============= ============= ============= ========================= +``param[0]` :math:`R` 2.5 0.1 20 Sphere radius (nm) +============== ============== ============= ============= ============= ========================= -Example using default parameters: +Example using start values: .. image:: ../images/model_dd_spheresurf.png :width: 650px diff --git a/docsrc/source/models/dd_triangle.rst b/docsrc/source/models/dd_triangle.rst index 10af8ad52..8cb14af33 100644 --- a/docsrc/source/models/dd_triangle.rst +++ b/docsrc/source/models/dd_triangle.rst @@ -14,16 +14,16 @@ Model This provides a simple triangular distribution. -============== ======================== ========= ============= ============= ======================== - Variable Symbol Default Lower bound Upper bound Description -============== ======================== ========= ============= ============= ======================== -``param(1)`` :math:`r_0` 3.5 1.0 20 mode -``param(2)`` :math:`w_\mathrm{L}` 0.3 0.1 5 left width -``param(3)`` :math:`w_\mathrm{R}` 0.3 0.1 5 right width -============== ======================== ========= ============= ============= ======================== +============== ======================== ============= ============= ============= ========================= + Variable Symbol Start Value Lower bound Upper bound Description +============== ======================== ============= ============= ============= ========================= +``param[0]`` :math:`r_0` 3.5 1.0 20 Mode (nm) +``param[1]`` :math:`w_\mathrm{L}` 0.3 0.1 5 Left width (nm) +``param[2]`` :math:`w_\mathrm{R}` 0.3 0.1 5 Right width (nm) +============== ======================== ============= ============= ============= ========================= -Example using default parameters: +Example using start values: .. image:: ../images/model_dd_triangle.png :width: 650px diff --git a/docsrc/source/models/dd_uniform.rst b/docsrc/source/models/dd_uniform.rst index 81a895417..7784da4de 100644 --- a/docsrc/source/models/dd_uniform.rst +++ b/docsrc/source/models/dd_uniform.rst @@ -15,15 +15,15 @@ Model This provides a simple uniform distribution. -============== ======================== ========= ============= ============= ======================== - Variable Symbol Default Lower bound Upper bound Description -============== ======================== ========= ============= ============= ======================== -``param(1)`` :math:`r_\mathrm{L}` 2.5 0.1 6 left edge -``param(2)`` :math:`r_\mathrm{R}` 3.0 0.2 20 right edge -============== ======================== ========= ============= ============= ======================== +============== ======================== ============= ============= ============= ========================= + Variable Symbol Start Value Lower bound Upper bound Description +============== ======================== ============= ============= ============= ========================= +``param[0]`` :math:`r_\mathrm{L}` 2.5 0.1 6 Left edge (nm) +``param[1]`` :math:`r_\mathrm{R}` 3.0 0.2 20 Right edge (nm) +============== ======================== ============= ============= ============= ========================= -Example using default parameters: +Example using start values: .. image:: ../images/model_dd_uniform.png :width: 650px diff --git a/docsrc/source/models/dd_wormchain.rst b/docsrc/source/models/dd_wormchain.rst index 920c50328..38db9f793 100644 --- a/docsrc/source/models/dd_wormchain.rst +++ b/docsrc/source/models/dd_wormchain.rst @@ -11,14 +11,14 @@ Model ========================================= -============== =========== ======== ======== ======== =============================== - Variable Symbol Default Lower Upper Description -============== =========== ======== ======== ======== =============================== -``param(1)`` :math:`L` 3.7 1.5 10 Contour length -``param(2)`` :math:`L_p` 10 2 100 Persistence length -============== =========== ======== ======== ======== =============================== - -Example: +============== ============ ============= ============= ============= ========================= + Variable Symbol Start Value Lower bound Upper bound Description +============== ============ ============= ============= ============= ========================= +``param[0]`` :math:`L` 3.7 1.5 10 Contour length (nm) +``param[1]`` :math:`L_p` 10 2 100 Persistence length (nm) +============== ============ ============= ============= ============= ========================= + +Example using start values: .. image:: ../images/model_dd_wormchain.png :width: 650px diff --git a/docsrc/source/models/dd_wormgauss.rst b/docsrc/source/models/dd_wormgauss.rst index 67443a63b..743d56aa5 100644 --- a/docsrc/source/models/dd_wormgauss.rst +++ b/docsrc/source/models/dd_wormgauss.rst @@ -12,15 +12,15 @@ Model ========================================= -============== =============== ======== ======== ======== =============================== - Variable Symbol Default Lower Upper Description -============== =============== ======== ======== ======== =============================== -``param(1)`` :math:`L` 3.7 1.5 10 Contour length -``param(2)`` :math:`L_p` 10 2 100 Persistence length -``param(3)`` :math:`\sigma` 0.2 0.01 2 Gaussian standard deviation -============== =============== ======== ======== ======== =============================== - -Example using default parameters: +============== ============== ============= ============= ============= ================================== + Variable Symbol Start Value Lower bound Upper bound Description +============== ============== ============= ============= ============= ================================== +``param[0]`` :math:`L` 3.7 1.5 10 Contour length (nm) +``param[1]`` :math:`L_p` 10 2 100 Persistence length (nm) +``param[2]`` :math:`\sigma` 0.2 0.01 2 Gaussian standard deviation (nm) +============== ============== ============= ============= ============= ================================== + +Example using start values: .. image:: ../images/model_dd_wormgauss.png :width: 650px diff --git a/docsrc/source/models/ex_4pdeer.rst b/docsrc/source/models/ex_4pdeer.rst index 0b19d51c8..b6c49029f 100644 --- a/docsrc/source/models/ex_4pdeer.rst +++ b/docsrc/source/models/ex_4pdeer.rst @@ -24,14 +24,14 @@ This experiment model has one modulated pathway and an unmodulated contribution. where :math:`T_0^{(1)}=0` is the refocusing time of the modulated dipolar pathway. -============== ================ ============ ============ ============ ================================================ - Variable Symbol Default Lower Upper Description -============== ================ ============ ============ ============ ================================================ -``param(1)`` :math:`\lambda` 0.3 0 1 modulated pathway amplitude (modulation depth) -============== ================ ============ ============ ============ ================================================ +============== ================ ============= ============ ============ ================================================ + Variable Symbol Start Values Lower Upper Description +============== ================ ============= ============ ============ ================================================ +``param[0]`` :math:`\lambda` 0.3 0 1 Modulated pathway amplitude (modulation depth) +============== ================ ============= ============ ============ ================================================ -Example of a simulated signal using default parameters: +Example of a simulated signal using start values: .. image:: ../images/model_ex_4pdeer.png :width: 550px diff --git a/docsrc/source/models/ex_5pdeer.rst b/docsrc/source/models/ex_5pdeer.rst index 8209b6d23..a3f64bd66 100644 --- a/docsrc/source/models/ex_5pdeer.rst +++ b/docsrc/source/models/ex_5pdeer.rst @@ -25,17 +25,17 @@ This experiment model has two modulated pathway and an unmodulated contribution. where :math:`T_0^{(1)}=0` and :math:`T_0^{(2)}` are the refocusing times of the two modulated dipolar pathways. -============== ======================== ================= ==================== ==================== ============================================= - Variable Symbol Default Lower Upper Description -============== ======================== ================= ==================== ==================== ============================================= -``param(1)`` :math:`\varLambda_0` 0.4 0 1 unmodulated pathways, amplitude -``param(2)`` :math:`\lambda_1` 0.4 0 1 1st modulated pathway, amplitude -``param(3)`` :math:`\lambda_2` 0.2 0 1 2nd modulated pathway, amplitude -``param(4)`` :math:`T_0^{(2)}` 5.0 0 20 2nd modulated pathway, refocusing time (us) -============== ======================== ================= ==================== ==================== ============================================= +============== ======================== ============= ============ ============ ================================================ + Variable Symbol Start Values Lower Upper Description +============== ======================== ============= ============ ============ ================================================ +``param[0]`` :math:`\varLambda_0` 0.4 0 1 Unmodulated pathways, amplitude +``param[1]`` :math:`\lambda_1` 0.4 0 1 1st modulated pathway, amplitude +``param[2]`` :math:`\lambda_2` 0.2 0 1 2nd modulated pathway, amplitude +``param[3]`` :math:`T_0^{(2)}` 5.0 0 20 2nd modulated pathway, refocusing time (μs) +============== ======================== ============= ============ ============ ================================================ -Example of a simulated signal using default parameters: +Example of a simulated signal using start values: .. image:: ../images/model_ex_5pdeer.png :width: 550px diff --git a/docsrc/source/models/ex_7pdeer.rst b/docsrc/source/models/ex_7pdeer.rst index 6a9909f37..89e30e33a 100644 --- a/docsrc/source/models/ex_7pdeer.rst +++ b/docsrc/source/models/ex_7pdeer.rst @@ -25,19 +25,19 @@ In order to reduce the parameter space, only the dipolar pathways refocusing at where :math:`T_0^{(1)}=0\;\mu s`, :math:`T_0^{(2)}`, and :math:`T_0^{(3)}` are the refocusing times of the three modulated dipolar pathways at positive evolution times. -============== ======================== ================= ==================== ==================== ============================================= - Variable Symbol Default Lower Upper Description -============== ======================== ================= ==================== ==================== ============================================= -``param(1)`` :math:`\varLambda_0` 0.4 0 1 unmodulated pathways, amplitude -``param(2)`` :math:`\lambda_1` 0.4 0 1 1st modulated pathway, amplitude -``param(3)`` :math:`\lambda_2` 0.2 0 1 2nd modulated pathway, amplitude -``param(4)`` :math:`\lambda_3` 0.2 0 1 3rd modulated pathway, amplitude -``param(5)`` :math:`T_0^{(2)}` 1.5 0 20 2nd modulated pathway, refocusing time (us) -``param(6)`` :math:`T_0^{(3)}` 3.5 0 20 3rd modulated pathway, refocusing time (us) -============== ======================== ================= ==================== ==================== ============================================= - - -Example of a simulated signal using default parameters: +============== ======================== ============= ============ ============ ================================================ + Variable Symbol Start Values Lower Upper Description +============== ======================== ============= ============ ============ ================================================ +``param[0]`` :math:`\varLambda_0` 0.4 0 1 Unmodulated pathways, amplitude +``param[1]`` :math:`\lambda_1` 0.4 0 1 1st modulated pathway, amplitude +``param[2]`` :math:`\lambda_2` 0.2 0 1 2nd modulated pathway, amplitude +``param[3]`` :math:`\lambda_3` 0.2 0 1 3rd modulated pathway, amplitude +``param[4]`` :math:`T_0^{(2)}` 1.5 0 20 2nd modulated pathway, refocusing time (μs) +``param[5]`` :math:`T_0^{(3)}` 3.5 0 20 3rd modulated pathway, refocusing time (μs) +============== ======================== ============= ============ ============ ================================================ + + +Example of a simulated signal using start values: .. image:: ../images/model_ex_7pdeer.png :width: 550px diff --git a/docsrc/source/models/ex_ovl4pdeer.rst b/docsrc/source/models/ex_ovl4pdeer.rst index c592475d0..e6d505adb 100644 --- a/docsrc/source/models/ex_ovl4pdeer.rst +++ b/docsrc/source/models/ex_ovl4pdeer.rst @@ -27,17 +27,17 @@ This experiment model has two modulated pathways and an unmodulated contribution where :math:`T_0^{(1)}=0` and :math:`T_0^{(2)}` are the refocusing times of the two modulated dipolar pathways. -============== ======================== ================= ==================== ==================== ============================================== - Variable Symbol Default Lower Upper Description -============== ======================== ================= ==================== ==================== ============================================== -``param(1)`` :math:`\varLambda_0` 0.1 0 1 unmodulated pathways, amplitude -``param(2)`` :math:`\lambda_1` 0.8 0 1 1st modulated pathway, amplitude -``param(3)`` :math:`\lambda_2` 0.1 0 1 2nd modulated pathway, amplitude -``param(4)`` :math:`T_0^{(2)}` 5.0 0 20 2nd modulated pathway, refocusing time (us) -============== ======================== ================= ==================== ==================== ============================================== +============== ======================== ============= ============ ============ ================================================ + Variable Symbol Start Values Lower Upper Description +============== ======================== ============= ============ ============ ================================================ +``param[0]`` :math:`\varLambda_0` 0.1 0 1 Unmodulated pathways, amplitude +``param[1]`` :math:`\lambda_1` 0.8 0 1 1st modulated pathway, amplitude +``param[2]`` :math:`\lambda_2` 0.1 0 1 2nd modulated pathway, amplitude +``param[3]`` :math:`T_0^{(2)}` 5.0 0 20 2nd modulated pathway, refocusing time (μs) +============== ======================== ============= ============ ============ ================================================ -Example of a simulated signal using default parameters: +Example of a simulated signal using start values: .. image:: ../images/model_ex_ovl4pdeer.png :width: 550px diff --git a/docsrc/source/uncertainty.rst b/docsrc/source/uncertainty.rst index b54db9671..3df8ccbeb 100644 --- a/docsrc/source/uncertainty.rst +++ b/docsrc/source/uncertainty.rst @@ -102,7 +102,7 @@ To plot the resulting 95% and 50% confidence interval for the non-parametric dis import matplotlib.pyplot as plt plt.plot(r,Pfit,'k') - plt.fill_between(r,Pci50[:,0]; Pci50[:,1],color='r',alpha=,0.5) + plt.fill_between(r,Pci50[:,0]; Pci50[:,1],color='r',alpha=0.5) plt.fill_between(r,Pci95[:,0]; Pci95[:,1],color='r',alpha=0.2) Assumptions: diff --git a/examples/plot_analyzing_pake_pattern.py b/examples/plot_analyzing_pake_pattern.py index 893fba997..4b7f29436 100644 --- a/examples/plot_analyzing_pake_pattern.py +++ b/examples/plot_analyzing_pake_pattern.py @@ -3,12 +3,12 @@ Analyzing the Pake pattern of a dipolar signal ============================================================================ -A very basic example for displaying the Pake pattern of a given dipolar signal. +A very basic example for displaying the dipolar spectrum (Pake pattern) of a given dipolar signal. """ # %% import numpy as np import matplotlib.pyplot as plt -from deerlab import * +import deerlab as dl # %% [markdown] # Generate a dipolar signal @@ -17,14 +17,13 @@ # %% # Prepare components -t = np.linspace(0,5,400) -r = np.linspace(2,5,100) -P = dd_gauss2(r,[3.5, 0.3, 0.2, 4, 0.2, 0.8]) -B = bg_exp(t,0.2) +t = np.linspace(0,5,400) # µs +r = np.linspace(2,5,100) # nm +P = dl.dd_gauss2(r,[3.5, 0.1, 0.2, 4, 0.05, 0.8]) +B = dl.bg_exp(t,0.2) lam = 0.3 -K = dipolarkernel(t,r,lam,B) -np.random.seed(0) -V = K@P + whitegaussnoise(t,0.005) +K = dl.dipolarkernel(t,r,lam,B) +V = K@P + dl.whitegaussnoise(t,0.005,seed=0) # Plot plt.plot(t,V,'k.') @@ -46,22 +45,22 @@ # %% -tstart = 3 # Time to start fitting background, in us +tstart = 3 # Time to start fitting background, in µs mask = t>tstart # Model for the background component (1-lambda)*B def Bmodel(par): lam,kappa = par - B = (1 - lam)*bg_exp(t[mask],kappa) + B = (1 - lam)*dl.bg_exp(t[mask],kappa) return B # Fit the background function -fit = fitsignal(V,t,r,'P',bg_exp,ex_4pdeer,uqanalysis=False) +fit = dl.fitsignal(V,t,r,'P',dl.bg_exp,dl.ex_4pdeer,uqanalysis=False) Bfit = fit.B lam = fit.exparam kappa = fit.bgparam # %% [markdown] -# Now we can use these fitted variables to isolate the dipolar evolution function +# Now we can use these fitted parameters to isolate the dipolar evolution function # from the primary data. Removal of the background via division leads to a noise # increase at later times and thus to an approximation ``Vcorr`` of the real dipolar # evolution function. @@ -78,16 +77,16 @@ def Bmodel(par): plt.ylabel('V(t)') # %% [markdown] -# Computing the Pake pattern -# --------------------------- +# Computing the dipolar spectrum +# -------------------------------- # # Now that the signal has the appropiate structure for Fourier transform it, -# we can call the ``fftspec`` function to obtained the Pake pattern. +# we can call the ``fftspec`` function to obtain the dipolar spectrum. # %% # Compute spectrum -nu,pake = fftspec(Vcorr,t,apodization=False) +nu,pake = dl.fftspec(Vcorr,t,apodization=False) # %% [markdown] # In order to avoid truncation ripples in the Fourier spectrum and at the same @@ -96,7 +95,7 @@ def Bmodel(par): # %% # Compute spectrum with apodization -nuapo,pakeapo = fftspec(Vcorr,t,apodization=False,mode='real') +nuapo,pakeapo = dl.fftspec(Vcorr,t,apodization=False,mode='real') # Plot results plt.plot(nu,pake,'k',nuapo,pakeapo,'b',linewidth=1.5) @@ -114,3 +113,5 @@ def Bmodel(par): # The improvement will only be visual as no further information can be gained # from additional zero-filling. + +# %% diff --git a/examples/plot_bootstrapped_parameter_distributions.py b/examples/plot_bootstrapped_parameter_distributions.py index cdc4450c1..4c886a2ef 100644 --- a/examples/plot_bootstrapped_parameter_distributions.py +++ b/examples/plot_bootstrapped_parameter_distributions.py @@ -4,28 +4,27 @@ ============================================ This example shows how to generate probability density functions of -values for fit parameters using bootstrapping, showcased for 5pDEER. +values for fit parameters using bootstrapping, using 5pDEER as an example. """ import numpy as np import matplotlib.pyplot as plt -from deerlab import * +import deerlab as dl # %% [markdown] # Generate data # ------------- -t = np.linspace(-0.1,6.5,100) # time axis, us -r = np.linspace(1.5,6,100) # distance axis, ns -param0 = [3, 0.3, 0.2, 3.5, 0.3, 0.65, 3.8, 0.2, 0.15] # parameters for three-Gaussian model -P = dd_gauss3(r,param0) # model distance distribution -B = lambda t,lam: bg_hom3d(t,300,lam) # background decay +t = np.linspace(-0.1,6.5,100) # time axis, µs +r = np.linspace(1.5,6,100) # distance axis, nm +param0 = [3, 0.1, 0.2, 3.5, 0.1, 0.65, 3.8, 0.05, 0.15] # parameters for three-Gaussian model +P = dl.dd_gauss3(r,param0) # model distance distribution +B = lambda t,lam: dl.bg_hom3d(t,300,lam) # background decay exparam = [0.6, 0.3, 0.1, 3.2] # parameters for 5pDEER experiment -pathinfo = ex_5pdeer(exparam) # pathways information +pathinfo = dl.ex_5pdeer(exparam) # pathways information -np.random.seed(0) -K = dipolarkernel(t,r,pathinfo,B) -Vexp = K@P + whitegaussnoise(t,0.01) +K = dl.dipolarkernel(t,r,pathinfo,B) +Vexp = K@P + dl.whitegaussnoise(t,0.01,seed=0) # %% [markdown] # Analysis @@ -44,7 +43,7 @@ def fitroutine(V): # When running the fit, since we are only interested in the parameters we'll ignore # the rest (otherwise the ``Bfit``,``Pfit``,etc. could be bootstrapped as well) # We need the Vfit to pass it to bootan as well, so we'll request that one too. - fit = fitsignal(V,t,r,'P',bg_hom3d,ex_5pdeer,par0,lb,ub,uqanalysis=False) + fit = dl.fitsignal(V,t,r,'P',dl.bg_hom3d,dl.ex_5pdeer,par0,lb,ub,uqanalysis=False) Vfit = fit.V exparam = fit.exparam exparam[0:3] /=sum(exparam[0:3]) @@ -57,7 +56,7 @@ def fitroutine(V): exparfit,bgparfit,Vfit = fitroutine(Vexp) # Bootstrapping with 100 samples -bootuq = bootan(fitroutine,Vexp,Vfit,100) +bootuq = dl.bootan(fitroutine,Vexp,Vfit,100) # Extract the uncertainty quantification for the parameters exparam_uq = bootuq[0] @@ -67,7 +66,7 @@ def fitroutine(V): Lam0_values,Lam0_pdf = exparam_uq.pardist(0) lam1_values,lam1_pdf = exparam_uq.pardist(1) lam2_values,lam2_pdf = exparam_uq.pardist(2) -T02_values,T02_pdf = exparam_uq.pardist(3) +T02_values,T02_pdf = exparam_uq.pardist(3) # Extract distributions for the background parameters conc_values,conc_pdf = bgparam_uq.pardist(0) @@ -114,6 +113,4 @@ def fitroutine(V): plt.xlabel('Spin conc. [$\mu M$]') plt.ylabel('PDF') - - # %% diff --git a/examples/plot_comparing_uncertainties_for_regularization.py b/examples/plot_comparing_uncertainties_for_regularization.py index c6867d092..8986488bd 100644 --- a/examples/plot_comparing_uncertainties_for_regularization.py +++ b/examples/plot_comparing_uncertainties_for_regularization.py @@ -10,7 +10,7 @@ import numpy as np import matplotlib.pyplot as plt -from deerlab import * +import deerlab as dl # %% [markdown] # Simulate the data @@ -18,21 +18,13 @@ # # Let's start by generating some data. - -# Prepare signal components -t = np.linspace(-0.4,3.5,200) - -# Use a distance-axis with less points to make analysis faster -r = np.linspace(2,5,200) - -P = dd_gauss2(r,[3, 0.3, 0.6, 3.5, 0.4, 0.4]) -B = bg_strexp(t,[0.04,1]) -lam = 0.32 - -# Generate the dipolar kernel -K = dipolarkernel(t,r,lam,B) -# Simulate signal -V = K@P + whitegaussnoise(t,0.01) +t = np.linspace(-0.4,3.5,160) # time axis, µs +r = np.linspace(2,5,160) # distance axis, nm +P = dl.dd_gauss2(r,[3, 0.1, 0.6, 3.5, 0.2, 0.4]) # model distribution +lam = 0.32 # modulatio depth +B = dl.bg_strexp(t,[0.04,1],lam) # background decay +K = dl.dipolarkernel(t,r,lam,B) # dipolar kernel matrix +V = K@P + dl.whitegaussnoise(t,0.01) # signal with added noise # %% [markdown] # For the sake of simplicity, in this examples we will assume that we know the @@ -41,60 +33,57 @@ # Covariance-based confidence intervals # ------------------------------------- # -# We now have all the elements required to fit our distance distribution via -# regularization. We will use the AIC to select the regularization parameter in -# the Tikhonov regularization. +# Fit a Tikhonov model to the data, using AIC to select the regularization parameter +fit = dl.fitregmodel(V,K,r,'tikhonov','aic') +Pfit = fit.P # fitted distribution +Vfit = fit.V # fitted DEER trace -# Fit data via regularization -fit = fitregmodel(V,K,r,'tikhonov','aic') -Pfit = fit.P -Puq = fit.uncertainty -# Obtain time-domain fit -Vfit = K@Pfit - -plt.plot(r,P,'k',r,Pfit,'r',linewidth=1) -Pci95 = Puq.ci(95) -Pci50 = Puq.ci(50) -plt.fill_between(r,Pci50[:,0],Pci50[:,1],color='r',linestyle='None',alpha=0.45) -plt.fill_between(r,Pci95[:,0],Pci95[:,1],color='r',linestyle='None',alpha=0.25) -plt.grid() -plt.xlabel('r [nm]') -plt.ylabel('P(r) [nm$^{-1}$]') -plt.title('Curvature Matrix CI') -plt.legend(['Truth','Fit','50%-CI','95%-CI']) +# curvature matrix confidence intervals for distribution +Pci95_cm = fit.uncertainty.ci(95) +Pci50_cm = fit.uncertainty.ci(50) # %% [markdown] # Bootstrapped confidence intervals # --------------------------------- # -# Now we are interested in the bootstrap confidence intervals. For this, we -# need to define a boot function e.g. ``mybootfcn()`` which takes a signal as -# output and returns the outputs of interest (``Pfit`` in our example). - +# Now we calculate the bootstrap confidence intervals. For this, we +# need to define a function that takes a signal as input and returns +# the outputs of interest (``Pfit`` in our example). def mybootfcn(V): - fit = fitregmodel(V,K,r,'tikhonov','aic') + fit = dl.fitregmodel(V,K,r,'tikhonov','aic') return fit.P # Launch bootstrapping Nsamples = 50 -booci = bootan(mybootfcn,V,Vfit,Nsamples) -Pci95 = booci.ci(95) -Pci50 = booci.ci(50) +booci = dl.bootan(mybootfcn,V,Vfit,Nsamples) +Pci95_bs = booci.ci(95) +Pci50_bs = booci.ci(50) # %% [markdown] # By plotting the results, one can see that the bootstrapped confidence intervals # are narrower in comparison to the ones obtained via the curvature -# matrices. This is due to the inherent accurate nature of bootstrapping. - -plt.plot(r,P,'k',r,Pfit,'b',linewidth=1) -plt.fill_between(r,Pci50[:,0],Pci50[:,1],color='b',linestyle='None',alpha=0.45) -plt.fill_between(r,Pci95[:,0],Pci95[:,1],color='b',linestyle='None',alpha=0.25) -plt.grid(alpha=0.3) -plt.xlabel('r [nm]') -plt.ylabel('P(r) [nm$^{-1}$]') -plt.title('Bootstrapped CI') -plt.legend(['Truth','Fit','50%-CI','95%-CI']) - +# matrices. This is because bootstrapping takes the nonnegativity constraint of P into +# account, whereas the curvature matrix CIs do not. + +fig, ax = plt.subplots(1,2,sharey=True) +ax[0].plot(r,P,'k',r,Pfit,'r',linewidth=1) +ax[0].fill_between(r,Pci50_cm[:,0],Pci50_cm[:,1],color='r',linestyle='None',alpha=0.45) +ax[0].fill_between(r,Pci95_cm[:,0],Pci95_cm[:,1],color='r',linestyle='None',alpha=0.25) + +ax[1].plot(r,P,'k',r,Pfit,'b',linewidth=1) +ax[1].fill_between(r,Pci50_bs[:,0],Pci50_bs[:,1],color='b',linestyle='None',alpha=0.45) +ax[1].fill_between(r,Pci95_bs[:,0],Pci95_bs[:,1],color='b',linestyle='None',alpha=0.25) + +ax[0].grid(alpha=0.5) +ax[0].set_xlabel('r [nm]') +ax[0].set_ylabel('P [nm$^{-1}$]') +ax[0].set_title('Curvature Matrix CI') +ax[0].legend(['Truth','Fit','50%-CI','95%-CI']) + +ax[1].grid(alpha=0.5) +ax[1].set_xlabel('r [nm]') +ax[1].set_title('Bootstrapped CI') +ax[1].legend(['Truth','Fit','50%-CI','95%-CI']) # %% diff --git a/examples/plot_emulating_deeranalysis_workflow.py b/examples/plot_emulating_deeranalysis_workflow.py index c2f07a2ec..aa87f0e3c 100644 --- a/examples/plot_emulating_deeranalysis_workflow.py +++ b/examples/plot_emulating_deeranalysis_workflow.py @@ -10,7 +10,7 @@ import numpy as np import matplotlib.pyplot as plt -from deerlab import * +import deerlab as dl # %% [markdown] # Generating a dataset @@ -21,23 +21,21 @@ # Parameters t = np.linspace(-0.1,3,250) rtrue = np.linspace(1,7,200) -Ptrue = dd_gauss3(rtrue,[4.5, 0.6, 0.4, 3, 0.4, 0.3, 4, 0.7, 0.5]) +Ptrue = dl.dd_gauss3(rtrue,[4.5, 0.35, 0.4, 3, 0.25, 0.3, 4, 0.4, 0.5]) lam = 0.3 conc = 180 #uM # Simulate an experimental signal with some a.u. and phase offset -Bmodel = lambda t, lam: bg_hom3d(t,conc,lam) -K = dipolarkernel(t,rtrue,lam,Bmodel) +Bmodel = lambda t, lam: dl.bg_hom3d(t,conc,lam) +K = dl.dipolarkernel(t,rtrue,lam,Bmodel) V = K@Ptrue*np.exp(1j*np.pi/16) # add a phase shift -np.random.seed(1) -rnoise = whitegaussnoise(t,0.01) # real-component noise -np.random.seed(2) -inoise = 1j*whitegaussnoise(t,0.01) # imaginary-component noise +rnoise = dl.whitegaussnoise(t,0.01,seed=1) # real-component noise +inoise = 1j*dl.whitegaussnoise(t,0.01,see=2) # imaginary-component noise V = V + rnoise + inoise # complex-valued noisy signal V = V*3e6 # add an arbitrary amplitude scale plt.plot(t,V.real,'.',t,V.imag,'.'), -plt.xlabel('t [$\mu s$]') +plt.xlabel('t [µs]') plt.ylabel('V(t)') plt.grid(alpha=0.3) plt.legend(['real','imag']) @@ -48,36 +46,36 @@ # # Pre-processing -V = correctphase(V) -t = correctzerotime(V,t) +V = dl.correctphase(V) +t = dl.correctzerotime(V,t) V = V/max(V) # Distance axis estimation -r = time2dist(t) +r = dl.time2dist(t) # Background fit -tstart = 1.0 # background fit start, in us +tstart = 1.0 # background fit start, in µs mask = t>tstart def Bmodel(par): lam,kappa,d = par # unpack parameters - B = (1-lam)*bg_strexp(t[mask],[kappa,d],lam) + B = (1-lam)*dl.bg_strexp(t[mask],[kappa,d],lam) return B # lam k d par0 = [0.5, 0.5, 3] lb = [0.1, 0.1, 1] ub = [1, 5, 6] -fit = fitparamodel(V[mask],Bmodel,par0,lb,ub,rescale=False) +fit = dl.fitparamodel(V[mask],Bmodel,par0,lb,ub,rescale=False) lamfit,kappa,d = fit.param -Bfit = bg_strexp(t,[kappa,d],lamfit) +Bfit = dl.bg_strexp(t,[kappa,d],lamfit) # Background "correction" by division Vcorr = (V/Bfit - 1 + lamfit)/lamfit # Tikhonov regularization using the L-curve criterion -K = dipolarkernel(t,r) -fit = fitregmodel(Vcorr,K,r,'tikhonov','lr',) +K = dl.dipolarkernel(t,r) +fit = dl.fitregmodel(Vcorr,K,r,'tikhonov','lr',) Pfit = fit.P # %% [markdown] diff --git a/examples/plot_extracting_gauss_constraints.py b/examples/plot_extracting_gauss_constraints.py index 8091cf9ee..648d45aee 100644 --- a/examples/plot_extracting_gauss_constraints.py +++ b/examples/plot_extracting_gauss_constraints.py @@ -1,58 +1,44 @@ # %% [markdown] """ -Extracting Gaussian constraints from a parameter-free distribution fit +Fitting Gaussians to a non-parametric distance distribution ======================================================================= -How to extract Gaussian constraints from a parameter-free fit. - -While parameter-free distance distributions are the most robust way to -analyze dipolar signals, many structural biology modelling programs -accept only estimators such as mean distances or Gaussian constraints. - -This example shows how to extract Gaussian constraints from a -parameter-free fit of a dipolar signal and how to calculate the -corresponding uncertainty. +This example shows how to fit Gaussians to a non-parametric distance +distribution obtained via Tikhonov regularization and how to calculate +the corresponding uncertainty. """ # %% import numpy as np import matplotlib.pyplot as plt -from deerlab import * +import deerlab as dl # %% -# Generating a dataset -# For this example we will simulate a simple 4pDEER signal +# Simulate a 4pDEER signal -# Parameters -t = np.linspace(0,5,250) -r = np.linspace(1,7,200) -P = dd_gauss3(r,[4.5, 0.6, 0.4, 3, 0.4, 0.3, 4, 0.7, 0.5]) -lam = 0.3 -conc = 80; #uM +t = np.linspace(0,5,250) # time axis, µs +r = np.linspace(1,7,200) # distance axis, nm +P = dl.dd_gauss3(r,[4.5, 0.35, 0.4, 3, 0.25, 0.3, 4, 0.4, 0.5]) # distance distribution +lam = 0.3 # modulation depth +conc = 80 # spin concentration, µM -# Simulate the signal -Bmodel = lambda t,lam: bg_hom3d(t,conc,lam) -K = dipolarkernel(t,r,lam,Bmodel) -np.random.seed(0) -V = K@P + whitegaussnoise(t,0.01) - -# %% +B = dl.bg_hom3d(t,conc,lam) # background +K = dl.dipolarkernel(t,r,lam,B) # kernel matrix +V = K@P + dl.whitegaussnoise(t,0.01,seed=0) # DEER trace, with added noise # %% [markdown] # Fit the dipolar signal # ---------------------- -# First, we need to fit the parameter-free distance distribution using ``fitsignal()`` -# We are only interested right now on the fitted distribution and the -# corresponding uncertainty quantification, so we will ignore the rest of -# the outputs. +# First, we fit the non-parametric distance distribution using ``fitsignal()`` # %% -fit = fitsignal(V,t,r,'P',bg_exp,ex_4pdeer,display=True) -Pfit = fit.P +fit = dl.fitsignal(V,t,r,'P',dl.bg_exp,dl.ex_4pdeer) +fit.plot() + # %% [markdown] -# Extract Gaussian constraints from the fit +# Fit Gaussians to the distance distribution # ----------------------------------------- -# Next, we will fit a multi-Gauss distribution to the fitted parameter-free +# Next, we fit a multi-Gauss distribution to the fitted non-parametric # distribution. We can do this by using the ``fitparamodel()`` function (in # this example, fitting a two-Gauss model). # @@ -62,9 +48,10 @@ # ``Pfit`` to the Gauss constraints that we then fit. # %% -# Extract the uncertainty quantification of the fitted distribution... + +# From the fit results, extract the distribution and the covariance matrix +Pfit = fit.P Pfit_uq = fit.Puncert -# ...specifically its covariance matrix Pfit_covmat = Pfit_uq.covmat # %% [markdown] @@ -74,35 +61,38 @@ # - ``PGauss```: the corresponding distribution # - ``paruq```: the uncertainty quantification of our constraints - # %% -Pmodel = lambda p: dd_gauss2(r,p) +Pmodel = lambda p: dl.dd_gauss2(r,p) # Get information on the model -info = dd_gauss2() +info = dl.dd_gauss2() par0 = info['Start'] lb = info['Lower'] ub = info['Upper'] -fit = fitparamodel(Pfit,Pmodel,par0,lb,ub,covmatrix=Pfit_covmat) + +# Fit the Gaussians +fit = dl.fitparamodel(Pfit,Pmodel,par0,lb,ub,covmatrix=Pfit_covmat,rescale=False) + +# Extract the fit results parfit = fit.param paruq = fit.uncertainty -PGauss = dd_gauss2(r,parfit) +PGauss = dl.dd_gauss2(r,parfit) -# Extract the 95#-confidence intervals... +# Extract the 95%-confidence intervals... par95 = paruq.ci(95) -# ... and print the results of the constraints -print('\nGaussian constraints:') -info = dd_gauss2() +# ... and print the results +print('\nGaussian components:') +info = dl.dd_gauss2() for i in range(len(parfit)): print(' parfit[{}] = {:2.2f} ({:2.2f}, {:2.2f}) {}'.format(i,parfit[i],par95[i,0],par95[i,1],info['Parameters'][i])) -# Now propagate the error of the constraints on the model -lb = np.zeros_like(r) # Non-negativity constraint -PGauss_uq = paruq.propagate(lambda par: dd_gauss2(r,par),lb) +# Now propagate the error of the Gaussian parameters to the distribution +lb = np.zeros_like(r) # non-negativity constraint +PGauss_uq = paruq.propagate(lambda par: dl.dd_gauss2(r,par),lb) PGauss95 = PGauss_uq.ci(95) # %% -# Plot the fitted constraints model on top of the parameter-free case +# Plot the fitted constraints model on top of the non-parametric case plt.plot(r,Pfit,'r',linewidth=1.5) plt.fill_between(r,Pfit_uq.ci(95)[:,0], Pfit_uq.ci(95)[:,1],facecolor='r',linestyle='None',alpha=0.2) @@ -113,8 +103,5 @@ plt.ylabel('P [nm$^{-1}$]') plt.tight_layout() plt.grid(alpha=0.3) -plt.legend(['Fit','95%-CI','2G-constraints','95%-CI']) +plt.legend(['Tikhonov fit','95%-CI','2G fit to Tikhonov','95%-CI']) plt.show() - - -# %% diff --git a/examples/plot_fitting_4pdeer.py b/examples/plot_fitting_4pdeer.py index 99904f1b7..eb1b8dd8d 100644 --- a/examples/plot_fitting_4pdeer.py +++ b/examples/plot_fitting_4pdeer.py @@ -3,29 +3,31 @@ Basic fitting of a 4-pulse DEER signal, non-parametric distribution ------------------------------------------------------------------- -How to fit a simple 4-pulse DEER signal with a parameter-free distribution. +Fit a simple 4-pulse DEER signal with a model with a non-parametric +distribution distribution and a homogeneous background, using Tikhonov regularization. """ import numpy as np import matplotlib.pyplot as plt -from deerlab import * +import deerlab as dl # %% [markdown] #Generate data #-------------- -t = np.linspace(-0.1,4,250) # time axis, us -r = np.linspace(1.5,6,len(t)) # distance axis, ns -param = [3, 0.3, 0.2, 3.5, 0.3, 0.65, 3.8, 0.2, 0.15] # parameters for three-Gaussian model -P = dd_gauss3(r,param) # model distance distribution -lam = 0.5 # modulation depth -B = bg_hom3d(t,300,lam) # background decay -K = dipolarkernel(t,r,lam,B) -np.random.seed(0) -Vexp = K@P + whitegaussnoise(t,0.01) - +t = np.linspace(-0.1,4,250) # time axis, µs +r = np.linspace(1.5,6,len(t)) # distance axis, ns +param = [3, 0.1, 0.2, 3.5, 0.1, 0.65, 3.8, 0.05, 0.15] # parameters for three-Gaussian model +P = dl.dd_gauss3(r,param) # model distance distribution +lam = 0.5 # modulation depth +B = dl.bg_hom3d(t,300,lam) # background decay +K = dl.dipolarkernel(t,r,lam,B) # kernel matrix +Vexp = K@P + dl.whitegaussnoise(t,0.01,seed=0) # %% [markdown] # Run fit #--------- -fit = fitsignal(Vexp,t,r,'P',bg_hom3d,ex_4pdeer,display=True) +fit = dl.fitsignal(Vexp,t,r,'P',dl.bg_hom3d,dl.ex_4pdeer) +fit.plot() + +# %% diff --git a/examples/plot_fitting_5pdeer.py b/examples/plot_fitting_5pdeer.py index 397cd63d9..bf680a460 100644 --- a/examples/plot_fitting_5pdeer.py +++ b/examples/plot_fitting_5pdeer.py @@ -9,22 +9,20 @@ # %% import numpy as np import matplotlib.pyplot as plt -from deerlab import * - +import deerlab as dl # %% # Generate data -np.random.seed(1) -t = np.linspace(-0.1,6.5,200) # time axis, us -r = np.linspace(1.5,6,100) # distance axis, ns -param0 = [3, 0.3, 0.2, 3.5, 0.3, 0.65, 3.8, 0.2, 0.15] # parameters for three-Gaussian model -P = dd_gauss3(r,param0) # model distance distribution -B = lambda t,lam: bg_hom3d(t,300,lam) # background decay +t = np.linspace(-0.1,6.5,200) # time axis, µs +r = np.linspace(1.5,6,100) # distance axis, nm +param0 = [3, 0.1, 0.2, 3.5, 0.1, 0.65, 3.8, 0.05, 0.15] # parameters for three-Gaussian model +P = dl.dd_gauss3(r,param0) # model distance distribution +B = lambda t,lam: dl.bg_hom3d(t,300,lam) # background decay exparam = [0.6, 0.3, 0.1, 3.2] # parameters for 5pDEER experiment -pathinfo = ex_5pdeer(exparam) # pathways information +pathinfo = dl.ex_5pdeer(exparam) # pathways information -K = dipolarkernel(t,r,pathinfo,B) -Vexp = K@P + whitegaussnoise(t,0.005) +K = dl.dipolarkernel(t,r,pathinfo,B) +Vexp = K@P + dl.whitegaussnoise(t,0.005,seed=1) # %% [markdown] # Now, 5pDEER data contain 3 additional parameters compared to 4pDEER (due @@ -57,5 +55,5 @@ # Run the fit with a 5-pulse DEER signal model # %% -fit = fitsignal(Vexp,t,r,'P',bg_hom3d,ex_5pdeer,par0,lb,ub,display=True) -plt.show() \ No newline at end of file +fit = dl.fitsignal(Vexp,t,r,'P',dl.bg_hom3d,dl.ex_5pdeer,par0,lb,ub) +fit.plot() diff --git a/examples/plot_fitting_custom_kernel_with_snlls.py b/examples/plot_fitting_custom_kernel_with_snlls.py index 9625f3670..b63a5d154 100644 --- a/examples/plot_fitting_custom_kernel_with_snlls.py +++ b/examples/plot_fitting_custom_kernel_with_snlls.py @@ -8,7 +8,7 @@ """ import numpy as np import matplotlib.pyplot as plt -from deerlab import * +import deerlab as dl # %% [markdown] # Generating a dataset @@ -19,12 +19,12 @@ r = np.linspace(2,6,200) # Generate ground truth and input signal -P = dd_gauss2(r,[3.5, 0.4, 0.4, 4.5, 0.7, 0.6]) +P = dl.dd_gauss2(r,[3.5, 0.25, 0.4, 4.5, 0.4, 0.6]) lam = 0.36 c0 = 250 #uM -B = bg_hom3d(t,c0,lam) -K = dipolarkernel(t,r,lam,B) -V = K@P + whitegaussnoise(t,0.01) +B = dl.bg_hom3d(t,c0,lam) +K = dl.dipolarkernel(t,r,lam,B) +V = K@P + dl.whitegaussnoise(t,0.01) # %% [markdown] # Fitting via SNLLS @@ -43,9 +43,9 @@ def Kmodel(p): # Unpack parameters lam,c0 = p # Get background - B = bg_hom3d(t,c0,lam) + B = dl.bg_hom3d(t,c0,lam) # Generate 4pDEER kernel - K = dipolarkernel(t,r,lam,B) + K = dl.dipolarkernel(t,r,lam,B) return K @@ -72,7 +72,7 @@ def Kmodel(p): ubl = [] # Unconstrained upper boundary # Run SNLLS optimization -fit = snlls(V,Kmodel,par0,lb,ub,lbl,ubl) +fit = dl.snlls(V,Kmodel,par0,lb,ub,lbl,ubl) parfit = fit.nonlin Pfit = fit.lin paruq = fit.uncertainty @@ -112,8 +112,4 @@ def Kmodel(p): plt.ylabel('P(r) [nm$^{-1}$]') plt.legend(['truth','fit','50%-CI','95%-CI']) - - - - # %% diff --git a/examples/plot_fitting_mixed_model.py b/examples/plot_fitting_mixed_model.py index a4f1130f4..bd3ffec4e 100644 --- a/examples/plot_fitting_mixed_model.py +++ b/examples/plot_fitting_mixed_model.py @@ -9,7 +9,7 @@ import numpy as np import matplotlib.pyplot as plt -from deerlab import * +import deerlab as dl # %% [markdown] # Simulate the data @@ -24,20 +24,19 @@ # Distribution parameters rmean = 4.5 -width = 0.3 +sigma = 0.2 chain = 4.3 pers = 10 amp = 0.35 # Generate distribution -P = dd_gauss(r,[rmean, width]) -P = amp*P + (1 - amp)*dd_wormchain(r,[chain, pers]) +P = dl.dd_gauss(r,[rmean, sigma]) +P = amp*P + (1 - amp)*dl.dd_wormchain(r,[chain, pers]) # Normalize distribution P = P/sum(P)/np.mean(np.diff(r)) # Generate dipolar evolution function -np.random.seed(0) -K = dipolarkernel(t,r) -V = K@P + whitegaussnoise(t,0.02) +K = dl.dipolarkernel(t,r) +V = K @ P + dl.whitegaussnoise(t,0.02,seed=0) # %% # Generating a mixed parametric model @@ -53,7 +52,7 @@ # parametric models as lambda functions. #Mix the models into new one -gausswlc = mixmodels(dd_gauss,dd_wormchain) +gausswlc = dl.mixmodels(dl.dd_gauss,dl.dd_wormchain) # %% [markdown] # Our new model ``gausswlc`` will now describe our sought linear combination of @@ -81,15 +80,15 @@ # a very basic dipolar kernel. #Generate the dipolar evolution function kernel -K = dipolarkernel(t,r) +K = dl.dipolarkernel(t,r) #Fit the model to the data -Vmodel = lambda par: K@gausswlc(r,par) +Vmodel = lambda par: K @ gausswlc(r,par) info = gausswlc() par0 = info['Start'] # built-in start values lb = info['Lower'] # built-in lower bounds ub = info['Upper'] # built-in upper bounds -fit = fitparamodel(V,Vmodel,par0,lb,ub,multistart=10) +fit = dl.fitparamodel(V,Vmodel,par0,lb,ub,multistart=10) fitpar = fit.param # %% [markdown] # From the fitted parameter set ``fitpar`` we can now generate our fitted distance @@ -116,5 +115,4 @@ plt.ylabel('P(r) [nm$^{-1}$]') plt.legend(['truth','fit']) - # %% diff --git a/examples/plot_globalfits_local_and_global_parameters.py b/examples/plot_globalfits_local_and_global_parameters.py index 70d7b6a7d..d1b951ae7 100644 --- a/examples/plot_globalfits_local_and_global_parameters.py +++ b/examples/plot_globalfits_local_and_global_parameters.py @@ -11,7 +11,7 @@ import numpy as np import matplotlib.pyplot as plt -from deerlab import * +import deerlab as dl # %% [markdown] # Generate two datasets @@ -29,26 +29,24 @@ # Parameters rmeanA = 3.45 # mean distance state A, in nm rmeanB = 5.05 # mean distance state B, in nm -fwhmA = 0.5 # FWHM state A, in nm -fwhmB = 0.3 # FWHM state B, in nm +sigmaA = 0.3 # standard deviation state A, in nm +sigmaB = 0.2 # standard deviation state B, in nm fracA1 = 0.8 # Molar fraction of state A under conditions 1 fracA2 = 0.2 # Molar fraction of state A under conditions 2 # The molar fraction of state B is not required as it follows fracB = 1 - fracA # Generate the two distributions for conditions 1 & 2 -P1 = dd_gauss2(r,[rmeanA, fwhmA, fracA1, rmeanB, fwhmB, 1-fracA1]) -P2 = dd_gauss2(r,[rmeanA, fwhmA, fracA2, rmeanB, fwhmB, 1-fracA2]) +P1 = dl.dd_gauss2(r,[rmeanA, sigmaA, fracA1, rmeanB, sigmaB, 1-fracA1]) +P2 = dl.dd_gauss2(r,[rmeanA, sigmaA, fracA2, rmeanB, sigmaB, 1-fracA2]) # Generate the corresponding dipolar kernels -K1 = dipolarkernel(t1,r) -K2 = dipolarkernel(t2,r) +K1 = dl.dipolarkernel(t1,r) +K2 = dl.dipolarkernel(t2,r) # ...and the two corresponding signals -np.random.seed(0) -V1 = K1@P1 + whitegaussnoise(t1,0.01) -np.random.seed(1) -V2 = K2@P2 + whitegaussnoise(t2,0.02) +V1 = K1@P1 + dl.whitegaussnoise(t1,0.01,seed=0) +V2 = K2@P2 + dl.whitegaussnoise(t2,0.02,seed=1) # (for the sake of simplicity no background and 100# modulation depth are assumed) # %% [markdown] @@ -60,7 +58,7 @@ # signal (local). # # In this examples we have the following parameters: -# - fixed: ``fwhmA``, ``fwhmB`` (known paramters) +# - fixed: ``sigmaA``, ``sigmaB`` (known paramters) # - global: ``rmeanA``, ``rmeanB`` (same for both signals) # - local: ``fracA1``, ``fracA2`` (different for both signals/conditions) # @@ -74,8 +72,8 @@ def myABmodel(par): #Fixed parameters - fwhmA = 0.5 - fwhmB = 0.3 + sigmaA = 0.5 + sigmaB = 0.3 #Global parameters rmeanA = par[0] rmeanB = par[1] @@ -84,13 +82,13 @@ def myABmodel(par): fracA2 = par[3] # Generate the signal-specific distribution - Pfit1 = dd_gauss2(r,[rmeanA, fwhmA, fracA1, rmeanB, fwhmB, max(1-fracA1,0)]) - Pfit2 = dd_gauss2(r,[rmeanA, fwhmA, fracA2, rmeanB, fwhmB, max(1-fracA2,0)]) + Pfit1 = dl.dd_gauss2(r,[rmeanA, sigmaA, fracA1, rmeanB, sigmaB, max(1-fracA1,0)]) + Pfit2 = dl.dd_gauss2(r,[rmeanA, sigmaA, fracA2, rmeanB, sigmaB, max(1-fracA2,0)]) # Generate signal #1 - V1fit = K1@Pfit1 + V1fit = K1 @ Pfit1 # Generate signal #2 - V2fit = K2@Pfit2 + V2fit = K2 @ Pfit2 # Return as a list Vfits = [V1fit,V2fit] @@ -119,7 +117,7 @@ def myABmodel(par): Vs = [V1,V2] # Fit the global parametric model to both signals -fit = fitparamodel(Vs,model,par0,lower,upper,multistart=40) +fit = dl.fitparamodel(Vs,model,par0,lower,upper,multistart=40) # The use of the option 'multistart' will help the solver to find the # global minimum and not to get stuck at local minima. @@ -162,8 +160,4 @@ def myABmodel(par): plt.tight_layout() - -# %% - - # %% diff --git a/examples/plot_multigauss_fitting_4pdeer.py b/examples/plot_multigauss_fitting_4pdeer.py index 012186653..315e85f7d 100644 --- a/examples/plot_multigauss_fitting_4pdeer.py +++ b/examples/plot_multigauss_fitting_4pdeer.py @@ -4,13 +4,12 @@ ======================================== This example showcases how to fit a simple 4-pulse DEER signal with -background using a multi-Gauss model, i.e automatically optimizing the -number of Gaussians in the model. +background using a multi-Gauss distance distribution model, determining +the optimal number of Gaussians. """ import numpy as np import matplotlib.pyplot as plt -from deerlab import * - +import deerlab as dl # %% [markdown] # Model function for a 4pDEER dipolar kernel @@ -23,9 +22,9 @@ def K4pdeer(par,t,r): # Unpack parameters lam,conc = par # Simualte background - B = bg_hom3d(t,conc,lam) + B = dl.bg_hom3d(t,conc,lam) # Generate dipolar kernel - K = dipolarkernel(t,r,lam,B) + K = dl.dipolarkernel(t,r,lam,B) return K @@ -35,24 +34,24 @@ def K4pdeer(par,t,r): t = np.linspace(-0.25,4,300) # time axis, us r = np.linspace(2.5,4.5,300) # distance axis, nm -param0 = [3, 0.3, 0.2, 3.5, 0.3, 0.45, 3.9, 0.2, 0.20] # parameters for three-Gaussian model -P = dd_gauss3(r,param0) # ground truth distance distribution +param0 = [3, 0.1, 0.2, 3.5, 0.1, 0.45, 3.9, 0.05, 0.20] # parameters for three-Gaussian model +P = dl.dd_gauss3(r,param0) # ground truth distance distribution lam = 0.3 # modulation depth conc = 250 # spin concentration, uM noiselvl = 0.005 # noise level # Generate 4pDEER dipolar signal with noise -np.random.seed(0) -V = K4pdeer([lam,conc],t,r)@P + whitegaussnoise(t,noiselvl) +K = K4pdeer([lam,conc],t,r) +V = K @ P + dl.whitegaussnoise(t,noiselvl,seed=0) # %% [markdown] # Multi-Gauss fitting # ------------------- # Parameter bounds: -# lambda conc rmean fwhm +# lambda conc rmean sigma lb = [1, 0.05] # distribution basis function lower bounds -ub = [20, 5] # distribution basis function upper bounds +ub = [20, 2.5] # distribution basis function upper bounds lbK = [0, 0.05] # kernel parameters lower bounds ubK = [1, 1500] # kernel parameters upper bounds @@ -61,7 +60,7 @@ def K4pdeer(par,t,r): NGauss = 5 # maximum number of Gaussians # Fit the kernel parameters with an optimized multi-Gauss distribution -fit = fitmultimodel(V,Kmodel,r,dd_gauss,NGauss,'aic',lb,ub,lbK,ubK) +fit = dl.fitmultimodel(V,Kmodel,r,dl.dd_gauss,NGauss,'aic',lb,ub,lbK,ubK) #Extract results Pfit = fit.P Kparfit = fit.Kparam @@ -96,11 +95,11 @@ def K4pdeer(par,t,r): plt.subplot(3,2,1) plt.plot(t,V,'k.') plt.plot(t,Vfit,'b',linewidth=1.5) -plt.plot(t,(1-Kparfit[0])*bg_hom3d(t,Kparfit[1],Kparfit[0]),'b--',linewidth=1.5) +plt.plot(t,(1-Kparfit[0])*dl.bg_hom3d(t,Kparfit[1],Kparfit[0]),'b--',linewidth=1.5) plt.tight_layout() plt.grid(alpha=0.3) plt.legend(['data','Vfit','Bfit']) -plt.xlabel('t [$\mu s$]') +plt.xlabel('t [µs]') plt.ylabel('V(t)') plt.subplot(322) diff --git a/examples/plot_param_uncertainty_propagation.py b/examples/plot_param_uncertainty_propagation.py index 6a94955a5..dd205c527 100644 --- a/examples/plot_param_uncertainty_propagation.py +++ b/examples/plot_param_uncertainty_propagation.py @@ -9,7 +9,7 @@ import numpy as np import matplotlib.pyplot as plt -from deerlab import * +import deerlab as dl # %% [markdown] # Generate data @@ -21,11 +21,10 @@ width = 0.3 # [nm] Rician width lam = 0.27 # Modulation depth conc = 150 # [uM] Spin concentration -P = dd_rice(r,[center, width]) -B = bg_hom3d(t,conc,lam) -K = dipolarkernel(t,r,lam,B) -np.random.seed(0) -V = K@P + whitegaussnoise(t,0.03) +P = dl.dd_rice(r,[center, width]) +B = dl.bg_hom3d(t,conc,lam) +K = dl.dipolarkernel(t,r,lam,B) +V = K@P + dl.whitegaussnoise(t,0.03,seed=0) # %% [markdown] # Fit the data @@ -47,10 +46,10 @@ # parameter set directly. # Pre-calculate the elemental dipolar kernel (for speed) -K0 = dipolarkernel(t,r) +K0 = dl.dipolarkernel(t,r) -Pmodel = lambda par: dd_rice(r,par[1:3]) -Bmodel = lambda par: bg_hom3d(t,par[3],par[0]) +Pmodel = lambda par: dl.dd_rice(r,par[1:3]) +Bmodel = lambda par: dl.bg_hom3d(t,par[3],par[0]) Vmodel = lambda par: (1 - par[0] + par[0]*K0@Pmodel(par))*Bmodel(par) # %% [markdown] @@ -63,7 +62,7 @@ upper = [0.50, 7.0, 0.5, 1500] # upper bounds # Finally we can run the fit and get the fitted parameters and their uncertainties -fit = fitparamodel(V,Vmodel,par0,lower,upper) +fit = dl.fitparamodel(V,Vmodel,par0,lower,upper) parfit = fit.param paruq = fit.uncertainty diff --git a/examples/plot_pseudotitration_parameter_free.py b/examples/plot_pseudotitration_parameter_free.py index 80cab364f..3e99fdfd6 100644 --- a/examples/plot_pseudotitration_parameter_free.py +++ b/examples/plot_pseudotitration_parameter_free.py @@ -10,7 +10,7 @@ import numpy as np import matplotlib.pyplot as plt -from deerlab import * +import deerlab as dl # %% [markdown] # Generating multiple datasets @@ -63,11 +63,11 @@ def Kmodel(par,ts,rA,rB,L): Ks = [[]]*Nsignals # General the dipolar kernels for i in range(Nsignals): - B = bg_exp(ts[i],k,lam) + B = dl.bg_exp(ts[i],k,lam) # Kernel for fraction A - KstateA = dipolarkernel(ts[i],rA,lam,B) + KstateA = dl.dipolarkernel(ts[i],rA,lam,B) # Kernel for fraction B - KstateB = dipolarkernel(ts[i],rB,lam,B) + KstateB = dl.dipolarkernel(ts[i],rB,lam,B) Ks[i] = np.concatenate((xA[i]*KstateA, xB[i]*KstateB),axis=1) return Ks @@ -92,8 +92,8 @@ def Kmodel(par,ts,rA,rB,L): rB = np.linspace(1,8,100) # Distributions for states A and B -PstateA = dd_gauss(rA,[5.5, 0.4]) -PstateB = dd_gauss2(rB,[4.5, 0.7, 0.4, 3.5, 0.6, 0.6]) +PstateA = dl.dd_gauss(rA,[5.5, 0.25]) +PstateB = dl.dd_gauss2(rB,[4.5, 0.4, 0.4, 3.5, 0.35, 0.6]) L = [0.3, 1, 3, 10, 30, 100, 300] # total ligand concentration, uM Kdis = 5.65 # dissociation constant, uM @@ -107,8 +107,7 @@ def Kmodel(par,ts,rA,rB,L): # Simulate dipolar signals Vs = [[]]*Nsignals for i in range(Nsignals): - np.random.seed(i) - Vs[i] = Ks[i]@np.concatenate((PstateA, PstateB)) + whitegaussnoise(ts[i],0.005) + Vs[i] = Ks[i]@np.concatenate((PstateA, PstateB)) + dl.whitegaussnoise(ts[i],0.005,seed=i) # %% [markdown] # Psuedotitration SNLLS Analysis @@ -130,7 +129,7 @@ def Kmodel(par,ts,rA,rB,L): ubl = [] # Unconstrained # Run SNLLS optimization -fit = snlls(Vs,lambda p: Kmodel(p,ts,rA,rB,L),par0,lb,ub,lbl,ubl) +fit = dl.snlls(Vs,lambda p: Kmodel(p,ts,rA,rB,L),par0,lb,ub,lbl,ubl) # Extract fit results parfit = fit.nonlin Pfit = fit.lin @@ -181,7 +180,3 @@ def Kmodel(par,ts,rA,rB,L): plt.ylabel('Fractions') plt.legend(['state A','state B']) plt.ylim([0,1]) - - - -# %% diff --git a/examples/plot_selecting_an_optimal_model.py b/examples/plot_selecting_an_optimal_model.py index afb1be8fa..45dcfd9a2 100644 --- a/examples/plot_selecting_an_optimal_model.py +++ b/examples/plot_selecting_an_optimal_model.py @@ -8,7 +8,7 @@ import numpy as np import matplotlib.pyplot as plt -from deerlab import * +import deerlab as dl # %% [markdown] # Data Generation @@ -18,13 +18,11 @@ # from a bimodal Gaussian distance distribution. # Prepare the signal components -t = np.linspace(-0.3,6,300) -r = np.linspace(2,6,200) -P = dd_gauss2(r,[3.8, 0.7, 0.7, 4.5, 0.3, 0.7]) - -# Prepare the dipolar kernel and get the signal -K = dipolarkernel(t,r) -V = K@P + whitegaussnoise(t,0.01) +t = np.linspace(-0.3,3.5,300) # time axis, µs +r = np.linspace(2,6,200) # distance axis, nm +P = dl.dd_gauss2(r,[3.8, 0.4, 0.7, 4.5, 0.2, 0.7]) # distance distribution +K = dl.dipolarkernel(t,r) # dipolar kernel matrix +V = K@P + dl.whitegaussnoise(t,0.02) # DEER signal, with added noise # %% [markdown] # Selecting an optimal model @@ -45,17 +43,17 @@ # The last model we can construct from built-in models using the ``mixmodels`` function. # Prepare the mixed model -dd_rice_gauss = mixmodels(dd_rice,dd_gauss) +dd_rice_gauss = dl.mixmodels(dl.dd_rice,dl.dd_gauss) # Prepare list of candidate parametric models -models = [dd_rice,dd_rice2,dd_rice3,dd_gauss,dd_gauss2,dd_gauss3,dd_rice_gauss] +models = [dl.dd_rice,dl.dd_rice2,dl.dd_rice3,dl.dd_gauss,dl.dd_gauss2,dl.dd_gauss3,dd_rice_gauss] # %% [markdown] # In order to make an appropiate choice, we need some liklihood estimator. All fit functions is DeerLab returns a stats # dictionary which contains (amongst other estimators) likelihood estimators such as the Akaike information criterion (AIC). # The model with the lowers AIC value can be considered to most likely to be the optimal model. # -# To do this, we jus have to evaluate the parametric models with ``fitparamodel`` while looping over all the distribution models +# To do this, we just have to evaluate the parametric models with ``fitparamodel`` while looping over all the distribution models # we listed above, and collecting the AIC-values for each model. aic = [] @@ -64,7 +62,7 @@ # Prepare the signal model with the new distance model Vmodel = lambda par: K@model(r,par) # Fit the signal - fit = fitparamodel(V,Vmodel,par0=info['Start'],lb=info['Lower'],ub=info['Upper']) + fit = dl.fitparamodel(V,Vmodel,par0=info['Start'],lb=info['Lower'],ub=info['Upper']) parfit = fit.param stats= fit.stats # Add current AIC value to the list @@ -93,8 +91,7 @@ plt.subplot(2,2,1) plt.plot(t,V,'k.') plt.grid(alpha=0.2) -plt.xlabel('t [$\mu s$]') -plt.ylabel('V(t)') +plt.xlabel('t [µs]') plt.legend(['data']) plt.subplot(2,2,2) diff --git a/setup.py b/setup.py index fddc8e451..389bbee84 100644 --- a/setup.py +++ b/setup.py @@ -9,7 +9,7 @@ def install_dependencies(): subprocess.run(['pip','install','memoization'],check=True) subprocess.run(['pip','install','matplotlib'],check=True) # Dependencies with OS-specific BLAS - if platform.system() is 'Windows': + if platform.system() == 'Windows': # Install Numpy,SciPy, CVXopt linked to MKL subprocess.run(['pip','install','pipwin'],check=True) subprocess.run(['pipwin','install','numpy'],check=False) diff --git a/test/test_bgmodels.py b/test/test_bgmodels.py index 43124ce1d..3cb33b3b3 100644 --- a/test/test_bgmodels.py +++ b/test/test_bgmodels.py @@ -1,5 +1,5 @@ import numpy as np -from deerlab.bg_models import * +import deerlab as dl def assert_bgmodel(model,Bref): "Check the correct behaviour of the core functionality of a background model" @@ -46,36 +46,36 @@ def assert_bgmodel(model,Bref): def test_bg_exp(): Bref = 0.416862019678508 # Reference from DeerLab 0.9.2 on MATLAB - assert_bgmodel(bg_exp,Bref) + assert_bgmodel(dl.bg_exp,Bref) def test_bg_hom3d(): Bref = 0.882785339742350 # Reference from DeerLab 0.9.2 on MATLAB - assert_bgmodel(bg_hom3d,Bref) + assert_bgmodel(dl.bg_hom3d,Bref) def test_bg_hom3dex(): Bref = 0.882896490000000 # Reference from DeerLab 0.9.2 on MATLAB - assert_bgmodel(bg_hom3dex,Bref) + assert_bgmodel(dl.bg_hom3dex,Bref) def test_bg_strexp(): Bref = 0.535261428518990 # Reference from DeerLab 0.9.2 on MATLAB - assert_bgmodel(bg_strexp,Bref) + assert_bgmodel(dl.bg_strexp,Bref) def test_bg_prodstrexp(): Bref = 0.286504796860190 # Reference from DeerLab 0.9.2 on MATLAB - assert_bgmodel(bg_prodstrexp,Bref) + assert_bgmodel(dl.bg_prodstrexp,Bref) def test_bg_sumstrexp(): Bref = 0.535261428518990 # Reference from DeerLab 0.9.2 on MATLAB - assert_bgmodel(bg_sumstrexp,Bref) + assert_bgmodel(dl.bg_sumstrexp,Bref) def test_bg_poly1(): Bref = -1.500000000000000 # Reference from DeerLab 0.9.2 on MATLAB - assert_bgmodel(bg_poly1,Bref) + assert_bgmodel(dl.bg_poly1,Bref) def test_bg_poly2(): Bref = -7.750000000000000 # Reference from DeerLab 0.9.2 on MATLAB - assert_bgmodel(bg_poly2,Bref) + assert_bgmodel(dl.bg_poly2,Bref) def test_bg_poly3(): Bref = -23.37500000000000 # Reference from DeerLab 0.9.2 on MATLAB - assert_bgmodel(bg_poly3,Bref) \ No newline at end of file + assert_bgmodel(dl.bg_poly3,Bref) diff --git a/test/test_bootan.py b/test/test_bootan.py index bb41c843f..0c0d5d8f1 100644 --- a/test/test_bootan.py +++ b/test/test_bootan.py @@ -54,7 +54,7 @@ def bootfcn(V): paruq2 = bootan(bootfcn,Vexp,Vfit,10,resampling='gaussian') - assert all(abs(paruq1.mean - paruq2.mean) < 1e-2) + assert all(abs(paruq1.mean - paruq2.mean) < 1.5e-2) # ====================================================================== diff --git a/test/test_correctscale.py b/test/test_correctscale.py index ca013327e..129692eee 100644 --- a/test/test_correctscale.py +++ b/test/test_correctscale.py @@ -5,7 +5,7 @@ def test_correction(): # ====================================================================== - "Check that scale correction works properly" + "Check that scale correction using DEER model works" t = np.linspace(-1,5,300) r = np.linspace(3,7,200) @@ -14,14 +14,14 @@ def test_correction(): scale = 1e8 V = dipolarkernel(t,r,0.25,B)@P - Vc = correctscale(V*scale,t) + Vc = correctscale(V*scale,t,model='deer') assert max(abs(Vc-V)) < 1.5e-1 # ====================================================================== def test_model_gauss(): # ====================================================================== - "Check scale correction using the gauss model works" + "Check scale correction using the Gauss model works" t = np.linspace(-1,5,300) r = np.linspace(3,7,200) diff --git a/test/test_correctzerotime.py b/test/test_correctzerotime.py index 3ba01aca0..98c43e438 100644 --- a/test/test_correctzerotime.py +++ b/test/test_correctzerotime.py @@ -32,3 +32,37 @@ def test_late_times(): assert max(abs(t_corr - t_truth)) < 1e-10 #======================================================================= + +def test_first_element(): +#======================================================================= + "Check that zero-time correction works when maximum is at later times" + + + t_truth = np.linspace(0,5,400) + y = 10 - t_truth**2 + y = y/max(y) + tshift = 1.2343 + + t = t_truth + tshift + + t_corr = correctzerotime(y,t) + + assert max(abs(t_corr - t_truth)) < 1e-10 +#======================================================================= + +def test_last_element(): +#======================================================================= + "Check that zero-time correction works when maximum is at later times" + + + t_truth = np.linspace(-5,0,400) + y = 10 - t_truth**2 + y = y/max(y) + tshift = 1.2343 + + t = t_truth + tshift + + t_corr = correctzerotime(y,t) + + assert max(abs(t_corr - t_truth)) < 1e-10 +#======================================================================= diff --git a/test/test_ddmodels.py b/test/test_ddmodels.py index fe51babba..b95a8956c 100644 --- a/test/test_ddmodels.py +++ b/test/test_ddmodels.py @@ -1,5 +1,5 @@ import numpy as np -from deerlab.dd_models import * +import deerlab as dl def assert_ddmodel(model): #============================================================================ @@ -45,70 +45,70 @@ def assert_ddmodel(model): #============================================================================= def test_dd_gauss(): - assert_ddmodel(dd_gauss) + assert_ddmodel(dl.dd_gauss) def test_dd_gauss2(): - assert_ddmodel(dd_gauss2) + assert_ddmodel(dl.dd_gauss2) def test_dd_gauss3(): - assert_ddmodel(dd_gauss3) + assert_ddmodel(dl.dd_gauss3) def test_dd_gengauss(): - assert_ddmodel(dd_gengauss) + assert_ddmodel(dl.dd_gengauss) def test_dd_skewgauss(): - assert_ddmodel(dd_skewgauss) + assert_ddmodel(dl.dd_skewgauss) def test_dd_rice(): - assert_ddmodel(dd_rice) + assert_ddmodel(dl.dd_rice) def test_dd_rice2(): - assert_ddmodel(dd_rice2) + assert_ddmodel(dl.dd_rice2) def test_dd_rice3(): - assert_ddmodel(dd_rice3) + assert_ddmodel(dl.dd_rice3) def test_dd_randcoil(): - assert_ddmodel(dd_randcoil) + assert_ddmodel(dl.dd_randcoil) def test_dd_circle(): - assert_ddmodel(dd_circle) + assert_ddmodel(dl.dd_circle) def test_dd_cos(): - assert_ddmodel(dd_cos) + assert_ddmodel(dl.dd_cos) def test_dd_sphere(): - assert_ddmodel(dd_sphere) + assert_ddmodel(dl.dd_sphere) def test_dd_shell(): - assert_ddmodel(dd_shell) + assert_ddmodel(dl.dd_shell) def test_dd_shellshell(): - assert_ddmodel(dd_shellshell) + assert_ddmodel(dl.dd_shellshell) def test_dd_shellsphere(): - assert_ddmodel(dd_shellsphere) + assert_ddmodel(dl.dd_shellsphere) def test_dd_shellvoidsphere(): - assert_ddmodel(dd_shellvoidsphere) + assert_ddmodel(dl.dd_shellvoidsphere) def test_dd_shellvoidshell(): - assert_ddmodel(dd_shellvoidshell) + assert_ddmodel(dl.dd_shellvoidshell) def test_dd_spherepoint(): - assert_ddmodel(dd_spherepoint) + assert_ddmodel(dl.dd_spherepoint) def test_dd_spheresurf(): - assert_ddmodel(dd_spheresurf) + assert_ddmodel(dl.dd_spheresurf) def test_dd_triangle(): - assert_ddmodel(dd_triangle) + assert_ddmodel(dl.dd_triangle) def test_dd_uniform(): - assert_ddmodel(dd_uniform) + assert_ddmodel(dl.dd_uniform) def test_dd_wormchain(): - assert_ddmodel(dd_wormchain) + assert_ddmodel(dl.dd_wormchain) def test_dd_wormgauss(): - assert_ddmodel(dd_wormgauss) \ No newline at end of file + assert_ddmodel(dl.dd_wormgauss) \ No newline at end of file diff --git a/test/test_dipolarkernel.py b/test/test_dipolarkernel.py index 21415ab2a..535078ec5 100644 --- a/test/test_dipolarkernel.py +++ b/test/test_dipolarkernel.py @@ -406,5 +406,5 @@ def test_nonuniform_r(): P = dd_gauss(r,[3,0.5]) K = dipolarkernel(t,r) V0 = K@P - assert V0 == 1 + assert np.round(V0,3) == 1 #======================================================================= diff --git a/test/test_fitmultimodel.py b/test/test_fitmultimodel.py index 3fcb66311..b8ea50925 100644 --- a/test/test_fitmultimodel.py +++ b/test/test_fitmultimodel.py @@ -12,7 +12,7 @@ def test_multigauss(): r = np.linspace(2,6,300) t = np.linspace(-0.5,6,500) K = dipolarkernel(t,r) - parin = [4, 0.2, 0.4, 4, 1, 0.4, 3, 0.4, 0.2] + parin = [4, 0.05, 0.4, 4, 0.4, 0.4, 3, 0.15, 0.2] P = dd_gauss3(r,parin) V = K@P @@ -20,7 +20,8 @@ def test_multigauss(): assert ovl(P,fit.P) > 0.95 # more than 99% overlap #======================================================================= -test_multigauss() + + def test_multirice(): #======================================================================= "Check that the fit of a multi-Rician model works" @@ -28,7 +29,7 @@ def test_multirice(): r = np.linspace(2,6,300) t = np.linspace(0,6,500) K = dipolarkernel(t,r) - parin = [4, 0.2, 0.4, 4, 1, 0.4, 3, 0.4, 0.4] + parin = [4, 0.05, 0.4, 4, 0.4, 0.4, 3, 0.15, 0.2] P = dd_rice3(r,parin) V = K@P @@ -45,7 +46,7 @@ def test_multigengauss(): r = np.linspace(2,6,300) t = np.linspace(0,6,500) K = dipolarkernel(t,r) - P = dd_gengauss(r,[2.5, 0.5, 5]) + 0.8*dd_gengauss(r,[3, 0.7, 2]) + P = dd_gengauss(r,[2.5, 0.2, 5]) + 0.8*dd_gengauss(r,[3, 0.7, 2]) P /= np.trapz(P,r) V = K@P @@ -61,11 +62,11 @@ def test_bounds(): r = np.linspace(2,6,300) t = np.linspace(0,6,500) K = dipolarkernel(t,r) - parin = [4, 0.2, 0.4, 4, 1, 0.4, 3, 0.4, 0.2] + parin = [4, 0.05, 0.4, 4, 0.4, 0.4, 3, 0.15, 0.2] P = dd_gauss3(r,parin) V = K@P - fit = fitmultimodel(V,K,r,dd_gauss,3,'aicc',lb = [2,0.1],ub=[5.5,1.5], uqanalysis=False) + fit = fitmultimodel(V,K,r,dd_gauss,3,'aicc',lb = [2,0.02],ub=[5.5,1.5], uqanalysis=False) assert ovl(P,fit.P) > 0.95 # more than 99% overlap #======================================================================= @@ -92,7 +93,7 @@ def test_global_multigauss(): "Check that the global fit of a multi-Gauss model works" r = np.linspace(2,6,300) - parin = [4, 0.2, 0.4, 4, 1, 0.4, 3, 0.4, 0.2] + parin = [4, 0.05, 0.4, 4, 0.4, 0.4, 3, 0.15, 0.2] P = dd_gauss3(r,parin) t1 = np.linspace(-0.5,6,500) @@ -113,7 +114,7 @@ def test_global_multirice(): "Check that the global fit of a multi-Rician model works" r = np.linspace(2,6,300) - parin = [4, 0.2, 0.4, 4, 1, 0.4, 3, 0.4, 0.4] + parin = [4, 0.05, 0.4, 4, 0.4, 0.4, 3, 0.15, 0.2] P = dd_rice3(r,parin) t1 = np.linspace(-0.5,6,500) @@ -136,7 +137,7 @@ def test_background_fit(): t = np.linspace(-0.3,4,100) r = np.linspace(3,6,200) - InputParam = [4, 0.2, 0.5, 4.3, 0.3, 0.4] + InputParam = [4, 0.15, 0.5, 4.3, 0.1, 0.4] P = dd_gauss2(r,InputParam) B = bg_exp(t,0.15) V = dipolarkernel(t,r,0.25,B)@P @@ -147,7 +148,7 @@ def Kmodel(par): K = dipolarkernel(t,r,lam,B) return K - fit = fitmultimodel(V,Kmodel,r,dd_gauss,2,'aicc',lb=[1,0.1],ub=[6,1],lbK=[0.2,0.01],ubK=[0.9,1],uqanalysis=False) + fit = fitmultimodel(V,Kmodel,r,dd_gauss,2,'aicc',lb=[1,0.02],ub=[6,1],lbK=[0.2,0.01],ubK=[0.9,1],uqanalysis=False) assert ovl(P,fit.P) > 0.95 # more than 99% overlap #======================================================================= @@ -182,12 +183,12 @@ def test_confinter_Pfit(): t = np.linspace(-0.3,4,100) r = np.linspace(3,6,200) - InputParam = [4, 0.2, 0.5, 4.3, 0.3, 0.4] + InputParam = [4, 0.05, 0.5, 4.3, 0.1, 0.4] P = dd_gauss2(r,InputParam) K = dipolarkernel(t,r) V = K@P - fit = fitmultimodel(V,K,r,dd_gauss,3,'aicc',lb=[1,0.1],ub=[6,1]) + fit = fitmultimodel(V,K,r,dd_gauss,3,'aicc',lb=[1,0.02],ub=[6,1]) lbP = np.zeros(len(r)) ubP = np.full(len(r), np.inf) @@ -200,14 +201,14 @@ def test_confinter_parfit(): t = np.linspace(-0.3,4,100) r = np.linspace(3,6,200) - InputParam = [4, 0.2, 0.5, 4.3, 0.3, 0.4] + InputParam = [4, 0.05, 0.5, 4.3, 0.1, 0.4] P = dd_gauss2(r,InputParam) K = dipolarkernel(t,r) V = K@P - lbPpar = [1,0.1] + lbPpar = [1,0.02] ubPpar = [6,1] - fit = fitmultimodel(V,K,r,dd_gauss,3,'aicc',lb=[1,0.1],ub=[6,1]) + fit = fitmultimodel(V,K,r,dd_gauss,3,'aicc',lb=lbPpar,ub=ubPpar) paruq = fit.paramUncert parfit = fit.Pparam assert_confidence_intervals(paruq.ci(50)[0:2,:],paruq.ci(95)[0:2,:],parfit,lbPpar,ubPpar) @@ -221,7 +222,7 @@ def test_goodness_of_fit(): r = np.linspace(2,6,300) t = np.linspace(-0.5,6,500) K = dipolarkernel(t,r) - parin = [4, 0.2, 0.4, 4, 1, 0.4, 3, 0.4, 0.2] + parin = [4, 0.05, 0.4, 4, 0.4, 0.4, 3, 0.15, 0.2] P = dd_gauss3(r,parin) V = K@P @@ -236,7 +237,7 @@ def test_globalfit_scales(): "Check that the global fit of a multi-Gauss model works" r = np.linspace(2,6,300) - parin = [4, 0.2, 0.4, 4, 1, 0.4, 3, 0.4, 0.2] + parin = [4, 0.05, 0.4, 4, 0.4, 0.4, 3, 0.15, 0.2] P = dd_gauss3(r,parin) scales = [1e3, 1e9] diff --git a/test/test_fitparamodel.py b/test/test_fitparamodel.py index 9bb22869d..74b21e851 100644 --- a/test/test_fitparamodel.py +++ b/test/test_fitparamodel.py @@ -10,7 +10,7 @@ def test_gaussian(): t = np.linspace(0,3,300) r = np.linspace(2,5,200) - par = np.array([3,0.5]) + par = np.array([3,0.2]) P = dd_gauss(r,par) K = dipolarkernel(t,r) @@ -32,7 +32,7 @@ def test_rice(): t = np.linspace(0,3,400) r = np.linspace(1,6,150) - par = np.array([4, 0.4]) + par = np.array([4, 0.17]) P = dd_rice(r,par) K = dipolarkernel(t,r) @@ -55,7 +55,7 @@ def test_unconstrained(): t = np.linspace(0,3,300) r = np.linspace(2,5,200) - par = np.array([3,0.5]) + par = np.array([3,0.2]) P = dd_gauss(r,par) K = dipolarkernel(t,r) @@ -75,7 +75,7 @@ def test_rescaling(): t = np.linspace(0,5,100) r = np.linspace(2,6,100) - P = dd_gauss(r,[4, 0.5]) + P = dd_gauss(r,[4, 0.2]) K = dipolarkernel(t,r) scale = 1e9 @@ -122,7 +122,7 @@ def test_confinter_Pfit(): t = np.linspace(0,3,300) r = np.linspace(2,5,200) - par = np.array([3,0.5]) + par = np.array([3,0.2]) P = dd_gauss(r,par) K = dipolarkernel(t,r) @@ -145,18 +145,17 @@ def test_manual_covmatrix(): t = np.linspace(0,3,300) r = np.linspace(2,5,200) - par = np.array([3,0.5]) + par = np.array([3,0.2]) P = dd_gauss(r,par) K = dipolarkernel(t,r) - np.random.seed(1) sig = 0.01 - V = K@P + whitegaussnoise(t,sig) + V = K@P + whitegaussnoise(t,sig,seed=1) covmat = sig**2*np.eye(len(t)) - par0 = [5, 0.5] - lb = [2, 0.3] - ub = [5, 0.7] + par0 = [5, 0.50] + lb = [2, 0.05] + ub = [5, 0.70] model = lambda p: K@dd_gauss(r,p) fitmanual = fitparamodel(V,model,par0,lb,ub, covmatrix = covmat) fitauto = fitparamodel(V,model,par0,lb,ub) @@ -175,7 +174,7 @@ def test_globalfit(): "Check that global fitting yields correct results" r = np.linspace(2,5,200) - par = np.array([3,0.5]) + par = np.array([3,0.2]) P = dd_gauss(r,par) t1 = np.linspace(0,3,300) @@ -206,7 +205,7 @@ def test_globalfit_scales(): t1 = np.linspace(0,5,300) t2 = np.linspace(0,2,300) r = np.linspace(3,5,100) - P = dd_gauss(r,[4,0.6]) + P = dd_gauss(r,[4,0.25]) K1 = dipolarkernel(t1,r) K2 = dipolarkernel(t2,r) scales = [1e3, 1e9] diff --git a/test/test_fitregmodel.py b/test/test_fitregmodel.py index c12ad184b..43e2356f7 100644 --- a/test/test_fitregmodel.py +++ b/test/test_fitregmodel.py @@ -10,7 +10,7 @@ def assert_solver(regtype,solver): np.random.seed(1) t = np.linspace(-2,4,300) r = np.linspace(2,6,100) - P = dd_gauss(r,[3,0.5]) + P = dd_gauss(r,[3,0.2]) K = dipolarkernel(t,r) V = K@P + whitegaussnoise(t,0.01) fit = fitregmodel(V,K,r,regtype,'aic',solver=solver) @@ -89,7 +89,7 @@ def test_tikh_with_noise(): np.random.seed(1) t = np.linspace(0,3,200) r = np.linspace(1,5,100) - P = dd_gauss(r,[3,0.2]) + P = dd_gauss(r,[3,0.08]) K = dipolarkernel(t,r) V = K@P + whitegaussnoise(t,0.01) @@ -105,7 +105,7 @@ def test_tv_with_noise(): np.random.seed(1) t = np.linspace(0,3,200) r = np.linspace(1,5,100) - P = dd_gauss(r,[3,0.2]) + P = dd_gauss(r,[3,0.08]) K = dipolarkernel(t,r) V = K@P + whitegaussnoise(t,0.01) @@ -121,7 +121,7 @@ def test_huber_with_noise(): np.random.seed(1) t = np.linspace(0,3,200) r = np.linspace(1,5,100) - P = dd_gauss(r,[3,0.2]) + P = dd_gauss(r,[3,0.08]) K = dipolarkernel(t,r) V = K@P + whitegaussnoise(t,0.01) @@ -137,12 +137,12 @@ def test_unconstrained(): t = np.linspace(-2,4,300) r = np.linspace(2,6,100) - P = dd_gauss(r,[3,0.5]) + P = dd_gauss(r,[3,0.2]) K = dipolarkernel(t,r) L = regoperator(r,2) V = K@P - fit = fitregmodel(V,K,r,'tikhonov',alpha=0.05,nonnegativity=False,renormalize=False) + fit = fitregmodel(V,K,r,'tikhonov',regparam=0.05,nonnegativity=False,renormalize=False) Pfit_ref = np.linalg.solve(K.T@K + 0.05**2*L.T@L, K.T@V) assert ovl(Pfit_ref,fit.P) > 0.99 # more than 99% overlap @@ -153,7 +153,7 @@ def generate_global_dataset(): t1 = np.linspace(0,3,100) t2 = np.linspace(0,4,200) r = np.linspace(2.5,5,80) - P = dd_gauss2(r,[3.7,0.5,0.5,4.3,0.3,0.5]) + P = dd_gauss2(r,[3.7,0.2,0.5,4.3,0.1,0.5]) K1 = dipolarkernel(t1,r) K2 = dipolarkernel(t2,r) np.random.seed(1) @@ -190,7 +190,7 @@ def test_confinter_tikh(): t = np.linspace(-2,4,300) r = np.linspace(2,6,100) - P = dd_gauss(r,[3,0.5]) + P = dd_gauss(r,[3,0.2]) K = dipolarkernel(t,r) V = K@P @@ -205,7 +205,7 @@ def test_confinter_tv(): t = np.linspace(-2,4,300) r = np.linspace(2,6,100) - P = dd_gauss(r,[3,0.5]) + P = dd_gauss(r,[3,0.2]) K = dipolarkernel(t,r) V = K@P @@ -221,7 +221,7 @@ def test_confinter_huber(): t = np.linspace(-2,4,300) r = np.linspace(2,6,100) - P = dd_gauss(r,[3,0.5]) + P = dd_gauss(r,[3,0.2]) K = dipolarkernel(t,r) V = K@P @@ -247,7 +247,7 @@ def test_renormalize(): t = np.linspace(-2,4,300) r = np.linspace(2,6,100) - P = dd_gauss(r,[3,0.5]) + P = dd_gauss(r,[3,0.2]) K = dipolarkernel(t,r) V = K@P @@ -263,7 +263,7 @@ def test_scale_agnostic(): t = np.linspace(-2,4,300) r = np.linspace(2,6,100) - P = dd_gauss(r,[3,0.5]) + P = dd_gauss(r,[3,0.2]) K = dipolarkernel(t,r) scale = 1e8 V = scale*K@P @@ -279,7 +279,7 @@ def test_scale_fit(): t = np.linspace(-2,4,300) r = np.linspace(2,6,100) - P = dd_gauss(r,[3,0.5]) + P = dd_gauss(r,[3,0.2]) K = dipolarkernel(t,r) scale = 1e8 V = scale*K@P @@ -295,7 +295,7 @@ def test_nonuniform_r(): t = np.linspace(0,4,300) r = np.sqrt(np.linspace(1,7**2,200)) - P = dd_gauss(r,[3,0.5]) + P = dd_gauss(r,[3,0.2]) K = dipolarkernel(t,r) V = K@P @@ -312,14 +312,14 @@ def test_obir(): np.random.seed(1) t = np.linspace(0,3,200) r = np.linspace(2,7,100) - P = dd_gauss2(r,[3.7,0.5,0.5,4.3,0.3,0.5]) + P = dd_gauss2(r,[3.7,0.2,0.5,4.3,0.1,0.5]) K = dipolarkernel(t,r) V = K@P + whitegaussnoise(t,0.05) fit = fitregmodel(V,K,r,'tikhonov','aic',obir=True) - assert ovl(P,fit.P) > 0.80 # more than 80% overlap + assert ovl(P,fit.P) > 0.75 # more than 80% overlap #============================================================ def test_obir_global(): @@ -330,7 +330,7 @@ def test_obir_global(): fit = fitregmodel([V1,V2],[K1,K2],r,'tikhonov','aic',obir=True) - assert ovl(P,fit.P) > 0.80 # more than 80% overlap + assert ovl(P,fit.P) > 0.75 # more than 80% overlap #============================================================ @@ -341,7 +341,7 @@ def test_goodnes_of_fit(): np.random.seed(1) t = np.linspace(0,8,500) r = np.linspace(2,5,300) - P = dd_gauss(r,[4.5,0.5]) + P = dd_gauss(r,[4.5,0.2]) K = dipolarkernel(t,r) V = K@P + whitegaussnoise(t,0.01) @@ -359,7 +359,7 @@ def test_goodnes_of_fit_global(): t1 = np.linspace(0,3,100) t2 = np.linspace(0,4,200) r = np.linspace(2.5,5,80) - P = dd_gauss2(r,[3.7,0.5,0.5,4.3,0.3,0.5]) + P = dd_gauss2(r,[3.7,0.2,0.5,4.3,0.1,0.5]) K1 = dipolarkernel(t1,r) K2 = dipolarkernel(t2,r) np.random.seed(1) @@ -379,7 +379,7 @@ def test_globalfit_scales(): t1 = np.linspace(0,5,300) t2 = np.linspace(0,2,300) r = np.linspace(3,5,100) - P = dd_gauss(r,[4,0.6]) + P = dd_gauss(r,[4,0.25]) K1 = dipolarkernel(t1,r) K2 = dipolarkernel(t2,r) scales = [1e3, 1e9] diff --git a/test/test_fitsignal.py b/test/test_fitsignal.py index 9f8fb67ec..13f3d2cbb 100644 --- a/test/test_fitsignal.py +++ b/test/test_fitsignal.py @@ -12,7 +12,7 @@ def assert_experiment_model(model): # -------------------------------------------------------------------- t = np.linspace(0,5,100) r = np.linspace(2,6,150) - P = dd_gauss(r,[4.5, 0.6]) + P = dd_gauss(r,[4.5, 0.255]) info = model() parIn = info['Start'] @@ -62,7 +62,7 @@ def test_dipevo_function(): t = np.linspace(0,5,100) r = np.linspace(2,6,150) - P = dd_gauss(r,[4.5, 0.6]) + P = dd_gauss(r,[4.5, 0.25]) K = dipolarkernel(t,r) V = K@P fit = fitsignal(V,t,r,'P',None,None,uqanalysis=False) @@ -75,7 +75,7 @@ def test_form_factor(): t = np.linspace(0,5,100) r = np.linspace(2,6,150) - P = dd_gauss(r,[4.5, 0.6]) + P = dd_gauss(r,[4.5, 0.25]) K = dipolarkernel(t,r,0.3) V = K@P fit = fitsignal(V,t,r,'P',None,ex_4pdeer,uqanalysis=False) @@ -88,7 +88,7 @@ def test_full_parametric(): t = np.linspace(0,5,100) r = np.linspace(2,6,200) - P = dd_gauss(r,[4.5, 0.6]) + P = dd_gauss(r,[4.5, 0.25]) lam = 0.3 kappa = 0.2 B = bg_exp(t,lam*kappa) @@ -118,7 +118,7 @@ def test_start_values(): t = np.linspace(0,5,100) r = np.linspace(2,6,150) - P = dd_gauss(r,[4.5, 0.6]) + P = dd_gauss(r,[4.5, 0.25]) lam = 0.4 Bmodel = lambda t,lam: bg_exp(t,0.4,lam) @@ -137,7 +137,7 @@ def test_boundaries(): t = np.linspace(0,5,100) r = np.linspace(2,6,150) - P = dd_gauss(r,[4.5, 0.6]) + P = dd_gauss(r,[4.5, 0.25]) lam = 0.4 Bmodel = lambda t,lam: bg_exp(t,0.4,lam) @@ -158,7 +158,7 @@ def test_global_4pdeer(): "Check the correct fit of two 4-DEER signals" r = np.linspace(2,6,90) - P = dd_gauss(r,[4.5, 0.6]) + P = dd_gauss(r,[4.5, 0.3]) info = ex_4pdeer() parIn = info['Start'] @@ -182,7 +182,7 @@ def test_global_full_parametric(): "Check global fitting with fully parametric models" r = np.linspace(2,6,200) - P = dd_gauss(r,[4.5, 0.6]) + P = dd_gauss(r,[4.5, 0.25]) lam = 0.3 kappa = 0.4 B = lambda t,lam: bg_exp(t,lam*kappa,lam) @@ -201,7 +201,7 @@ def test_global_mixed_backgrounds(): "Check global fitting with different background models" r = np.linspace(2,6,200) - P = dd_gauss(r,[4.5, 0.6]) + P = dd_gauss(r,[4.5, 0.25]) lam = 0.3 kappa = 0.4 B = lambda t,lam: bg_exp(t,lam*kappa,lam) @@ -220,7 +220,7 @@ def test_global_mixed_experiments(): "Check global fitting with different experiment models" r = np.linspace(2,6,200) - P = dd_gauss(r,[4.5, 0.6]) + P = dd_gauss(r,[4.5, 0.25]) lam = 0.3 t1 = np.linspace(0,5,100) t2 = np.linspace(0,3,200) @@ -262,7 +262,7 @@ def assert_confinter_param(subset): bgmodel = bg_exp r = np.linspace(2,6,90) - P = ddmodel(r,[4.5, 0.6]) + P = ddmodel(r,[4.5, 0.25]) info = exmodel() parIn = info['Start'] @@ -277,17 +277,17 @@ def assert_confinter_param(subset): fit = fitsignal(V,t,r,ddmodel,bgmodel,exmodel,uqanalysis=True) - if subset is 'ex': + if subset == 'ex': info = exmodel() pfit = fit.exparam ci50 = fit.exparamUncert.ci(50) ci95 = fit.exparamUncert.ci(95) - if subset is 'dd': + if subset == 'dd': info = ddmodel() pfit = fit.ddparam ci50 = fit.ddparamUncert.ci(50) ci95 = fit.ddparamUncert.ci(95) - if subset is 'bg': + if subset == 'bg': info = bgmodel() pfit = fit.bgparam ci50 = fit.bgparamUncert.ci(50) @@ -320,7 +320,7 @@ def test_confinter_ddparam(): def assert_confinter_models(subset): #---------------------------------------------------------------------- exmodel = ex_4pdeer - if subset is 'Pfitfree': + if subset == 'Pfitfree': ddmodel = 'P' subset = 'Pfit' else: @@ -328,7 +328,7 @@ def assert_confinter_models(subset): bgmodel = bg_exp r = np.linspace(2,6,90) - P = dd_gauss(r,[4.5, 0.6]) + P = dd_gauss(r,[4.5, 0.25]) info = exmodel() parIn = info['Start'] @@ -343,19 +343,19 @@ def assert_confinter_models(subset): fit = fitsignal(V,t,r,ddmodel,bgmodel,exmodel,uqanalysis=True) - if subset is 'Pfit': + if subset == 'Pfit': modelfit = fit.P lb = np.zeros_like(r) ub = np.full_like(r,inf) ci50 = fit.Puncert.ci(50) ci95 = fit.Puncert.ci(95) - if subset is 'Bfit': + if subset == 'Bfit': modelfit = fit.B lb = np.full_like(t,-inf) ub = np.full_like(t,inf) ci50 = fit.Buncert.ci(50) ci95 = fit.Buncert.ci(95) - if subset is 'Vfit': + if subset == 'Vfit': modelfit = fit.V lb = np.full_like(t,-inf) ub = np.full_like(t,inf) @@ -396,7 +396,7 @@ def assert_confinter_noforeground(): bgmodel = bg_exp t = np.linspace(0,5,100) r = np.linspace(2,6,90) - P = dd_gauss(r,[4.5, 0.6]) + P = dd_gauss(r,[4.5, 0.25]) kappa = 0.4 lam = 0.3 @@ -422,7 +422,7 @@ def assert_confinter_dipevofun(): "Check that the confidence inervals for a dipolar evolution function fit are correct" r = np.linspace(2,6,90) - P = dd_gauss(r,[4.5, 0.6]) + P = dd_gauss(r,[4.5, 0.25]) t = np.linspace(0,5,100) np.random.seed(0) @@ -444,7 +444,7 @@ def test_global_scale_4pdeer(): # ====================================================================== "Check the correct fit of two 4-DEER signals" r = np.linspace(2,6,90) - P = dd_gauss(r,[4.5, 0.6]) + P = dd_gauss(r,[4.5, 0.25]) info = ex_4pdeer() parIn = info['Start'] diff --git a/test/test_selregparam.py b/test/test_selregparam.py index 690d5bf8c..82a8b35d0 100644 --- a/test/test_selregparam.py +++ b/test/test_selregparam.py @@ -8,7 +8,7 @@ def test_compensate_condition(): "Check that alpha compensates for larger condition numbers" r = np.linspace(2,6,100) - P = dd_gauss(r,[3,0.5]) + P = dd_gauss(r,[3,0.2]) # Lower condition number t1 = np.linspace(0,3,200) @@ -29,7 +29,7 @@ def get_alpha_from_method(method): t = np.linspace(0,5,500) r = np.linspace(2,5,80) - P = dd_gauss(r,[3,0.4]) + P = dd_gauss(r,[3,0.16986436005760383]) K = dipolarkernel(t,r) V = K@P @@ -182,7 +182,7 @@ def test_algorithms(): t = np.linspace(0,5,80) r = np.linspace(2,5,80) - P = dd_gauss(r,[3,0.4]) + P = dd_gauss(r,[3,0.16986436005760383]) K = dipolarkernel(t,r) V = K@P @@ -199,7 +199,7 @@ def test_nonuniform_r(): t = np.linspace(0,3,200) r = np.sqrt(np.linspace(1,7**2,200)) - P = dd_gauss(r,[3,0.5]) + P = dd_gauss(r,[3,0.2]) K = dipolarkernel(t,r) V = K@P @@ -219,7 +219,7 @@ def test_tikh_global(): t3 = np.linspace(0,4,70) r = np.linspace(2,5,80) - P = dd_gauss2(r,[3,0.4,0.3,3.5,0.4,0.7]) + P = dd_gauss2(r,[3,0.15,0.3,3.5,0.15,0.7]) K1 = dipolarkernel(t1,r) S1 = K1@P + whitegaussnoise(t1,0.03) @@ -242,7 +242,7 @@ def test_tv_global(): t3 = np.linspace(0,4,70) r = np.linspace(2,5,80) - P = dd_gauss2(r,[3,0.4,0.3,3.5,0.4,0.7]) + P = dd_gauss2(r,[3,0.15,0.3,3.5,0.15,0.7]) K1 = dipolarkernel(t1,r) S1 = K1@P + whitegaussnoise(t1,0.03) @@ -265,7 +265,7 @@ def test_huber_global(): t3 = np.linspace(0,4,70) r = np.linspace(2,5,80) - P = dd_gauss2(r,[3,0.4,0.3,3.5,0.4,0.7]) + P = dd_gauss2(r,[3,0.15,0.3,3.5,0.15,0.7]) K1 = dipolarkernel(t1,r) S1 = K1@P + whitegaussnoise(t1,0.03) @@ -319,7 +319,7 @@ def test_unconstrained(): t = np.linspace(0,5,80) r = np.linspace(2,5,80) - P = dd_gauss(r,[3,0.4]) + P = dd_gauss(r,[3,0.15]) K = dipolarkernel(t,r) V = K@P @@ -335,7 +335,7 @@ def test_manual_candidates(): t = np.linspace(0,5,80) r = np.linspace(2,5,80) - P = dd_gauss(r,[3,0.4]) + P = dd_gauss(r,[3,0.15]) K = dipolarkernel(t,r) L = regoperator(r,2) alphas = regparamrange(K,L) @@ -352,7 +352,7 @@ def get_alpha_from_regtype(regtype): np.random.seed(1) t = np.linspace(0,5,500) r = np.linspace(2,5,80) - P = dd_gauss(r,[3,0.4]) + P = dd_gauss(r,[3,0.15]) K = dipolarkernel(t,r) V = K@P + whitegaussnoise(t,0.01) diff --git a/test/test_snlls.py b/test/test_snlls.py index 3d60a3f3f..252aab03c 100644 --- a/test/test_snlls.py +++ b/test/test_snlls.py @@ -44,7 +44,7 @@ def Kmodel(p,t,r): ubl = [] # Separable LSQ fit - fit = snlls(V,lambda p: Kmodel(p,t,r),nlpar0,lb,ub,lbl,ubl, penalty=False, uqanalysis=False) + fit = snlls(V,lambda p: Kmodel(p,t,r),nlpar0,lb,ub,lbl,ubl, reg=False, uqanalysis=False) nonlinfit = fit.nonlin linfit = fit.lin parout = [nonlinfit[0], nonlinfit[1], linfit[0], nonlinfit[2], nonlinfit[3], linfit[1]] @@ -290,7 +290,7 @@ def test_goodnes_of_fit(): t = np.linspace(0,4,200) lam = 0.25 K = dipolarkernel(t,r,lam) - parin = [3.5, 0.4, 0.6, 4.5, 0.5, 0.4] + parin = [3.5, 0.15, 0.6, 4.5, 0.2, 0.4] P = dd_gauss2(r,parin) V = K@P diff --git a/test/test_whitegaussnoise.py b/test/test_whitegaussnoise.py new file mode 100644 index 000000000..513510796 --- /dev/null +++ b/test/test_whitegaussnoise.py @@ -0,0 +1,55 @@ +import numpy as np +from deerlab import whitegaussnoise + +def test_length(): +# ====================================================================== + "Check that the array length is correct" + + N = 151 + t = np.linspace(0,3,N) + noise = whitegaussnoise(t,0.2) + + assert len(noise) == N +# ====================================================================== + + +def test_rescale(): +# ====================================================================== + "Check that the rescale keyword is working" + + N = 10 + t = np.linspace(0,3,N) + sig = 1.1234 + noise = whitegaussnoise(t,sig,rescale=True) + + assert np.isclose(np.std(noise),sig) +# ====================================================================== + + +def test_seed(): +# ====================================================================== + "Check that the seed keyword is working" + + N = 151 + t = np.linspace(0,3,N) + sig = 0.2 + s = 432 + noise1 = whitegaussnoise(t,sig,seed=s) + noise2 = whitegaussnoise(t,sig,seed=s) + + assert np.array_equal(noise1,noise2) +# ====================================================================== + + +def test_noseed(): +# ====================================================================== + "Check that the seed keyword is working if set to None" + + N = 151 + t = np.linspace(0,3,N) + sig = 0.2 + noise1 = whitegaussnoise(t,sig,seed=None) + noise2 = whitegaussnoise(t,sig,seed=None) + + assert not np.array_equal(noise1,noise2) +# ======================================================================