deconvolution.py のソースコード解析に基づき、Sphinx (MyST) でビルド可能なMarkdownドキュメントを作成します。
deconvolution.py ドキュメント
概要
deconvolution.py は、1次元スペクトルデータに対し、ガウス型またはローレンツ型の装置関数（フィルタ関数）を用いて逆畳み込み演算を実行するPythonプログラムです。主に、分光データなどの測定結果から装置による広がり効果を除去し、真のスペクトル形状を回復することを目的としています。
複数の逆畳み込み手法をサポートしており、コマンドライン引数を通じて選択できます。また、結果のグラフ表示とExcelファイルへの出力機能も備えています。
機能詳細
このプログラムは、以下の逆畳み込みおよび関連する処理モードをサポートしています。
逆畳み込み手法
プログラムは以下の逆畳み込みアルゴリズムを実装または利用しています。
Jacobi法 (
jacobi):反復法による逆畳み込み手法の一つ。
参考文献：「科学計測のための波形データ処理」、南茂夫 編著 (CQ出版社, 1986) に基づいています。
各反復ステップで、前の反復の結果のみを用いて次の結果を計算します。
Gauss-Seidel法 (
gs):反復法による逆畳み込み手法の一つ。
参考文献：「科学計測のための波形データ処理」、南茂夫 編著 (CQ出版社, 1986) に基づいています。
Jacobi法と同様ですが、新しい反復値を計算する際に、すでに更新された値があればそれを利用します。これにより、Jacobi法よりも速く収束する傾向があります。
Smoothing Penalty法 (
sp):平滑化ペナルティ項を導入した逆畳み込み手法。
リッジ回帰と同様に正則化を利用しますが、係数の1次差分をペナルティとして課すことで、より滑らかな解を導きます。
Ridge回帰 (
ridge):リッジ正則化を伴う回帰ベースの逆畳み込み手法。
係数の大きさにペナルティを課すことで、過学習を防ぎ、安定した解を求めます。
非負最小二乗法 (NNLS) (
nnls):scipy.optimize.nnlsを利用した逆畳み込み手法。結果の係数（逆畳み込みされたスペクトル）が非負であることを制約として課します。
必要に応じてリッジまたはスムージングペナルティによる正則化を適用できます。
動作確認・補助機能
以下のモードは、逆畳み込みの結果の検証や、単純な畳み込み・既知のライブラリ関数による逆畳み込みを目的としています。
numpy.convolveによる畳み込み (convolve):入力データと装置関数を畳み込むことで、装置の影響を受けたスペクトルをシミュレートします。
np.convolve()のmode引数 (same,full,valid) をサポートします。
scipy.signal.deconvolveによる逆畳み込み (deconvolve):scipyライブラリの標準的な逆畳み込み関数を使用します。入力データの品質に敏感なため、注意が必要です。
FFTによる逆畳み込み (
fft):高速フーリエ変換 (FFT) を用いて周波数空間で逆畳み込みを実行します。
FFTの性質上、データ長が2のべき乗に調整されます。
入力データのみプロット (
plot):選択された入力データのみをプロットし、逆畳み込み処理は実行しません。
サポートされる装置関数タイプ
以下の装置関数（フィルタ関数）タイプがサポートされています。
ガウス関数 (
gauss)ローレンツ関数 (
lorentzまたはlorentzian)
データ前処理
以下のデータ前処理オプションを組み合わせることができます。
スムージング (
average/convolve):SmoothingBySimpleAverage()(単純移動平均) またはSmoothingByPolynomialFit()(Savitzky-Golay風の多項式フィッティング) によるデータ平滑化。smoothmode引数で指定します。
データ拡張 (
extend):extend_smooth()関数により、データの両端にゼロ領域や線形補間領域を追加します。特に畳み込みやFFTによる逆畳み込みにおいて、エッジ効果を低減するために使用されます。
使用ライブラリ
このプログラムは、以下のPythonライブラリに依存しています。
標準ライブラリ
sys: システム固有のパラメータと関数を提供します。os.path: パス名を操作するための関数を提供します。argparse: コマンドライン引数を解析するためのモジュールです。csv: CSV (Comma Separated Values) ファイルを読み書きするためのモジュールです。re: 正規表現操作を提供します。math: 数学関数を提供します (exp,log,sqrt)。
非標準ライブラリ
numpy(数値計算ライブラリ): 大規模な多次元配列と行列の操作、およびこれらに対する高度な数学関数を提供します。pandas(データ分析ライブラリ): データ構造 (DataFrameなど) とデータ分析ツールを提供します。scipy(科学計算ライブラリ): 科学的および技術的な計算のためのアルゴリズムとツールを提供します。特に、信号処理 (scipy.signal) や最適化 (scipy.optimize) の機能を使用しています。matplotlib(プロットライブラリ): 静的、アニメーション、対話的な可視化を作成します。tklib(カスタムライブラリ): このプログラムは、ファイルI/O (tkVariousData)、アプリケーション管理 (tkApplication)、グラフィックイベント処理 (tkPlotEvent) など、特定の機能のためにtklibライブラリに依存しています。このライブラリは、このプログラムの実行環境で利用可能である必要があります。
インストール方法
deconvolution.py を実行するには、非標準ライブラリをインストールする必要があります。
pip install numpy pandas scipy matplotlib
tklib はカスタムライブラリであり、一般的な pip コマンドではインストールできません。別途提供されるか、環境に配置する必要があります。
コマンドライン引数
プログラムは argparse モジュールを使用してコマンドライン引数を処理します。
python deconvolution.py --help
上記コマンドで表示されるヘルプメッセージは以下の通りです。
usage: deconvolution.py [-h] [--mode {plot,gs,jacobi,sp,ridge,nnls,convolve,fft,deconvolve}] [--xmin XMIN] [--xmax XMAX] [--xlabel XLABEL] [--ylabel YLABEL] [--flip-x {0,1}] [--func-type {gauss,lorentz,lorentzian}] [--wa WA] [--grange GRANGE] [--alpha ALPHA] [--dump DUMP] [--nmaxiter NMAXITER] [--eps EPS] [--nsmooth NSMOOTH] [--zero-correction {0,1}] [--convmode {same,full,valid}] [--smoothmode SMOOTHMODE] [--kzero KZERO] [--klin KLIN] [--nnls-maxiter NNLS_MAXITER] [--nnls-penalty {ridge,smooth}] [--template TEMPLATE] [--figsize WIDTH HEIGHT] [infile]
Deconvolute one-dimensional spectral data with a Gaussian or Lorentzian instrumental function.
positional arguments:
infile Input data file readable by tkVariousData. (default: pes.csv)
options:
-h, --help show this help message and exit
--mode {plot,gs,jacobi,sp,ridge,nnls,convolve,fft,deconvolve}
Calculation mode. nnls is non-negative least squares using the instrumental function as basis functions. (default: gs)
--xmin XMIN Minimum x value. Use '*' or 'none' for no lower limit. (default: -100000000.0)
--xmax XMAX Maximum x value. Use '*' or 'none' for no upper limit. (default: 100000000.0)
--xlabel XLABEL X column index or label. (default: 0)
--ylabel YLABEL Y column index or label. (default: 1)
--flip-x {0,1} If 1, multiply input x values by -1 before range filtering, plotting, and calculation. (default: 0)
--func-type {gauss,lorentz,lorentzian}
Instrumental/window-function type. (default: gauss)
--wa WA Half width at half maximum of the instrumental function. (default: 0.12)
--grange GRANGE Window-function range in units of wa for sampled filters. (default: 2.0)
--alpha ALPHA Regularization coefficient for ridge/sp modes. (default: 0.1)
--dump DUMP Damping factor for Jacobi/Gauss-Seidel iterations. (default: 1.0)
--nmaxiter NMAXITER Maximum iteration count for Jacobi/Gauss-Seidel modes. (default: 300)
--eps EPS Relative convergence criterion for iterative modes. (default: 0.0001)
--nsmooth NSMOOTH Smoothing width used by polynomial smoothing. (default: 5)
--zero-correction {0,1}
Clip negative values to zero after each iterative update. (default: 0)
--convmode {same,full,valid}
np.convolve mode used when smoothmode contains convolve. (default: same)
--smoothmode SMOOTHMODE
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. (default: None)
--kzero KZERO Number of half-window widths used for zero extension. (default: 5)
--klin KLIN Number of half-window widths used for linear edge shaping. (default: 5)
--nnls-maxiter NNLS_MAXITER
Maximum iteration count for nnls / bounded least-squares solvers. (default: None)
--nnls-penalty {ridge,smooth}
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. (default: ridge)
--template TEMPLATE Excel template used by tkVariousData.to_excel. (default: StandardGraph.xlsm)
--figsize WIDTH HEIGHT
Matplotlib figure size. (default: [6.0, 6.0])
引数詳細
infile説明:
tkVariousDataで読み取り可能な入力データファイル。デフォルト値:
pes.csv
--mode説明: 実行する計算モード。
nnlsは、装置関数を基底関数とする非負最小二乗法です。選択肢:
plot,gs,jacobi,sp,ridge,nnls,convolve,fft,deconvolveデフォルト値:
gs
--xmin説明: X軸の最小値。下限を設定しない場合は
*またはnoneを使用します。デフォルト値:
-1.0e8
--xmax説明: X軸の最大値。上限を設定しない場合は
*またはnoneを使用します。デフォルト値:
1.0e8
--xlabel説明: X軸のデータ列のインデックスまたはラベル。
デフォルト値:
0
--ylabel説明: Y軸のデータ列のインデックスまたはラベル。
デフォルト値:
1
--flip-x説明:
1を指定すると、範囲フィルタリング、プロット、計算の前にX入力値を-1倍します。選択肢:
0,1デフォルト値:
0
--func-type説明: 装置関数（窓関数）のタイプ。
選択肢:
gauss,lorentz,lorentzianデフォルト値:
gauss
--wa説明: 装置関数の半値半幅 (HWHM)。
デフォルト値:
0.12
--grange説明: サンプリングされたフィルタの窓関数の範囲を
waの単位で指定。デフォルト値:
2.0
--alpha説明: リッジ/SPモードの正則化係数。
デフォルト値:
0.1
--dump説明: Jacobi/Gauss-Seidel反復の減衰因子。
デフォルト値:
1.0
--nmaxiter説明: Jacobi/Gauss-Seidelモードの最大反復回数。
デフォルト値:
300
--eps説明: 反復モードの相対収束基準。
デフォルト値:
1.0e-4
--nsmooth説明: 多項式スムージングで使用されるスムージング幅。
デフォルト値:
5
--zero-correction説明: 各反復更新後に負の値をゼロにクリップするかどうか。
選択肢:
0,1デフォルト値:
0
--convmode説明:
smoothmodeにconvolveが含まれる場合に使用されるnp.convolve()のモード。選択肢:
same,full,validデフォルト値:
same
--smoothmode説明:
average,extend, およびconvolveの組み合わせ。前処理を行わない場合はnoneまたは空文字列を使用します。デフォルトはモード依存で、convolve/fftの場合はconvolve+extend、それ以外は前処理なしです。デフォルト値: None (モード依存)
--kzero説明: ゼロ拡張に使用される半窓幅の数。
デフォルト値:
5
--klin説明: 線形エッジ整形に使用される半窓幅の数。
デフォルト値:
5
--nnls-maxiter説明: NNLS / 有界最小二乗ソルバーの最大反復回数。
デフォルト値: None
--nnls-penalty説明:
--mode nnlsおよび--alpha > 0の場合に使用される正則化ペナルティ。ridgeは係数のサイズにペナルティを課し、smoothは1次差分にペナルティを課します。正則化なしの純粋なNNLSには--alpha 0を使用します。選択肢:
ridge,smoothデフォルト値:
ridge
--template説明:
tkVariousData.to_excel()で使用されるExcelテンプレート。デフォルト値:
StandardGraph.xlsm
--figsize説明: Matplotlibの図のサイズを
WIDTHとHEIGHTで指定します。デフォルト値:
[6.0, 6.0]
入出力ファイル形式
入力
ファイル形式:
infileで指定される入力データファイルは、tklib.tkvariousdata.tkVariousDataオブジェクトが読み取り可能な形式である必要があります (例: CSV)。データ構造: X軸およびY軸のデータは、コマンドライン引数
--xlabelおよび--ylabelで指定された列インデックスまたはラベルに従って読み込まれます。
出力
グラフ表示:
実行時には、Matplotlibウィンドウに処理結果のグラフが表示されます。
反復計算モード (Jacobi, Gauss-Seidel) では、反復ごとにグラフが更新されます。
Excelファイル:
処理されたデータは、入力ファイル名を基にしたExcelファイル (
.xlsx) に保存されます。デフォルトの出力ファイル名は{infilebody}-deconvoluted.xlsxです。出力されるデータフレームには、通常、
x,y(input),y(deconvoluted)の列が含まれます。--template引数で指定されたExcelテンプレートを使用してファイルが生成されます。
ログファイル:
標準出力の内容は、入力ファイル名を基にしたログファイルにリダイレクトされ保存されます。このログファイルには、プログラムの実行情報、パラメータ、およびデバッグ情報が含まれます。デフォルトのログファイル名は入力ファイル名と同じです。
主要な関数
プログラム内で定義されている主要な関数の概要は以下の通りです。
parse_xlimit(value)X軸範囲の制限値を解析します。
*またはnoneは制限なしを意味します。
parse_label(value)データ列セレクタを整数インデックスまたはラベル文字列として解析します。
build_arg_parser()コマンドライン引数パーサーを構築し、プログラムの説明と引数を定義します。
usage(app=None)ヘルプメッセージを表示します。
apply_args_to_globals(args)解析されたコマンドライン引数をグローバル変数に適用します。
Gaussian(x, x0, whalf)ガウス関数を計算します。
Lorentzian(x, x0, whalf)ローレンツ関数を計算します。
convolve_func(x, x0, whalf=None)func_typeに応じてガウス関数またはローレンツ関数を返します。
Hij(xstep, Wa, Grange, i, j)畳み込みカーネルの要素
H(i, j)を計算します。
Ridge(x, y, m, lsqfunc, w, nGmax, alpha = 0.0, is_print = False)リッジ回帰を実行し、係数を計算します。
Smoothing_Penalty_Regression(x, y, m, lsqfunc, w, nGmax, alpha = 0.0, nsmooth = 1, is_print = False)スムージングペナルティ回帰を実行し、係数を計算します。
make_wf(Wa, Grange, xstep)サンプリングされた窓関数（フィルタ関数）を生成します。
convolve(xraw, yraw, ywf, **kwargs)numpy.convolve()を使用してデータを畳み込みます。
extend_smooth(x, y, nzero, nlin, xstep = 0.0)データにゼロおよび線形のエッジ拡張を適用します。
SmoothingBySimpleAverage(y, n)単純移動平均法でデータを平滑化します。
SmoothingByPolynomialFit(y, n)多項式フィッティング（Savitzky-Golayフィルターに類似）でデータを平滑化します。
deconvolute_fft(xRaw, yRaw, xG, yG)FFTを用いて逆畳み込みを実行します。
deconvolute_deconvolve(xRaw, yRaw, xG, yG)scipy.signal.deconvolve()を用いて逆畳み込みを実行します。
make_edge_correction_weights(xRaw, whalf, nGmax, w_floor=1.0e-6)リッジ回帰などで使用されるエッジ補正重みを計算します。
edge_weight_at(w, n, nGmax, j, w_floor=1.0e-6)特定の係数インデックス
jの対称的なエッジ補正重みを返します。
make_basis_matrix(xRaw, whalf, w=None, nGmax=None)装置関数を基にした密な基底行列
Kを構築します。
make_first_difference_matrix(n)スムージングペナルティ用の1次差分行列
Dを返します。
deconvolute_nnls(xRaw, yRaw, xG, yG)非負最小二乗法 (NNLS) を用いて逆畳み込みを実行します。
deconvolute_ridge(xRaw, yRaw, xG, yG)リッジ回帰を用いて逆畳み込みを実行します。
deconvolute_sp(xRaw, yRaw, xG, yG, nsmooth)スムージングペナルティ法を用いて逆畳み込みを実行します。
_set_iteration_ylim(ax0, yRaw, y)反復プロットのY軸範囲を安定して設定する内部ヘルパー関数。
_init_iteration_plot(fig, ax0, xRaw, yRaw, y, xgmin, xgmax)反復更新グラフを一度だけ初期化する内部ヘルパー関数。
_update_iteration_plot(fig, ax0, line_updated, yRaw, y)反復逆畳み込み中にYデータのみを更新する内部ヘルパー関数。
deconvolute_jacobi(xRaw, yRaw, xG, yG, fig, ax)Jacobi法を用いて逆畳み込みを実行し、結果をプロットします。
deconvolute_gauss_seidel(xRaw, yRaw, xG, yG, fig, ax)Gauss-Seidel法を用いて逆畳み込みを実行し、結果をプロットします。
plot_input_data(xRaw, yRaw, label_x='', label_y='')選択された入力データのみをプロットします。
main()プログラムのエントリポイントです。引数解析、データ読み込み、前処理、逆畳み込み実行、結果のプロットと保存を調整します。
注意事項
numpy.convolve()の使用:畳み込み関数を規格化するため、フィルタ関数
ywfの和で割る必要があります。戻り値のリストは、フィルタ関数の点数の1/2だけ、前後に拡張される場合があります。
mode="same"を指定すると、戻り値のリストの範囲は入力関数の範囲に一致させられます。元のデータに逆畳み込みで戻すには、拡張部分のデータも必要となる可能性があります。
scipy.signal.deconvolve()の使用:scipy.signal.deconvolve()は入力データの質に非常に敏感です。以下の条件を満たさないと、適切な結果が得られない可能性があります。フィルタ関数の幅は入力データよりも短い必要があります。
フィルタ関数の値はある程度の値よりも大きい必要があります（ゼロに近い値があると問題が発生する場合があります）。
入力データの前には、フィルタ関数の幅よりも大きいゼロ領域が必要です。
入力データにノイズがある場合は、適切に平滑化処理が必要です。