Develop code review

modules/L1L2.py

@@ -0,0 +1,67 @@
+""" Module to implement L1 and/or L2 regularization with
+SciPy linear_models module
+"""
+
+def L1L2(t, F, bound, Nz, alpha1, alpha2, iterations = 10000):
+    """
+    Returns solution using mixed L1 and/or L2 regularization
+    with simple gradient descent
+
+    F(t) = ∫f(s)*exp(-s*t)ds
+
+    or
+
+    min = ||C*f - F||2 + alpha1*||I*f||1 + alpha2*||I*f||2
+
+    Parameters:
+    ------------
+    t : array 
+        Time domain data from experiment
+    F : array,
+        Transient data from experiment F(t)
+    bound : list
+        [lowerbound, upperbound] of s domain points
+    Nz : int
+        Number of points z to compute, must be smaller than len(F)
+    alpha1, alpha 2 : floats
+        Regularization parameters for L1 and L2 regularizers
+    iterations : int 
+        Maximum number of iterations. Optional
+
+
+    Returns:
+    ------------
+    s : array
+        Emission rates domain points (evenly spaced on log scale)
+    f : array
+        Inverse Laplace transform f(s)
+    F_restored : array 
+        Reconstructed transient from C@f + intercept
+    """
+
+    import numpy as np
+    from scipy.sparse import diags
+    from sklearn.linear_model import ElasticNet
+
+    # set up grid points (# = Nz)
+    h = np.log(bound[1]/bound[0])/(Nz - 1)      # equally spaced on logscale
+    s = bound[0]*np.exp(np.arange(Nz)*h)        # z (Nz by 1)
+
+    # construct C matrix from [1]
+    s_mesh, t_mesh = np.meshgrid(s, t)
+    C = np.exp(-t_mesh*s_mesh)       
+    C[:, 0] /= 2.
+    C[:, -1] /= 2.
+    C *= h
+
+    alpha = alpha1 + alpha2
+    l1_ratio = alpha1/alpha
+
+    model = ElasticNet(alpha=alpha, l1_ratio=l1_ratio, tol = 1e-12,
+                       fit_intercept = True, max_iter = iterations)
+    model.fit(C, F)
+
+    f = model.coef_
+
+    F_restored = C@f + model.intercept_
+    return s, f, F_restored

modules/L1_FISTA.py

@@ -0,0 +1,21 @@
+from fista import Fista
+import numpy as np
+
+def l1_fista(t, F, bound, Nz, alpha, iterations = 50000, penalty = 'l11'):
+
+    "Function to implement FISTA algorithm"
+
+    h = np.log(bound[1]/bound[0])/(Nz - 1)      # equally spaced on logscale
+    z = bound[0]*np.exp(np.arange(Nz)*h)        # z (Nz by 1)
+
+    z_mesh, t_mesh = np.meshgrid(z, t)
+    C = np.exp(-t_mesh*z_mesh)                   # specify integral kernel
+    C[:, 0] /= 2.
+    C[:, -1] /= 2.
+    C *= h
+
+    fista = Fista(loss='least-square', penalty = penalty, lambda_=alpha, n_iter = iterations)
+    fista.fit(C, F)
+
+
+    return z, fista.coefs_, fista.predict(C)

modules/L2.py

