import sys
import os.path
import csv
import re
from math import exp, log, sqrt
import numpy as np
import pandas as pd
import scipy.signal
import matplotlib.pyplot as plt
import matplotlib.widgets as wg


from tklib.tkutils import replace_path, getarg, getintarg, getfloatarg, pint, pfloat
from tklib.tksci.tksci import kB, e
from tklib.tkvariousdata import tkVariousData
from tklib.tkapplication import tkApplication
from tklib.tkgraphic.tkplotevent import tkPlotEvent


"""
逆畳み込み
1. Jacobi法による逆畳み込み
2. Gauss-Seidel法による逆畳み込み
参考文献：「科学計測のための波形データ処理」、南茂夫 編著 (CQ出版社, 1986)

テスト中
1. Smoothing penalty法による逆畳み込み
2. Ridge回帰による逆畳み込み

動作確認用
3. numpy.convolveによる畳み込み
4. scipy.signal.deconvolveによる逆畳み込み
5. fftによる逆畳み込み
注意点:
yconv = np.convolve(yraw, ywf, **kwargs) / sum(ywf)
　畳み込み関数を規格化するため、フィルター関数ywfの和で割る必要がある。
　戻り値のリストは、フィルター関数の点数の１／２だけ、前後に拡張される。
　　mode = "same"を指定すると、戻り値のリストの範囲は入力関数の範囲に一致させられる。
　　逆畳み込みで元のデータに戻すには、拡張部分のデータも必要。
IDec, remainder = scipy.signal.deconvolve(Iconv, ywf) * sum(ywf)
　逆畳み込みは入力データの質に非常に敏感。以下の条件を満足しないと、ちゃんとした結果が得られない。
・フィルター関数の幅は入力データよりも短い
・フィルター関数の値はある程度の値よりも大きい（０に近い値があるとだめ？）
・入力データの前には、フィルター関数の幅よりも大きいゼロ領域が必要
・入力データにノイズがある場合は、適切に平滑化処理が必要
"""


#===================
# Global parametrs
#===================
app = tkApplication()

pi = 3.1415926535

# Input file
infile = "pes.csv"
template = "StandardGraph.xlsm"


# Analysis range
xlabel = 0
ylabel = 1
xmin = -1.0e8  # -4.5 # eV
xmax =  1.0e8  # 2.0 # eV

# mode = [gs|jacobi|ridge|convolve|fft]
mode = 'gs'

# Filter parameters for window function
func_type = 'gauss'
Wa     = 0.12  # eV
Grange = 2.0

# for update graph
sleeptime = 0.001

# for Ridge regression
alpha = 0.1

# for Jacobi/Gauss-Seidel method
dump     = 1.0
nmaxiter = 300
eps      = 1.0e-4
nsmooth  = 5

zero_correction = 0

# only for 'convolve': convmode = [same|full|]
convmode = 'same'
# smoothmode = [convolve|extend|average]
smoothmode = 'convolve+extend'

# only for 'convolve' and 'fft': 'exnted data' parameters
kzero = 5
klin  = 5

figsize = (6, 6)

#===================
# 起動時引数
#===================
argv = sys.argv
scriptname = argv[0]
def usage(app):
    print("")
    print("usage: python {} file mode xmin xmax Wa func_type dump nmaxiter eps nsmooth zeroc xlable ylabel".format(scriptname))
    print("       mode = [gs|jacobi|sp|ridge]   gs=Gauss-Seidel method   jacobi=Jacobi method   sp=Smoothing penalty method")
    print("       zeroc = [0|1]  zero correction after a Jacobi/Gauss-Seidel cycle]")
    print("  ex.: python {} pes.csv gs -6.0 2.0 0.12 gauss 1.0 300 1.0e-4 5 0 0 1".format(scriptname))
    print("usage: python {} file mode convmode smoothmode xmin xmax Wa func_type Grange kzero klin".format(scriptname))
    print("       mode = [convolve|fft]")
    print("       convmode = [''|full|same]")
    print("       smoothmode = combination of [none|convolve|extend|average]")
    print("  ex.: python {} SnSe-DOS.csv convolve same convolve+extend -4.5 2.0 0.12 gauss 2.0 1 5".format(scriptname))
    print("  ex.: python {} SnSe-DOS.csv fft full convolve+extend -4.5 2.0 0.12 gauss 2.0 5 5".format(scriptname))
    print("")


