deconvolution.py ダウンロード/コピー
deconvolution.py
deconvolution.py
1import sys
2import os.path
3import argparse
4import csv
5import re
6from math import exp, log, sqrt
7import numpy as np
8import pandas as pd
9import scipy.signal
10import scipy.optimize
11import matplotlib.pyplot as plt
12import matplotlib.widgets as wg
13
14
15from tklib.tkutils import replace_path, getarg, getintarg, getfloatarg, pint, pfloat
16from tklib.tksci.tksci import kB, e
17from tklib.tkvariousdata import tkVariousData
18from tklib.tkapplication import tkApplication
19from tklib.tkgraphic.tkplotevent import tkPlotEvent
20
21
22"""
23概要:
24 スペクトルデータの逆畳み込みを実行するスクリプト。
25詳細説明:
26 このスクリプトは、1次元スペクトルデータに対してガウス関数またはローレンツ関数を
27 用いた逆畳み込み処理を提供します。Jacobi法、Gauss-Seidel法、Smoothing penalty法、
28 Ridge回帰、非負最小二乗法（NNLS）、FFT法、scipy.signal.deconvolve
29 など、複数の逆畳み込みアルゴリズムをサポートしています。
30
31主な機能:
32 Jacobi法による反復逆畳み込み、Gauss-Seidel法による反復逆畳み込み、
33 Smoothing penalty法による正則化逆畳み込み、Ridge回帰による正則化逆畳み込み、
34 非負最小二乗法（NNLS）による制約付き逆畳み込み、FFT（高速フーリエ変換）を用いた逆畳み込み、
35 scipy.signal.deconvolve を使用した直接逆畳み込み、
36 numpy.convolve を用いた畳み込みによる動作確認、入力データの平滑化、拡張、整形オプション。
37
38注意点:
39 畳み込み関数を規格化するため、フィルター関数（ywf）の和で割る必要があります。
40 numpy.convolve の戻り値はフィルター関数の点数の約1/2だけ前後に拡張されます。
41 mode='same' を指定すると、戻り値の範囲は入力関数の範囲に一致します。
42 逆畳み込みで元のデータに戻すには、拡張部分のデータも必要となる場合があります。
43 scipy.signal.deconvolve は入力データの質に非常に敏感です。以下の条件を満足しないと、
44 期待する結果が得られない場合があります。フィルター関数の幅は入力データよりも短いこと。
45 フィルター関数の値がある程度の値よりも大きいこと（0に近い値があると問題になる可能性）。
46 入力データの前には、フィルター関数の幅よりも大きいゼロ領域が必要な場合があること。
47 入力データにノイズがある場合は、適切に平滑化処理が必要となること。
48
49参考文献:
50 「科学計測のための波形データ処理」、南茂夫 編著 (CQ出版社, 1986)
51関連リンク:
52 deconvolution_usage
53"""
54
55
56#===================
57# Global parametrs
58#===================
59app = tkApplication()
60
61pi = 3.1415926535
62
63# Input file
64infile = "pes.csv"
65template = "StandardGraph.xlsm"
66
67
68# Analysis range
69xlabel = 0
70ylabel = 1
71xmin = -1.0e8 # -4.5 # eV
72xmax = 1.0e8 # 2.0 # eV
73
74# mode = [plot|gs|jacobi|sp|ridge|nnls|convolve|fft|deconvolve]
75mode = 'gs'
76
77# Filter parameters for window function
78func_type = 'gauss'
79Wa = 0.12 # eV
80Grange = 2.0
81
82# for update graph
83sleeptime = 0.05
84
85# for Ridge regression
86alpha = 0.1
87
88# for NNLS / non-negative regularized deconvolution
89nnls_maxiter = None
90nnls_penalty = 'ridge'
91
92# for Jacobi/Gauss-Seidel method
93dump = 1.0
94nmaxiter = 300
95eps = 1.0e-4
96nsmooth = 5
97
98zero_correction = 0
99
100# only for 'convolve': convmode = [same|full|]
101convmode = 'same'
102# smoothmode = [convolve|extend|average]
103smoothmode = 'convolve+extend'
104
105# only for 'convolve' and 'fft': 'exnted data' parameters
106kzero = 5
107klin = 5
108
109figsize = (6, 6)
110
111# x-axis sign conversion. Useful for spectra stored as binding energy
112# when plotting/analyzing on the opposite energy axis.
113flip_x = 0
114
115#===================
116# 起動時引数
117#===================
118scriptname = sys.argv[0]
119
120
121def parse_xlimit(value):
122 """
123 概要:
124 x軸の範囲制限を解析します。
125 詳細説明:
126 '*'または'none'は無制限を示し、数値文字列は浮動小数点数に変換されます。
127 これは、xmin/xmaxがオープンに設定できる以前の動作を維持します。
128 引数:
129 :param value: 解析する値。数値、文字列、またはNone。
130 :type value: any
131 戻り値:
132 :returns: 解析されたx軸の制限値、またはNone。
133 :rtype: float or None
134 """
135 if value is None:
136 return None
137 if isinstance(value, (int, float)):
138 return float(value)
139 text = str(value).strip()
140 if text == '' or text == '*' or text.lower() in ('none', 'null'):
141 return None
142 return float(text)
143
144
145def parse_label(value):
146 """
147 概要:
148 データ列セレクタを解析します。
149 詳細説明:
150 整数インデックスまたはラベル文字列として解析します。
151 引数:
152 :param value: 解析する値。整数または文字列。
153 :type value: int or str
154 戻り値:
155 :returns: 解析された列インデックス（整数）またはラベル（文字列）。
156 :rtype: int or str
157 """
158 if isinstance(value, int):
159 return value
160 text = str(value).strip()
161 try:
162 return int(text)
163 except ValueError:
164 return text
165
166
167def build_arg_parser():
168 """
169 概要:
170 コマンドライン引数パーサーを構築します。
171 詳細説明:
172 1次元スペクトルデータをガウス関数またはローレンツ関数型の装置関数で
173 逆畳み込みするための引数パーサーを作成します。
174 戻り値:
175 :returns: 設定済みのArgument Parserオブジェクト。
176 :rtype: argparse.ArgumentParser
177 """
178 parser = argparse.ArgumentParser(
179 description=(
180 "Deconvolute one-dimensional spectral data with a Gaussian or "
181 "Lorentzian instrumental function."
182 ),
183 formatter_class=argparse.ArgumentDefaultsHelpFormatter,
184 )
185
186 parser.add_argument(
187 'infile',
188 nargs='?',
189 default=infile,
190 help='Input data file readable by tkVariousData.',
191 )
192 parser.add_argument(
193 '--mode',
194 choices=['plot', 'gs', 'jacobi', 'sp', 'ridge', 'nnls', 'convolve', 'fft', 'deconvolve'],
195 default=mode,
196 help=(
197 'Calculation mode. nnls is non-negative least squares using the '
198 'instrumental function as basis functions.'
199 ),
200 )
201 parser.add_argument('--xmin', type=parse_xlimit, default=xmin,
202 help="Minimum x value. Use '*' or 'none' for no lower limit.")
203 parser.add_argument('--xmax', type=parse_xlimit, default=xmax,
204 help="Maximum x value. Use '*' or 'none' for no upper limit.")
205 parser.add_argument('--xlabel', type=parse_label, default=xlabel,
206 help='X column index or label.')
207 parser.add_argument('--ylabel', type=parse_label, default=ylabel,
208 help='Y column index or label.')
209 parser.add_argument('--flip-x', '--flip_x', dest='flip_x', type=int, choices=[0, 1], default=flip_x,
210 help='If 1, multiply input x values by -1 before range filtering, plotting, and calculation.')
211
212 parser.add_argument('--func-type', choices=['gauss', 'lorentz', 'lorentzian'],
213 default=func_type,
214 help='Instrumental/window-function type.')
215 parser.add_argument('--wa', type=float, default=Wa,
216 help='Half width at half maximum of the instrumental function.')
217 parser.add_argument('--grange', type=float, default=Grange,
218 help='Window-function range in units of wa for sampled filters.')
219
220 parser.add_argument('--alpha', type=float, default=alpha,
221 help='Regularization coefficient for ridge/sp modes.')
222 parser.add_argument('--dump', type=float, default=dump,
223 help='Damping factor for Jacobi/Gauss-Seidel iterations.')
224 parser.add_argument('--nmaxiter', type=int, default=nmaxiter,
225 help='Maximum iteration count for Jacobi/Gauss-Seidel modes.')
226 parser.add_argument('--eps', type=float, default=eps,
227 help='Relative convergence criterion for iterative modes.')
228 parser.add_argument('--nsmooth', type=int, default=nsmooth,
229 help='Smoothing width used by polynomial smoothing.')
230 parser.add_argument('--zero-correction', type=int, choices=[0, 1], default=zero_correction,
231 help='Clip negative values to zero after each iterative update.')
232
233 parser.add_argument('--convmode', choices=['same', 'full', 'valid'], default=convmode,
234 help='np.convolve mode used when smoothmode contains convolve.')
235 parser.add_argument('--smoothmode', default=None,
236 help=(
237 "Combination of 'average', 'extend', and 'convolve'. "
238 "Use 'none' or an empty string for no preprocessing. "
239 "Default is mode-dependent: 'convolve+extend' for "
240 "convolve/fft, otherwise no preprocessing."
241 ))
242 parser.add_argument('--kzero', type=int, default=kzero,
243 help='Number of half-window widths used for zero extension.')
244 parser.add_argument('--klin', type=int, default=klin,
245 help='Number of half-window widths used for linear edge shaping.')
246
247 parser.add_argument('--nnls-maxiter', type=int, default=nnls_maxiter,
248 help='Maximum iteration count for nnls / bounded least-squares solvers.')
249 parser.add_argument('--nnls-penalty', choices=['ridge', 'smooth'], default=nnls_penalty,
250 help=(
251 'Regularization penalty used when --mode nnls and --alpha > 0. '
252 'ridge penalizes coefficient size; smooth penalizes first differences. '
253 'Use --alpha 0 for pure NNLS without regularization.'
254 ))
255 parser.add_argument('--template', default=template,
256 help='Excel template used by tkVariousData.to_excel.')
257 parser.add_argument('--figsize', type=float, nargs=2, default=figsize,
258 metavar=('WIDTH', 'HEIGHT'),
259 help='Matplotlib figure size.')
260 return parser
261
262
263def usage(app=None):
264 """
265 概要:
266 スクリプトの使用方法を表示します。
267 詳細説明:
268 build_arg_parser() によって生成されたヘルプメッセージを出力します。
269 引数:
270 :param app: Tkinterアプリケーションオブジェクト。現在は使用されません。
271 :type app: tkApplication or None
272 """
273 print(build_arg_parser().format_help())
274
275
276def apply_args_to_globals(args):
277 """
278 概要:
279 解析されたコマンドライン引数をグローバル変数に適用します。
280 詳細説明:
281 argsオブジェクトから取得した値で、スクリプト全体の動作を制御するグローバル変数を更新します。
282 引数:
283 :param args: コマンドライン引数を格納したNamespaceオブジェクト。
284 :type args: argparse.Namespace
285 """
286 global infile, template, xlabel, ylabel, xmin, xmax, mode, flip_x
287 global func_type, Wa, Grange, sleeptime, alpha, nnls_maxiter, nnls_penalty
288 global dump, nmaxiter, eps, nsmooth, zero_correction
289 global convmode, smoothmode, kzero, klin, figsize
290
291 infile = args.infile
292 template = args.template
293 xlabel = args.xlabel
294 ylabel = args.ylabel
295 xmin = args.xmin
296 xmax = args.xmax
297 mode = args.mode
298 flip_x = args.flip_x
299
300 func_type = 'lorentz' if args.func_type == 'lorentzian' else args.func_type
301 Wa = args.wa
302 Grange = args.grange
303 alpha = args.alpha
304 nnls_maxiter = args.nnls_maxiter
305 nnls_penalty = args.nnls_penalty
306
307 dump = args.dump
308 nmaxiter = args.nmaxiter
309 eps = args.eps
310 nsmooth = args.nsmooth
311 zero_correction = args.zero_correction
312
313 convmode = args.convmode
314 if args.smoothmode is None:
315 smoothmode = 'convolve+extend' if mode in ('convolve', 'fft') else ''
316 elif args.smoothmode.lower() == 'none':
317 smoothmode = ''
318 else:
319 smoothmode = args.smoothmode
320
321 kzero = args.kzero
322 klin = args.klin
323 figsize = tuple(args.figsize)
324
325
326#header, ext = os.path.splitext(infile)
327#outfile = header + '-deconvoluted.xlsx'
328
329def Gaussian(x, x0, whalf):
330 """
331 概要:
332 ガウス関数を計算します。
333 詳細説明:
334 中心 x0、半値半幅 whalf のガウス関数を返します。
335 引数:
336 :param x: 評価点。
337 :type x: float
338 :param x0: 関数の中心。
339 :type x0: float
340 :param whalf: 半値半幅（HWHM）。
341 :type whalf: float
342 戻り値:
343 :returns: ガウス関数の値。
344 :rtype: float
345 """
346#A = 1/whalf * sqrt(ln2 / pi)
347 A = 0.469718639 / whalf
348#a = whalf / sqrt(ln2)
349 a = whalf / 0.832554611
350 X = (x - x0) / a
351 return A * exp(-X*X)
352
353def Lorentzian(x, x0, whalf):
354 """
355 概要:
356 ローレンツ関数を計算します。
357 詳細説明:
358 中心 x0、半値半幅 whalf のローレンツ関数を返します。
359 引数:
360 :param x: 評価点。
361 :type x: float
362 :param x0: 関数の中心。
363 :type x0: float
364 :param whalf: 半値半幅（HWHM）。
365 :type whalf: float
366 戻り値:
367 :returns: ローレンツ関数の値。
368 :rtype: float
369 """
370#A = 1/whalf/pi
371 A = 1.0 / whalf / pi
372 X = (x - x0) / whalf
373 return A / (1.0 + X * X)
374
375def convolve_func(x, x0, whalf=None):
376 """
377 概要:
378 設定された関数タイプに基づいて畳み込み関数を計算します。
379 詳細説明:
380 func_type グローバル変数に応じてガウス関数またはローレンツ関数を呼び出します。
381 whalf が指定されない場合は、グローバル変数 Wa が使用されます。
382 引数:
383 :param x: 評価点。
384 :type x: float
385 :param x0: 関数の中心。
386 :type x0: float
387 :param whalf: 半値半幅（HWHM）。指定されない場合はグローバルなWaが使用されます。
388 :type whalf: float or None
389 戻り値:
390 :returns: 畳み込み関数の値。
391 :rtype: float
392 """
393 if whalf is None:
394 whalf = Wa
395 if func_type == 'gauss':
396 return Gaussian(x, x0, whalf)
397 else:
398 return Lorentzian(x, x0, whalf)
399
400def Hij(xstep, Wa, Grange, i, j):
401 """
402 概要:
403 畳み込み行列の要素 Hij を計算します。
404 詳細説明:
405 ステップサイズ xstep、半値半幅 Wa、範囲 Grange を用いて、
406 インデックス i と j に対応する畳み込み関数の値を計算します。
407 引数:
408 :param xstep: x軸のステップサイズ。
409 :type xstep: float
410 :param Wa: 半値半幅（HWHM）。
411 :type Wa: float
412 :param Grange: 畳み込み関数の範囲を Wa の単位で指定。
413 :type Grange: float
414 :param i: 行インデックス。
415 :type i: int
416 :param j: 列インデックス。
417 :type j: int
418 戻り値:
419 :returns: Hij 行列の要素の値。
420 :rtype: float
421 """
422 return convolve_func(xstep * i, xstep * j, Wa)
423
424def Ridge(x, y, m, lsqfunc, w, nGmax, alpha = 0.0, is_print = False):
425 """
426 概要:
427 Ridge回帰により係数を計算します。
428 詳細説明:
429 最小二乗関数 lsqfunc とエッジ補正重み w を用いて、
430 観測データ y から係数 ci を推定します。正則化パラメータ alpha を使用します。
431 引数:
432 :param x: x軸のデータ点配列。
433 :type x: numpy.ndarray or list
434 :param y: 観測データ配列。
435 :type y: numpy.ndarray or list
436 :param m: 係数の数（データ点の数）。
437 :type m: int
438 :param lsqfunc: 最小二乗法で使用する関数。2つの引数 (x1, x2) を取ります。
439 :type lsqfunc: callable
440 :param w: エッジ補正重み配列。
441 :type w: list
442 :param nGmax: 畳み込み関数の最大インデックス範囲。
443 :type nGmax: int
444 :param alpha: 正則化パラメータ（リッジパラメータ）。デフォルトは0.0。
445 :type alpha: float
446 :param is_print: 中間結果を出力するかどうかのフラグ。デフォルトはFalse。
447 :type is_print: bool
448 戻り値:
449 :returns: 推定された係数のリスト。
450 :rtype: list
451 """
452 n = len(x)
453 Si = np.empty([m, 1])
454 Sij = np.empty([m, m])
455
456 def cal_w(j):
457 if j <= nGmax:
458 wj = w[j]
459 elif n - 1 - j < nGmax:
460 wj = w[n-1-j]
461 else:
462 wj = w[nGmax - 1]
463 return wj
464
465 for l in range(0, m):
466 v = 0.0
467 for i in range(0, n):
468 v += y[i] * lsqfunc(x[l], x[i])
469# Si[l, 0] = v
470 Si[l, 0] = v / cal_w(l)
471
472 for j in range(0, m):
473 for l in range(j, m):
474 v = 0.0
475 for i in range(0, n):
476 v += lsqfunc(x[j], x[i]) * lsqfunc(x[l], x[i])
477# Sij[j, l] = Sij[l, j] = v
478 Sij[j, l] = Sij[l, j] = v / cal_w(j) / cal_w(l)
479 Sij[j, j] += alpha
480
481 if is_print:
482 print("Vector and Matrix:")
483 print("Si=")
484 print(Si)
485 print("Sij=")
486 print(Sij)
487 print("")
488
489 ci = np.linalg.solve(Sij, Si)
490 ci = ci.transpose().tolist()
491 return ci[0]
492
493def Smoothing_Penalty_Regression(x, y, m, lsqfunc, w, nGmax, alpha = 0.0, nsmooth = 1, is_print = False):
494 """
495 概要:
496 Smoothing Penalty回帰により係数を計算します。
497 詳細説明:
498 最小二乗関数 lsqfunc とエッジ補正重み w を用いて、
499 観測データ y から係数 ci を推定します。正則化パラメータ alpha と平滑化幅 nsmooth を使用します。
500 引数:
501 :param x: x軸のデータ点配列。
502 :type x: numpy.ndarray or list
503 :param y: 観測データ配列。
504 :type y: numpy.ndarray or list
505 :param m: 係数の数（データ点の数）。
506 :type m: int
507 :param lsqfunc: 最小二乗法で使用する関数。2つの引数 (x1, x2) を取ります。
508 :type lsqfunc: callable
509 :param w: エッジ補正重み配列。
510 :type w: list
511 :param nGmax: 畳み込み関数の最大インデックス範囲。
512 :type nGmax: int
513 :param alpha: 正則化パラメータ（ペナルティパラメータ）。デフォルトは0.0。
514 :type alpha: float
515 :param nsmooth: 平滑化の幅。デフォルトは1。
516 :type nsmooth: int
517 :param is_print: 中間結果を出力するかどうかのフラグ。デフォルトはFalse。
518 :type is_print: bool
519 戻り値:
520 :returns: 推定された係数のリスト。
521 :rtype: list
522 """
523 n = len(x)
524 Si = np.empty([m, 1])
525 Sij = np.empty([m, m])
526
527 def cal_w(j):
528 if j <= nGmax:
529 wj = w[j]
530 elif n - 1 - j < nGmax:
531 wj = w[n-1-j]
532 else:
533 wj = w[nGmax - 1]
534 return max([wj, 1.0e-6])
535# return wj
536
537 for l in range(0, m):
538 v = 0.0
539 for i in range(0, n):
540 v += y[i] * lsqfunc(x[l], x[i])
541# Si[l, 0] = v
542 Si[l, 0] = v / cal_w(l)
543
544 for j in range(0, m):
545 for l in range(j, m):
546 v = 0.0
547 for i in range(0, n):
548 v += lsqfunc(x[j], x[i]) * lsqfunc(x[l], x[i])
549# Sij[j, l] = Sij[l, j] = v
550 Sij[j, l] = Sij[l, j] = v / cal_w(j) / cal_w(l)
551
552 nband = int(nsmooth / 2 + 1.0001)
553 use_ridge = 1
554 if use_ridge > 1:
555 for j in range(0, m):
556 Sij[j, j] += alpha
557 elif use_ridge > 0:
558 for j in range(0, m):
559 j0 = j - nband
560 j1 = j + nband
561 if j0 < 0:
562 j0 = 0
563 if j1 >= m:
564 j1 = m
565 for k in range(j0, j):
566 Sij[j, k] -= alpha
567 for k in range(j+1, j1):
568 Sij[j, k] -= alpha
569
570 if is_print:
571 print("Vector and Matrix:")
572 print("Si=")
573 print(Si)
574 print("Sij=")
575 print(Sij)
576 print("")
577
578 ci = np.linalg.solve(Sij, Si)
579 ci = ci.transpose().tolist()
580 return ci[0]
581
582def make_wf(Wa, Grange, xstep):
583 """
584 概要:
585 サンプリングされたウィンドウ関数（フィルター関数）を作成します。
586 詳細説明:
587 半値半幅 Wa と範囲 Grange を用いて、xstep 間隔でサンプリングされた
588 ウィンドウ関数 yG を生成し、その積算値 SG で正規化します。
589 x軸が降順であっても、ウィンドウ関数が空にならないように
590 xstep の絶対値が使用されます。
591 引数:
592 :param Wa: 半値半幅（HWHM）。
593 :type Wa: float
594 :param Grange: 畳み込み関数の範囲を Wa の単位で指定。
595 :type Grange: float
596 :param xstep: x軸のステップサイズ。
597 :type xstep: float
598 戻り値:
599 :returns: x軸のサンプリング点リストと、正規化されたウィンドウ関数値リストのタプル。
600 :rtype: tuple of (list, list)
601 例外:
602 :raises ValueError: xstepが非ゼロでない場合、または正規化積算値SGがゼロ以下の場合。
603 """
604 # The sampled window function assumes a positive grid spacing.
605 # Experimental spectra, especially XPS/HAXPES binding-energy axes,
606 # are often stored in descending x order. Use abs(xstep) here so the
607 # window does not become empty when the input axis is descending.
608 dx = abs(xstep)
609 if dx <= 0.0:
610 raise ValueError("xstep must be non-zero in make_wf().")
611
612 ixG0 = int(Grange * Wa / dx + 1.0001)
613 ixGmax = 2 * ixG0
614 nxGmax = ixG0 + 1
615 xG0 = ixG0 * dx
616
617 xG = []
618 yG = []
619 for i in range(ixGmax+1):
620 x = i * dx
621 xG.append(x)
622 yG.append(convolve_func(x, xG0, Wa))
623
624 SG = 0.0
625 for i in range(len(yG)-1):
626 SG += (yG[i] + yG[i+1]) / 2.0 * (xG[i+1] - xG[i])
627
628 if SG <= 0.0:
629 raise ValueError(
630 f"Invalid sampled window function: SG={SG}. "
631 f"Check wa={Wa}, grange={Grange}, xstep={xstep}."
632 )
633
634 for i in range(ixGmax+1):
635 yG[i] /= SG
636
637 print(" Range: {} in width".format(Grange * Wa))
638 print(" i range: {} - {} at center {}".format(0,ixGmax, xG0))
639 print(" ixGmax = ", ixGmax)
640 print(" SG = ", SG)
641
642 return xG, yG
643
644
645def convolve(xraw, yraw, ywf, **kwargs):
646 """
647 概要:
648 データを畳み込みます。
649 詳細説明:
650 numpy.convolve を使用して生データ yraw をウィンドウ関数 ywf で畳み込みます。
651 畳み込み関数の和で結果を正規化し、必要に応じてx軸を拡張します。
652 引数:
653 :param xraw: 生データのx軸配列。
654 :type xraw: numpy.ndarray or list
655 :param yraw: 生データのy軸配列。
656 :type yraw: numpy.ndarray or list
657 :param ywf: ウィンドウ関数（フィルター関数）のy値配列。
658 :type ywf: numpy.ndarray or list
659 :param kwargs: numpy.convolve に渡される追加のキーワード引数。
660 :type kwargs: dict
661 戻り値:
662 :returns: 畳み込み後のx軸配列とy軸配列のタプル。
663 :rtype: tuple of (numpy.ndarray, numpy.ndarray)
664 """
665 yconv = np.convolve(yraw, ywf, **kwargs) / sum(ywf)
666 n_new = len(yconv)
667 dn = n_new - len(yraw)
668 if dn > 0:
669 offset = int(dn / 2)
670 xmin = xraw[0]
671 xstep = xraw[1] - xmin
672 xmin_new = xmin - offset * xstep
673 print("convolve: the length of the output data has been changed "
674 + "from {} to {}".format(len(yraw), n_new))
675 print(" Add offset = {}".format(offset))
676 print(" xmin changes: {} => {}".format(xmin, xmin_new))
677 x = np.array([xmin_new + i * xstep for i in range(n_new)])
678 return x, yconv
679 return xraw, yconv
680
681def extend_smooth(x, y, nzero, nlin, xstep = 0.0):
682 """
683 概要:
684 データの端をゼロ埋めと線形補間で拡張・平滑化します。
685 詳細説明:
686 データ y の先頭に nzero 個のゼロを追加し、最初の nlin 点を線形に整形して
687 滑らかな接続を実現します。xstepは内部で計算されます。
688 引数:
689 :param x: 元のx軸配列。
690 :type x: numpy.ndarray or list
691 :param y: 元のy軸配列。
692 :type y: numpy.ndarray or list
693 :param nzero: データ先頭に追加するゼロの数。
694 :type nzero: int
695 :param nlin: 線形整形を適用するデータ点の数。
696 :type nlin: int
697 :param xstep: x軸のステップサイズ。デフォルトは0.0ですが、xから計算されます。
698 :type xstep: float
699 戻り値:
700 :returns: 拡張・平滑化されたx軸配列とy軸配列のタプル。
701 :rtype: tuple of (numpy.ndarray, numpy.ndarray)
702 """
703 xmin = x[0]
704 xstep = x[1] - x[0]
705 xmin_new = x[0] - nzero * xstep
706 n_new = nzero + len(x)
707 print("extend_smooth:")
708 print(" Add {} zeros at top of the data".format(nzero))
709 print(" xmin changes: {} => {}".format(xmin, xmin_new))
710 print(" Reshape {} input data with a linear filter".format(nlin))
711
712 xx = np.array([xmin_new + i * xstep for i in range(n_new)])
713 yy = np.zeros(n_new)
714 for i in range(nlin):
715 k = i / (nlin - 1)
716 yy[i+nzero] = k * y[i]
717 for i in range(len(x) - nlin):
718 yy[i+nzero+nlin] = y[i+nzero]
719 return xx, yy
720
721def SmoothingBySimpleAverage(y, n):
722 """
723 概要:
724 移動平均でデータを平滑化します。
725 詳細説明:
726 各データ点について、幅 n の範囲内の平均値で置き換えることでデータを平滑化します。
727 データ端では利用可能な点のみを使用します。
728 引数:
729 :param y: 平滑化するデータ配列。
730 :type y: list
731 :param n: 移動平均の窓幅。
732 :type n: int
733 戻り値:
734 :returns: 平滑化されたデータ配列。
735 :rtype: list
736 """
737 n2 = int(n / 2);
738 ndata = len(y);
739 ys = []
740 for i in range(0, ndata):
741 c = 0;
742 ys.append(0.0);
743 for k in range(i - n2, i + n2 + 1):
744 if k < 0 or k >= ndata:
745 continue
746 ys[i] += y[k]
747 c += 1
748 if c > 0:
749 ys[i] /= c;
750 else:
751 ys[i] = y[i]
752 return ys;
753
754def SmoothingByPolynomialFit(y, n):
755 """
756 概要:
757 多項式フィッティング（Savitzky-Golay風）でデータを平滑化します。
758 詳細説明:
759 各データ点について、幅 n の範囲内で多項式フィッティングを行い、
760 その中心点の値で置き換えることでデータを平滑化します。
761 Savitzky-Golayフィルターに似た重み付けが使用されます。
762 引数:
763 :param y: 平滑化するデータ配列。
764 :type y: list
765 :param n: 平滑化の窓幅。
766 :type n: int
767 戻り値:
768 :returns: 平滑化されたデータ配列。
769 :rtype: list
770 """
771 m = int(n / 2);
772 W23 = (4.0 * m * m - 1.0) * (2.0 * m + 3.0) / 3.0;
773 w23j = [0.0]*(2*m+1)
774 for j in range(-m, m+1):
775 w23j[j + m] = (3.0 * m * (m+1.0) - 1.0 - 5.0 * j * j) / W23
776
777 ndata = len(y)
778 ys = []
779 for i in range(0, ndata):
780 if(i-m < 0 or ndata <= i+m):
781 ys.append(y[i])
782 else:
783 c = 0.0;
784 ys.append(0.0);
785 for j in range(-m, m+1):
786 k = i + j
787 if k < 0 or k >= ndata:
788 continue
789 ys[i] += w23j[j+m] * y[k]
790 c += w23j[j+m]
791
792 if c > 0:
793 ys[i] /= c
794 else:
795 ys[i] = y[i]
796
797 return ys;
798
799
800def deconvolute_fft(xRaw, yRaw, xG, yG):
801 """
802 概要:
803 FFT（高速フーリエ変換）を用いて逆畳み込みを実行します。
804 詳細説明:
805 生データとウィンドウ関数をFFTし、周波数領域で除算を行うことで逆畳み込みを実行します。
806 データ長は2のべき乗に調整されます。
807 引数:
808 :param xRaw: 生データのx軸配列。
809 :type xRaw: list
810 :param yRaw: 生データのy軸配列。
811 :type yRaw: list
812 :param xG: ウィンドウ関数のx軸配列。
813 :type xG: list
814 :param yG: ウィンドウ関数のy軸配列。
815 :type yG: list
816 戻り値:
817 :returns: 逆畳み込み後のx軸配列、y軸配列、FFT処理された生データのx軸配列、y軸配列、FFT処理されたウィンドウ関数のx軸配列、y軸配列のタプル。
818 :rtype: tuple of (list, list, list, list, list, list)
819 """
820 k = sum(yG)
821
822 print("Deconvolution by FFT")
823 n = len(xRaw)
824 nlog = int(log(n) / log(2) + 1.0 - 1.0e-5)
825 nfft = pow(2, nlog)
826 print(" Data number is changed from {} to 2^{} = {} for FFT".format(n, nlog, nfft))
827
828 xmin = xRaw[0]
829 xstep = xRaw[1] - xmin
830 xRawFFT = [xmin + i * xstep for i in range(nfft)]
831 yRawFFT = np.insert(yRaw, len(yRaw), np.zeros(nfft - n))
832# filterの中心位置の原点からのずれによって、iFFT後の原点がずれる
833 yGFFT = np.insert(yG, len(yG), np.zeros(nfft - len(yG)))
834 xminG = xmin + len(xG) / 2 * xstep
835 print(" xmin: ", xmin, xminG)
836 xGFFT = [xminG + i * xstep for i in range(nfft)]
837# nadd = int((nfft - len(yG)) / 2)
838# yGFFT = np.insert(yG, len(yG), np.zeros(nadd))
839# yGFFT = np.insert(yGFFT, 0, np.zeros(nfft - len(yGFFT)))
840
841 yRawFFTed = np.fft.fft(yRawFFT)
842 yGFFTed = np.fft.fft(yGFFT)
843 ycFFTed = yRawFFTed / yGFFTed
844 ydeconv = np.fft.ifft(ycFFTed)
845 ydeconv = [float(ydeconv[i]) for i in range(len(ydeconv))]
846
847 print("")
848
849 return xGFFT, ydeconv, xRawFFT, yRawFFT, xRawFFT, yGFFT
850
851def deconvolute_deconvolve(xRaw, yRaw, xG, yG):
852 """
853 概要:
854 scipy.signal.deconvolve を用いて逆畳み込みを実行します。
855 詳細説明:
856 scipy.signal.deconvolve を直接使用して、観測データ yRaw からウィンドウ関数 yG を逆畳み込みします。
857 引数:
858 :param xRaw: 生データのx軸配列。
859 :type xRaw: list
860 :param yRaw: 生データのy軸配列。
861 :type yRaw: list
862 :param xG: ウィンドウ関数のx軸配列。
863 :type xG: list
864 :param yG: ウィンドウ関数のy軸配列。
865 :type yG: list
866 戻り値:
867 :returns: 逆畳み込み後のx軸配列とy軸配列のタプル。
868 :rtype: tuple of (list, numpy.ndarray)
869 """
870 k = sum(yG)
871
872 print("Deconvolution by scipy.signal.deconvove")
873 print("")
874 IDec, remainder = scipy.signal.deconvolve(yRaw, yG)
875 IDec *= k
876 ndata = len(xRaw)
877 nGhalf = int(len(xG) / 2)
878
879 return xRaw[nGhalf:ndata-nGhalf], IDec
880
881
882def make_edge_correction_weights(xRaw, whalf, nGmax, w_floor=1.0e-6):
883 """
884 概要:
885 Ridge回帰のような手法で使用されるエッジ補正重みを作成します。
886 詳細説明:
887 データ範囲の端で係数が不自然に大きくなるのを防ぐための重みリストを生成します。
888 重みは畳み込み関数の累積和に基づいて計算されます。
889 引数:
890 :param xRaw: 生データのx軸配列。
891 :type xRaw: numpy.ndarray
892 :param whalf: 半値半幅（HWHM）。
893 :type whalf: float
894 :param nGmax: 畳み込み関数の最大インデックス範囲。
895 :type nGmax: int
896 :param w_floor: 重みの最小値。これより小さくならないようにクリップされます。デフォルトは1.0e-6。
897 :type w_floor: float
898 戻り値:
899 :returns: エッジ補正重みのリスト。
900 :rtype: list
901 """
902 n = len(xRaw)
903 w = [0.0]
904 for i in range(1, n):
905 w.append(w[i - 1] + convolve_func(xRaw[i], 0.0, whalf))
906 wGmax = w[n - 1]
907 if abs(wGmax) < w_floor:
908 return [1.0] * n
909 return [max(1.0 - w[i] / wGmax, w_floor) for i in range(n)]
910
911
912def edge_weight_at(w, n, nGmax, j, w_floor=1.0e-6):
913 """
914 概要:
915 係数インデックス j に対する対称的なエッジ補正重みを返します。
916 詳細説明:
917 与えられた重みリスト w とインデックス j に基づいて、対称的なエッジ補正重みを計算します。
918 引数:
919 :param w: エッジ補正重み配列。
920 :type w: list
921 :param n: データ点の総数。
922 :type n: int
923 :param nGmax: 畳み込み関数の最大インデックス範囲。
924 :type nGmax: int
925 :param j: 係数のインデックス。
926 :type j: int
927 :param w_floor: 重みの最小値。これより小さくならないようにクリップされます。デフォルトは1.0e-6。
928 :type w_floor: float
929 戻り値:
930 :returns: 係数インデックス j におけるエッジ補正重み。
931 :rtype: float
932 """
933 if j <= nGmax:
934 wj = w[j]
935 elif n - 1 - j < nGmax:
936 wj = w[n - 1 - j]
937 else:
938 wj = w[nGmax - 1]
939 return max(wj, w_floor)
940
941
942def make_basis_matrix(xRaw, whalf, w=None, nGmax=None):
943 """
944 概要:
945 装置関数の稠密な基底行列を構築します。
946 詳細説明:
947 K[i, j] は係数 j を観測点 i にマッピングします。
948 エッジ重みが提供された場合、各基底列は同じ補正で使用されます。
949 引数:
950 :param xRaw: 生データのx軸配列。
951 :type xRaw: numpy.ndarray or list
952 :param whalf: 半値半幅（HWHM）。
953 :type whalf: float
954 :param w: エッジ補正重み配列。デフォルトはNone。
955 :type w: list or None
956 :param nGmax: 畳み込み関数の最大インデックス範囲。デフォルトはNone。
957 :type nGmax: int or None
958 戻り値:
959 :returns: 装置関数の基底行列。
960 :rtype: numpy.ndarray
961 """
962 x = np.asarray(xRaw, dtype=float)
963 n = len(x)
964 K = np.empty((n, n), dtype=float)
965 for i in range(n):
966 xi = x[i]
967 for j in range(n):
968 v = convolve_func(xi, x[j], whalf)
969 if w is not None and nGmax is not None:
970 v /= edge_weight_at(w, n, nGmax, j)
971 K[i, j] = v
972 return K
973
974
975def make_first_difference_matrix(n):
976 """
977 概要:
978 平滑化ペナルティのための1階差分行列 D を返します。
979 詳細説明:
980 D @ c は [c[1]-c[0], c[2]-c[1], ...] のようなベクトルを生成します。
981 このベクトルにペナルティを課すことで、係数の急激な振動を抑制します。
982 引数:
983 :param n: 係数の数。
984 :type n: int
985 戻り値:
986 :returns: 1階差分行列。
987 :rtype: numpy.ndarray
988 """
989 if n < 2:
990 return np.zeros((0, n), dtype=float)
991 D = np.zeros((n - 1, n), dtype=float)
992 for i in range(n - 1):
993 D[i, i] = -1.0
994 D[i, i + 1] = 1.0
995 return D
996
997
998def deconvolute_nnls(xRaw, yRaw, xG, yG):
999 """
1000 概要:
1001 非負最小二乗法（NNLS）を用いて逆畳み込みを実行します。
1002 詳細説明:
1003 観測データ yRaw と装置関数の基底 K を用いて、係数が非負であるという制約の下で
1004 最小二乗問題を解きます。任意でRidgeまたはSmoothな正則化を適用できます。
1005 引数:
1006 :param xRaw: 生データのx軸配列。
1007 :type xRaw: numpy.ndarray or list
1008 :param yRaw: 生データのy軸配列。
1009 :type yRaw: numpy.ndarray or list
1010 :param xG: ウィンドウ関数のx軸配列。
1011 :type xG: list
1012 :param yG: ウィンドウ関数のy軸配列。
1013 :type yG: list
1014 戻り値:
1015 :returns: 逆畳み込み後のx軸配列とy軸配列のタプル。
1016 :rtype: tuple of (numpy.ndarray, list)
1017 """
1018 global Wa, Grange, alpha, nnls_maxiter, nnls_penalty
1019
1020 dx = abs(xRaw[1] - xRaw[0])
1021
1022 print("Deconvolution by non-negative least squares (NNLS)")
1023 print(" Constraint: all Gaussian/Lorentzian basis coefficients >= 0")
1024 if alpha > 0.0:
1025 print(f" Regularization: {nnls_penalty} penalty with alpha={alpha:g}")
1026 else:
1027 print(" Regularization: none; pure NNLS")
1028 print("")
1029
1030 n = len(xRaw)
1031 nGmax = int(10.0 * Wa / dx + 1.0001)
1032 w = make_edge_correction_weights(xRaw, Wa, nGmax)
1033 K = make_basis_matrix(xRaw, Wa, w=w, nGmax=nGmax)
1034
1035 yvec = np.asarray(yRaw, dtype=float)
1036
1037 if alpha > 0.0:
1038 if nnls_penalty == 'smooth':
1039 R = make_first_difference_matrix(n)
1040 else:
1041 R = np.eye(n, dtype=float)
1042
1043 A = np.vstack([K, sqrt(alpha) * R])
1044 b = np.concatenate([yvec, np.zeros(R.shape[0], dtype=float)])
1045
1046 result = scipy.optimize.lsq_linear(
1047 A,
1048 b,
1049 bounds=(0.0, np.inf),
1050 max_iter=nnls_maxiter,
1051 )
1052 coeff = result.x
1053 rnorm = float(np.linalg.norm(K @ coeff - yvec))
1054 regnorm = float(np.linalg.norm(R @ coeff))
1055 print(f" bounded LS status = {result.status}: {result.message}")
1056 print(f" data residual norm = {rnorm:g}")
1057 print(f" penalty norm = {regnorm:g}")
1058 else:
1059 coeff, rnorm = scipy.optimize.nnls(K, yvec, maxiter=nnls_maxiter)
1060 print(f" NNLS residual norm = {rnorm:g}")
1061
1062 y = coeff / dx
1063
1064 n_active = int(np.count_nonzero(coeff > 0.0))
1065 print(f" active coefficients = {n_active} / {n}")
1066 print(f" coefficient range = {coeff.min():g} - {coeff.max():g}")
1067 print("")
1068
1069 return xRaw, y.tolist()
1070
1071def deconvolute_ridge(xRaw, yRaw, xG, yG):
1072 """
1073 概要:
1074 Ridge回帰を用いて逆畳み込みを実行します。
1075 詳細説明:
1076 観測データ yRaw からRidge回帰により逆畳み込み係数を推定します。
1077 エッジ補正重みと正則化パラメータ alpha を使用します。
1078 引数:
1079 :param xRaw: 生データのx軸配列。
1080 :type xRaw: list
1081 :param yRaw: 生データのy軸配列。
1082 :type yRaw: list
1083 :param xG: ウィンドウ関数のx軸配列。
1084 :type xG: list
1085 :param yG: ウィンドウ関数のy軸配列。
1086 :type yG: list
1087 戻り値:
1088 :returns: 畳み込み後のx軸配列とy軸配列のタプル。
1089 :rtype: tuple of (list, list)
1090 """
1091 global Wa, Grange
1092
1093 k = sum(yG)
1094 dx = abs(xRaw[1] - xRaw[0])
1095
1096 print("Deconvolution by Ridge regression")
1097 print("")
1098
1099 nGmax = int(10.0 * Wa / dx + 1.0001)
1100# print("nGmax=", nGmax)
1101
1102 w = [0.0]
1103 n = len(xRaw)
1104 for i in range(1, len(xRaw)):
1105 w.append(w[i-1] + convolve_func(xRaw[i], 0, Wa))
1106 wGmax = w[n-1]
1107 w_floor = 1.0e-6
1108 w = [max(1.0 - w[i] / wGmax, w_floor) for i in range(n)]
1109# w = [1.0 - w[i] / wGmax for i in range(n)]
1110
1111 y = Ridge(xRaw, yRaw, len(xRaw), lambda x1, x2: convolve_func(x1, x2, Wa), w, nGmax, alpha)
1112 y = [v / dx for v in y]
1113
1114 return xRaw, y
1115
1116def deconvolute_sp(xRaw, yRaw, xG, yG, nsmooth):
1117 """
1118 概要:
1119 Smoothing Penalty法を用いて逆畳み込みを実行します。
1120 詳細説明:
1121 観測データ yRaw からSmoothing Penalty回帰により逆畳み込み係数を推定します。
1122 エッジ補正重み、正則化パラメータ alpha、および平滑化幅 nsmooth を使用します。
1123 引数:
1124 :param xRaw: 生データのx軸配列。
1125 :type xRaw: list
1126 :param yRaw: 生データのy軸配列。
1127 :type yRaw: list
1128 :param xG: ウィンドウ関数のx軸配列。
1129 :type xG: list
1130 :param yG: ウィンドウ関数のy軸配列。
1131 :type yG: list
1132 :param nsmooth: 平滑化の幅。
1133 :type nsmooth: int
1134 戻り値:
1135 :returns: 逆畳み込み後のx軸配列とy軸配列のタプル。
1136 :rtype: tuple of (list, list)
1137 """
1138 global Wa, Grange
1139
1140 k = sum(yG)
1141 dx = abs(xRaw[1] - xRaw[0])
1142
1143 print("Deconvolution by smoothing penalty method")
1144 print("")
1145
1146 nGmax = int(10.0 * Wa / dx + 1.0001)
1147# print("nGmax=", nGmax)
1148
1149 w_floor = 1.0e-6
1150 w = [w_floor]
1151# w = [0.0]
1152 n = len(xRaw)
1153 for i in range(1, len(xRaw)):
1154 w.append(w[i-1] + convolve_func(xRaw[i], 0, Wa))
1155 wGmax = w[n-1]
1156 w = [max(1.0 - w[i] / wGmax, w_floor) for i in range(n)]
1157# w = [1.0 - w[i] / wGmax for i in range(n)]
1158
1159 y = Smoothing_Penalty_Regression(xRaw, yRaw, len(xRaw),
1160 lambda x1, x2: convolve_func(x1, x2, Wa), w, nGmax, alpha, nsmooth)
1161 y = [v / dx for v in y]
1162
1163 return xRaw, y
1164
1165#window = None
1166def _set_iteration_ylim(ax0, yRaw, y):
1167 """
1168 概要:
1169 反復プロットのy軸範囲を安定して設定します。
1170 詳細説明:
1171 生データ yRaw と更新されたデータ y に基づいて、Matplotlibの軸 ax0 のy軸範囲を
1172 動的に調整します。これにより、反復処理中のグラフ表示が安定します。
1173 引数:
1174 :param ax0: Matplotlibの軸オブジェクト。
1175 :type ax0: matplotlib.axes.Axes
1176 :param yRaw: 生データのy軸配列。
1177 :type yRaw: list or numpy.ndarray
1178 :param y: 更新されたデータ（反復結果）のy軸配列。
1179 :type y: list or numpy.ndarray
1180 """
1181 ymin = min([0.0, min(yRaw), min(y)])
1182 ymax = max([0.0, max(yRaw), max(y)])
1183 if ymax <= ymin:
1184 ymax = ymin + 1.0
1185 margin = 0.05 * (ymax - ymin)
1186 ax0.set_ylim([ymin - margin, ymax + margin])
1187
1188
1189def _init_iteration_plot(fig, ax0, xRaw, yRaw, y, xgmin, xgmax):
1190 """
1191 概要:
1192 反復更新グラフを一度だけ初期化します。
1193 詳細説明:
1194 反復処理の開始時にグラフをセットアップし、rawデータと初期の更新データをプロットします。
1195 Line2Dオブジェクトを再利用することで、反復ごとの描画コストを削減します。
1196 引数:
1197 :param fig: MatplotlibのFigureオブジェクト。
1198 :type fig: matplotlib.pyplot.Figure
1199 :param ax0: Matplotlibの軸オブジェクト。
1200 :type ax0: matplotlib.axes.Axes
1201 :param xRaw: 生データのx軸配列。
1202 :type xRaw: list or numpy.ndarray
1203 :param yRaw: 生データのy軸配列。
1204 :type yRaw: list or numpy.ndarray
1205 :param y: 現在の更新データ（反復結果）のy軸配列。
1206 :type y: list or numpy.ndarray
1207 :param xgmin: x軸の最小表示範囲。
1208 :type xgmin: float
1209 :param xgmax: x軸の最大表示範囲。
1210 :type xgmax: float
1211 戻り値:
1212 :returns: 更新されるプロットのLine2Dオブジェクト。
1213 :rtype: matplotlib.lines.Line2D
1214 """
1215 ax0.cla()
1216 ax0.plot(xRaw, yRaw, label = 'raw/initial')
1217 line_updated, = ax0.plot(xRaw, y, label = 'updated')
1218 ax0.set_xlim([xgmin, xgmax])
1219 _set_iteration_ylim(ax0, yRaw, y)
1220 ax0.legend()
1221
1222 # Layout adjustment is expensive; do it once before the loop.
1223 app.plotevent.layout()
1224 fig.canvas.draw_idle()
1225 fig.canvas.flush_events()
1226 plt.pause(max(sleeptime, 0.05))
1227 return line_updated
1228
1229
1230def _update_iteration_plot(fig, ax0, line_updated, yRaw, y):
1231 """
1232 概要:
1233 反復逆畳み込み中にyデータのみを更新します。
1234 詳細説明:
1235 既存のLine2Dオブジェクトのyデータを更新し、グラフのy軸範囲を調整して再描画します。
1236 これにより、高速な更新が可能になります。
1237 引数:
1238 :param fig: MatplotlibのFigureオブジェクト。
1239 :type fig: matplotlib.pyplot.Figure
1240 :param ax0: Matplotlibの軸オブジェクト。
1241 :type ax0: matplotlib.axes.Axes
1242 :param line_updated: 更新されるプロットのLine2Dオブジェクト。
1243 :type line_updated: matplotlib.lines.Line2D
1244 :param yRaw: 生データのy軸配列。y軸範囲調整のために使用されます。
1245 :type yRaw: list or numpy.ndarray
1246 :param y: 更新されたデータ（反復結果）のy軸配列。
1247 :type y: list or numpy.ndarray
1248 """
1249 line_updated.set_ydata(y)
1250 _set_iteration_ylim(ax0, yRaw, y)
1251 fig.canvas.draw_idle()
1252 fig.canvas.flush_events()
1253 plt.pause(max(sleeptime, 0.05))
1254
1255
1256def deconvolute_jacobi(xRaw, yRaw, xG, yG, fig, ax):
1257 """
1258 概要:
1259 Jacobi法を用いて逆畳み込みを実行します。
1260 詳細説明:
1261 反復計算により、畳み込み方程式をJacobi法で解き、データを逆畳み込みします。
1262 途中で平滑化や負値クリッピングを行い、グラフをリアルタイムで更新します。
1263 収束条件または最大反復回数に達するまで処理を続行します。
1264 引数:
1265 :param xRaw: 生データのx軸配列。
1266 :type xRaw: list
1267 :param yRaw: 生データのy軸配列。
1268 :type yRaw: list
1269 :param xG: ウィンドウ関数のx軸配列。
1270 :type xG: list
1271 :param yG: ウィンドウ関数のy軸配列。
1272 :type yG: list
1273 :param fig: MatplotlibのFigureオブジェクト。
1274 :type fig: matplotlib.pyplot.Figure
1275 :param ax: Matplotlibの軸オブジェクトのリスト。最初の要素が使用されます。
1276 :type ax: list of matplotlib.axes.Axes
1277 戻り値:
1278 :returns: 逆畳み込み後のx軸配列とy軸配列のタプル。
1279 :rtype: tuple of (list, list)
1280 """
1281 global Wa, Grange
1282# global window
1283
1284 k = sum(yG)
1285
1286 print("Deconvolution by Jacobi method")
1287 print("")
1288
1289 xstep = xRaw[1] - xRaw[0]
1290
1291 xgmin = min(xRaw)
1292 xgmax = max(xRaw)
1293
1294 n = len(xRaw)
1295 Sg = np.zeros(n)
1296 for i in range(n):
1297 for j in range(n):
1298 Sg[i] += Hij(xstep, Wa, Grange, i, j)
1299 print("Filter area w.r.t. i: Sg=", Sg[int(n/2)])
1300
1301 ymax = max([abs(yRaw[i]) for i in range(n)])
1302 y = yRaw.copy()
1303 yPrev = yRaw.copy()
1304
1305 line_updated = _init_iteration_plot(fig, ax[0], xRaw, yRaw, y, xgmin, xgmax)
1306
1307 for it in range(nmaxiter):
1308 Hx = np.zeros(n)
1309 print("iter=", it)
1310
1311 for i in range(n):
1312 for j in range(n):
1313 Hx[i] += Hij(xstep, Wa, Grange, i, j) * yPrev[j]
1314 h = Hij(xstep, Wa, Grange, i, i)
1315 y[i] = yPrev[i] + (yRaw[i] - Hx[i] / Sg[i]) / h * dump
1316
1317# y = SmoothingBySimpleAverage(y, nsmooth)
1318 y = SmoothingByPolynomialFit(y, nsmooth)
1319 if zero_correction:
1320 for i in range(n):
1321 if y[i] < 0.0:
1322 y[i] = 0.0
1323
1324# w = plt.get_current_fig_manager().window
1325# if w != window:
1326 if app.plotevent.stop_flag:
1327 print("Terminated by user")
1328 break
1329
1330 _update_iteration_plot(fig, ax[0], line_updated, yRaw, y)
1331
1332 max_err = max([abs(y[i] - yPrev[i]) for i in range(n)])
1333 rel_err = max_err / ymax
1334 print(" max error: ", max_err, " relative error: ", rel_err, " eps=", eps)
1335 if max_err / ymax < eps:
1336 print("Converged at max_err={} ({} relative) < {}".format(max_err, rel_err, eps))
1337 break
1338
1339 yPrev = y.copy()
1340 else:
1341 print("Not converged")
1342
1343 return xRaw, y
1344
1345
1346def deconvolute_gauss_seidel(xRaw, yRaw, xG, yG, fig, ax):
1347 """
1348 概要:
1349 Gauss-Seidel法を用いて逆畳み込みを実行します。
1350 詳細説明:
1351 反復計算により、畳み込み方程式をGauss-Seidel法で解き、データを逆畳み込みします。
1352 途中で平滑化や負値クリッピングを行い、グラフをリアルタイムで更新します。
1353 収束条件または最大反復回数に達するまで処理を続行します。
1354 引数:
1355 :param xRaw: 生データのx軸配列。
1356 :type xRaw: list
1357 :param yRaw: 生データのy軸配列。
1358 :type yRaw: list
1359 :param xG: ウィンドウ関数のx軸配列。
1360 :type xG: list
1361 :param yG: ウィンドウ関数のy軸配列。
1362 :type yG: list
1363 :param fig: MatplotlibのFigureオブジェクト。
1364 :type fig: matplotlib.pyplot.Figure
1365 :param ax: Matplotlibの軸オブジェクトのリスト。最初の要素が使用されます。
1366 :type ax: list of matplotlib.axes.Axes
1367 戻り値:
1368 :returns: 逆畳み込み後のx軸配列とy軸配列のタプル。
1369 :rtype: tuple of (list, list)
1370 """
1371 global Wa, Grange
1372# global window
1373
1374 k = sum(yG)
1375
1376 print("Deconvolution by Gauss-Seidel method")
1377 print("")
1378
1379 xstep = xRaw[1] - xRaw[0]
1380
1381 xgmin = min(xRaw)
1382 xgmax = max(xRaw)
1383
1384 n = len(xRaw)
1385 Sg = np.zeros(n)
1386 for i in range(n):
1387 for j in range(n):
1388 Sg[i] += Hij(xstep, Wa, Grange, i, j)
1389 print("Filter area w.r.t. i: Sg=", Sg[int(n/2)])
1390
1391 ymax = max([abs(yRaw[i]) for i in range(n)])
1392 y = yRaw.copy()
1393 yPrev = yRaw.copy()
1394
1395 line_updated = _init_iteration_plot(fig, ax[0], xRaw, yRaw, y, xgmin, xgmax)
1396
1397 for it in range(nmaxiter):
1398 Hx = np.zeros(n)
1399 print(f"iter={it:4} ", end = '')
1400
1401 for i in range(n):
1402 for j in range(i):
1403 Hx[i] += Hij(xstep, Wa, Grange, i, j) * y[j]
1404 for j in range(i, n):
1405 Hx[i] += Hij(xstep, Wa, Grange, i, j) * yPrev[j]
1406 h = Hij(xstep, Wa, Grange, i, i)
1407 y[i] = yPrev[i] + (yRaw[i] - Hx[i] / Sg[i]) / h * dump
1408
1409# y = SmoothingBySimpleAverage(y, nsmooth)
1410 y = SmoothingByPolynomialFit(y, nsmooth)
1411 if zero_correction:
1412 for i in range(n):
1413 if y[i] < 0.0:
1414 y[i] = 0.0
1415
1416# w = plt.get_current_fig_manager().window
1417# if w != window:
1418 if app.plotevent.stop_flag:
1419 print("Terminated by user")
1420 break
1421
1422 _update_iteration_plot(fig, ax[0], line_updated, yRaw, y)
1423
1424 max_err = max([abs(y[i] - yPrev[i]) for i in range(n)])
1425 rel_err = max_err / ymax
1426 print(f" max error: {max_err:10.4g} relative error: {rel_err:10.4g} eps={eps:10.4g}")
1427 if max_err / ymax < eps:
1428 print("Converged at max_err={} ({} relative) < {}".format(max_err, rel_err, eps))
1429 break
1430
1431 yPrev = y.copy()
1432 else:
1433 print("Not converged")
1434
1435 return xRaw, y
1436
1437
1438
1439
1440def plot_input_data(xRaw, yRaw, label_x='', label_y=''):
1441 """
1442 概要:
1443 選択された入力データのみをプロットします。
1444 詳細説明:
1445 逆畳み込みを実行せず、指定された入力データをグラフに表示するだけのモードです。
1446 逆畳み込みの結果はExcelファイルに出力されません。
1447 引数:
1448 :param xRaw: 入力データのx軸配列。
1449 :type xRaw: list
1450 :param yRaw: 入力データのy軸配列。
1451 :type yRaw: list
1452 :param label_x: x軸のラベル文字列。
1453 :type label_x: str
1454 :param label_y: y軸のラベル文字列。
1455 :type label_y: str
1456 """
1457 print("Plot input data only")
1458 print(" No deconvolution is performed.")
1459 print(" No deconvoluted Excel output is written.")
1460 print("")
1461
1462 fig = plt.figure(figsize=figsize)
1463 ax0 = fig.add_subplot(1, 1, 1)
1464
1465 app.plotevent = tkPlotEvent(plt)
1466 app.plotevent.add_button()
1467
1468 ax0.plot(xRaw, yRaw, label='input')
1469 ax0.set_xlabel(str(label_x))
1470 ax0.set_ylabel(str(label_y))
1471 ax0.set_xlim([min(xRaw), max(xRaw)])
1472
1473 ymin = min(yRaw)
1474 ymax = max(yRaw)
1475 if ymax <= ymin:
1476 ymax = ymin + 1.0
1477 margin = 0.05 * (ymax - ymin)
1478 ax0.set_ylim([ymin - margin, ymax + margin])
1479 ax0.legend()
1480
1481 app.plotevent.layout()
1482 fig.canvas.draw_idle()
1483 fig.canvas.flush_events()
1484 plt.pause(0.1)
1485
1486 app.plotevent.finalize_button()
1487 app.terminate("Press ENTER to exit>>", usage=lambda app: usage(app), pause=True)
1488
1489
1490#======================
1491# main
1492#======================
1493def main():
1494 """
1495 概要:
1496 メイン処理を実行します。
1497 詳細説明:
1498 コマンドライン引数を解析し、入力データを読み込み、
1499 選択されたモードに基づいて逆畳み込み処理を実行します。
1500 結果はプロットされ、Excelファイルに保存されます。
1501 """
1502 global xmin, xmax
1503
1504 args = build_arg_parser().parse_args()
1505 apply_args_to_globals(args)
1506
1507 logfile = app.replace_path(infile)
1508 outfile = app.replace_path(infile, template = ["{dirname}", "{filebody}-deconvoluted.xlsx"])
1509 print(f"Open logfile [{logfile}]")
1510 app.redirect(targets = ["stdout", logfile], mode = 'w')
1511
1512 print("")
1513 print(f"Process data in [{infile}] with {func_type} function of width={Wa}")
1514 print("infile : ", infile)
1515 print("xlabel : ", xlabel)
1516 print("ylabel : ", ylabel)
1517 print("outfile : ", outfile)
1518 print("mode : ", mode)
1519 print("flip_x : ", flip_x)
1520 if mode == 'convolve' or mode == 'fft':
1521 print("For mode = 'convolve' or 'fft':")
1522 print(" convmode: ", convmode)
1523 print(" x range : ", xmin, xmax)
1524 if mode == 'gs' or mode == 'jacobi':
1525 print(" x range : ", xmin, xmax)
1526 print("For mode = 'gs' or 'jacobi':")
1527 print(" dump=", dump)
1528 print(" nmaxiter=", nmaxiter)
1529 print(" eps=", eps)
1530 print(" nsmooth=", nsmooth)
1531 print(" zero correction=", zero_correction)
1532 if mode in ('ridge', 'sp'):
1533 print(" alpha=", alpha)
1534 if mode == 'nnls':
1535 print("For mode = 'nnls':")
1536 print(" non-negative coefficients: enabled")
1537 print(" alpha=", alpha)
1538 if alpha > 0.0:
1539 print(" regularization penalty=", nnls_penalty)
1540 else:
1541 print(" regularization penalty= none (pure NNLS)")
1542 print(" nnls_maxiter=", nnls_maxiter)
1543 if mode == 'plot':
1544 print("For mode = 'plot':")
1545 print(" plot selected input data only")
1546 print("")
1547
1548 if isinstance(xlabel, str) and '*** choose ***' in xlabel:
1549 app.terminate(f"\nError: Choose x label\nTerminated.\n", usage = lambda app: usage(app), pause = True)
1550 if isinstance(ylabel, str) and '*** choose ***' in ylabel:
1551 app.terminate(f"\nError: Choose y label\nTerminated.\n", usage = lambda app: usage(app), pause = True)
1552
1553 print("Input data")
1554 print("X range to be calculated:", xmin, "-", xmax)
1555 datafile = tkVariousData(infile)
1556 header, datalist = datafile.Read_minimum_matrix(close_fp = True, usage = lambda app: usage(app))
1557 label_x, _xin = datafile.FindDataArray(xlabel, flag = 'i')
1558 label_y, _yin = datafile.FindDataArray(ylabel, flag = 'i')
1559 if _xin is None:
1560 app.terminate(f"Error: Can not find the data with [{xlabel}]", usage = lambda app: usage(app), pause = True)
1561 if _yin is None:
1562 app.terminate(f"Error: Can not find the data with [{xlabel}]", usage = lambda app: usage(app), pause = True)
1563
1564 xin_all = []
1565 yin_all = []
1566 for i in range(len(_xin)):
1567 xval = -_xin[i] if flip_x else _xin[i]
1568 xin_all.append(xval)
1569 yin_all.append(_yin[i])
1570
1571 if flip_x:
1572 print(" flip_x=1: multiply input x values by -1 before filtering/calculation.")
1573
1574 xin = []
1575 yin = []
1576 for xval, yval in zip(xin_all, yin_all):
1577 if not (xmin is None or xmin == '' or xmin == '*') and xval < xmin:
1578 continue
1579 if not (xmax is None or xmax == '' or xmax == '*') and xmax < xval:
1580 continue
1581 xin.append(xval)
1582 yin.append(yval)
1583# print("x=", xin)
1584# print("y=", yin)
1585
1586 if len(xin) < 2:
1587 app.terminate(
1588 f"\nError in main(): Need at least two data points after x-range filtering.\n",
1589 usage = lambda app: usage(app), pause = True,
1590 )
1591
1592 # Internally use ascending x order. This avoids negative xstep problems in
1593 # sampled filters and keeps coefficient-to-spectrum scaling positive.
1594 # Descending binding-energy axes are common in photoemission data.
1595 if xin[1] < xin[0]:
1596 print(" Input x axis is descending; reverse data internally for calculation.")
1597 xin = list(reversed(xin))
1598 yin = list(reversed(yin))
1599 elif any(xin[i + 1] < xin[i] for i in range(len(xin) - 1)):
1600 print(" Input x axis is not monotonic; sort data by x for calculation.")
1601 xy = sorted(zip(xin, yin), key=lambda t: t[0])
1602 xin = [v[0] for v in xy]
1603 yin = [v[1] for v in xy]
1604
1605 nindata = len(xin)
1606 print(" ndata = ", nindata)
1607 if xmin is None or xmin == '' or xmin == '*':
1608 xmin = min(xin)
1609 if xmax is None or xmax == '' or xmax == '*':
1610 xmax = max(xin)
1611 print("x range=", xmin, xmax)
1612
1613 xinmin = xin[0]
1614 xinstep = xin[1] - xin[0]
1615 if xinstep == 0.0:
1616 app.terminate(f"\nError in main(): xin[0] and xin[1] are the same. Check input file [{infile}].\n",
1617 usage = lambda app: usage(app), pause = True)
1618
1619 xinmax = xinmin + xinstep * (nindata - 1)
1620 print(" x range: {} - {} at {} step".format(xinmin, xinmax, xinstep))
1621 print("Smooth mode: ", smoothmode)
1622 print("")
1623
1624 xRaw = xin
1625 yRaw = yin
1626 xstep = xinstep
1627
1628 if mode == 'plot':
1629 plot_input_data(xRaw, yRaw, label_x=label_x, label_y=label_y)
1630 return
1631
1632 print("Filter: Wa = {} Grange = {}".format(Wa, Grange))
1633 xG, yG = make_wf(Wa, Grange, xstep)
1634 print("")
1635
1636 fig = plt.figure(figsize = figsize)
1637 ax = []
1638 ax.append(fig.add_subplot(3, 1, 1))
1639 ax.append(fig.add_subplot(3, 1, 2))
1640 ax.append(fig.add_subplot(3, 1, 3))
1641# ax.append(fig.add_subplot(4, 1, 4))
1642
1643 app.plotevent = tkPlotEvent(plt)
1644 app.plotevent.add_button()
1645 app.plotevent.layout()
1646
1647 fig.canvas.draw_idle()
1648 fig.canvas.flush_events()
1649 plt.pause(0.1)
1650# window = plt.get_current_fig_manager().window
1651
1652# make data to be deconvolved
1653 nw = int(len(xG) / 2)
1654 if 'average' in smoothmode:
1655 yRaw = SmoothingBySimpleAverage(yRaw, 31)
1656 if 'extend' in smoothmode:
1657 xRaw, yRaw = extend_smooth(xRaw, yRaw, nw*kzero, nw*klin, xstep)
1658 if 'convolve' in smoothmode:
1659 xconv, yconv = convolve(xRaw, yRaw, yG, mode = convmode)
1660 else:
1661 xconv, yconv = xRaw, yRaw
1662
1663# deconvolution
1664 if mode == 'fft':
1665 xDec, yDec, xRawFFT, yRawFFT, xGFFT, yGFFT = deconvolute_fft(xconv, yconv, xG, yG)
1666 xconv = xRawFFT
1667 yconv = yRawFFT
1668 datafft = ax[2].plot(xGFFT, yGFFT, label = 'WF for FFT')
1669 elif mode == 'jacobi':
1670 xDec, yDec = deconvolute_jacobi(xconv, yconv, xG, yG, fig, ax)
1671 yG = [convolve_func(xRaw[i], 0.0, Wa) for i in range(len(xRaw))]
1672 datafft = ax[2].plot(xRaw, yG, label = 'WF for Jacobi method')
1673 elif mode == 'gs':
1674 xDec, yDec = deconvolute_gauss_seidel(xconv, yconv, xG, yG, fig, ax)
1675 yG = [convolve_func(xRaw[i], 0.0, Wa) for i in range(len(xRaw))]
1676 datafft = ax[2].plot(xRaw, yG, label = 'WF for Gauss-Seidel method')
1677 elif mode == 'sp':
1678 yconv = SmoothingByPolynomialFit(yconv, nsmooth)
1679 xDec, yDec = deconvolute_sp(xconv, yconv, xG, yG, nsmooth)
1680# yconv = SmoothingByPolynomialFit(yconv, nsmooth)
1681 datafft = ax[2].plot(xG, yG, label = 'WF for Gauss kernel regression')
1682 elif mode == 'ridge':
1683 yconv = SmoothingByPolynomialFit(yconv, nsmooth)
1684 xDec, yDec = deconvolute_ridge(xconv, yconv, xG, yG)
1685 datafft = ax[2].plot(xG, yG, label = 'WF for Ridge regression')
1686 elif mode == 'nnls':
1687 yconv = SmoothingByPolynomialFit(yconv, nsmooth)
1688 xDec, yDec = deconvolute_nnls(xconv, yconv, xG, yG)
1689 datafft = ax[2].plot(xG, yG, label = 'WF for NNLS')
1690 elif mode == 'convolve' or mode == 'deconvolve':
1691 xDec, yDec = deconvolute_deconvolve(xconv, yconv, xG, yG)
1692 datafft = ax[2].plot(xG, yG, label = 'WF for scipy.signal.deconvolve')
1693 else:
1694 raise ValueError(f"Unsupported mode: {mode}")
1695
1696 app.plotevent.finalize_button()
1697 if app.plotevent.stop_flag:
1698 print("")
1699 print("Iteration terminated by user")
1700 else:
1701 xgmin = min([min(xconv), min(xDec)])
1702 xgmax = max([max(xconv), max(xDec)])
1703
1704 ax[0].cla()
1705 data1 = ax[0].plot(xconv, yconv, label = 'input(convoluted)')
1706 data3 = ax[0].plot(xDec, yDec, label = 'deconvoluted')
1707 data1 = ax[1].plot(xRaw, yRaw, label = 'input(raw)')
1708 data2 = ax[1].plot(xconv, yconv, label = 'input(convoluted)')
1709# data5 = ax[2].plot(xDec, yDec, label = 'deconvoluted')
1710# data6 = ax[2].plot([xgmin, xgmax], [0.0, 0.0], color = 'red', linestyle = 'dashed')
1711 ax[0].set_xlim([xgmin, xgmax])
1712 ax[0].set_ylim([0.0, max([max(yRaw), max(yDec)])])
1713 ax[1].set_xlim([xgmin, xgmax])
1714 ax[1].set_ylim([0.0, max(yRaw)])
1715# ax[2].set_xlim([xgmin, xgmax])
1716# ax[2].set_ylim([-100, 10000])
1717 ax[2].set_xlim([xgmin, xgmax])
1718 ax[2].set_ylim([0.0, max(yG) if len(yG) > 0 else 1.0])
1719# ax[1].set_ylim([0.0, max(yRaw)])
1720
1721 ax[0].legend()
1722 ax[1].legend()
1723 ax[2].legend()
1724# ax[3].legend()
1725
1726# plt.tight_layout()
1727 app.plotevent.layout()
1728
1729 fig.canvas.draw_idle()
1730 fig.canvas.flush_events()
1731 plt.pause(0.1)
1732
1733 print("")
1734 print("Save deconvoluted data to [{}]".format(outfile))
1735 df = pd.DataFrame(np.array([xconv, yconv, yDec]).T, columns = ['x', 'y(input)', 'y(deconvoluted)'])
1736# df.to_excel(outfile, index = False)
1737 data_list = df.values.tolist()
1738 tkVariousData().to_excel(outfile, list(df.columns), list(map(list, zip(*data_list))), template = template)
1739
1740 app.terminate("Press ENTER to exit>>", usage = lambda app: usage(app), pause = True)
1741
1742
1743if __name__ == '__main__':
1744 main()