import sys
import os.path
import argparse
import csv
import re
from math import exp, log, sqrt
import numpy as np
import pandas as pd
import scipy.signal
import scipy.optimize
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法、Gauss-Seidel法、Smoothing penalty法、
    Ridge回帰、非負最小二乗法（NNLS）、FFT法、scipy.signal.deconvolve
    など、複数の逆畳み込みアルゴリズムをサポートしています。

主な機能:
    Jacobi法による反復逆畳み込み、Gauss-Seidel法による反復逆畳み込み、
    Smoothing penalty法による正則化逆畳み込み、Ridge回帰による正則化逆畳み込み、
    非負最小二乗法（NNLS）による制約付き逆畳み込み、FFT（高速フーリエ変換）を用いた逆畳み込み、
    scipy.signal.deconvolve を使用した直接逆畳み込み、
    numpy.convolve を用いた畳み込みによる動作確認、入力データの平滑化、拡張、整形オプション。

注意点:
    畳み込み関数を規格化するため、フィルター関数（ywf）の和で割る必要があります。
    numpy.convolve の戻り値はフィルター関数の点数の約1/2だけ前後に拡張されます。
    mode='same' を指定すると、戻り値の範囲は入力関数の範囲に一致します。
    逆畳み込みで元のデータに戻すには、拡張部分のデータも必要となる場合があります。
    scipy.signal.deconvolve は入力データの質に非常に敏感です。以下の条件を満足しないと、
    期待する結果が得られない場合があります。フィルター関数の幅は入力データよりも短いこと。
    フィルター関数の値がある程度の値よりも大きいこと（0に近い値があると問題になる可能性）。
    入力データの前には、フィルター関数の幅よりも大きいゼロ領域が必要な場合があること。
    入力データにノイズがある場合は、適切に平滑化処理が必要となること。

参考文献:
    「科学計測のための波形データ処理」、南茂夫 編著 (CQ出版社, 1986)
関連リンク:
    deconvolution_usage