if len(argv) == 1:
    app.terminate("", usage = lambda app: usage(app), pause = True)
    exit()
    
if len(argv) >= 2:
    infile = argv[1]
if len(argv) >= 3:
    mode = argv[2]

if mode == 'convolve' or mode == 'fft':
    if len(argv) >= 4:
        convmode = argv[3]
    if len(argv) >= 5:
        smoothmode = argv[4]
    if len(argv) >= 6:
        xmin = float(argv[5])
    if len(argv) >= 7:
        xmax = float(argv[6])
    if len(argv) >= 8:
        Wa = float(argv[7])
    if len(argv) >= 9:
        func_type = argv[8]
    if len(argv) >= 10:
        Grange = float(argv[9])
    if len(argv) >= 11:
        kzero = float(argv[10])
    if len(argv) >= 12:
        klin = float(argv[11])
elif mode == 'gs' or mode == 'jacobi' or mode == 'ridge' or mode == 'sp':
    convmode   = ''
    smoothmode = ''
    if len(argv) >= 4:
        xmin = pfloat(argv[3], defval = argv[3])
    if len(argv) >= 5:
        xmax = pfloat(argv[4], defval = argv[4])
    if len(argv) >= 6:
        Wa = float(argv[5])
    if len(argv) >= 7:
        func_type = argv[6]
    if len(argv) >= 8:
        dump  = float(argv[7])
        alpha = float(argv[7])
    if len(argv) >= 9:
        nmaxiter = int(argv[8])
    if len(argv) >= 10:
        eps = float(argv[9])
    if len(argv) >= 11:
        nsmooth = int(argv[10])
    if len(argv) >= 12:
        zero_correction = int(argv[11])
    if len(argv) >= 13:
        xlabel = pint(argv[12], defval = argv[12])
    if len(argv) >= 14:
        ylabel = pint(argv[12], defval = argv[13])
else:
    app.terminate(f"\nError: Invalid mode=[{mode}]\nPress ENTER to exit>>", usage = lambda app: usage(app), pause = True)

#header, ext = os.path.splitext(infile)
#outfile = header + '-deconvoluted.xlsx'

def Gaussian(x, x0, whalf):
#A = 1/whalf * sqrt(ln2 / pi)
    A = 0.469718639 / whalf
#a = whalf / sqrt(ln2)
    a = whalf / 0.832554611
    X = (x - x0) / a
    return A * exp(-X*X)

def Lorentzian(x, x0, whalf):
#A = 1/whalf/pi
    A = 1.0 / whalf / pi
    X = (x - x0) / whalf
    return A / (1.0 + X * X)

def convolve_func(x, x0, whalf):
    if func_type == 'gauss':
        return Gaussian(x, x0, Wa)
    else:
        return Lorentzian(x, x0, Wa)

def Hij(xstep, Wa, Grange, i, j):
    return convolve_func(xstep * i, xstep * j, Wa)
    
def Ridge(x, y, m, lsqfunc, w, nGmax, alpha = 0.0, is_print = False):
    n = len(x)
    Si  = np.empty([m, 1])
    Sij = np.empty([m, m])

    def cal_w(j):
        if j <= nGmax:
            wj = w[j]
        elif n - 1 - j < nGmax:
            wj = w[n-1-j]
        else:
            wj = w[nGmax - 1]
        return wj
        
    for l in range(0, m):
        v = 0.0
        for i in range(0, n):
            v += y[i] * lsqfunc(x[l], x[i])
#        Si[l, 0] = v
        Si[l, 0] = v / cal_w(l)

    for j in range(0, m):
        for l in range(j, m):
            v = 0.0
            for i in range(0, n):
                v += lsqfunc(x[j], x[i]) * lsqfunc(x[l], x[i])