@@ -0,0 +1,67 @@
+"""Module to implement routines of numerical inverse 
+Laplace tranform using L2 regularization algorithm
+"""
+
+import numpy as np
+from scipy.sparse import diags
+
+def L2(t, F, bound, Nz, alpha):
+    """
+    Returns solution for problem imposing L2 regularization
+
+    F(t) = ∫f(s)*exp(-s*t)ds
+
+    or
+
+    min = ||C*f - F||2 + alpha*||I*f||2
+
+
+    Parameters:
+    ------------
+    t : array 
+        Time domain data from experiment
+    F : array,
+        Transient data from experiment F(t)
+    bound : list
+        [lowerbound, upperbound] of s domain points
+    Nz : int
+        Number of points z to compute, must be smaller than len(F)
+    alpha : float
+        Regularization parameters for L2 regularizers
+    iterations : int 
+        Maximum number of iterations. Optional
+
+
+    Returns:
+    ------------
+    s : array
+        Emission rates domain points (evenly spaced on log scale)
+    f : array
+        Inverse Laplace transform f(s)
+    F_restored : array 
+        Reconstructed transient from C@f + intercept
+    """
+
+    # set up grid points (# = Nz)
+    h = np.log(bound[1]/bound[0])/(Nz - 1)      # equally spaced on logscale
+    s = bound[0]*np.exp(np.arange(Nz)*h)        # z (Nz by 1)
+
+    # construct C matrix from [1]
+    s_mesh, t_mesh = np.meshgrid(s, t)
+    C = np.exp(-t_mesh*s_mesh)       
+    C[:, 0] /= 2.
+    C[:, -1] /= 2.
+    C *= h
+
+    l2 = alpha
+
+    data = [-2*np.ones(Nz), 1*np.ones(Nz), 1*np.ones(Nz)]
+    positions = [-1, -2, 0]
+    I = diags(data, positions, (Nz+2, Nz)).toarray()
+    #I      = np.identity(Nz)
+
+    f   = np.linalg.solve(l2*np.dot(I.T,I) + np.dot(C.T,C), np.dot(C.T,F))
+
+    F_restored = C@f
+
+    return s, f, F_restored#, res_norm, sol_norm

modules/contin.py

@@ -0,0 +1,131 @@
+"""Module to implement routines of numerical inverse 
+Laplace tranform using Contin  algorithm [1]
+
+[1] Provencher, S. (1982)
+"""
+
+from __future__ import division
+import numpy as np
+
+def Contin(t, F, bound, Nz, alpha):
+    """
+    Function returns processed data f(z) from the equation
+
+    F(t) = ∫f(z)*exp(-z*t)
+
+    Parameters:
+    --------------
+    t : array 
+        Time domain data from experiment
+    F : array,
+        Transient data from experiment F(t)
+    bound : list
+        [lowerbound, upperbound] of z domain points
+    Nz : int
+        Number of points z to compute, must be smaller than len(F)
+    alpha : float
+        Regularization parameter
+
+    Returns:
+    --------------
+    z : array 
+        Emission rates domain points (evenly spaced on log scale)
+    f : array
+        Inverse Laplace transform f(z)
+    F_restored : array 
+        Reconstructed transient from C@f(z)
+    """
+
+
+    # set up grid points (# = Nz)
+    h = np.log(bound[1]/bound[0])/(Nz - 1)      # equally spaced on logscale
+    z = bound[0]*np.exp(np.arange(Nz)*h)        # z (Nz by 1)
+
+    # construct C matrix from [1]
+    z_mesh, t_mesh = np.meshgrid(z, t)
+    C = np.exp(-t_mesh*z_mesh)       
+    C[:, 0] /= 2.
+    C[:, -1] /= 2.
+    C *= h
+
+    # construct regularization matrix R to impose gaussian-like peaks in f(z)
+    # R - tridiagonal matrix (1,-2,1)
+    Nreg = Nz + 2
+    R = np.zeros([Nreg, Nz])
+    R[0, 0] = 1.
+    R[-1, -1] = 1.
+    R[1:-1, :] = -2*np.diag(np.ones(Nz)) + np.diag(np.ones(Nz-1), 1) \
+        + np.diag(np.ones(Nz-1), -1)
+
+    #R = U*H*Z^T 
+    U, H, Z = np.linalg.svd(R, full_matrices=False)     # H diagonal
+    Z = Z.T
+    H = np.diag(H)
+
+    #C*Z*inv(H) = Q*S*W^T
+    Hinv = np.diag(1.0/np.diag(H))
+    Q, S, W = np.linalg.svd(C.dot(Z).dot(Hinv), full_matrices=False)  # S diag
+    W = W.T
+    S = np.diag(S)
+
+    # construct GammaTilde & Stilde
+    # ||GammaTilde - Stilde*f5||^2 = ||Xi||^2
+    Gamma = np.dot(Q.T, F)
+    Sdiag = np.diag(S)
+    Salpha = np.sqrt(Sdiag**2 + alpha**2)
+    GammaTilde = Gamma*Sdiag/Salpha
+    Stilde = np.diag(Salpha)
+
+    # construct LDP matrices G = Z*inv(H)*W*inv(Stilde), B = -G*GammaTilde
+    # LDP: ||Xi||^2 = min, with constraint G*Xi >= B
+    Stilde_inv = np.diag(1.0/np.diag(Stilde))
+    G = Z @ Hinv @ W @ Stilde_inv
+    B = -G @ GammaTilde
+
+    # call LDP solver
+    Xi = ldp(G, B)
+
+    # final solution
+    zf = np.dot(G, Xi + GammaTilde)
+    f = zf/z
+
+    F_restored = C@zf
+
+    return z, f, F_restored
+
+
+def ldp(G, h):
+    """
+    Helper for Contin() for solving NNLS [1]
+
+    [1] - Lawson and Hanson’s (1974)
+    Parameters:
+    -------------
+    G : matrix
+        Z*inv(H)*W*inv(Stilde)
+    h : array
+        -G*GammaTilde
+
+    Returns:
+    -------------
+    x : array
+        Solution of argmin_x || Ax - b ||_2 
+    """
+
+    from scipy.optimize import nnls
+
+    m, n = G.shape
+    A = np.concatenate((G.T, h.reshape(1, m)))
+    b = np.zeros(n+1)
+    b[n] = 1.
+
+    # Solving for argmin_x || Ax - b ||_2 
+    x, resnorm = nnls(A, b)
+
+    r = A@x - b
+
+    if np.linalg.norm(r) == 0:
+        print('\n No solution found, try different input!')
+    else:
+        x = -r[0:-1]/r[-1]
+    return x

