Fix behavior of weights in global fitting (#171)

* fitregmodel: fix behaviour of weights in global fitting

* fitparamodel: add test for global weights behaviour

* lsqcomponents: fix usage of global weights

* fitmultimodel: fix behaviour of weights in global fitting

* fitmultimodel: return list of fitted signals instead of single array, fix scaling of output fit signals and UQ, fix plotting

* fitmultimodel: fix scaling issue in global fitting

* snlls: fix behaviour of weights in global fitting

* fitmodel: add tests to check behaviour of weights in global fitting

* fitmultimodel: fix error when plotting single dataset fits

* fitmultimodel: fix bug in plottin of single dataset fits

* fitmultimodel: fix error in previous commit

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}]')

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')
 

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

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':

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)

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
+# ======================================================================

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
-#=======================================================================
+#=======================================================================
+
+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
+# ======================================================================