#            Sij[j, l] = Sij[l, j] = v
            Sij[j, l] = Sij[l, j] = v / cal_w(j) / cal_w(l)
        Sij[j, j] += alpha

    if is_print:
        print("Vector and Matrix:")
        print("Si=")
        print(Si)
        print("Sij=")
        print(Sij)
        print("")

    ci = np.linalg.inv(Sij) @ Si
    ci = ci.transpose().tolist()
    return ci[0]

def Smoothing_Penalty_Regression(x, y, m, lsqfunc, w, nGmax, alpha = 0.0, nsmooth = 1, is_print = False):
    n = len(x)
    Si  = np.empty([m, 1])
    Sij = np.empty([m, m])

    def cal_w(j):
        if j <= nGmax:
            wj = w[j]
        elif n - 1 - j < nGmax:
            wj = w[n-1-j]
        else:
            wj = w[nGmax - 1]
        return wj

    for l in range(0, m):
        v = 0.0
        for i in range(0, n):
            v += y[i] * lsqfunc(x[l], x[i])
#        Si[l, 0] = v
        Si[l, 0] = v / cal_w(l)

    for j in range(0, m):
        for l in range(j, m):
            v = 0.0
            for i in range(0, n):
                v += lsqfunc(x[j], x[i]) * lsqfunc(x[l], x[i])
#            Sij[j, l] = Sij[l, j] = v
            Sij[j, l] = Sij[l, j] = v / cal_w(j) / cal_w(l)

    nband = int(nsmooth / 2 + 1.0001)
    use_ridge = 1
    if use_ridge > 1:
        for j in range(0, m):
            Sij[j, j] += alpha
    elif use_ridge > 0:
        for j in range(0, m):
            j0 = j - nband
            j1 = j + nband
            if j0 < 0:
                j0 = 0
            if j1 >= m:
                j1 = m
            for k in range(j0, j):
                Sij[j, k] -= alpha
            for k in range(j+1, j1):
                Sij[j, k] -= alpha

    if is_print:
        print("Vector and Matrix:")
        print("Si=")
        print(Si)
        print("Sij=")
        print(Sij)
        print("")

    ci = np.linalg.inv(Sij) @ Si
    ci = ci.transpose().tolist()
    return ci[0]

def make_wf(Wa, Grange, xstep):    
    ixG0   = int(Grange * Wa / xstep + 1.0001)
    ixGmax = 2 * ixG0
    nxGmax = ixG0 + 1
    xG0   = ixG0 * xstep

    xG = []
    yG = []
    for i in range(ixGmax+1):
        x = i * xstep
        xG.append(x)
        yG.append(convolve_func(x, xG0, Wa))

    SG = 0.0
    for i in range(len(yG)-1):
        SG += (yG[i] + yG[i+1]) / 2.0 * (xG[i+1] - xG[i])

    for i in range(ixGmax+1):
        yG[i] /= SG

    print("   Range: {} in width".format(Grange * Wa))
    print("   i range: {} - {} at center {}".format(0,ixGmax, xG0))
    print("   ixGmax = ", ixGmax)
    print("   SG = ", SG)

    return xG, yG


def convolve(xraw, yraw, ywf, **kwargs):
    yconv = np.convolve(yraw, ywf, **kwargs) / sum(ywf)
    n_new = len(yconv)
    dn = n_new - len(yraw)
    if dn > 0:
        offset = int(dn / 2)
        xmin = xraw[0]
        xstep = xraw[1] - xmin
        xmin_new = xmin - offset * xstep
        print("convolve: the length of the output data has been changed "
                + "from {} to {}".format(len(yraw), n_new))
        print("  Add offset = {}".format(offset))
        print("  xmin changes: {} => {}".format(xmin, xmin_new))
        x = np.array([xmin_new + i * xstep for i in range(n_new)])
        return x, yconv
    return xraw, yconv