"""


#===================
# 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 = [plot|gs|jacobi|sp|ridge|nnls|convolve|fft|deconvolve]
mode = 'gs'

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

# for update graph
sleeptime = 0.05

# for Ridge regression
alpha = 0.1

# for NNLS / non-negative regularized deconvolution
nnls_maxiter = None
nnls_penalty = 'ridge'

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

# x-axis sign conversion.  Useful for spectra stored as binding energy
# when plotting/analyzing on the opposite energy axis.
flip_x = 0

#===================
# 起動時引数
#===================
scriptname = sys.argv[0]


def parse_xlimit(value):
    """
    概要:
        x軸の範囲制限を解析します。
    詳細説明:
        '*'または'none'は無制限を示し、数値文字列は浮動小数点数に変換されます。
        これは、xmin/xmaxがオープンに設定できる以前の動作を維持します。
    引数:
        :param value: 解析する値。数値、文字列、またはNone。
        :type value: any
    戻り値:
        :returns: 解析されたx軸の制限値、またはNone。
        :rtype: float or None
    """
    if value is None:
        return None
    if isinstance(value, (int, float)):
        return float(value)
    text = str(value).strip()
    if text == '' or text == '*' or text.lower() in ('none', 'null'):
        return None
    return float(text)


def parse_label(value):
    """
    概要:
        データ列セレクタを解析します。
    詳細説明:
        整数インデックスまたはラベル文字列として解析します。
    引数:
        :param value: 解析する値。整数または文字列。
        :type value: int or str
    戻り値:
        :returns: 解析された列インデックス（整数）またはラベル（文字列）。
        :rtype: int or str
    """
    if isinstance(value, int):
        return value
    text = str(value).strip()
    try:
        return int(text)
    except ValueError:
        return text


def build_arg_parser():
    """
    概要:
        コマンドライン引数パーサーを構築します。
    詳細説明:
        1次元スペクトルデータをガウス関数またはローレンツ関数型の装置関数で
        逆畳み込みするための引数パーサーを作成します。
    戻り値:
        :returns: 設定済みのArgument Parserオブジェクト。
        :rtype: argparse.ArgumentParser
    """
    parser = argparse.ArgumentParser(
        description=(
            "Deconvolute one-dimensional spectral data with a Gaussian or "
            "Lorentzian instrumental function."
        ),
        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
    )

    parser.add_argument(
        'infile',
        nargs='?',
        default=infile,
        help='Input data file readable by tkVariousData.',
    )
    parser.add_argument(
        '--mode',
        choices=['plot', 'gs', 'jacobi', 'sp', 'ridge', 'nnls', 'convolve', 'fft', 'deconvolve'],
        default=mode,
        help=(
            'Calculation mode. nnls is non-negative least squares using the '
            'instrumental function as basis functions.'
        ),
    )
    parser.add_argument('--xmin', type=parse_xlimit, default=xmin,
                        help="Minimum x value. Use '*' or 'none' for no lower limit.")
    parser.add_argument('--xmax', type=parse_xlimit, default=xmax,
                        help="Maximum x value. Use '*' or 'none' for no upper limit.")
    parser.add_argument('--xlabel', type=parse_label, default=xlabel,
                        help='X column index or label.')
    parser.add_argument('--ylabel', type=parse_label, default=ylabel,
                        help='Y column index or label.')
    parser.add_argument('--flip-x', '--flip_x', dest='flip_x', type=int, choices=[0, 1], default=flip_x,
                        help='If 1, multiply input x values by -1 before range filtering, plotting, and calculation.')

    parser.add_argument('--func-type', choices=['gauss', 'lorentz', 'lorentzian'],
                        default=func_type,
                        help='Instrumental/window-function type.')
    parser.add_argument('--wa', type=float, default=Wa,
                        help='Half width at half maximum of the instrumental function.')
    parser.add_argument('--grange', type=float, default=Grange,
                        help='Window-function range in units of wa for sampled filters.')

    parser.add_argument('--alpha', type=float, default=alpha,
                        help='Regularization coefficient for ridge/sp modes.')
    parser.add_argument('--dump', type=float, default=dump,
                        help='Damping factor for Jacobi/Gauss-Seidel iterations.')
    parser.add_argument('--nmaxiter', type=int, default=nmaxiter,
                        help='Maximum iteration count for Jacobi/Gauss-Seidel modes.')
    parser.add_argument('--eps', type=float, default=eps,
                        help='Relative convergence criterion for iterative modes.')
    parser.add_argument('--nsmooth', type=int, default=nsmooth,
                        help='Smoothing width used by polynomial smoothing.')
    parser.add_argument('--zero-correction', type=int, choices=[0, 1], default=zero_correction,
                        help='Clip negative values to zero after each iterative update.')

    parser.add_argument('--convmode', choices=['same', 'full', 'valid'], default=convmode,
                        help='np.convolve mode used when smoothmode contains convolve.')
    parser.add_argument('--smoothmode', default=None,
                        help=(
                            "Combination of 'average', 'extend', and 'convolve'. "
                            "Use 'none' or an empty string for no preprocessing. "
                            "Default is mode-dependent: 'convolve+extend' for "
                            "convolve/fft, otherwise no preprocessing."
                        ))
    parser.add_argument('--kzero', type=int, default=kzero,
                        help='Number of half-window widths used for zero extension.')
    parser.add_argument('--klin', type=int, default=klin,
                        help='Number of half-window widths used for linear edge shaping.')

    parser.add_argument('--nnls-maxiter', type=int, default=nnls_maxiter,
                        help='Maximum iteration count for nnls / bounded least-squares solvers.')
    parser.add_argument('--nnls-penalty', choices=['ridge', 'smooth'], default=nnls_penalty,
                        help=(
                            'Regularization penalty used when --mode nnls and --alpha > 0. '
                            'ridge penalizes coefficient size; smooth penalizes first differences. '
                            'Use --alpha 0 for pure NNLS without regularization.'
                        ))
    parser.add_argument('--template', default=template,
                        help='Excel template used by tkVariousData.to_excel.')
    parser.add_argument('--figsize', type=float, nargs=2, default=figsize,
                        metavar=('WIDTH', 'HEIGHT'),
                        help='Matplotlib figure size.')
    return parser


def usage(app=None):
    """
    概要:
        スクリプトの使用方法を表示します。
    詳細説明:
        build_arg_parser() によって生成されたヘルプメッセージを出力します。
    引数:
        :param app: Tkinterアプリケーションオブジェクト。現在は使用されません。
        :type app: tkApplication or None
    """
    print(build_arg_parser().format_help())


def apply_args_to_globals(args):
    """
    概要:
        解析されたコマンドライン引数をグローバル変数に適用します。
    詳細説明:
        argsオブジェクトから取得した値で、スクリプト全体の動作を制御するグローバル変数を更新します。
    引数:
        :param args: コマンドライン引数を格納したNamespaceオブジェクト。
        :type args: argparse.Namespace
    """
    global infile, template, xlabel, ylabel, xmin, xmax, mode, flip_x
    global func_type, Wa, Grange, sleeptime, alpha, nnls_maxiter, nnls_penalty
    global dump, nmaxiter, eps, nsmooth, zero_correction
    global convmode, smoothmode, kzero, klin, figsize

    infile = args.infile
    template = args.template
    xlabel = args.xlabel
    ylabel = args.ylabel
    xmin = args.xmin
    xmax = args.xmax
    mode = args.mode
    flip_x = args.flip_x

    func_type = 'lorentz' if args.func_type == 'lorentzian' else args.func_type
    Wa = args.wa
    Grange = args.grange
    alpha = args.alpha
    nnls_maxiter = args.nnls_maxiter
    nnls_penalty = args.nnls_penalty

    dump = args.dump
    nmaxiter = args.nmaxiter
    eps = args.eps
    nsmooth = args.nsmooth
    zero_correction = args.zero_correction

    convmode = args.convmode
    if args.smoothmode is None:
        smoothmode = 'convolve+extend' if mode in ('convolve', 'fft') else ''
    elif args.smoothmode.lower() == 'none':
        smoothmode = ''
    else:
        smoothmode = args.smoothmode

    kzero = args.kzero
    klin = args.klin
    figsize = tuple(args.figsize)


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

def Gaussian(x, x0, whalf):
    """
    概要:
        ガウス関数を計算します。
    詳細説明:
        中心 x0、半値半幅 whalf のガウス関数を返します。
    引数:
        :param x: 評価点。
        :type x: float
        :param x0: 関数の中心。
        :type x0: float
        :param whalf: 半値半幅（HWHM）。
        :type whalf: float
    戻り値:
        :returns: ガウス関数の値。
        :rtype: float
    """
#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):
    """
    概要:
        ローレンツ関数を計算します。
    詳細説明:
        中心 x0、半値半幅 whalf のローレンツ関数を返します。
    引数:
        :param x: 評価点。
        :type x: float
        :param x0: 関数の中心。
        :type x0: float
        :param whalf: 半値半幅（HWHM）。
        :type whalf: float
    戻り値:
        :returns: ローレンツ関数の値。
        :rtype: float
    """
#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=None):
    """
    概要:
        設定された関数タイプに基づいて畳み込み関数を計算します。
    詳細説明:
        func_type グローバル変数に応じてガウス関数またはローレンツ関数を呼び出します。
        whalf が指定されない場合は、グローバル変数 Wa が使用されます。
    引数:
        :param x: 評価点。
        :type x: float
        :param x0: 関数の中心。
        :type x0: float
        :param whalf: 半値半幅（HWHM）。指定されない場合はグローバルなWaが使用されます。
        :type whalf: float or None
    戻り値:
        :returns: 畳み込み関数の値。
        :rtype: float
    """
    if whalf is None:
        whalf = Wa
    if func_type == 'gauss':
        return Gaussian(x, x0, whalf)
    else:
        return Lorentzian(x, x0, whalf)

def Hij(xstep, Wa, Grange, i, j):
    """
    概要:
        畳み込み行列の要素 Hij を計算します。
    詳細説明:
        ステップサイズ xstep、半値半幅 Wa、範囲 Grange を用いて、
        インデックス i と j に対応する畳み込み関数の値を計算します。
    引数:
        :param xstep: x軸のステップサイズ。
        :type xstep: float
        :param Wa: 半値半幅（HWHM）。
        :type Wa: float
        :param Grange: 畳み込み関数の範囲を Wa の単位で指定。
        :type Grange: float
        :param i: 行インデックス。
        :type i: int
        :param j: 列インデックス。
        :type j: int
    戻り値:
        :returns: Hij 行列の要素の値。
        :rtype: float
    """
    return convolve_func(xstep * i, xstep * j, Wa)
    
def Ridge(x, y, m, lsqfunc, w, nGmax, alpha = 0.0, is_print = False):
    """
    概要:
        Ridge回帰により係数を計算します。
    詳細説明:
        最小二乗関数 lsqfunc とエッジ補正重み w を用いて、
        観測データ y から係数 ci を推定します。正則化パラメータ alpha を使用します。
    引数:
        :param x: x軸のデータ点配列。
        :type x: numpy.ndarray or list
        :param y: 観測データ配列。
        :type y: numpy.ndarray or list
        :param m: 係数の数（データ点の数）。
        :type m: int
        :param lsqfunc: 最小二乗法で使用する関数。2つの引数 (x1, x2) を取ります。
        :type lsqfunc: callable
        :param w: エッジ補正重み配列。
        :type w: list
        :param nGmax: 畳み込み関数の最大インデックス範囲。
        :type nGmax: int
        :param alpha: 正則化パラメータ（リッジパラメータ）。デフォルトは0.0。
        :type alpha: float
        :param is_print: 中間結果を出力するかどうかのフラグ。デフォルトはFalse。
        :type is_print: bool
    戻り値:
        :returns: 推定された係数のリスト。
        :rtype: list
    """
    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.solve(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):
    """
    概要:
        Smoothing Penalty回帰により係数を計算します。
    詳細説明:
        最小二乗関数 lsqfunc とエッジ補正重み w を用いて、
        観測データ y から係数 ci を推定します。正則化パラメータ alpha と平滑化幅 nsmooth を使用します。
    引数:
        :param x: x軸のデータ点配列。
        :type x: numpy.ndarray or list
        :param y: 観測データ配列。
        :type y: numpy.ndarray or list
        :param m: 係数の数（データ点の数）。
        :type m: int
        :param lsqfunc: 最小二乗法で使用する関数。2つの引数 (x1, x2) を取ります。
        :type lsqfunc: callable
        :param w: エッジ補正重み配列。
        :type w: list
        :param nGmax: 畳み込み関数の最大インデックス範囲。
        :type nGmax: int
        :param alpha: 正則化パラメータ（ペナルティパラメータ）。デフォルトは0.0。
        :type alpha: float
        :param nsmooth: 平滑化の幅。デフォルトは1。
        :type nsmooth: int
        :param is_print: 中間結果を出力するかどうかのフラグ。デフォルトはFalse。
        :type is_print: bool
    戻り値:
        :returns: 推定された係数のリスト。
        :rtype: list
    """
    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 max([wj, 1.0e-6])
#        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.solve(Sij, Si)
    ci = ci.transpose().tolist()
    return ci[0]

def make_wf(Wa, Grange, xstep):    
    """
    概要:
        サンプリングされたウィンドウ関数（フィルター関数）を作成します。
    詳細説明:
        半値半幅 Wa と範囲 Grange を用いて、xstep 間隔でサンプリングされた
        ウィンドウ関数 yG を生成し、その積算値 SG で正規化します。
        x軸が降順であっても、ウィンドウ関数が空にならないように
        xstep の絶対値が使用されます。
    引数:
        :param Wa: 半値半幅（HWHM）。
        :type Wa: float
        :param Grange: 畳み込み関数の範囲を Wa の単位で指定。
        :type Grange: float
        :param xstep: x軸のステップサイズ。
        :type xstep: float
    戻り値:
        :returns: x軸のサンプリング点リストと、正規化されたウィンドウ関数値リストのタプル。
        :rtype: tuple of (list, list)
    例外:
        :raises ValueError: xstepが非ゼロでない場合、または正規化積算値SGがゼロ以下の場合。
    """
    # The sampled window function assumes a positive grid spacing.
    # Experimental spectra, especially XPS/HAXPES binding-energy axes,
    # are often stored in descending x order.  Use abs(xstep) here so the
    # window does not become empty when the input axis is descending.
    dx = abs(xstep)
    if dx <= 0.0:
        raise ValueError("xstep must be non-zero in make_wf().")

    ixG0   = int(Grange * Wa / dx + 1.0001)
    ixGmax = 2 * ixG0
    nxGmax = ixG0 + 1
    xG0   = ixG0 * dx

    xG = []
    yG = []
    for i in range(ixGmax+1):
        x = i * dx
        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])

    if SG <= 0.0:
        raise ValueError(
            f"Invalid sampled window function: SG={SG}. "
            f"Check wa={Wa}, grange={Grange}, xstep={xstep}."
        )

    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):
    """
    概要:
        データを畳み込みます。
    詳細説明:
        numpy.convolve を使用して生データ yraw をウィンドウ関数 ywf で畳み込みます。
        畳み込み関数の和で結果を正規化し、必要に応じてx軸を拡張します。
    引数:
        :param xraw: 生データのx軸配列。
        :type xraw: numpy.ndarray or list
        :param yraw: 生データのy軸配列。
        :type yraw: numpy.ndarray or list
        :param ywf: ウィンドウ関数（フィルター関数）のy値配列。
        :type ywf: numpy.ndarray or list
        :param kwargs: numpy.convolve に渡される追加のキーワード引数。
        :type kwargs: dict
    戻り値:
        :returns: 畳み込み後のx軸配列とy軸配列のタプル。
        :rtype: tuple of (numpy.ndarray, numpy.ndarray)
    """
    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):
    """
    概要:
        データの端をゼロ埋めと線形補間で拡張・平滑化します。
    詳細説明:
        データ y の先頭に nzero 個のゼロを追加し、最初の nlin 点を線形に整形して
        滑らかな接続を実現します。xstepは内部で計算されます。
    引数:
        :param x: 元のx軸配列。
        :type x: numpy.ndarray or list
        :param y: 元のy軸配列。
        :type y: numpy.ndarray or list
        :param nzero: データ先頭に追加するゼロの数。
        :type nzero: int
        :param nlin: 線形整形を適用するデータ点の数。
        :type nlin: int
        :param xstep: x軸のステップサイズ。デフォルトは0.0ですが、xから計算されます。
        :type xstep: float
    戻り値:
        :returns: 拡張・平滑化されたx軸配列とy軸配列のタプル。
        :rtype: tuple of (numpy.ndarray, numpy.ndarray)
    """
    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):
    """
    概要:
        移動平均でデータを平滑化します。
    詳細説明:
        各データ点について、幅 n の範囲内の平均値で置き換えることでデータを平滑化します。
        データ端では利用可能な点のみを使用します。
    引数:
        :param y: 平滑化するデータ配列。
        :type y: list
        :param n: 移動平均の窓幅。
        :type n: int
    戻り値:
        :returns: 平滑化されたデータ配列。
        :rtype: list
    """
    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):
    """
    概要:
        多項式フィッティング（Savitzky-Golay風）でデータを平滑化します。
    詳細説明:
        各データ点について、幅 n の範囲内で多項式フィッティングを行い、
        その中心点の値で置き換えることでデータを平滑化します。
        Savitzky-Golayフィルターに似た重み付けが使用されます。
    引数:
        :param y: 平滑化するデータ配列。
        :type y: list
        :param n: 平滑化の窓幅。
        :type n: int
    戻り値:
        :returns: 平滑化されたデータ配列。
        :rtype: list
    """
    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):
    """
    概要:
        FFT（高速フーリエ変換）を用いて逆畳み込みを実行します。
    詳細説明:
        生データとウィンドウ関数をFFTし、周波数領域で除算を行うことで逆畳み込みを実行します。
        データ長は2のべき乗に調整されます。
    引数:
        :param xRaw: 生データのx軸配列。
        :type xRaw: list
        :param yRaw: 生データのy軸配列。
        :type yRaw: list
        :param xG: ウィンドウ関数のx軸配列。
        :type xG: list
        :param yG: ウィンドウ関数のy軸配列。
        :type yG: list
    戻り値:
        :returns: 逆畳み込み後のx軸配列、y軸配列、FFT処理された生データのx軸配列、y軸配列、FFT処理されたウィンドウ関数のx軸配列、y軸配列のタプル。
        :rtype: tuple of (list, list, list, list, list, list)
    """
    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):
    """
    概要:
        scipy.signal.deconvolve を用いて逆畳み込みを実行します。
    詳細説明:
        scipy.signal.deconvolve を直接使用して、観測データ yRaw からウィンドウ関数 yG を逆畳み込みします。
    引数:
        :param xRaw: 生データのx軸配列。
        :type xRaw: list
        :param yRaw: 生データのy軸配列。
        :type yRaw: list
        :param xG: ウィンドウ関数のx軸配列。
        :type xG: list
        :param yG: ウィンドウ関数のy軸配列。
        :type yG: list
    戻り値:
        :returns: 逆畳み込み後のx軸配列とy軸配列のタプル。
        :rtype: tuple of (list, numpy.ndarray)
    """
    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 make_edge_correction_weights(xRaw, whalf, nGmax, w_floor=1.0e-6):
    """
    概要:
        Ridge回帰のような手法で使用されるエッジ補正重みを作成します。
    詳細説明:
        データ範囲の端で係数が不自然に大きくなるのを防ぐための重みリストを生成します。
        重みは畳み込み関数の累積和に基づいて計算されます。
    引数:
        :param xRaw: 生データのx軸配列。
        :type xRaw: numpy.ndarray
        :param whalf: 半値半幅（HWHM）。
        :type whalf: float
        :param nGmax: 畳み込み関数の最大インデックス範囲。
        :type nGmax: int
        :param w_floor: 重みの最小値。これより小さくならないようにクリップされます。デフォルトは1.0e-6。
        :type w_floor: float
    戻り値:
        :returns: エッジ補正重みのリスト。
        :rtype: list
    """
    n = len(xRaw)
    w = [0.0]
    for i in range(1, n):
        w.append(w[i - 1] + convolve_func(xRaw[i], 0.0, whalf))
    wGmax = w[n - 1]
    if abs(wGmax) < w_floor:
        return [1.0] * n
    return [max(1.0 - w[i] / wGmax, w_floor) for i in range(n)]


def edge_weight_at(w, n, nGmax, j, w_floor=1.0e-6):
    """
    概要:
        係数インデックス j に対する対称的なエッジ補正重みを返します。
    詳細説明:
        与えられた重みリスト w とインデックス j に基づいて、対称的なエッジ補正重みを計算します。
    引数:
        :param w: エッジ補正重み配列。
        :type w: list
        :param n: データ点の総数。
        :type n: int
        :param nGmax: 畳み込み関数の最大インデックス範囲。
        :type nGmax: int
        :param j: 係数のインデックス。
        :type j: int
        :param w_floor: 重みの最小値。これより小さくならないようにクリップされます。デフォルトは1.0e-6。
        :type w_floor: float
    戻り値:
        :returns: 係数インデックス j におけるエッジ補正重み。
        :rtype: float
    """
    if j <= nGmax:
        wj = w[j]
    elif n - 1 - j < nGmax:
        wj = w[n - 1 - j]
    else:
        wj = w[nGmax - 1]
    return max(wj, w_floor)


def make_basis_matrix(xRaw, whalf, w=None, nGmax=None):
    """
    概要:
        装置関数の稠密な基底行列を構築します。
    詳細説明:
        K[i, j] は係数 j を観測点 i にマッピングします。
        エッジ重みが提供された場合、各基底列は同じ補正で使用されます。
    引数:
        :param xRaw: 生データのx軸配列。
        :type xRaw: numpy.ndarray or list
        :param whalf: 半値半幅（HWHM）。
        :type whalf: float
        :param w: エッジ補正重み配列。デフォルトはNone。
        :type w: list or None
        :param nGmax: 畳み込み関数の最大インデックス範囲。デフォルトはNone。
        :type nGmax: int or None
    戻り値:
        :returns: 装置関数の基底行列。
        :rtype: numpy.ndarray
    """
    x = np.asarray(xRaw, dtype=float)
    n = len(x)
    K = np.empty((n, n), dtype=float)
    for i in range(n):
        xi = x[i]
        for j in range(n):
            v = convolve_func(xi, x[j], whalf)
            if w is not None and nGmax is not None:
                v /= edge_weight_at(w, n, nGmax, j)
            K[i, j] = v
    return K


def make_first_difference_matrix(n):
    """
    概要:
        平滑化ペナルティのための1階差分行列 D を返します。
    詳細説明:
        D @ c は [c[1]-c[0], c[2]-c[1], ...] のようなベクトルを生成します。
        このベクトルにペナルティを課すことで、係数の急激な振動を抑制します。
    引数:
        :param n: 係数の数。
        :type n: int
    戻り値:
        :returns: 1階差分行列。
        :rtype: numpy.ndarray
    """
    if n < 2:
        return np.zeros((0, n), dtype=float)
    D = np.zeros((n - 1, n), dtype=float)
    for i in range(n - 1):
        D[i, i] = -1.0
        D[i, i + 1] = 1.0
    return D


def deconvolute_nnls(xRaw, yRaw, xG, yG):
    """
    概要:
        非負最小二乗法（NNLS）を用いて逆畳み込みを実行します。
    詳細説明:
        観測データ yRaw と装置関数の基底 K を用いて、係数が非負であるという制約の下で
        最小二乗問題を解きます。任意でRidgeまたはSmoothな正則化を適用できます。
    引数:
        :param xRaw: 生データのx軸配列。
        :type xRaw: numpy.ndarray or list
        :param yRaw: 生データのy軸配列。
        :type yRaw: numpy.ndarray or list
        :param xG: ウィンドウ関数のx軸配列。
        :type xG: list
        :param yG: ウィンドウ関数のy軸配列。
        :type yG: list
    戻り値:
        :returns: 逆畳み込み後のx軸配列とy軸配列のタプル。
        :rtype: tuple of (numpy.ndarray, list)
    """
    global Wa, Grange, alpha, nnls_maxiter, nnls_penalty

    dx = abs(xRaw[1] - xRaw[0])

    print("Deconvolution by non-negative least squares (NNLS)")
    print("  Constraint: all Gaussian/Lorentzian basis coefficients >= 0")
    if alpha > 0.0:
        print(f"  Regularization: {nnls_penalty} penalty with alpha={alpha:g}")
    else:
        print("  Regularization: none; pure NNLS")
    print("")

    n = len(xRaw)
    nGmax = int(10.0 * Wa / dx + 1.0001)
    w = make_edge_correction_weights(xRaw, Wa, nGmax)
    K = make_basis_matrix(xRaw, Wa, w=w, nGmax=nGmax)

    yvec = np.asarray(yRaw, dtype=float)

    if alpha > 0.0:
        if nnls_penalty == 'smooth':
            R = make_first_difference_matrix(n)
        else:
            R = np.eye(n, dtype=float)

        A = np.vstack([K, sqrt(alpha) * R])
        b = np.concatenate([yvec, np.zeros(R.shape[0], dtype=float)])

        result = scipy.optimize.lsq_linear(
            A,
            b,
            bounds=(0.0, np.inf),
            max_iter=nnls_maxiter,
        )
        coeff = result.x
        rnorm = float(np.linalg.norm(K @ coeff - yvec))
        regnorm = float(np.linalg.norm(R @ coeff))
        print(f"  bounded LS status  = {result.status}: {result.message}")
        print(f"  data residual norm = {rnorm:g}")
        print(f"  penalty norm       = {regnorm:g}")
    else:
        coeff, rnorm = scipy.optimize.nnls(K, yvec, maxiter=nnls_maxiter)
        print(f"  NNLS residual norm = {rnorm:g}")

    y = coeff / dx

    n_active = int(np.count_nonzero(coeff > 0.0))
    print(f"  active coefficients = {n_active} / {n}")
    print(f"  coefficient range   = {coeff.min():g} - {coeff.max():g}")
    print("")

    return xRaw, y.tolist()

def deconvolute_ridge(xRaw, yRaw, xG, yG):
    """
    概要:
        Ridge回帰を用いて逆畳み込みを実行します。
    詳細説明:
        観測データ yRaw からRidge回帰により逆畳み込み係数を推定します。
        エッジ補正重みと正則化パラメータ alpha を使用します。
    引数:
        :param xRaw: 生データのx軸配列。
        :type xRaw: list
        :param yRaw: 生データのy軸配列。
        :type yRaw: list
        :param xG: ウィンドウ関数のx軸配列。
        :type xG: list
        :param yG: ウィンドウ関数のy軸配列。
        :type yG: list
    戻り値:
        :returns: 畳み込み後のx軸配列とy軸配列のタプル。
        :rtype: tuple of (list, list)
    """
    global Wa, Grange

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

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

    nGmax = int(10.0 * Wa / dx + 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_floor = 1.0e-6
    w = [max(1.0 - w[i] / wGmax, w_floor) for i in range(n)]
#    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):
    """
    概要:
        Smoothing Penalty法を用いて逆畳み込みを実行します。
    詳細説明:
        観測データ yRaw からSmoothing Penalty回帰により逆畳み込み係数を推定します。
        エッジ補正重み、正則化パラメータ alpha、および平滑化幅 nsmooth を使用します。
    引数:
        :param xRaw: 生データのx軸配列。
        :type xRaw: list
        :param yRaw: 生データのy軸配列。
        :type yRaw: list
        :param xG: ウィンドウ関数のx軸配列。
        :type xG: list
        :param yG: ウィンドウ関数のy軸配列。
        :type yG: list
        :param nsmooth: 平滑化の幅。
        :type nsmooth: int
    戻り値:
        :returns: 逆畳み込み後のx軸配列とy軸配列のタプル。
        :rtype: tuple of (list, list)
    """
    global Wa, Grange

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

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

    nGmax = int(10.0 * Wa / dx + 1.0001)
#    print("nGmax=", nGmax)

    w_floor = 1.0e-6
    w = [w_floor]
#    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 = [max(1.0 - w[i] / wGmax, w_floor) for i in range(n)]
#    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 _set_iteration_ylim(ax0, yRaw, y):
    """
    概要:
        反復プロットのy軸範囲を安定して設定します。
    詳細説明:
        生データ yRaw と更新されたデータ y に基づいて、Matplotlibの軸 ax0 のy軸範囲を
        動的に調整します。これにより、反復処理中のグラフ表示が安定します。
    引数:
        :param ax0: Matplotlibの軸オブジェクト。
        :type ax0: matplotlib.axes.Axes
        :param yRaw: 生データのy軸配列。
        :type yRaw: list or numpy.ndarray
        :param y: 更新されたデータ（反復結果）のy軸配列。
        :type y: list or numpy.ndarray
    """
    ymin = min([0.0, min(yRaw), min(y)])
    ymax = max([0.0, max(yRaw), max(y)])
    if ymax <= ymin:
        ymax = ymin + 1.0
    margin = 0.05 * (ymax - ymin)
    ax0.set_ylim([ymin - margin, ymax + margin])


def _init_iteration_plot(fig, ax0, xRaw, yRaw, y, xgmin, xgmax):
    """
    概要:
        反復更新グラフを一度だけ初期化します。
    詳細説明:
        反復処理の開始時にグラフをセットアップし、rawデータと初期の更新データをプロットします。
        Line2Dオブジェクトを再利用することで、反復ごとの描画コストを削減します。
    引数:
        :param fig: MatplotlibのFigureオブジェクト。
        :type fig: matplotlib.pyplot.Figure
        :param ax0: Matplotlibの軸オブジェクト。
        :type ax0: matplotlib.axes.Axes
        :param xRaw: 生データのx軸配列。
        :type xRaw: list or numpy.ndarray
        :param yRaw: 生データのy軸配列。
        :type yRaw: list or numpy.ndarray
        :param y: 現在の更新データ（反復結果）のy軸配列。
        :type y: list or numpy.ndarray
        :param xgmin: x軸の最小表示範囲。
        :type xgmin: float
        :param xgmax: x軸の最大表示範囲。
        :type xgmax: float
    戻り値:
        :returns: 更新されるプロットのLine2Dオブジェクト。
        :rtype: matplotlib.lines.Line2D
    """
    ax0.cla()
    ax0.plot(xRaw, yRaw, label = 'raw/initial')
    line_updated, = ax0.plot(xRaw, y, label = 'updated')
    ax0.set_xlim([xgmin, xgmax])
    _set_iteration_ylim(ax0, yRaw, y)
    ax0.legend()

    # Layout adjustment is expensive; do it once before the loop.
    app.plotevent.layout()
    fig.canvas.draw_idle()
    fig.canvas.flush_events()
    plt.pause(max(sleeptime, 0.05))
    return line_updated


def _update_iteration_plot(fig, ax0, line_updated, yRaw, y):
    """
    概要:
        反復逆畳み込み中にyデータのみを更新します。
    詳細説明:
        既存のLine2Dオブジェクトのyデータを更新し、グラフのy軸範囲を調整して再描画します。
        これにより、高速な更新が可能になります。
    引数:
        :param fig: MatplotlibのFigureオブジェクト。
        :type fig: matplotlib.pyplot.Figure
        :param ax0: Matplotlibの軸オブジェクト。
        :type ax0: matplotlib.axes.Axes
        :param line_updated: 更新されるプロットのLine2Dオブジェクト。
        :type line_updated: matplotlib.lines.Line2D
        :param yRaw: 生データのy軸配列。y軸範囲調整のために使用されます。
        :type yRaw: list or numpy.ndarray
        :param y: 更新されたデータ（反復結果）のy軸配列。
        :type y: list or numpy.ndarray
    """
    line_updated.set_ydata(y)
    _set_iteration_ylim(ax0, yRaw, y)
    fig.canvas.draw_idle()
    fig.canvas.flush_events()
    plt.pause(max(sleeptime, 0.05))


def deconvolute_jacobi(xRaw, yRaw, xG, yG, fig, ax):
    """
    概要:
        Jacobi法を用いて逆畳み込みを実行します。
    詳細説明:
        反復計算により、畳み込み方程式をJacobi法で解き、データを逆畳み込みします。
        途中で平滑化や負値クリッピングを行い、グラフをリアルタイムで更新します。
        収束条件または最大反復回数に達するまで処理を続行します。
    引数:
        :param xRaw: 生データのx軸配列。
        :type xRaw: list
        :param yRaw: 生データのy軸配列。
        :type yRaw: list
        :param xG: ウィンドウ関数のx軸配列。
        :type xG: list
        :param yG: ウィンドウ関数のy軸配列。
        :type yG: list
        :param fig: MatplotlibのFigureオブジェクト。
        :type fig: matplotlib.pyplot.Figure
        :param ax: Matplotlibの軸オブジェクトのリスト。最初の要素が使用されます。
        :type ax: list of matplotlib.axes.Axes
    戻り値:
        :returns: 逆畳み込み後のx軸配列とy軸配列のタプル。
        :rtype: tuple of (list, list)
    """
    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()

    line_updated = _init_iteration_plot(fig, ax[0], xRaw, yRaw, y, xgmin, xgmax)

    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

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

        _update_iteration_plot(fig, ax[0], line_updated, yRaw, y)

        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):
    """
    概要:
        Gauss-Seidel法を用いて逆畳み込みを実行します。
    詳細説明:
        反復計算により、畳み込み方程式をGauss-Seidel法で解き、データを逆畳み込みします。
        途中で平滑化や負値クリッピングを行い、グラフをリアルタイムで更新します。
        収束条件または最大反復回数に達するまで処理を続行します。
    引数:
        :param xRaw: 生データのx軸配列。
        :type xRaw: list
        :param yRaw: 生データのy軸配列。
        :type yRaw: list
        :param xG: ウィンドウ関数のx軸配列。
        :type xG: list
        :param yG: ウィンドウ関数のy軸配列。
        :type yG: list
        :param fig: MatplotlibのFigureオブジェクト。
        :type fig: matplotlib.pyplot.Figure
        :param ax: Matplotlibの軸オブジェクトのリスト。最初の要素が使用されます。
        :type ax: list of matplotlib.axes.Axes
    戻り値:
        :returns: 逆畳み込み後のx軸配列とy軸配列のタプル。
        :rtype: tuple of (list, list)
    """
    global Wa, Grange
#    global window

    k = sum(yG)

    print("Deconvolution by Gauss-Seidel 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()

    line_updated = _init_iteration_plot(fig, ax[0], xRaw, yRaw, y, xgmin, xgmax)

    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

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

        _update_iteration_plot(fig, ax[0], line_updated, yRaw, y)

        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




def plot_input_data(xRaw, yRaw, label_x='', label_y=''):
    """
    概要:
        選択された入力データのみをプロットします。
    詳細説明:
        逆畳み込みを実行せず、指定された入力データをグラフに表示するだけのモードです。
        逆畳み込みの結果はExcelファイルに出力されません。
    引数:
        :param xRaw: 入力データのx軸配列。
        :type xRaw: list
        :param yRaw: 入力データのy軸配列。
        :type yRaw: list
        :param label_x: x軸のラベル文字列。
        :type label_x: str
        :param label_y: y軸のラベル文字列。
        :type label_y: str
    """
    print("Plot input data only")
    print("  No deconvolution is performed.")
    print("  No deconvoluted Excel output is written.")
    print("")

    fig = plt.figure(figsize=figsize)
    ax0 = fig.add_subplot(1, 1, 1)

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

    ax0.plot(xRaw, yRaw, label='input')
    ax0.set_xlabel(str(label_x))
    ax0.set_ylabel(str(label_y))
    ax0.set_xlim([min(xRaw), max(xRaw)])

    ymin = min(yRaw)
    ymax = max(yRaw)
    if ymax <= ymin:
        ymax = ymin + 1.0
    margin = 0.05 * (ymax - ymin)
    ax0.set_ylim([ymin - margin, ymax + margin])
    ax0.legend()

    app.plotevent.layout()
    fig.canvas.draw_idle()
    fig.canvas.flush_events()
    plt.pause(0.1)

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


#======================
# main
#======================
def main():
    """
    概要:
        メイン処理を実行します。
    詳細説明:
        コマンドライン引数を解析し、入力データを読み込み、
        選択されたモードに基づいて逆畳み込み処理を実行します。
        結果はプロットされ、Excelファイルに保存されます。
    """
    global xmin, xmax

    args = build_arg_parser().parse_args()
    apply_args_to_globals(args)

    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"Process 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)
    print("flip_x  : ", flip_x)
    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 in ('ridge', 'sp'):
        print("  alpha=", alpha)
    if mode == 'nnls':
        print("For mode = 'nnls':")
        print("  non-negative coefficients: enabled")
        print("  alpha=", alpha)
        if alpha > 0.0:
            print("  regularization penalty=", nnls_penalty)
        else:
            print("  regularization penalty= none (pure NNLS)")
        print("  nnls_maxiter=", nnls_maxiter)
    if mode == 'plot':
        print("For mode = 'plot':")
        print("  plot selected input data only")
    print("")

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

    print("Input data")
    print("X range to be calculated:", xmin, "-", xmax)
    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_all = []
    yin_all = []
    for i in range(len(_xin)):
        xval = -_xin[i] if flip_x else _xin[i]
        xin_all.append(xval)
        yin_all.append(_yin[i])

    if flip_x:
        print("  flip_x=1: multiply input x values by -1 before filtering/calculation.")

    xin = []
    yin = []
    for xval, yval in zip(xin_all, yin_all):
        if not (xmin is None or xmin == '' or xmin == '*') and xval < xmin:
            continue
        if not (xmax is None or xmax == '' or xmax == '*') and xmax < xval:
            continue
        xin.append(xval)
        yin.append(yval)
#    print("x=", xin)
#    print("y=", yin)

    if len(xin) < 2:
        app.terminate(
            f"\nError in main(): Need at least two data points after x-range filtering.\n",
            usage = lambda app: usage(app), pause = True,
        )

    # Internally use ascending x order.  This avoids negative xstep problems in
    # sampled filters and keeps coefficient-to-spectrum scaling positive.
    # Descending binding-energy axes are common in photoemission data.
    if xin[1] < xin[0]:
        print("  Input x axis is descending; reverse data internally for calculation.")
        xin = list(reversed(xin))
        yin = list(reversed(yin))
    elif any(xin[i + 1] < xin[i] for i in range(len(xin) - 1)):
        print("  Input x axis is not monotonic; sort data by x for calculation.")
        xy = sorted(zip(xin, yin), key=lambda t: t[0])
        xin = [v[0] for v in xy]
        yin = [v[1] for v in xy]

    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

    if mode == 'plot':
        plot_input_data(xRaw, yRaw, label_x=label_x, label_y=label_y)
        return

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

    fig.canvas.draw_idle()
    fig.canvas.flush_events()
    plt.pause(0.1)
#    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')
    elif mode == 'nnls':
        yconv = SmoothingByPolynomialFit(yconv, nsmooth)
        xDec, yDec = deconvolute_nnls(xconv, yconv, xG, yG)
        datafft = ax[2].plot(xG, yG, label = 'WF for NNLS')
    elif mode == 'convolve' or mode == 'deconvolve':
        xDec, yDec = deconvolute_deconvolve(xconv, yconv, xG, yG)
        datafft = ax[2].plot(xG, yG, label = 'WF for scipy.signal.deconvolve')
    else:
        raise ValueError(f"Unsupported mode: {mode}")

    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) if len(yG) > 0 else 1.0])
#    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()

        fig.canvas.draw_idle()
        fig.canvas.flush_events()
        plt.pause(0.1)

    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__':
    main()