From 32c97ffdeb0abca0b978658781382f226ac5d42d Mon Sep 17 00:00:00 2001 From: Hugo Karas Date: Wed, 11 Sep 2024 19:50:20 +0100 Subject: [PATCH 01/17] Added modelUncert for bootstrapped --- deerlab/fit.py | 9 +++++++-- deerlab/model.py | 2 +- 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/deerlab/fit.py b/deerlab/fit.py index 7c757649..09523c15 100644 --- a/deerlab/fit.py +++ b/deerlab/fit.py @@ -491,14 +491,19 @@ def bootstrap_fcn(ysim): else: bootstrap_verbose = False - param_uq = bootstrap_analysis(bootstrap_fcn,ysplit,fitresults.model,samples=bootstrap,noiselvl=noiselvl,cores=bootcores, verbose=bootstrap_verbose) + param_uq = bootstrap_analysis(bootstrap_fcn,ysplit,fitresults.model,samples=bootstrap,noiselvl=noiselvl,cores=bootcores, verbose=bootstrap_verbose, full_output=True) # Include information on the boundaries for better uncertainty estimates paramlb = model._vecsort(model._getvector('lb'))[np.concatenate(param_idx)] paramub = model._vecsort(model._getvector('ub'))[np.concatenate(param_idx)] fitresults.paramUncert = UQResult('bootstrap',data=param_uq[0].samples,lb=paramlb,ub=paramub) fitresults.param = fitresults.paramUncert.median + # Get the uncertainty estimates for the model response + modellb = np.min(param_uq[1].samples,axis=0) + modelub = np.max(param_uq[1].samples,axis=0) + fitresults.model = [param_uq[n].median for n in range(1,len(param_uq))] + fitresults.modelUncert = UQResult('bootstrap',data=param_uq[1].samples,lb=modellb,ub=modelub) if len(fitresults.model)==1: fitresults.model = fitresults.model[0] # Get some basic information on the parameter vector @@ -509,7 +514,7 @@ def bootstrap_fcn(ysim): # Dictionary of parameter names and fit uncertainties FitResult_paramuq = {f'{key}Uncert': model._getparamuq(fitresults.paramUncert,idx) for key,idx in zip(keys,param_idx)} # Dictionary of other fit quantities of interest - FitResult_dict = {key: getattr(fitresults,key) for key in ['y','mask','param','paramUncert','model','cost','plot','residuals','stats','regparam','regparam_stats','__plot_inputs']} + FitResult_dict = {key: getattr(fitresults,key) for key in ['y','mask','param','paramUncert','model','modelUncert','cost','plot','residuals','stats','regparam','regparam_stats','__plot_inputs']} _paramlist = model._parameter_list('vector') param_idx = [[] for _ in _paramlist] diff --git a/deerlab/model.py b/deerlab/model.py index 7ba709dd..f033da81 100644 --- a/deerlab/model.py +++ b/deerlab/model.py @@ -867,7 +867,7 @@ def __call__(self,*args,**kargs): # Check that all parameters have been passed if len(θ)!=self.Nparam: - raise SyntaxError(f'The model requires {self.Nparam} parameters, but {len(args_list)} were specified.') + raise SyntaxError(f'The model requires {self.Nparam} parameters, but {len(θ)} were specified.') # Determine which parameters are linear and which nonlinear θlin, θnonlin = self._split_linear(θ) From 6b6d6290bd30b6f024bcf9dcca0a5200f46837f9 Mon Sep 17 00:00:00 2001 From: Hugo Karas Date: Thu, 12 Sep 2024 06:59:52 +0100 Subject: [PATCH 02/17] Remove unused input --- deerlab/fit.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/deerlab/fit.py b/deerlab/fit.py index 09523c15..646e48a8 100644 --- a/deerlab/fit.py +++ b/deerlab/fit.py @@ -491,7 +491,7 @@ def bootstrap_fcn(ysim): else: bootstrap_verbose = False - param_uq = bootstrap_analysis(bootstrap_fcn,ysplit,fitresults.model,samples=bootstrap,noiselvl=noiselvl,cores=bootcores, verbose=bootstrap_verbose, full_output=True) + param_uq = bootstrap_analysis(bootstrap_fcn,ysplit,fitresults.model,samples=bootstrap,noiselvl=noiselvl,cores=bootcores, verbose=bootstrap_verbose) # Include information on the boundaries for better uncertainty estimates paramlb = model._vecsort(model._getvector('lb'))[np.concatenate(param_idx)] paramub = model._vecsort(model._getvector('ub'))[np.concatenate(param_idx)] From 9a0d715df2d6e0dcd4608084caec591bc50d4250 Mon Sep 17 00:00:00 2001 From: Hugo Karas Date: Fri, 13 Sep 2024 14:32:06 +0100 Subject: [PATCH 03/17] Always create modelUncert quantification --- deerlab/fit.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/deerlab/fit.py b/deerlab/fit.py index 646e48a8..1ffd6edb 100644 --- a/deerlab/fit.py +++ b/deerlab/fit.py @@ -361,7 +361,9 @@ def fit(model_, y, *constants, par0=None, penalties=None, bootstrap=0, noiselvl= paramlist : list List of the fitted parameter names ordered according to the model parameter indices. model : ndarray - Fitted model response. + Fitted model response. + modelUncert : :ref:`UQResult` + Uncertainty quantification of the fitted model response. regparam : scalar Regularization parameter value used for the regularization of the linear parameters. penweights : scalar or list thereof @@ -470,7 +472,7 @@ def fit(model_, y, *constants, par0=None, penalties=None, bootstrap=0, noiselvl= fitfcn = lambda y,penweights: snlls(y, Amodel_fcn, par0, lb=lb, ub=ub, lbl=lbl, ubl=ubl, mask=mask, weights=weights, subsets=ysubsets, lin_frozen=linfrozen, nonlin_frozen=nonlinfrozen, regparam=regparam, reg=reg, regparamrange=regparamrange, noiselvl=noiselvl, - extrapenalty=extrapenalties(penweights), **kwargs) + extrapenalty=extrapenalties(penweights), modeluq=True, **kwargs) # Prepare outer optimization of the penalty weights, if necessary fitfcn = _outerOptimization(fitfcn,penalties,sigmas) From 86eeecf10f6dff5ef5e5e50669a4f5f86238a97f Mon Sep 17 00:00:00 2001 From: Hugo Karas Date: Fri, 13 Sep 2024 14:32:27 +0100 Subject: [PATCH 04/17] Allow bootstrap resampling to be variable in propagation --- deerlab/fitresult.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/deerlab/fitresult.py b/deerlab/fitresult.py index 4841f9e6..6980c626 100644 --- a/deerlab/fitresult.py +++ b/deerlab/fitresult.py @@ -199,7 +199,7 @@ def evaluate(self, model, *constants): response = model(*constants,**parameters) return response - def propagate(self, model, *constants, lb=None, ub=None): + def propagate(self, model, *constants, lb=None, ub=None,samples=None): """ Propagate the uncertainty in the fit results to a model's response. @@ -223,7 +223,10 @@ def propagate(self, model, *constants, lb=None, ub=None): lb : array_like, optional Lower bounds of the model response. ub : array_like, optional - Upper bounds of the model response. + Upper bounds of the model response. + samples : int, optional + Number of samples to use when propagating a bootstraped uncertainty. If not provided, default value is 1000. + Returns ------- @@ -241,7 +244,7 @@ def propagate(self, model, *constants, lb=None, ub=None): # Propagate the uncertainty from that subset to the model - modeluq = self.paramUncert.propagate(lambda param: model(*constants,*[param[s] for s in fitparam_idx]),lb,ub) + modeluq = self.paramUncert.propagate(lambda param: model(*constants,*[param[s] for s in fitparam_idx]),lb,ub,samples) return modeluq From cb0368361cf4a97bedb58a362bf68762a41e128c Mon Sep 17 00:00:00 2001 From: Hugo Karas Date: Fri, 13 Sep 2024 14:32:43 +0100 Subject: [PATCH 05/17] Update examples --- examples/advanced/ex_profileanalysis.py | 2 +- examples/basic/ex_bootstrapping.py | 3 ++- examples/basic/ex_fitting_5pdeer.py | 2 +- examples/intermediate/ex_fitting_5pdeer_pathways.py | 2 +- 4 files changed, 5 insertions(+), 4 deletions(-) diff --git a/examples/advanced/ex_profileanalysis.py b/examples/advanced/ex_profileanalysis.py index 82d5fc42..6a195bd5 100644 --- a/examples/advanced/ex_profileanalysis.py +++ b/examples/advanced/ex_profileanalysis.py @@ -48,7 +48,7 @@ #%% # Extract fitted dipolar signal Vfit = results.model -Vci = results.propagate(Vmodel).ci(95) +Vci = results.modelUncert.ci(95) # Extract fitted distance distribution Pfit = results.P diff --git a/examples/basic/ex_bootstrapping.py b/examples/basic/ex_bootstrapping.py index de9a558f..196bd772 100644 --- a/examples/basic/ex_bootstrapping.py +++ b/examples/basic/ex_bootstrapping.py @@ -59,7 +59,7 @@ # Extract fitted dipolar signal Vfit = results.model -Vci = results.propagate(Vmodel).ci(95) +Vci = results.modelUncert.ci(95) # Extract fitted distance distribution Pfit = results.P @@ -78,6 +78,7 @@ plt.plot(t,Vexp,'.',color='grey',label='Data') # Plot the fitted signal plt.plot(t,Vfit,linewidth=3,label='Bootstrap median',color=violet) +plt.fill_between(t,Vci[:,0],Vci[:,1],linewidth=0.1,label='Bootstrap median',color=violet,alpha=0.3) plt.plot(t,Bfit,'--',linewidth=3,color=violet,label='Unmodulated contribution') plt.legend(frameon=False,loc='best') plt.xlabel('Time $t$ (μs)') diff --git a/examples/basic/ex_fitting_5pdeer.py b/examples/basic/ex_fitting_5pdeer.py index 83d7519d..432d3129 100644 --- a/examples/basic/ex_fitting_5pdeer.py +++ b/examples/basic/ex_fitting_5pdeer.py @@ -58,7 +58,7 @@ # Extract fitted dipolar signal Vfit = results.model -Vci = results.propagate(Vmodel).ci(95) +Vci = results.modelUncert.ci(95) # Extract fitted distance distribution Pfit = results.P diff --git a/examples/intermediate/ex_fitting_5pdeer_pathways.py b/examples/intermediate/ex_fitting_5pdeer_pathways.py index 370a7d93..516e50dc 100644 --- a/examples/intermediate/ex_fitting_5pdeer_pathways.py +++ b/examples/intermediate/ex_fitting_5pdeer_pathways.py @@ -47,7 +47,7 @@ # Extract fitted dipolar signal Vfit = results.model -Vci = results.propagate(Vmodel).ci(95) +Vci = results.modelUncert.ci(95) # Extract fitted distance distribution Pfit = results.P From 4ae39640c5947c2c21d5f635fd35cd75826dc20b Mon Sep 17 00:00:00 2001 From: Hugo Karas Date: Fri, 13 Sep 2024 14:33:02 +0100 Subject: [PATCH 06/17] Add test for modelUncert output --- test/test_model_class.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/test/test_model_class.py b/test/test_model_class.py index b8f6b864..473ecd3d 100644 --- a/test/test_model_class.py +++ b/test/test_model_class.py @@ -546,6 +546,23 @@ def test_fit_evaluate_model(mock_data,mock_x,model_type): response *= results.scale assert np.allclose(response,mock_data) + +# ================================================================ +@pytest.mark.parametrize('method', ['bootstrap','moment']) +@pytest.mark.parametrize('model_type', model_types) +def test_fit_modelUncert(mock_data,mock_x,model_type,method): + model = _generate_model(model_type, fixed_axis=False) + + if method=='bootstrap': + results = fit(model,mock_data,mock_x, bootstrap=3) + else: + results = fit(model,mock_data,mock_x) + + assert hasattr(results,'modelUncert') + ci_lower = results.modelUncert.ci(95)[:,0] + ci_upper = results.modelUncert.ci(95)[:,1] + assert np.less_equal(ci_lower,ci_upper).all() + # ================================================================ # ================================================================ From b8ac0eccb3aca94f5fd0956a953644f9be808c5f Mon Sep 17 00:00:00 2001 From: Hugo Karas Date: Fri, 13 Sep 2024 14:33:11 +0100 Subject: [PATCH 07/17] Minor docstring update --- deerlab/classes.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/deerlab/classes.py b/deerlab/classes.py index 39063cbf..d41e4482 100644 --- a/deerlab/classes.py +++ b/deerlab/classes.py @@ -522,7 +522,7 @@ def propagate(self,model,lb=None,ub=None,samples=None): ubm : ndarray Upper bounds of the values returned by ``model``, by default assumed unconstrained. samples : int, optional - Number of samples to use when propagating uncertainty. If not provided, default value is 1000. + Number of samples to use when propagating a bootstraped uncertainty. If not provided, default value is 1000. Returns ------- From 14668464a63904187da84047f937ee6524c69bb0 Mon Sep 17 00:00:00 2001 From: Hugo Karas Date: Fri, 13 Sep 2024 15:01:44 +0100 Subject: [PATCH 08/17] Bootrstap Uncertainty sampling reduction Reduces uncertainty calculation for bootstrapped uncertainties. Previously it was always 1000 samples. This saves around 4s, when using <100 samples. --- deerlab/fit.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/deerlab/fit.py b/deerlab/fit.py index 1ffd6edb..0c516c15 100644 --- a/deerlab/fit.py +++ b/deerlab/fit.py @@ -543,8 +543,8 @@ def _scale(x): FitResult_param_[f'{key}_scale'] = _scale(FitResult_param_[key]) # Normalization factor FitResult_param_[key] = param.normalization(FitResult_param_[key]) # Normalized value - FitResult_paramuq_[f'{key}_scaleUncert'] = FitResult_paramuq_[f'{key}Uncert'].propagate(_scale) - FitResult_paramuq_[f'{key}Uncert'] = FitResult_paramuq_[f'{key}Uncert'].propagate(lambda x: x/FitResult_param_[f'{key}_scale'], lb=param.lb, ub=param.ub) # Normalization of the uncertainty + FitResult_paramuq_[f'{key}_scaleUncert'] = FitResult_paramuq_[f'{key}Uncert'].propagate(_scale,samples=bootstrap) + FitResult_paramuq_[f'{key}Uncert'] = FitResult_paramuq_[f'{key}Uncert'].propagate(lambda x: x/FitResult_param_[f'{key}_scale'], lb=param.lb, ub=param.ub,samples=bootstrap) # Normalization of the uncertainty if len(noiselvl)==1: noiselvl = noiselvl[0] From 691cdf0bd8cfb925545e3c8b2546af2f4f6de192 Mon Sep 17 00:00:00 2001 From: Hugo Karas Date: Mon, 16 Sep 2024 12:55:37 +0100 Subject: [PATCH 09/17] Add model copying --- deerlab/model.py | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/deerlab/model.py b/deerlab/model.py index f033da81..5e4866a3 100644 --- a/deerlab/model.py +++ b/deerlab/model.py @@ -218,6 +218,13 @@ def unfreeze(self): self.frozen = False self.value = None #--------------------------------------------------------------------------------------- + def copy(self): + """ + Return a deep copy of the parameter + """ + + return deepcopy(self) + #=================================================================================== @@ -973,6 +980,13 @@ def __str__(self): """ return self._parameter_table() #--------------------------------------------------------------------------------------- + def copy(self): + """ + Return a deep copy of the model + """ + + return deepcopy(self) + #=================================================================================== #============================================================================== From 2b7e527c50e34f97455e49a598b6fdc46eb78c58 Mon Sep 17 00:00:00 2001 From: Hugo Karas Date: Mon, 16 Sep 2024 12:56:06 +0100 Subject: [PATCH 10/17] General improvemements --- deerlab/classes.py | 2 +- deerlab/fit.py | 4 ++-- deerlab/fitresult.py | 31 +++++++++++++++++-------------- 3 files changed, 20 insertions(+), 17 deletions(-) diff --git a/deerlab/classes.py b/deerlab/classes.py index d41e4482..49a6ec46 100644 --- a/deerlab/classes.py +++ b/deerlab/classes.py @@ -583,7 +583,7 @@ def propagate(self,model,lb=None,ub=None,samples=None): # Get the parameter uncertainty distribution values,pdf = self.pardist(n) # Random sampling form the uncertainty distribution - sampled_parameters[n] = [np.random.choice(values, p=pdf/sum(pdf)) for _ in range(Nsamples)] + sampled_parameters[n] = np.random.choice(values, p=pdf/sum(pdf),size=Nsamples) # Convert to matrix sampled_parameters = np.atleast_2d(sampled_parameters) diff --git a/deerlab/fit.py b/deerlab/fit.py index 0c516c15..a1b5e8bc 100644 --- a/deerlab/fit.py +++ b/deerlab/fit.py @@ -543,8 +543,8 @@ def _scale(x): FitResult_param_[f'{key}_scale'] = _scale(FitResult_param_[key]) # Normalization factor FitResult_param_[key] = param.normalization(FitResult_param_[key]) # Normalized value - FitResult_paramuq_[f'{key}_scaleUncert'] = FitResult_paramuq_[f'{key}Uncert'].propagate(_scale,samples=bootstrap) - FitResult_paramuq_[f'{key}Uncert'] = FitResult_paramuq_[f'{key}Uncert'].propagate(lambda x: x/FitResult_param_[f'{key}_scale'], lb=param.lb, ub=param.ub,samples=bootstrap) # Normalization of the uncertainty + FitResult_paramuq_[f'{key}_scaleUncert'] = FitResult_paramuq_[f'{key}Uncert'].propagate(_scale,samples=bootstrap+1) + FitResult_paramuq_[f'{key}Uncert'] = FitResult_paramuq_[f'{key}Uncert'].propagate(lambda x: x/FitResult_param_[f'{key}_scale'], lb=param.lb, ub=param.ub,samples=bootstrap+1) # Normalization of the uncertainty if len(noiselvl)==1: noiselvl = noiselvl[0] diff --git a/deerlab/fitresult.py b/deerlab/fitresult.py index 6980c626..c1a552b0 100644 --- a/deerlab/fitresult.py +++ b/deerlab/fitresult.py @@ -97,20 +97,21 @@ def _extarct_params_from_model(self, model): if not hasattr(self,'param'): raise ValueError('The fit object does not contain any fitted parameters.') - # # Enforce model normalization - # normfactor_keys = [] - # for key in modelparam: - # param = getattr(model,key) - # if np.all(param.linear): - # if param.normalization is not None: - # normfactor_key = f'{key}_scale' - # normfactor_keys.append(normfactor_key) - # try: - # model.addnonlinear(normfactor_key,lb=-np.inf,ub=np.inf,par0=1,description=f'Normalization factor of {key}') - # getattr(model,normfactor_key).freeze(1) - # except KeyError: - # pass - + # Enforce model normalization + normfactor_keys = [] + for key in modelparam: + param = getattr(model,key) + if np.all(param.linear): + if param.normalization is not None: + normfactor_key = f'{key}_scale' + normfactor_keys.append(normfactor_key) + try: + model.addnonlinear(normfactor_key,lb=-np.inf,ub=np.inf,par0=1,description=f'Normalization factor of {key}') + getattr(model,normfactor_key).freeze(1) + except KeyError: + pass + modelparam += normfactor_keys + # # Get some basic information on the parameter vector # modelparam = model._parameter_list(order='vector') @@ -187,6 +188,7 @@ def evaluate(self, model, *constants): Model response at the fitted parameter values. """ try: + model = model.copy() modelparam = model._parameter_list('vector') modelparam, fitparams, fitparam_idx = self._extarct_params_from_model(model) except AttributeError: @@ -234,6 +236,7 @@ def propagate(self, model, *constants, lb=None, ub=None,samples=None): responseUncert : :ref:`UQResult` Uncertainty quantification of the model's response. """ + model = model.copy() try: modelparam = model._parameter_list('vector') From 76cfed180d63ec17a3d5913837cd96484444294b Mon Sep 17 00:00:00 2001 From: Hugo Karas Date: Tue, 26 Nov 2024 17:07:46 +0100 Subject: [PATCH 11/17] Fixes issues with functions being deep copyied --- deerlab/fitresult.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/deerlab/fitresult.py b/deerlab/fitresult.py index c1a552b0..27d33aa9 100644 --- a/deerlab/fitresult.py +++ b/deerlab/fitresult.py @@ -236,9 +236,10 @@ def propagate(self, model, *constants, lb=None, ub=None,samples=None): responseUncert : :ref:`UQResult` Uncertainty quantification of the model's response. """ - model = model.copy() + try: + model = model.copy() modelparam = model._parameter_list('vector') modelparam, fitparams, fitparam_idx = self._extarct_params_from_model(model) From 6e41ffea67f68b434af5b4b87699a014a6c84afc Mon Sep 17 00:00:00 2001 From: Hugo Karas Date: Wed, 27 Nov 2024 10:41:42 +0100 Subject: [PATCH 12/17] Test update -`test_fit_model_confidence_intervals` increased the number of bootstraps so it passes -`test_fit_modelUncert` expanded to test the bootstrapped --- test/test_model_class.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/test/test_model_class.py b/test/test_model_class.py index 473ecd3d..8fe63611 100644 --- a/test/test_model_class.py +++ b/test/test_model_class.py @@ -42,7 +42,9 @@ def mock_x(): @pytest.fixture(scope='module') def mock_data(mock_x): - return bigauss(mock_x,mean1=3,mean2=4,std1=0.5,std2=0.2,amp1=0.5,amp2=0.6) + data = bigauss(mock_x,mean1=3,mean2=4,std1=0.5,std2=0.2,amp1=0.5,amp2=0.6) + data += whitegaussnoise(mock_x,0.01,seed=1) + return data model_types = ['parametric','semiparametric','semiparametric_vec', 'nonparametric','nonparametric_vec','nonparametric_vec_normalized'] @@ -511,7 +513,7 @@ def test_fit_model_confidence_intervals(mock_data,model_type,method): model = _generate_model(model_type, fixed_axis=True) if method=='bootstrap': - results = fit(model,mock_data + whitegaussnoise(x,0.01,seed=1), bootstrap=3) + results = fit(model,mock_data + whitegaussnoise(x,0.01,seed=1), bootstrap=100) else: results = fit(model,mock_data + whitegaussnoise(x,0.01,seed=1)) @@ -551,6 +553,7 @@ def test_fit_evaluate_model(mock_data,mock_x,model_type): @pytest.mark.parametrize('method', ['bootstrap','moment']) @pytest.mark.parametrize('model_type', model_types) def test_fit_modelUncert(mock_data,mock_x,model_type,method): + "Check that the uncertainty of fit results can be calculated and is the uncertainty of the model is non zero for all but nonparametric models" model = _generate_model(model_type, fixed_axis=False) if method=='bootstrap': @@ -562,6 +565,9 @@ def test_fit_modelUncert(mock_data,mock_x,model_type,method): ci_lower = results.modelUncert.ci(95)[:,0] ci_upper = results.modelUncert.ci(95)[:,1] assert np.less_equal(ci_lower,ci_upper).all() + if model_type != 'nonparametric' and model_type != 'nonparametric_vec' and model_type != 'nonparametric_vec_normalized': + assert np.all(np.round(ci_lower) <= np.round(results.model)) and np.less(ci_lower.sum(),results.model.sum()) + assert np.all(np.round(ci_upper,5) >= np.round(results.model,5)) and np.greater(ci_upper.sum(),results.model.sum()) # ================================================================ From d70b721d2100cb058eea63a911cc20e671276dfc Mon Sep 17 00:00:00 2001 From: Hugo Karas Date: Wed, 27 Nov 2024 10:42:07 +0100 Subject: [PATCH 13/17] Correct number of bootstrap samples --- deerlab/fit.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/deerlab/fit.py b/deerlab/fit.py index a1b5e8bc..71d26231 100644 --- a/deerlab/fit.py +++ b/deerlab/fit.py @@ -493,7 +493,7 @@ def bootstrap_fcn(ysim): else: bootstrap_verbose = False - param_uq = bootstrap_analysis(bootstrap_fcn,ysplit,fitresults.model,samples=bootstrap,noiselvl=noiselvl,cores=bootcores, verbose=bootstrap_verbose) + param_uq = bootstrap_analysis(bootstrap_fcn,ysplit,fitresults.model,samples=bootstrap-1,noiselvl=noiselvl,cores=bootcores, verbose=bootstrap_verbose) # Include information on the boundaries for better uncertainty estimates paramlb = model._vecsort(model._getvector('lb'))[np.concatenate(param_idx)] paramub = model._vecsort(model._getvector('ub'))[np.concatenate(param_idx)] @@ -543,8 +543,8 @@ def _scale(x): FitResult_param_[f'{key}_scale'] = _scale(FitResult_param_[key]) # Normalization factor FitResult_param_[key] = param.normalization(FitResult_param_[key]) # Normalized value - FitResult_paramuq_[f'{key}_scaleUncert'] = FitResult_paramuq_[f'{key}Uncert'].propagate(_scale,samples=bootstrap+1) - FitResult_paramuq_[f'{key}Uncert'] = FitResult_paramuq_[f'{key}Uncert'].propagate(lambda x: x/FitResult_param_[f'{key}_scale'], lb=param.lb, ub=param.ub,samples=bootstrap+1) # Normalization of the uncertainty + FitResult_paramuq_[f'{key}_scaleUncert'] = FitResult_paramuq_[f'{key}Uncert'].propagate(_scale,samples=bootstrap) + FitResult_paramuq_[f'{key}Uncert'] = FitResult_paramuq_[f'{key}Uncert'].propagate(lambda x: x/FitResult_param_[f'{key}_scale'], lb=param.lb, ub=param.ub,samples=bootstrap) # Normalization of the uncertainty if len(noiselvl)==1: noiselvl = noiselvl[0] From 744baf88350155ff1eea78b37c173b6e61d1f2a9 Mon Sep 17 00:00:00 2001 From: Hugo Karas Date: Wed, 27 Nov 2024 11:27:51 +0100 Subject: [PATCH 14/17] Fix gaussian normalisation issues --- deerlab/dd_models.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/deerlab/dd_models.py b/deerlab/dd_models.py index 1c4d6d3f..c8f4f132 100644 --- a/deerlab/dd_models.py +++ b/deerlab/dd_models.py @@ -99,7 +99,7 @@ def docstr_example(fcnstr): # ================================================================= def _normalize(r,P): if not all(P==0): - P = P/np.trapz(P,r) + P = P/np.trapezoid(P,r) return P # ================================================================= @@ -110,7 +110,7 @@ def _multigaussfun(r,r0,sig): P = np.sqrt(1/(2*np.pi))*1/sig*np.exp(-0.5*((r.T-r0)/sig)**2) if not np.all(P==0): # Normalization - P = np.squeeze(P)/np.sum([np.trapezoid(c,r) for c in P.T]) + P = np.squeeze(P)/np.array([np.trapezoid(c,r) for c in P.T]).reshape(1, -1) else: P = np.squeeze(P) return P From 80a3b86f52792e4ecbb7e7f9d0a561b11bb03a5c Mon Sep 17 00:00:00 2001 From: Hugo Karas Date: Wed, 27 Nov 2024 15:25:57 +0100 Subject: [PATCH 15/17] Bug fix in test --- test/test_model_class.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/test/test_model_class.py b/test/test_model_class.py index 8fe63611..1a4c8a98 100644 --- a/test/test_model_class.py +++ b/test/test_model_class.py @@ -42,6 +42,11 @@ def mock_x(): @pytest.fixture(scope='module') def mock_data(mock_x): + data = bigauss(mock_x,mean1=3,mean2=4,std1=0.5,std2=0.2,amp1=0.5,amp2=0.6) + return data + +@pytest.fixture(scope='module') +def mock_data_noise(mock_x): data = bigauss(mock_x,mean1=3,mean2=4,std1=0.5,std2=0.2,amp1=0.5,amp2=0.6) data += whitegaussnoise(mock_x,0.01,seed=1) return data @@ -552,14 +557,14 @@ def test_fit_evaluate_model(mock_data,mock_x,model_type): # ================================================================ @pytest.mark.parametrize('method', ['bootstrap','moment']) @pytest.mark.parametrize('model_type', model_types) -def test_fit_modelUncert(mock_data,mock_x,model_type,method): +def test_fit_modelUncert(mock_data_noise,mock_x,model_type,method): "Check that the uncertainty of fit results can be calculated and is the uncertainty of the model is non zero for all but nonparametric models" model = _generate_model(model_type, fixed_axis=False) if method=='bootstrap': - results = fit(model,mock_data,mock_x, bootstrap=3) + results = fit(model,mock_data_noise,mock_x, bootstrap=3) else: - results = fit(model,mock_data,mock_x) + results = fit(model,mock_data_noise,mock_x) assert hasattr(results,'modelUncert') ci_lower = results.modelUncert.ci(95)[:,0] From 93c6c8ebe0bcb8f8ee90a00e001230dde7daaac5 Mon Sep 17 00:00:00 2001 From: Hugo Karas Date: Thu, 9 Jan 2025 14:00:30 +0100 Subject: [PATCH 16/17] multi-guass_background fix removal --- deerlab/dd_models.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/deerlab/dd_models.py b/deerlab/dd_models.py index c8f4f132..e590f7c2 100644 --- a/deerlab/dd_models.py +++ b/deerlab/dd_models.py @@ -110,7 +110,7 @@ def _multigaussfun(r,r0,sig): P = np.sqrt(1/(2*np.pi))*1/sig*np.exp(-0.5*((r.T-r0)/sig)**2) if not np.all(P==0): # Normalization - P = np.squeeze(P)/np.array([np.trapezoid(c,r) for c in P.T]).reshape(1, -1) + P = np.squeeze(P)/np.array([np.trapezoid(c,r) for c in P.T]).sum() else: P = np.squeeze(P) return P From 69a370c3da994d5996e101028e9a0b17d8d0cf80 Mon Sep 17 00:00:00 2001 From: Hugo Karas Date: Thu, 9 Jan 2025 15:03:42 +0100 Subject: [PATCH 17/17] Fix example for latest matplotlib --- examples/advanced/ex_dipolarpathways_selection.py | 5 +---- setup.py | 2 +- 2 files changed, 2 insertions(+), 5 deletions(-) diff --git a/examples/advanced/ex_dipolarpathways_selection.py b/examples/advanced/ex_dipolarpathways_selection.py index 57fb0e18..bf255174 100644 --- a/examples/advanced/ex_dipolarpathways_selection.py +++ b/examples/advanced/ex_dipolarpathways_selection.py @@ -11,7 +11,6 @@ import numpy as np import matplotlib.pyplot as plt from matplotlib.gridspec import GridSpec -from mpl_toolkits.axes_grid1.inset_locator import inset_axes, InsetPosition import deerlab as dl @@ -79,9 +78,7 @@ Pfit = fits[n].P Pci = fits[n].PUncert.ci(95) # Setup the inset plot - axins = inset_axes(ax1,width="30%", height="30%", loc='upper left') - ip = InsetPosition(ax1,[0.35, 0.17+0.24*n, 0.6, 0.1]) - axins.set_axes_locator(ip) + axins = ax1.inset_axes([0.35, 0.17+0.24*n, 0.6, 0.1]) axins.yaxis.set_ticklabels([]) axins.yaxis.set_visible(False) # Plot the distance distributions and their confidence bands diff --git a/setup.py b/setup.py index a0747967..5d237fc5 100644 --- a/setup.py +++ b/setup.py @@ -25,7 +25,7 @@ 'joblib>=1.0.0', 'dill>=0.3.0', 'tqdm>=4.51.0', - 'matplotlib>=3.3.4', + 'matplotlib>=3.6.0', 'memoization>=0.3.1', 'pytest>=6.2.2', 'setuptools>=53.0.0',