def extend_smooth(x, y, nzero, nlin, xstep = 0.0):
    xmin = x[0]
    xstep = x[1] - x[0]
    xmin_new = x[0] - nzero * xstep
    n_new = nzero + len(x)
    print("extend_smooth:")
    print("  Add {} zeros at top of the data".format(nzero))
    print("    xmin changes: {} => {}".format(xmin, xmin_new))
    print("  Reshape {} input data with a linear filter".format(nlin))

    xx = np.array([xmin_new + i * xstep for i in range(n_new)])
    yy = np.zeros(n_new)
    for i in range(nlin):
        k = i / (nlin - 1)
        yy[i+nzero] = k * y[i]
    for i in range(len(x) - nlin):
        yy[i+nzero+nlin] = y[i+nzero]
    return xx, yy

def SmoothingBySimpleAverage(y, n):
    n2 = int(n / 2);
    ndata = len(y);
    ys = []
    for i in range(0, ndata):
        c = 0;
        ys.append(0.0);
        for k in range(i - n2, i + n2 + 1):
            if k < 0 or k >= ndata:
                continue
            ys[i] += y[k]
            c += 1
        if c > 0:
            ys[i] /= c;
        else:
            ys[i] = y[i]
    return ys;

def SmoothingByPolynomialFit(y, n):
    m = int(n / 2);
    W23 = (4.0 * m * m - 1.0) * (2.0 * m + 3.0) / 3.0;
    w23j = [0.0]*(2*m+1)
    for j in range(-m, m+1):
        w23j[j + m] = (3.0 * m * (m+1.0) - 1.0 - 5.0 * j * j) / W23

    ndata = len(y)
    ys = []
    for i in range(0, ndata):
        if(i-m < 0 or ndata <= i+m):
            ys.append(y[i])
        else:
            c = 0.0;
            ys.append(0.0);
            for j in range(-m, m+1):
                k = i + j
                if k < 0 or k >= ndata:
                    continue
                ys[i] += w23j[j+m] * y[k]
                c += w23j[j+m]

            if c > 0:
                ys[i] /= c
            else:
                ys[i] = y[i]

    return ys;


def deconvolute_fft(xRaw, yRaw, xG, yG):
    k = sum(yG)

    print("Deconvolution by FFT")
    n = len(xRaw)
    nlog = int(log(n) / log(2) + 1.0 - 1.0e-5)
    nfft = pow(2, nlog)
    print("  Data number is changed from {} to 2^{} = {} for FFT".format(n, nlog, nfft))
    
    xmin = xRaw[0]
    xstep = xRaw[1] - xmin
    xRawFFT = [xmin + i * xstep for i in range(nfft)]
    yRawFFT = np.insert(yRaw, len(yRaw), np.zeros(nfft - n))
# filterの中心位置の原点からのずれによって、iFFT後の原点がずれる
    yGFFT   = np.insert(yG, len(yG), np.zeros(nfft - len(yG)))
    xminG = xmin + len(xG) / 2 * xstep
    print("  xmin: ", xmin, xminG) 
    xGFFT = [xminG + i * xstep for i in range(nfft)]
#    nadd = int((nfft - len(yG)) / 2)
#    yGFFT   = np.insert(yG, len(yG), np.zeros(nadd))
#    yGFFT   = np.insert(yGFFT, 0, np.zeros(nfft - len(yGFFT)))

    yRawFFTed = np.fft.fft(yRawFFT)
    yGFFTed   = np.fft.fft(yGFFT)
    ycFFTed = yRawFFTed / yGFFTed
    ydeconv = np.fft.ifft(ycFFTed)
    ydeconv = [float(ydeconv[i]) for i in range(len(ydeconv))]

    print("")

    return xGFFT, ydeconv, xRawFFT, yRawFFT, xRawFFT, yGFFT

def deconvolute_deconvolve(xRaw, yRaw, xG, yG):
    k = sum(yG)

    print("Deconvolution by scipy.signal.deconvove")
    print("")
    IDec, remainder = scipy.signal.deconvolve(yRaw, yG)
    IDec *= k
    ndata = len(xRaw)
    nGhalf = int(len(xG) / 2)

    return xRaw[nGhalf:ndata-nGhalf], IDec

