diff --git a/deerlab/fitmultimodel.py b/deerlab/fitmultimodel.py index b777aaf7e..ec302f2fb 100644 --- a/deerlab/fitmultimodel.py +++ b/deerlab/fitmultimodel.py @@ -207,14 +207,14 @@ def Kmodel(Kpar): try: Kmodel(np.random.uniform(size=nKparam)) notEnoughParam = False - except ValueError: + except: notEnoughParam = True else: # If the kernel is just a matrix make it a callable without parameters nKparam = 0 K = copy.deepcopy(Kmodel) # need a copy to avoid infite recursion on next step Kmodel = lambda _: K - + # Extract information about the model nparam = len(model.start) if lb is None: @@ -231,7 +231,7 @@ def Kmodel(Kpar): # Ensure that all arrays are numpy.nparray lb,ub,lbK,ubK = np.atleast_1d(lb,ub,lbK,ubK) - if len(lbK) is not nKparam or len(ubK) is not nKparam: + if len(lbK)!=nKparam or len(ubK)!=nKparam: raise ValueError('The upper/lower bounds of the kernel parameters must be ',nKparam,'-element arrays') areLocations = [str in ['Mean','Location'] for str in paramnames] @@ -264,6 +264,7 @@ def nonlinmodel(par,Nmodels): Pbasis = model(r,par[subset]) # Combine all non-linear functions into one Knonlin[:,iModel] = K@Pbasis + Knonlin = weights[:,np.newaxis]*Knonlin return Knonlin #=============================================================================== @@ -423,16 +424,16 @@ def initialize_nonlin_param(par,Ncomp,lb,ub,strategy): lin_ub = np.full(Ncomp,np.inf) # Separable non-linear least-squares (SNLLS) fit - scale = 1e2 - fit = dl.snlls(V*scale,Knonlin,par0,nlin_lb,nlin_ub,lin_lb,lin_ub, + upscale = 1e2 + fit = dl.snlls(V*upscale,Knonlin,par0,nlin_lb,nlin_ub,lin_lb,lin_ub, weights=weights, reg=False, nonlin_tol=tol, nonlin_maxiter=maxiter) pnonlin = fit.nonlin plin = fit.lin par_prev = pnonlin - plin /= scale - fit.model /= scale - fit.modelUncert = fit.modelUncert.propagate(lambda x: x/scale) + plin /= upscale + fit.model /= upscale + fit.modelUncert = fit.modelUncert.propagate(lambda x: x/upscale) # Store the fitted parameters pnonlin_.append(pnonlin) @@ -501,50 +502,76 @@ def initialize_nonlin_param(par,Ncomp,lb,ub,strategy): Puq = dl.UQResult('void') paramuq = dl.UQResult('void') - # Goodness of fit - stats = [] - for subset in Vsubsets: - Ndof = len(V[subset]) - (nKparam + nparam + Nopt) - stats.append(goodness_of_fit(V[subset],Vfit[subset],Ndof)) - if len(stats)==1: - stats = stats[0] - # If requested re-normalize the distribution postscale = np.trapz(Pfit,r) if renormalize: Pfit = Pfit/postscale - fitparam_amp = fitparam_amp/sum(fitparam_amp) if uq: Puq_ = copy.deepcopy(Puq) # need a copy to avoid infite recursion on next step Puq = Puq_.propagate(lambda P: P/postscale, lbm=np.zeros_like(Pfit)) + postscale /= sum(fitparam_amp) + fitparam_amp = fitparam_amp/sum(fitparam_amp) # Dataset scales scales = [] for i in range(len(Vsubsets)): scales.append(prescales[i]*postscale) - if len(scales)==1: - scales = scales[0] + + # Get fitted signals and their uncertainty + modelfit, modelfituq = [],[] + for i,subset in enumerate(Vsubsets): + V[subset] = V[subset]*prescales[i] + modelfit.append( scales[i]*fit.model[subset] ) + if uq: + modelfituq.append( fit.modelUncert.propagate(lambda V: scales[i]*V[subset]) ) + else: + modelfituq.append(dl.UQResult('void')) + + # Goodness of fit + stats = [] + for i,subset in enumerate(Vsubsets): + Ndof = len(V[subset]) - (nKparam + nparam + Nopt) + stats.append(goodness_of_fit(V[subset],modelfit[i],Ndof)) # Results display function - def plotfcn(show=False): - fig = _plot(Vsubsets,V,Vfit,r,Pfit,Puq,fcnals,maxModels,method,uq,show) + def plotfcn(show=True): + fig = _plot(Vsubsets,V,modelfit,modelfituq,r,Pfit,Puq,fcnals,maxModels,method,uq,show) return fig - return FitResult(P=Pfit, Pparam=fitparam_P, Kparam=fitparam_K, amps=fitparam_amp, V=fit.model, Puncert=Puq, - paramUncert=paramuq, Vuncert=fit.modelUncert, selfun=fcnals, Nopt=Nopt, Pn=Peval, scale=scales, plot=plotfcn, + # If just one dataset, return vector instead of list + if len(Vsubsets)==1: + stats = stats[0] + scales = scales[0] + modelfit = modelfit[0] + modelfituq = modelfituq[0] + + + return FitResult(P=Pfit, Pparam=fitparam_P, Kparam=fitparam_K, amps=fitparam_amp, V=modelfit, Puncert=Puq, + paramUncert=paramuq, Vuncert=modelfituq, selfun=fcnals, Nopt=Nopt, Pn=Peval, scale=scales, plot=plotfcn, stats=stats, cost=fit.cost, residuals=fit.residuals, success=fit.success) # ========================================================================= -def _plot(Vsubsets,V,Vfit,r,Pfit,Puq,fcnals,maxModels,method,uq,show): +def _plot(Vsubsets,V,Vfit,Vuq,r,Pfit,Puq,fcnals,maxModels,method,uq,show): # ========================================================================= + if not isinstance(Vfit, list): + Vuq = [Vuq] + Vfit = [Vfit] + nSignals = len(Vsubsets) fig,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].plot(Vfit[i],'tab:blue') + if uq: + # Confidence intervals of the fitted datasets + Vci95 = Vuq[i].ci(95) # 95#-confidence interval + Vci50 = Vuq[i].ci(50) # 50#-confidence interval + tax = np.arange(len(subset)) + axs[i].fill_between(tax,Vci50[:,0],Vci50[:,1],color='tab:blue',linestyle='None',alpha=0.45) + axs[i].fill_between(tax,Vci95[:,0],Vci95[:,1],color='tab:blue',linestyle='None',alpha=0.25) axs[i].grid(alpha=0.3) axs[i].set_xlabel('Array Elements') axs[i].set_ylabel(f'V[{i}]') diff --git a/deerlab/fitregmodel.py b/deerlab/fitregmodel.py index 0f89baacd..e3171a124 100644 --- a/deerlab/fitregmodel.py +++ b/deerlab/fitregmodel.py @@ -159,6 +159,8 @@ def fitregmodel(V, K, r, regtype='tikhonov', regparam='aic', regorder=2, solver= else: prescales = [1] + + # Determine an optimal value of the regularization parameter if requested if type(regparam) is str: alpha = dl.selregparam(V,K,r,regtype,regparam,regorder=regorder,weights=weights, nonnegativity=nonnegativity,huberparam=huberparam) @@ -200,19 +202,20 @@ def fitregmodel(V, K, r, regtype='tikhonov', regparam='aic', regorder=2, solver= # Get fit final status Vfit = K@Pfit success = ~np.all(Pfit==0) - res = V - Vfit - fval = np.linalg.norm(V - Vfit)**2 + alpha**2*np.linalg.norm(L@Pfit)**2 + + # Construct residual parts for for the residual and regularization terms + res = weights*(V - K@Pfit) + + # Construct Jacobians for the residual and penalty terms + Jres = K*weights[:,np.newaxis] + res,J = _augment(res,Jres,regtype,alpha,L,Pfit,huberparam) + + # Get objective function value + fval = np.linalg.norm(res)**2 # Uncertainty quantification # ---------------------------------------------------------------- if uq: - # Construct residual parts for for the residual and regularization terms - res = weights*(V - K@Pfit) - - # Construct Jacobians for the residual and penalty terms - Jres = K*weights[:,np.newaxis] - res,J = _augment(res,Jres,regtype,alpha,L,Pfit,huberparam) - # Calculate the heteroscedasticity consistent covariance matrix covmat = hccm(J,res,'HC1') diff --git a/deerlab/lsqcomponents.py b/deerlab/lsqcomponents.py index ab25cfa36..b9179d22f 100644 --- a/deerlab/lsqcomponents.py +++ b/deerlab/lsqcomponents.py @@ -1,14 +1,19 @@ import numpy as np -def lsqcomponents(V, K, L=None, alpha=0, weights=1, regtype='tikhonov', huberparam=1.35): +def lsqcomponents(V, K, L=None, alpha=0, weights=None, regtype='tikhonov', huberparam=1.35): # ============================================================================================== """ Calculate the components needed for the least-squares (LSQ) solvers """ + if weights is None: + weights = np.ones_like(V) + else: + weights = np.atleast_1d(weights) + # Compute components of the LSQ normal equations - KtK = weights*K.T@K - KtV = weights*K.T@V + KtK = K.T@(weights[:,np.newaxis]*K) + KtV = K.T@(weights*V) # No regularization term -> done if L is None: @@ -39,7 +44,7 @@ def optimizeregterm(regfun, threshold, maxIter): # Compute the regularization term if regtype.lower() == 'tikhonov': regterm = L.T@L - + elif regtype.lower() == 'tv': maxIter = 500 changeThreshold = 1e-1 diff --git a/deerlab/selregparam.py b/deerlab/selregparam.py index 1fcf051e7..6ade73522 100644 --- a/deerlab/selregparam.py +++ b/deerlab/selregparam.py @@ -204,7 +204,8 @@ def _evalalpha(alpha,V,K,L,selmethod,nonneg,noiselvl,regtype,weights,HuberParame H = K@pK # Residual term - Residual = np.linalg.norm(K@P - V) + residuals = weights*(K@P - V) + Residual = np.linalg.norm(residuals) # Regularization penalty term if regtype.lower() == 'tikhonov': Penalty = np.linalg.norm(L@P) @@ -223,7 +224,7 @@ def _evalalpha(alpha,V,K,L,selmethod,nonneg,noiselvl,regtype,weights,HuberParame # Cross validation (CV) if selmethod =='cv': - f_ = sum(abs((V - K@P)/(np.ones(N) - np.diag(H)))**2) + f_ = sum(abs(residuals/(np.ones(N) - np.diag(H)))**2) # Generalized Cross Validation (GCV) elif selmethod =='gcv': @@ -261,11 +262,11 @@ def _evalalpha(alpha,V,K,L,selmethod,nonneg,noiselvl,regtype,weights,HuberParame # Extrapolated Error (EE) elif selmethod =='ee': - f_ = Residual**2/np.linalg.norm(K.T@(K@P - V)) + f_ = Residual**2/np.linalg.norm(K.T@(residuals)) # Normalized Cumulative Periodogram (NCP) elif selmethod == 'ncp': - resPeriodogram = abs(np.fft.fft(K@P - V))**2 + resPeriodogram = abs(np.fft.fft(residuals))**2 wnoisePeriodogram = np.zeros(len(resPeriodogram)) respowSpectrum = np.zeros(len(resPeriodogram)) for j in range(len(resPeriodogram)-1): @@ -279,12 +280,12 @@ def _evalalpha(alpha,V,K,L,selmethod,nonneg,noiselvl,regtype,weights,HuberParame eigs,_ = np.linalg.eig(np.eye(np.shape(H)[0],np.shape(H)[1]) - H) eigs[eigs < Treshold] = 0 nzeigs = np.real(eigs[eigs!=0]) - f_ = V.T@(V - K@P)/np.prod(nzeigs)**(1/len(nzeigs)) + f_ = V.T@(-residuals)/np.prod(nzeigs)**(1/len(nzeigs)) # Mallows' C_L (MCL) elif selmethod == 'mcl': if noiselvl==-1: - noiselvl = np.std(V - K@P) + noiselvl = np.std(residuals) f_ = Residual**2 + 2*noiselvl**2*np.trace(H) - 2*N*noiselvl**2 elif selmethod == 'lr' or selmethod == 'lc': diff --git a/deerlab/snlls.py b/deerlab/snlls.py index b1d3498b1..3628a1b2d 100644 --- a/deerlab/snlls.py +++ b/deerlab/snlls.py @@ -311,10 +311,10 @@ def linear_problem(A,optimize_alpha,alpha): # Optimiza the regularization parameter only if needed if optimize_alpha: - alpha = dl.selregparam(y, A, ax, regtype, regparam, regorder=regorder) + alpha = dl.selregparam(y, A, ax, regtype, regparam, weights=weights, regorder=regorder) # Components for linear least-squares - AtA, Aty = dl.lsqcomponents(y, A, L, alpha, weights, regtype=regtype) + AtA, Aty = dl.lsqcomponents(y, A, L, alpha, weights=weights, regtype=regtype) # Solve the linear least-squares problem result = linSolver(AtA, Aty) diff --git a/test/test_fitmodel.py b/test/test_fitmodel.py index 45fa1da32..ec9f2204a 100644 --- a/test/test_fitmodel.py +++ b/test/test_fitmodel.py @@ -785,3 +785,60 @@ def test_convergence_criteria(): assert ovl(P,fit.P) > 0.90 # ====================================================================== + +def test_global_weights(): +# ====================================================================== + "Check that the global weights properly work when specified" + + t = np.linspace(-0.3,5,300) + r = np.linspace(2,6,80) + + P1 = dl.dd_gauss(r,[3,0.2]) + P2 = dl.dd_gauss(r,[5,0.2]) + + K = dl.dipolarkernel(t,r,mod=0.2) + + scales = [1e3, 1e9] + sigma1 = 0.001 + V1 = K@P1 + dl.whitegaussnoise(t,sigma1,seed=1) + sigma2 = 0.001 + V2 = K@P2 + dl.whitegaussnoise(t,sigma2,seed=1) + + V1 = scales[0]*V1 + V2 = scales[1]*V2 + + Kmodel= lambda lam: [dl.dipolarkernel(t,r,mod=lam)]*2 + fit1 = dl.fitmodel([V1,V2],[t,t],r,'P',None,dl.ex_4pdeer,weights=[1,0]) + fit2 = dl.fitmodel([V1,V2],[t,t],r,'P',None,dl.ex_4pdeer,weights=[0,1]) + + assert ovl(P1,fit1.P) > 0.95 and ovl(P2,fit2.P) > 0.95 +# ====================================================================== + + +def test_global_weights_param(): +# ====================================================================== + "Check that the global weights properly work when specified" + + t = np.linspace(-0.3,5,300) + r = np.linspace(2,6,80) + + P1 = dl.dd_gauss(r,[3,0.2]) + P2 = dl.dd_gauss(r,[5,0.2]) + + K = dl.dipolarkernel(t,r,mod=0.2) + + scales = [1e3, 1e9] + sigma1 = 0.001 + V1 = K@P1 + dl.whitegaussnoise(t,sigma1,seed=1) + sigma2 = 0.001 + V2 = K@P2 + dl.whitegaussnoise(t,sigma2,seed=1) + + V1 = scales[0]*V1 + V2 = scales[1]*V2 + + Kmodel= lambda lam: [dl.dipolarkernel(t,r,mod=lam)]*2 + fit1 = dl.fitmodel([V1,V2],[t,t],r,dd_gauss,None,dl.ex_4pdeer,weights=[1,0]) + fit2 = dl.fitmodel([V1,V2],[t,t],r,dd_gauss,None,dl.ex_4pdeer,weights=[0,1]) + + assert ovl(P1,fit1.P) > 0.95 and ovl(P2,fit2.P) > 0.95 +# ====================================================================== \ No newline at end of file diff --git a/test/test_fitmultimodel.py b/test/test_fitmultimodel.py index d0b866070..d12394f22 100644 --- a/test/test_fitmultimodel.py +++ b/test/test_fitmultimodel.py @@ -430,4 +430,25 @@ def test_convergence_criteria(): fit = fitmultimodel(V,K,r,dd_gauss,3,'aicc', uq=False, tol=1e-3, maxiter=200) assert ovl(P,fit.P) > 0.95 # more than 99% overlap -#======================================================================= \ No newline at end of file +#======================================================================= + +def test_global_weights(): +# ====================================================================== + "Check that the global weights properly work when specified" + + t = np.linspace(0,5,300) + r = np.linspace(2,8,400) + K = dipolarkernel(t,r) + + param1 = [3,0.2] + param2 = [5,0.2] + P1 = dd_gauss(r,param1) + P2 = dd_gauss(r,param2) + V1 = K@P1 + whitegaussnoise(t,0.01,seed=1) + V2 = K@P2 + whitegaussnoise(t,0.01,seed=1) + + fit1 = fitmultimodel([V1,V2],[K,K],r,dd_gauss,3,'aic', weights=[1,0]) + fit2 = fitmultimodel([V1,V2],[K,K],r,dd_gauss,3,'aic', weights=[0,1]) + + assert ovl(P1,fit1.P) > 0.95 and ovl(P2,fit2.P) > 0.95 +# ====================================================================== \ No newline at end of file diff --git a/test/test_fitparamodel.py b/test/test_fitparamodel.py index 19bfd038b..19629f104 100644 --- a/test/test_fitparamodel.py +++ b/test/test_fitparamodel.py @@ -293,3 +293,28 @@ def test_confinter_values(): ci_match = lambda ci,ci_ref,truth:np.max(abs(np.array(ci) - np.array(ci_ref)))/truth < 0.01 assert ci_match(a_ci,a_ci_ref,p[0]) & ci_match(b_ci,b_ci_ref,p[1]) # ====================================================================== + +def test_global_weights(): +# ====================================================================== + "Check that the global weights properly work when specified" + + t = np.linspace(0,5,300) + r = np.linspace(2,8,600) + K = dipolarkernel(t,r) + + param1 = [3,0.2] + param2 = [5,0.2] + P1 = dd_gauss(r,param1) + P2 = dd_gauss(r,param2) + V1 = K@P1 + whitegaussnoise(t,0.01,seed=1) + V2 = K@P2 + whitegaussnoise(t,0.01,seed=1) + + par0 = [5, 0.5] + lb = [1, 0.1] + ub = [20, 1] + model = lambda p: [K@dd_gauss(r,p)]*2 + fit1 = fitparamodel([V1,V2],model,par0,lb,ub,weights=[1,0]) + fit2 = fitparamodel([V1,V2],model,par0,lb,ub,weights=[0,1]) + + assert all(abs(fit1.param/param1-1) < 0.03) and all(abs(fit2.param/param2-1) < 0.03) +# ====================================================================== diff --git a/test/test_fitregmodel.py b/test/test_fitregmodel.py index fb88f1ec2..401ba4319 100644 --- a/test/test_fitregmodel.py +++ b/test/test_fitregmodel.py @@ -465,3 +465,24 @@ def test_confinter_values(): ci_match = lambda ci,ci_ref,truth:np.max(abs(np.array(ci) - np.array(ci_ref)))/truth < 0.05 assert ci_match(a_ci,a_ci_ref,p[0]) & ci_match(b_ci,b_ci_ref,p[1]) # ====================================================================== + +def test_global_weights(): +# ====================================================================== + "Check that the global weights properly work when specified" + + t = np.linspace(0,5,300) + r = np.linspace(2,8,150) + K = dipolarkernel(t,r) + + param1 = [3,0.2] + param2 = [5,0.2] + P1 = dd_gauss(r,param1) + P2 = dd_gauss(r,param2) + V1 = K@P1 + whitegaussnoise(t,0.01,seed=1) + V2 = K@P2 + whitegaussnoise(t,0.01,seed=1) + + fit1 = fitregmodel([V1,V2],[K,K],r,weights=[1,0]) + fit2 = fitregmodel([V1,V2],[K,K],r,weights=[0,1]) + + assert ovl(P1,fit1.P) > 0.95 and ovl(P2,fit2.P) > 0.95 +# ====================================================================== \ No newline at end of file diff --git a/test/test_snlls.py b/test/test_snlls.py index ed7253535..ef1387ede 100644 --- a/test/test_snlls.py +++ b/test/test_snlls.py @@ -489,3 +489,32 @@ def test_confinter_scaling(): assert np.max(abs(1-ci2/ci1)) < 0.001 # Allow up to 0.1% error #============================================================ + + +def test_global_weights(): +# ====================================================================== + "Check that the global weights properly work when specified" + + t = np.linspace(-0.3,5,300) + r = np.linspace(2,6,150) + + P1 = dd_gauss(r,[3,0.2]) + P2 = dd_gauss(r,[5,0.2]) + + K = dipolarkernel(t,r,mod=0.2) + + scales = [1e3, 1e9] + sigma1 = 0.001 + V1 = K@P1 + whitegaussnoise(t,sigma1,seed=1) + sigma2 = 0.001 + V2 = K@P2 + whitegaussnoise(t,sigma2,seed=1) + + V1 = scales[0]*V1 + V2 = scales[1]*V2 + + Kmodel= lambda lam: [dipolarkernel(t,r,mod=lam)]*2 + fit1 = snlls([V1,V2],Kmodel,par0=[0.2],lb=0,ub=1,lbl=np.zeros_like(r),weights=[1,0]) + fit2 = snlls([V1,V2],Kmodel,par0=[0.2],lb=0,ub=1,lbl=np.zeros_like(r),weights=[0,1]) + + assert ovl(P1,fit1.lin) > 0.95 and ovl(P2,fit2.lin) > 0.95 +# ====================================================================== \ No newline at end of file