modules/demo.py

@@ -0,0 +1,63 @@
+def demo(Index, Nz, Reg_L1, Reg_L2, Reg_C, Reg_S, Bounds, dt, Methods, Plot, Residuals, LCurve, Arrhenius, Heatplot):
+    """Gets data from widgets initialized with interface(), 
+    calls laplace() to process data and calls plot_data(), hp()
+    to plot results
+
+    Parameters:
+    -------------
+    Index : int
+        Index of transient in dataset
+    Nz : int
+        Value which is length of calculated vector
+    Reg_L1, Reg_L2 : floats
+        Reg. parameters for FISTA(L1) and L2 regularization
+    Reg_C, Reg_S : floats
+        Reg. parameters for CONTIN and reSpect algorithms
+    Bounds : list
+        [lowerbound, upperbound ] bounds of emission rates domain points
+    dt : int
+        Time step of transient data points in ms
+    Methods : list 
+        Methods to process data
+    Plot : boolean
+        Calls plot_data() if True
+    Residuals : boolean
+        Calls regopt() and plots L-curve to control 
+        regularization if True
+    LCurve : boolean, 
+        Automatically picks optimal reg. parameter from 
+        L-curve if True
+    Heatplot : boolean
+        Plots heatplot for all dataset and saves data 
+        in .LDLTS if True
+    """
+
+    import numpy as np
+    import matplotlib.pyplot as plt
+
+    from read_file import read_file
+    from laplace import laplace
+    from plot_data import plot_data
+    from hp import hp
+    from read_file import read_file
+    from regopt import regopt
+    from interface import interface
+
+    Bounds = 10.0**np.asarray(Bounds)
+
+    t, C, T = read_file(interface.path, dt, proc = True)# time, transients, temperatures 
+
+    data = laplace(t, C[Index], Nz, Reg_L1, Reg_L2, Reg_C, Reg_S, Bounds, Methods)
+    #print(data)
+
+    if Plot:
+        ay, aph = plot_data(t, C[Index], data, T, Index)
+
+        if Heatplot:
+            hp(t, C, T, aph, Methods, Index, Reg_L1, Reg_L2, Reg_C, Reg_S, Bounds, Nz, LCurve, Arrhenius)
+
+    if Residuals:
+        regopt(t, C[Index], ay, Methods, Reg_L1, Reg_L2, Reg_C, Reg_S, Bounds, Nz)
+
+    plt.show()
+