def deconvolute_ridge(xRaw, yRaw, xG, yG):
    global Wa, Grange

    k  = sum(yG)
    dx = xRaw[1] - xRaw[0]

    print("Deconvolution by Ridge regression")
    print("")

    nGmax = int(10.0 * Wa / (xRaw[1] - xRaw[0]) + 1.0001)
#    print("nGmax=", nGmax)

    w = [0.0]
    n = len(xRaw)
    for i in range(1, len(xRaw)):
        w.append(w[i-1] + convolve_func(xRaw[i], 0, Wa))
    wGmax = w[n-1]
    w = [1.0 - w[i] / wGmax for i in range(n)]

    y = Ridge(xRaw, yRaw, len(xRaw), lambda x1, x2: convolve_func(x1, x2, Wa), w, nGmax, alpha)
    y = [v / dx for v in y]

    return xRaw, y

def deconvolute_sp(xRaw, yRaw, xG, yG, nsmooth):
    global Wa, Grange

    k  = sum(yG)
    dx = xRaw[1] - xRaw[0]

    print("Deconvolution by smoothing penalty method")
    print("")

    nGmax = int(10.0 * Wa / (xRaw[1] - xRaw[0]) + 1.0001)
#    print("nGmax=", nGmax)

    w = [0.0]
    n = len(xRaw)
    for i in range(1, len(xRaw)):
        w.append(w[i-1] + convolve_func(xRaw[i], 0, Wa))
    wGmax = w[n-1]
    w = [1.0 - w[i] / wGmax for i in range(n)]

    y = Smoothing_Penalty_Regression(xRaw, yRaw, len(xRaw), 
                lambda x1, x2: convolve_func(x1, x2, Wa), w, nGmax, alpha, nsmooth)
    y = [v / dx for v in y]

    return xRaw, y

#window = None
def deconvolute_jacobi(xRaw, yRaw, xG, yG, fig, ax):
    global Wa, Grange
#    global window

    k = sum(yG)

    print("Deconvolution by Jacobi method")
    print("")

    xstep = xRaw[1] - xRaw[0]

    xgmin = min(xRaw)
    xgmax = max(xRaw)

    n = len(xRaw)
    Sg = np.zeros(n)
    for i in range(n):
        for j in range(n):
            Sg[i] += Hij(xstep, Wa, Grange, i, j)
    print("Filter area w.r.t. i: Sg=", Sg[int(n/2)])

    ymax = max([abs(yRaw[i]) for i in range(n)])
    y     = yRaw.copy()
    yPrev = yRaw.copy()

    for it in range(nmaxiter):
        Hx = np.zeros(n)
        print("iter=", it)

        for i in range(n):
            for j in range(n):
                Hx[i] += Hij(xstep, Wa, Grange, i, j) * yPrev[j]
            h = Hij(xstep, Wa, Grange, i, i)
            y[i] = yPrev[i] + (yRaw[i] - Hx[i] / Sg[i]) / h * dump

#        y = SmoothingBySimpleAverage(y, nsmooth)
        y = SmoothingByPolynomialFit(y, nsmooth)
        if zero_correction:
            for i in range(n):
                if y[i] < 0.0:
                    y[i] = 0.0

        ax[0].cla()

        data1 = ax[0].plot(xRaw, yRaw, label = 'raw/initial')
        data1 = ax[0].plot(xRaw, y, label = 'updated')
#        data4 = ax[2].plot(xG, yG, label = 'filter')
        ax[0].set_xlim([xgmin, xgmax])
        ygmax = max([max(xRaw), max(y)])
        ax[0].set_ylim([0.0, ygmax])
#        ax[1].set_xlim([xgmin, xgmax])
#        ax[1].set_ylim([0.0, max(yRaw)])
#        ax[2].set_xlim([xgmin, xgmax])
#        ax[2].set_ylim([0.0, max(yG)])

#        w = plt.get_current_fig_manager().window
#        if w != window:
        if app.plotevent.stop_flag:
            print("Terminated by user")
            break

        ax[0].legend()
#        ax[2].legend()
#        plt.tight_layout()
        app.plotevent.layout()
        plt.pause(sleeptime)

        max_err = max([abs(y[i] - yPrev[i]) for i in range(n)])
        rel_err = max_err / ymax
        print("  max error: ", max_err, "  relative error: ", rel_err, "  eps=", eps)
        if max_err / ymax < eps:
            print("Converged at max_err={} ({} relative) < {}".format(max_err, rel_err, eps))
            break
        
        yPrev = y.copy()
    else:
        print("Not converged")

    return xRaw, y

def deconvolute_gauss_seidel(xRaw, yRaw, xG, yG, fig, ax):
    global Wa, Grange
#    global window

    k = sum(yG)

    print("Deconvolution by Jacobi method")
    print("")

    xstep = xRaw[1] - xRaw[0]

    xgmin = min(xRaw)
    xgmax = max(xRaw)

    n = len(xRaw)
    Sg = np.zeros(n)
    for i in range(n):
        for j in range(n):
            Sg[i] += Hij(xstep, Wa, Grange, i, j)
    print("Filter area w.r.t. i: Sg=", Sg[int(n/2)])

    ymax = max([abs(yRaw[i]) for i in range(n)])
    y     = yRaw.copy()
    yPrev = yRaw.copy()

    for it in range(nmaxiter):
        Hx    = np.zeros(n)
        print(f"iter={it:4}  ", end = '')

        for i in range(n):
            for j in range(i):
                Hx[i] += Hij(xstep, Wa, Grange, i, j) * y[j]
            for j in range(i, n):
                Hx[i] += Hij(xstep, Wa, Grange, i, j) * yPrev[j]
            h = Hij(xstep, Wa, Grange, i, i)
            y[i] = yPrev[i] + (yRaw[i] - Hx[i] / Sg[i]) / h * dump

#        y = SmoothingBySimpleAverage(y, nsmooth)
        y = SmoothingByPolynomialFit(y, nsmooth)
        if zero_correction:
            for i in range(n):
                if y[i] < 0.0:
                    y[i] = 0.0

        ax[0].cla()

        data1 = ax[0].plot(xRaw, yRaw, label = 'raw/initial')
        data1 = ax[0].plot(xRaw, y, label = 'updated')
#        data4 = ax[2].plot(xG, yG, label = 'filter')
        ax[0].set_xlim([xgmin, xgmax])
        ygmax = max([max(xRaw), max(y)])
        ax[0].set_ylim([0.0, ygmax])
#        ax[1].set_xlim([xgmin, xgmax])
#        ax[1].set_ylim([0.0, max(yRaw)])
#        ax[2].set_xlim([xgmin, xgmax])
#        ax[2].set_ylim([0.0, max(yG)])

#        w = plt.get_current_fig_manager().window
#        if w != window:
        if app.plotevent.stop_flag:
            print("Terminated by user")
            break

        ax[0].legend()
#        ax[2].legend()
#        plt.tight_layout()
        app.plotevent.layout()
        plt.pause(sleeptime)

        max_err = max([abs(y[i] - yPrev[i]) for i in range(n)])
        rel_err = max_err / ymax
        print(f"  max error: {max_err:10.4g}  relative error: {rel_err:10.4g}  eps={eps:10.4g}")
        if max_err / ymax < eps:
            print("Converged at max_err={} ({} relative) < {}".format(max_err, rel_err, eps))
            break
        
        yPrev = y.copy()
    else:
        print("Not converged")

    return xRaw, y


#======================
# main
#======================
def main():
    global xmin, xmax

    logfile = app.replace_path(infile)
    outfile = app.replace_path(infile, template = ["{dirname}", "{filebody}-deconvoluted.xlsx"])
    print(f"Open logfile [{logfile}]")
    app.redirect(targets = ["stdout", logfile], mode = 'w')

    print("")
    print(f"Deconvolute data in [{infile}] with {func_type} function of width={Wa}")
    print("infile  : ", infile)
    print("xlabel  : ", xlabel)
    print("ylabel  : ", ylabel)
    print("outfile : ", outfile)
    print("mode    : ", mode)
    if mode == 'convolve' or mode == 'fft':
        print("For mode = 'convolve' or 'fft':")
        print("  convmode: ", convmode)
        print("  x range : ", xmin, xmax)
    if mode == 'gs' or mode == 'jacobi':
        print("  x range : ", xmin, xmax)
        print("For mode = 'gs' or 'jacobi':")
        print("  dump=", dump)
        print("  nmaxiter=", nmaxiter)
        print("  eps=", eps)
        print("  nsmooth=", nsmooth)
        print("  zero correction=", zero_correction)
    if mode == 'ridge' or mode == 'sp':
        print("  alpha=", alpha)
    print("")

    if '*** choose ***' in xlabel:
        app.terminate(f"\nError: Choose x label\nTerminated.\n", usage = lambda app: usage(app), pause = True)
    if '*** choose ***' in ylabel:
        app.terminate(f"\nError: Choose y label\nTerminated.\n", usage = lambda app: usage(app), pause = True)

    print("Input data")
    print(f"X range to be calculated: {xmin:12} - {xmax:12}")
    datafile = tkVariousData(infile)
    header, datalist = datafile.Read_minimum_matrix(close_fp = True, usage = lambda app: usage(app))
    label_x, _xin = datafile.FindDataArray(xlabel, flag = 'i')
    label_y, _yin = datafile.FindDataArray(ylabel, flag = 'i')
    if _xin is None:
        app.terminate(f"Error: Can not find the data with [{xlabel}]", usage = lambda app: usage(app), pause = True)
    if _yin is None:
        app.terminate(f"Error: Can not find the data with [{xlabel}]", usage = lambda app: usage(app), pause = True)

    xin = []
    yin = []
    for i in range(len(_xin)):
        if not (xmin is None or xmin == '' or xmin == '*') and _xin[i] < xmin:
            continue
        if not (xmax is None or xmax == '' or xmax == '*') and xmax < _xin[i]:
            continue
        xin.append(_xin[i])
        yin.append(_yin[i])
#    print("x=", xin)
#    print("y=", yin)

    nindata = len(xin)
    print("  ndata = ", nindata)
    if xmin is None or xmin == '' or xmin == '*':
        xmin = min(xin)
    if xmax is None or xmax == '' or xmax == '*':
        xmax = max(xin)
    print("x range=", xmin, xmax)

    xinmin  = xin[0]
    xinstep = xin[1] - xin[0]
    if xinstep == 0.0:
        app.terminate(f"\nError in main(): xin[0] and xin[1] are the same. Check input file [{infile}].\n",
                    usage = lambda app: usage(app), pause = True)

    xinmax  = xinmin + xinstep * (nindata - 1)
    print("  x range: {} - {} at {} step".format(xinmin, xinmax, xinstep))
    print("Smooth mode: ", smoothmode)
    print("")

    xRaw  = xin
    yRaw  = yin
    xstep = xinstep
    print("Filter: Wa = {}  Grange = {}".format(Wa, Grange))
    xG, yG = make_wf(Wa, Grange, xstep)
    print("")

    fig = plt.figure(figsize = figsize)
    ax = []
    ax.append(fig.add_subplot(3, 1, 1))
    ax.append(fig.add_subplot(3, 1, 2))
    ax.append(fig.add_subplot(3, 1, 3))
#    ax.append(fig.add_subplot(4, 1, 4))

    app.plotevent = tkPlotEvent(plt)
    app.plotevent.add_button()
    app.plotevent.layout()

    plt.pause(0.001)
#    window = plt.get_current_fig_manager().window

# make data to be deconvolved
    nw = int(len(xG) / 2)
    if 'average' in smoothmode:
        yRaw = SmoothingBySimpleAverage(yRaw, 31)
    if 'extend' in smoothmode:
        xRaw, yRaw = extend_smooth(xRaw, yRaw, nw*kzero, nw*klin, xstep)
    if 'convolve' in smoothmode:
        xconv, yconv = convolve(xRaw, yRaw, yG, mode = convmode)
    else:
        xconv, yconv = xRaw, yRaw

# deconvolution
    if mode == 'fft':
        xDec, yDec, xRawFFT, yRawFFT, xGFFT, yGFFT = deconvolute_fft(xconv, yconv, xG, yG)
        xconv = xRawFFT
        yconv = yRawFFT
        datafft = ax[2].plot(xGFFT, yGFFT, label = 'WF for FFT')
    elif mode == 'jacobi':
        xDec, yDec = deconvolute_jacobi(xconv, yconv, xG, yG, fig, ax)
        yG = [convolve_func(xRaw[i], 0.0, Wa) for i in range(len(xRaw))]
        datafft = ax[2].plot(xRaw, yG, label = 'WF for Jacobi method')
    elif mode == 'gs':
        xDec, yDec = deconvolute_gauss_seidel(xconv, yconv, xG, yG, fig, ax)
        yG = [convolve_func(xRaw[i], 0.0, Wa) for i in range(len(xRaw))]
        datafft = ax[2].plot(xRaw, yG, label = 'WF for Gauss-Seidel method')
    elif mode == 'sp':
        yconv = SmoothingByPolynomialFit(yconv, nsmooth)
        xDec, yDec = deconvolute_sp(xconv, yconv, xG, yG, nsmooth)
#        yconv = SmoothingByPolynomialFit(yconv, nsmooth)
        datafft = ax[2].plot(xG, yG, label = 'WF for Gauss kernel regression')
    elif mode == 'ridge':
        yconv = SmoothingByPolynomialFit(yconv, nsmooth)
        xDec, yDec = deconvolute_ridge(xconv, yconv, xG, yG)
        datafft = ax[2].plot(xG, yG, label = 'WF for Ridge regression')
    else:
        xDec, yDec = deconvolute_deconvolve(xconv, yconv, xG, yG)
        datafft = ax[2].plot(xG, yG, label = 'WF for convolve')

    app.plotevent.finalize_button()
    if app.plotevent.stop_flag:
        print("")
        print("Iteration terminated by user")
    else:
        xgmin = min([min(xconv), min(xDec)])
        xgmax = max([max(xconv), max(xDec)])

        ax[0].cla()
        data1 = ax[0].plot(xconv, yconv, label = 'input(convoluted)')
        data3 = ax[0].plot(xDec, yDec, label = 'deconvoluted')
        data1 = ax[1].plot(xRaw, yRaw, label = 'input(raw)')
        data2 = ax[1].plot(xconv, yconv, label = 'input(convoluted)')
#    data5 = ax[2].plot(xDec, yDec, label = 'deconvoluted')
#    data6 = ax[2].plot([xgmin, xgmax], [0.0, 0.0], color = 'red', linestyle = 'dashed')
        ax[0].set_xlim([xgmin, xgmax])
        ax[0].set_ylim([0.0, max([max(yRaw), max(yDec)])])
        ax[1].set_xlim([xgmin, xgmax])
        ax[1].set_ylim([0.0, max(yRaw)])
#    ax[2].set_xlim([xgmin, xgmax])
#    ax[2].set_ylim([-100, 10000])
        ax[2].set_xlim([xgmin, xgmax])
        ax[2].set_ylim([0.0, max(yG)])
#    ax[1].set_ylim([0.0, max(yRaw)])

        ax[0].legend()
        ax[1].legend()
        ax[2].legend()
#    ax[3].legend()

#        plt.tight_layout()
        app.plotevent.layout()

        plt.pause(0.001)

    print("")
    print("Save deconvoluted data to [{}]".format(outfile))
    df = pd.DataFrame(np.array([xconv, yconv, yDec]).T, columns = ['x', 'y(input)', 'y(deconvoluted)'])
#    df.to_excel(outfile, index = False)
    data_list = df.values.tolist()
    tkVariousData().to_excel(outfile, list(df.columns), list(map(list, zip(*data_list))), template = template)

    app.terminate("Press ENTER to exit>>", usage = lambda app: usage(app), pause = True)


if __name__ == '__main__':
    usage(app)
    main()