"""
X線回折（XRD）データ解析スクリプト
概要:
このスクリプトは、実験的に測定されたX線回折（XRD）パターンと、結晶構造情報ファイル（CIF）から計算された理論的なXRDパターンを比較・分析するためのツールです。
主な機能として、パターンプロット、相関分析、オーバーラップチェック、およびLASSO回帰による混合相の定量フィッティングを提供します。
詳細説明:
本スクリプトは、様々なモードでXRDデータ解析を実行します。
- `plot_all`: 実験XRDパターンと各CIF由来の理論XRDパターンを個別のグラフに表示します。
- `plot_one`: 実験XRDパターンと各CIF由来の回折ピーク位置を一つのグラフに表示します。
- `plot_one2`: 実験XRDパターンと各CIF由来の理論XRDパターンを一つのグラフに表示します。
- `overwrap`: 異なるCIFパターン間のオーバーラップ（重なり）を評価します。
- `CIFcorrelation`: CIFパターン間の相関係数を評価します。
- `correlation`: 実験パターンとCIFパターン間の相関係数を評価します。
- `fit`: LASSO回帰を用いて実験パターンを背景とCIFパターンでフィッティングし、各相の寄与を定量化します。
スマアリング（ブロードニング）効果や背景処理もサポートしています。
関連リンク:
:doc:`xrd_fit_usage`
"""
import os
import sys
import glob
import numpy as np
from numpy import exp, log, sin, cos, tan, arcsin, arccos, arctan, pi
from scipy.special import legendre
from scipy.interpolate import interp1d
from matplotlib import pyplot as plt
import pandas as pd
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import Lasso
from sklearn.metrics import mean_absolute_error, mean_squared_error
from tklib.tkapplication import tkApplication
from tklib.tkfile import tkFile
from tklib.tkutils import getarg, getintarg, getfloatarg, pint, pfloat, split_file_path, replace_path
from tklib.tkvariousdata import tkVariousData
from tklib.tksci.tksci import Gaussian, Lorentzian, GaussLorentz
from tklib.tksci.tkmatrix import make_matrix1, make_matrix2, make_matrix3
from tklib.tkcrystal.tkxrd import Xray_wavelengths
from tklib.tkgraphic.tkplotevent import tkPlotEvent
# Xray_wavelengths
# URL: https://github.com/materialsproject/pymatgen/blob/v2023.5.10/pymatgen/analysis/diffraction/xrd.py
# 'CuKa', 'CuKa2', 'CuKa1', 'CuKb1', 'MoKa', 'MoKa2', 'MoKa1', 'MoKb1',
# 'CrKa', 'CrKa2', 'CrKa1', 'CrKb1', 'FeKa', 'FeKa2', 'FeKa1', 'FeKb1', 'CoKa', 'CoKa2', 'CoKa1', 'CoKb1',
# 'AgKa', 'AgKa2', 'AgKa1', 'AgKb1'
usage_str = '''
" (i) usage: python {} mode input_path cif_dir wavelength Q2min Q2max Q2step fwhm yscale BG_order LASSO_alpha".format(sys.argv[0])
" ex: python {} plot_all {} {} {} {} {} {} {} {}".format(sys.argv[0], cparams.infile, cparams.wavelength, cparams.xmin, cparams.xmax, cparams.xstep, cparams.fwhm, cparams.yscale, cparams.BGorder, cparams.alpha)
" ex: python {} plot_one {} {} {} {} {} {} {} {}".format(sys.argv[0], cparams.infile, cparams.wavelength, cparams.xmin, cparams.xmax, cparams.xstep, cparams.fwhm, cparams.yscale, cparams.BGorder, cparams.alpha)
" ex: python {} plot_one2 {} {} {} {} {} {} {} {} {}".format(sys.argv[0], cparams.infile, cparams.wavelength, cparams.xmin, cparams.xmax, cparams.xstep, cparams.fwhm, cparams.yscale, cparams.BGorder, cparams.alpha)
" ex: python {} fit {} {} {} {} {} {} {} {} {}".format(sys.argv[0], cparams.infile, cparams.wavelength, cparams.xmin, cparams.xmax, cparams.xstep, cparams.fwhm, cparams.yscale, cparams.BGorder, cparams.alpha)
'''[1:-1]
#================================
# global parameters
#================================
module_names = []
modules = []
markers = ['o', 's', '+', 'x', 'D', 'v', '^', '<', '>', '8', 'h', 'H']
colors = ['black', 'red', 'blue', 'darkgreen', 'darkorange', 'hotpink', 'lightgreen', 'cyan', 'yellow', 'magenta', 'chocolate',
'navy', 'slategray', 'olive' ]
figsize = (12, 8)
figsize_one = (10, 6)
figsize_coeff = (8, 6)
fontsize = 12
legend_fontsize = 8
legend_fontsize_one = 12
#=============================
# Treat argments
#=============================
[ドキュメント]
def usage(app = None, cparams = None):
"""
コマンドライン引数の使い方を表示する。
詳細説明:
アプリケーションとパラメータオブジェクトから利用方法文字列を取得し、整形して標準出力に出力する。
:param app: tkApplicationのインスタンス。
:param cparams: 設定パラメータを保持するオブジェクト。
:returns: なし
"""
cparams = app.get_params()
# app.usage(infile = cparams.infile)
for s in app.usage_str.split('\n'):
cmd = 'print({})'.format(s.rstrip())
eval(cmd)
[ドキュメント]
def initialize(app = None, cparams = None):
"""
アプリケーションのパラメータを初期値で設定する。
詳細説明:
デバッグフラグ、プラグインディレクトリ、モード、入力ファイルパス、CIFファイルパス、X線波長、
2θ範囲、FWHM、ガウス比率、Yスケール、背景次数、LASSOのalphaなどの初期値を設定する。
:param app: tkApplicationのインスタンス。
:param cparams: 設定パラメータを保持するオブジェクト。
:returns: なし
"""
cparams.debug = 0
cparams.findvalidstructure = True
cparams.plugin_dir = 'filter'
cparams.mode = 'plot'
cparams.infile = '*.txt'
cparams.cif_files = 'data/*.*'
cparams.beam = 'X-ray'
cparams.wavelength = "CuKa"
cparams.xmin = 20.0
cparams.xmax = 120.0
cparams.xstep = 0.02
cparams.fwhm = 0.2
cparams.Gfraction = 0.5
cparams.fwhm_smear = 0.0
cparams.Gfraction_smear = 0.0
cparams.yscale = 'linear'
cparams.BGorder = 3
cparams.alpha = 0.1
[ドキュメント]
def update_vars(app = None, cparams = None):
"""
コマンドライン引数に基づいて設定パラメータを更新する。
詳細説明:
`sys.argv` から引数を取得し、`cparams` オブジェクトの対応する属性に設定する。
数値型引数は適切な型に変換される。
:param app: tkApplicationのインスタンス。
:param cparams: 設定パラメータを保持するオブジェクト。
:returns: なし
"""
# if getarg(2, None) is None:
# app.terminate(usage = usage)
cparams.mode = getarg(1, cparams.mode)
cparams.plugin_dir = getarg(2, cparams.plugin_dir)
cparams.infile = getarg(3, cparams.infile)
cparams.cif_files = getarg(4, cparams.cif_files)
cparams.wavelength = getarg(5, cparams.wavelength)
cparams.xmin = getfloatarg(6, cparams.xmin)
cparams.xmax = getfloatarg(7, cparams.xmax)
cparams.xstep = getfloatarg(8, cparams.xstep)
cparams.fwhm = getfloatarg(9, cparams.fwhm)
cparams.Gfraction = getfloatarg(10, cparams.Gfraction)
cparams.fwhm_smear = getfloatarg(11, cparams.fwhm_smear)
cparams.Gfraction_smear = getfloatarg(12, cparams.Gfraction_smear)
cparams.yscale = getarg (13, cparams.yscale)
cparams.BGorder = getintarg (14, cparams.BGorder)
cparams.alpha = getfloatarg(15, cparams.alpha)
[ドキュメント]
def read_file(path, app, cparams):
"""
指定されたパスのファイルを読み込み、その内容を処理する。
詳細説明:
登録されているモジュールを順に試行し、ファイルのタイプを判別する。
ファイルタイプが一致したモジュールを使用してデータを読み込み、必要に応じて変換する。
:param path: str. 読み込むファイルのパス。
:param app: tkApplicationのインスタンス。
:param cparams: 設定パラメータを保持するオブジェクト。
:returns: tuple. (`module`, `inf`). `module` はファイルを読み込んだモジュールオブジェクト、
`inf` は読み込まれたデータを含む辞書。ファイルが読み込めない場合は `(None, None)`.
"""
module = None
for i in range(len(modules)):
name = module_names[i]
m = modules[i]
file_type = m.check_file_type(path, app = app, cparams = cparams)
# file_type = app.call(m, "check_file_type", path, app = app, cparams = cparams)
print(f"try [{name}] for [{path}]: file_type={file_type}")
if file_type is not None and 'Error' not in file_type:
print(" type matched.")
module = m
break
if module is None:
return None, None
# inf = app.call(module, "read_data", path, app = app, cparams = cparams)
inf = module.read_data(path, app = app, cparams = cparams)
return module, inf
[ドキュメント]
def read_all_files(app, cparams, input_only = False):
"""
入力XRDファイルとすべてのCIFファイルを読み込む。
詳細説明:
`cparams.infile` で指定された入力XRDファイルを読み込み、その2θ範囲に基づいて解析範囲を調整する。
`cparams.cif_files` で指定されたマスクに一致するCIFファイルをすべて読み込み、変換処理を行う。
`input_only` が`True`の場合、入力ファイルのみを読み込む。
:param app: tkApplicationのインスタンス。
:param cparams: 設定パラメータを保持するオブジェクト。
:param input_only: bool. Trueの場合、入力ファイルのみを読み込む。
:returns: dict. 読み込まれたモジュールとデータを含む辞書。
キーは "module_input", "module_cif", "inf_input", "inf_cif_list"。
"""
print("")
print(f"read_all_files(): Read input file [{cparams.infile}]")
module_input, inf_input = read_file(cparams.infile, app, cparams)
if module_input:
# module_input.print_data(inf_input)
inf_input = module_input.convert(inf_input, cparams = cparams)
# save_data([cparams.outfile], inf_input, cparams = cparams)
# app.call(module_input, "plot_data", inf_input, cparams = cparams)
else:
app.terminate(f"Error in read_all_files(): Could not read [{cparams.infile}]", pause = True)
xQ2_infile = inf_input["data_list"][0]
xmin = min(xQ2_infile)
xmax = max(xQ2_infile)
xstep = xQ2_infile[1] - xQ2_infile[0]
print(f" 2Theta range: {xmin} - {xmax}, {xstep} step")
print(f" fwhm: {cparams.fwhm}")
print(f" Gaussian fraction: {cparams.Gfraction}")
cparams.xmin = max([cparams.xmin, xmin])
cparams.xmax = min([cparams.xmax, xmax])
cparams.xstep = xstep
print(f"2Theta range to be calculated: {cparams.xmin} - {cparams.xmax} degrees, {cparams.xstep} step")
if input_only:
inf = {
"module_input": module_input,
"inf_input" : inf_input,
}
return inf
cif_mask = cparams.cif_files
files = glob.glob(cif_mask)
print("")
print(f"Read cif and xlsx files from [{cif_mask}]")
inf_cif_list = []
module_cif = None
for f in files:
print("")
print(f" Read [{f}]")
dirname, basename, filebody, ext = split_file_path(f)
if len(filebody) == 0 or filebody[0] == '~':
print(f" [{basename}] has '~' at the first character. may be a temprary file. skip")
continue
if '-out.' in basename.lower():
print(f" [{basename}] include '-out.'. maybe an output file of some program. skip")
continue
module_cif, inf_cif = read_file(f, app, cparams)
if module_cif:
print(f" File [{f}] is red by [{module_cif.name}] module")
# app.call(module_cif, "print_data", inf_cif)
inf = app.call(module_cif, "convert", inf_cif, app = app, cparams = cparams)
# app.call(module_cif, "plot_data", inf_cif, cparams = cparams)
inf_cif_list.append(inf_cif)
inf = {
"module_input": module_input,
"module_cif" : module_cif,
"inf_input" : inf_input,
"inf_cif_list": inf_cif_list,
}
return inf
[ドキュメント]
def max_none(x):
"""
None値を含むリストまたはイテラブルの最大値を返す。
詳細説明:
入力されたリスト `x` の要素からNoneを除外し、残りの数値の最大値を計算する。
数値要素が一つもない場合は、非常に小さい負の値を返す。
:param x: list or iterable. 数値とNoneを含むリストまたはイテラブル。
:returns: float. 最大値。
"""
m = -1.0e100
for v in x:
if v is not None and m < v:
m = v
return m
[ドキュメント]
def min_none(x):
"""
None値を含むリストまたはイテラブルの最小値を返す。
詳細説明:
入力されたリスト `x` の要素からNoneを除外し、残りの数値の最小値を計算する。
数値要素が一つもない場合は、非常に大きい正の値を返す。
:param x: list or iterable. 数値とNoneを含むリストまたはイテラブル。
:returns: float. 最小値。
"""
m = 1.0e100
for v in x:
if v is not None and m > v:
m = v
return m
[ドキュメント]
def normalize_none(l, Amin = 0.0, Amax = 1.0, vmin = None, vmax = None):
"""
None値を含むリストの要素を指定された範囲に正規化する。
詳細説明:
リスト `l` の数値要素を `[vmin, vmax]` の範囲から `[Amin, Amax]` の範囲に線形変換して正規化する。
`vmin` または `vmax` が指定されない場合は、リスト内のNone以外の要素から自動的に計算される。
None値は変更されない。
:param l: list. 数値とNoneを含むリスト。
:param Amin: float. 正規化後の最小値。デフォルトは 0.0。
:param Amax: float. 正規化後の最大値。デフォルトは 1.0。
:param vmin: float, optional. 元データの最小値。指定しない場合は自動計算。
:param vmax: float, optional. 元データの最大値。指定しない場合は自動計算。
:returns: list. 正規化されたリスト。
"""
if vmax is None:
vmax = max_none(l)
if vmin is None:
vmin = min_none(l)
if vmax - vmin == 0.0:
vmax = vmin + 1.0
for i in range(len(l)):
if l[i] is None:
continue
l[i] = (l[i] - vmin) / (vmax - vmin) * (Amax - Amin) + Amin
return l
[ドキュメント]
def normalize(l, Amin = 0.0, Amax = 1.0, vmin = None, vmax = None):
"""
数値リストの要素を指定された範囲に正規化する。
詳細説明:
リスト `l` の要素を `[vmin, vmax]` の範囲から `[Amin, Amax]` の範囲に線形変換して正規化する。
`vmin` または `vmax` が指定されない場合は、リスト内の要素から自動的に計算される。
入力リストがNoneの場合はNoneを返す。
:param l: list. 数値のリスト。
:param Amin: float. 正規化後の最小値。デフォルトは 0.0。
:param Amax: float. 正規化後の最大値。デフォルトは 1.0。
:param vmin: float, optional. 元データの最小値。指定しない場合は自動計算。
:param vmax: float, optional. 元データの最大値。指定しない場合は自動計算。
:returns: list or None. 正規化されたリスト、または入力がNoneの場合はNone。
"""
if l is None:
return None
if vmax is None:
vmax = max(l)
if vmin is None:
vmin = min(l)
if vmax - vmin == 0.0:
vmax = vmin + 1.0
for i in range(len(l)):
l[i] = (l[i] - vmin) / (vmax - vmin) * (Amax - Amin) + Amin
return l
[ドキュメント]
def plot_all(inf, app, cparams):
"""
実験XRDパターンとCIFからの理論パターンを個別のサブプロットに表示する。
詳細説明:
入力された実験XRDパターンを最上段のサブプロットに表示し、
その後、読み込まれた各CIFファイルから計算されたXRDパターンをそれぞれ異なるサブプロットに表示する。
各理論パターンには回折ピーク位置もプロットされる。すべてのプロットは共通の2θ軸を共有する。
:param inf: dict. 読み込まれたデータを含む辞書。
:param app: tkApplicationのインスタンス。
:param cparams: 設定パラメータを保持するオブジェクト。
:returns: なし
"""
print("")
print("#########################################")
print(" Plot one XRD curve in a different box")
print("#########################################")
module_input = inf["module_input"]
module_cif = inf["module_cif"]
inf_input = inf["inf_input"]
inf_cif_list = inf["inf_cif_list"]
if module_cif is None:
app.terminate(f"Error in plot_all(): CIFファイルが見つかりませんでした", pause = True)
ncif = len(inf_cif_list)
print("")
print("plot")
print(f"yscale: {cparams.yscale}")
print("# of cif data:", ncif)
sample_infile = inf_input["sample_name"]
xQ2_infile = inf_input["data_list"][0]
yobs_infile = inf_input["data_list"][1]
if len(inf_input["data_list"]) >= 3:
ysim_infile = inf_input["data_list"][2]
else:
ysim_infile = None
vmax = max(yobs_infile)
vmin = min(yobs_infile)
yobs_infile = normalize(yobs_infile, Amin = 0.0, Amax = 1.0, vmin = vmin, vmax = vmax)
if ysim_infile is not None:
ysim_infile = normalize(ysim_infile, Amin = 0.0, Amax = 1.0, vmin = vmin, vmax = vmax)
fig, axes = plt.subplots(ncif + 1, 1, figsize = figsize, sharex=True, gridspec_kw={'hspace': 0})
plot_event = tkPlotEvent(plt, distance = 'x')
ncolors = len(colors)
nmarkers = len(markers)
ax0 = axes[0]
ax0.tick_params(labelsize = fontsize)
ax0.set_yticks([])
ax0.set_xlabel(None)
ax_bottom = axes[ncif]
for i in range(1, ncif + 1):
ax = axes[i]
ax.tick_params(labelsize = fontsize)
if i < ncif:
# ax.set_xticks([])
ax.set_xlabel(None)
ax.set_yticks([])
ax = axes[0]
ax.set_title(f"{sample_infile}")
ax.tick_params(labelsize = fontsize)
ax.plot(xQ2_infile, yobs_infile, label = "obs", color = colors[0], linewidth = 0.5)
if ysim_infile is not None:
ax.plot(xQ2_infile, ysim_infile, label = "sim", color = colors[1], linewidth = 0.3)
ax.axhline(0.0, linestyle = 'dashed', color = 'red', linewidth = 0.5)
if cparams.yscale == 'log':
ax.set_yscale('log')
if ysim_infile is None:
ymax = max(yobs_infile)
else:
ymax = max([max(ysim_infile), max(yobs_infile)])
ax.set_ylim([1.0e-4 * ymax, ax.get_ylim()[1]])
ax.legend(fontsize = legend_fontsize)
for i in range(1, ncif + 1):
inf_cif = inf_cif_list[i - 1]
xQ2, xrd_cal = inf_cif["conv_data"]
filename = inf_cif["filename"]
dirname, basename, filebody, ext = split_file_path(filename)
src = inf_cif["diffractions"]["source"]
Q2 = inf_cif["diffractions"]["Q2"]
dhkl = inf_cif["diffractions"]["dhkl"]
hkl = inf_cif["diffractions"]["hkl"]
mul = inf_cif["diffractions"].get("mul", None)
if mul is None:
mul = [0 for i in range(len(src))]
Int = inf_cif["diffractions"]["intensity"]
ndiffractions = len(Q2)
vmax = max(xrd_cal)
vmin = min(xrd_cal)
xrd_cal = normalize(xrd_cal, Amin = 0.0, Amax = 1.0, vmin = vmin, vmax = vmax)
Int = normalize(Int, Amin = 0.0, Amax = 1.0)
ax = axes[i]
color = colors[(i - 1) % ncolors]
marker = markers[(i - 1) % nmarkers]
phase = [filebody] * ndiffractions
ax.plot(xQ2, xrd_cal, label = filebody, color = 'black', linewidth = 0.5)
data = ax.plot(Q2, Int, linestyle = '', marker = marker, markerfacecolor = color, markeredgecolor = color, markersize = 2.0)
data0 = ax0.plot(Q2, Int, linestyle = '', marker = marker, markerfacecolor = color, markeredgecolor = color, markersize = 2.0)
for j in range(ndiffractions):
# data0 = ax0.plot([Q2[j], Q2[j]], [0.0, Int[j]], linestyle = 'dashed', color = color, linewidth = 0.5)
ax.plot([Q2[j], Q2[j]], [0.0, Int[j]], linestyle = 'dashed', color = color, linewidth = 0.5)
ax0.plot([Q2[j], Q2[j]], [0.0, Int[j]], linestyle = 'dashed', color = color, linewidth = 0.5)
#plot_event.add_data({"label": "diffraction", "plot_type": "2D", "axis": ax, "data": [Q2, dhkl, hkl, mul, Int]})
hkl_str = [f"{hkl[i][0]} {hkl[i][1]} {hkl[i][2]}" for i in range(ndiffractions)]
plot_event.add_data({"label": "diffraction", "plot_type": "2D", "axis": ax, "data": data,
"x": Q2, 'y': Int, "xlabel": "2Theta", "ylabel": "Intensity",
"xlist": [src, phase, Q2, dhkl, hkl_str, mul, Int],
"xlabels": ['X-ray', 'phase', 'Q2', 'dhkl', 'hkl', 'multiplicity', 'Intensity']
})
plot_event.add_data({"label": "diffraction", "plot_type": "2D", "axis": ax0, "data": data0,
"x": Q2, 'y': Int, "xlabel": "2Theta", "ylabel": "Intensity",
"xlist": [src, phase, Q2, dhkl, hkl_str, mul, Int],
"xlabels": ['X-ray', 'phase', 'Q2', 'dhkl', 'hkl', 'multiplicity', 'Intensity']
})
ax.axhline(0.0, linestyle = 'dashed', color = 'red', linewidth = 0.5)
if cparams.yscale == 'log':
ax.set_yscale('log')
ymax = max(xrd_cal)
ax.set_ylim([1.0e-4 * ymax, ax.get_ylim()[1]])
ax.legend(fontsize = legend_fontsize)
xmin = max([inf_input["xmin"], cparams.xmin])
xmax = min([inf_input["xmax"], cparams.xmax])
for ax in axes:
for q2 in range(int(xmin), int(xmax)):
if q2 % 5 == 0:
ax.axvline(q2, linestyle = 'dotted', linewidth = 0.5, color = 'black')
else:
ax.axvline(q2, linestyle = 'dotted', linewidth = 0.3, color = 'gray')
ax.set_xlim(xmin, xmax)
plot_event.register_click(fig) #callback = lambda event: plot_event.onclick(event))
# plot_event.register_event(fig, event = "button_press_event",
# callback = lambda event: plot_event.onclick(event))
ax_bottom.set_xlabel(r'2$\theta$', fontsize = fontsize)
axes[int(ncif/2)].set_ylabel('Intensity', fontsize = fontsize)
ax.legend(fontsize = legend_fontsize)
plt.tight_layout()
plt.pause(0.1)
input("Press ENTER to terminate>>")
[ドキュメント]
def plot_one(inf, app, cparams):
"""
実験XRDパターンとCIFからの回折ピーク位置を同じグラフに表示する。
詳細説明:
入力された実験XRDパターンと、読み込まれた各CIFファイルから計算された回折ピーク位置（強度は正規化）を
同一のグラフにプロットする。これにより、実験パターンと各結晶相のピーク位置の関係を一度に確認できる。
:param inf: dict. 読み込まれたデータを含む辞書。
:param app: tkApplicationのインスタンス。
:param cparams: 設定パラメータを保持するオブジェクト。
:returns: なし
"""
print("")
print("#########################################")
print(" Plot input XRD curve and cif diffraction angles in one graph")
print("#########################################")
module_input = inf["module_input"]
module_cif = inf["module_cif"]
inf_input = inf["inf_input"]
inf_cif_list = inf["inf_cif_list"]
if module_cif is None:
app.terminate(f"Error in plot_all(): CIFファイルが見つかりませんでした", pause = True)
ncif = len(inf_cif_list)
print("")
print("plot")
print(f"yscale: {cparams.yscale}")
print("# of cif data:", ncif)
sample_infile = inf_input["sample_name"]
xQ2_infile = inf_input["data_list"][0]
yobs_infile = inf_input["data_list"][1]
if len(inf_input["data_list"]) >= 3:
ysim_infile = inf_input["data_list"][2]
else:
ysim_infile = None
vmax = max(yobs_infile)
vmin = min(yobs_infile)
yobs_infile = normalize(yobs_infile, Amin = 0.0, Amax = 1.0, vmin = vmin, vmax = vmax)
if ysim_infile is not None:
ysim_infile = normalize(ysim_infile, Amin = 0.0, Amax = 1.0, vmin = vmin, vmax = vmax)
fig = plt.figure(figsize = figsize_one)
plot_event = tkPlotEvent(plt, distance = 'x')
ncolors = len(colors)
nmarkers = len(markers)
ax = fig.add_subplot(1, 1, 1)
# ax.set_xticks([])
ax.set_yticks([])
ax.set_xlabel(None)
ax.set_title(f"{sample_infile}")
ax.tick_params(labelsize = fontsize)
ax.plot(xQ2_infile, yobs_infile, label = "obs", color = colors[0], linewidth = 0.5)
if ysim_infile is not None:
ax.plot(xQ2_infile, ysim_infile, label = "sim", color = colors[1], linewidth = 0.3)
ax.legend(fontsize = legend_fontsize)
if ysim_infile is None:
ymax = max(yobs_infile)
else:
ymax = max([max(yobs_infile), max(ysim_infile)])
for i in range(1, ncif + 1):
inf_cif = inf_cif_list[i - 1]
xQ2, xrd_cal = inf_cif["conv_data"]
filename = inf_cif["filename"]
dirname, basename, filebody, ext = split_file_path(filename)
src = inf_cif["diffractions"]["source"]
Q2 = inf_cif["diffractions"]["Q2"]
dhkl = inf_cif["diffractions"]["dhkl"]
hkl = inf_cif["diffractions"]["hkl"]
mul = inf_cif["diffractions"].get("mul", None)
if mul is None:
mul = [0 for i in range(len(src))]
Int = inf_cif["diffractions"]["intensity"]
ndiffractions = len(Q2)
vmax = max(xrd_cal)
vmin = min(xrd_cal)
xrd_cal = normalize(xrd_cal, Amin = 0.0, Amax = 1.0, vmin = vmin, vmax = vmax)
Int = normalize(Int, Amin = 0.0, Amax = 1.0, vmin = vmin, vmax = vmax)
color = colors[(i - 1) % ncolors]
marker = markers[(i - 1) % nmarkers]
phase = [filebody] * ndiffractions
data0 = ax.plot(Q2, Int, linestyle = '', marker = marker, markerfacecolor = color, markeredgecolor = color, markersize = 2.0)
ymax = max([max(Int), ymax])
for j in range(ndiffractions):
# data0 = ax.plot([Q2[j], Q2[j]], [0.0, Int[j]], linestyle = 'dashed', color = color, linewidth = 0.5)
ax.plot([Q2[j], Q2[j]], [0.0, Int[j]], linestyle = 'dashed', color = color, linewidth = 0.5)
if len(hkl[0]) == 3:
hkl_str = [f"{hkl[i][0]} {hkl[i][1]} {hkl[i][2]}" for i in range(ndiffractions)]
else:
hkl_str = [f"{hkl[i][0]} {hkl[i][1]} ({hkl[i][2]}) {hkl[i][3]}" for i in range(ndiffractions)]
plot_event.add_data({"label": "diffraction", "plot_type": "2D", "axis": ax, "data": data0,
"x": Q2, 'y': Int, "xlabel": "2Theta", "ylabel": "Intensity",
"xlist": [src, phase, Q2, dhkl, hkl_str, mul, Int],
"xlabels": ["X-ray", 'phase', 'Q2', 'dhkl', 'hkl', 'multiplicity', 'Intensity']
})
ymax = max([max(Int), ymax])
ax.legend(fontsize = legend_fontsize)
plot_event.register_click(fig) #callback = lambda event: plot_event.onclick(event))
# plot_event.register_event(fig, event = "button_press_event",
# callback = lambda event: plot_event.onclick(event))
Qmin = max([inf_input["xmin"], cparams.xmin])
Qmax = min([inf_input["xmax"], cparams.xmax])
ax.set_xlim(Qmin, Qmax)
ax.set_ylabel('Intensity', fontsize = fontsize)
if cparams.yscale == 'log':
ax.set_yscale('log')
ax.set_ylim([1.0e-4 * ymax, ax.get_ylim()[1]])
plt.tight_layout()
plt.pause(0.1)
input("Press ENTER to terminate>>")
[ドキュメント]
def plot_one2(inf, app, cparams):
"""
実験XRDパターンとCIFからの理論XRDパターンを同じグラフに表示する。
詳細説明:
入力された実験XRDパターンと、読み込まれた各CIFファイルから計算された理論XRDパターン（ブロードニング適用済み）を
同一のグラフにプロットする。これにより、実験パターンと各結晶相の回折プロファイルの形状を一度に比較できる。
:param inf: dict. 読み込まれたデータを含む辞書。
:param app: tkApplicationのインスタンス。
:param cparams: 設定パラメータを保持するオブジェクト。
:returns: なし
"""
print("")
print("#########################################")
print(" Plot XRD curves in one graph")
print("#########################################")
module_input = inf["module_input"]
module_cif = inf["module_cif"]
inf_input = inf["inf_input"]
inf_cif_list = inf["inf_cif_list"]
if module_cif is None:
app.terminate(f"Error in plot_all(): CIFファイルが見つかりませんでした", pause = True)
ncif = len(inf_cif_list)
print("")
print("plot")
print(f"yscale: {cparams.yscale}")
print("# of cif data:", ncif)
sample_infile = inf_input["sample_name"]
xQ2_infile = inf_input["data_list"][0]
yobs_infile = inf_input["data_list"][1]
if len(inf_input["data_list"]) >= 3:
ysim_infile = inf_input["data_list"][2]
else:
ysim_infile = None
vmax = max(yobs_infile)
vmin = min(yobs_infile)
yobs_infile = normalize(yobs_infile, Amin = 0.0, Amax = 1.0, vmin = vmin, vmax = vmax)
if ysim_infile is not None:
ysim_infile = normalize(inf_input["data_list"][2], Amin = 0.0, Amax = 1.0, vmin = vmin, vmax = vmax)
fig = plt.figure(figsize = figsize_one)
plot_event = tkPlotEvent(plt, distance = 'x')
ncolors = len(colors)
nmarkers = len(markers)
ax = fig.add_subplot(1, 1, 1)
# ax.set_xticks([])
ax.set_yticks([])
ax.set_xlabel(None)
ax.set_title(f"{sample_infile}")
ax.tick_params(labelsize = fontsize)
ax.plot(xQ2_infile, yobs_infile, label = "obs", color = colors[0], linewidth = 1.0)
# ax.plot(xQ2_infile, ysim_infile, label = "sim", color = colors[1], linewidth = 0.3)
ymax = max(yobs_infile)
for i in range(1, ncif + 1):
inf_cif = inf_cif_list[i - 1]
xQ2, xrd_cal = inf_cif["conv_data"]
filename = inf_cif["filename"]
dirname, basename, filebody, ext = split_file_path(filename)
vmax = max(xrd_cal)
vmin = min(xrd_cal)
xrd_cal = normalize(xrd_cal, Amin = 0.0, Amax = 1.0, vmin = vmin, vmax = vmax)
ax.plot(xQ2, xrd_cal, label = filebody, linestyle = '-', color = colors[i % ncolors], linewidth = 0.5)
ymax = max([ymax, max(xrd_cal)])
Qmin = max([inf_input["xmin"], cparams.xmin])
Qmax = min([inf_input["xmax"], cparams.xmax])
ax.set_xlim(Qmin, Qmax)
ax.set_ylabel('Intensity', fontsize = fontsize)
if cparams.yscale == 'log':
ax.set_yscale('log')
ax.set_ylim([1.0e-4 * ymax, ax.get_ylim()[1]])
ax.legend(fontsize = legend_fontsize)
plt.tight_layout()
plt.pause(0.1)
input("Press ENTER to terminate>>")
[ドキュメント]
def convolution(x, y, filter, nskip = 1):
"""
データ `y` を指定された `filter` で畳み込む。
詳細説明:
離散的な畳み込み演算を実行し、`y` の各点において `filter` を適用して
平滑化された（またはブロードニングされた）データ `yc` を生成する。
`nskip` パラメータで畳み込みの計算間隔を制御し、結果のデータ点数を減らすことができる。
:param x: list or array. x座標のリスト。
:param y: list or array. 畳み込み対象のy座標のリスト。
:param filter: list or array. 畳み込みに使用するフィルタの配列。
:param nskip: int. 畳み込み計算をスキップする間隔。1の場合、すべての点で計算する。
:returns: tuple. (`_x`, `_y`). 畳み込み後のx座標とy座標のリスト。`y` がNoneの場合はNoneを返す。
"""
if y is None:
return None
ny = len(y)
nconv = len(filter) // 2
yc = np.zeros(ny)
for i in range(0, ny, nskip):
for idf in range(-nconv, nconv + 1):
i1 = i + idf
if i1 % nskip == 0 and 0 <= i1 < ny:
yc[i] += filter[idf] * y[i1]
_x = [x[i] for i in range(0, ny, nskip)]
_y = [yc[i] for i in range(0, ny, nskip)]
return _x, _y
# return x, yc
[ドキュメント]
def collect_data(inf, app, cparams, is_print = True):
"""
フィッティングや相関分析のために、実験および理論XRDデータを収集・前処理する。
詳細説明:
入力された実験XRDデータと各CIFファイルから計算された理論XRDデータを統一された2θ範囲に補間し、
スマアリング（畳み込み）と正規化を適用する。さらに、背景多項式のための基底関数も準備する。
:param inf: dict. 読み込まれたデータを含む辞書。
:param app: tkApplicationのインスタンス。
:param cparams: 設定パラメータを保持するオブジェクト。
:param is_print: bool. 処理の詳細をコンソールに出力するかどうか。
:returns: tuple. (Q2_list, yobs_infile, ysim_infile, bg_names, bg_list, sample_names, ycif_list, labels_train, x_train)。
- `Q2_list`: list. 共通の2θ座標のリスト。
- `yobs_infile`: list. 前処理された実験XRD強度データ。
- `ysim_infile`: list or None. 前処理されたシミュレーションXRD強度データ（存在する場合）。
- `bg_names`: list. 背景多項式のラベルのリスト。
- `bg_list`: list. 背景多項式の基底関数（レジェンドル多項式）のリスト。
- `sample_names`: list. CIFファイルから取得した試料名またはファイルボディ名のリスト。
- `ycif_list`: list. 各CIFファイルから計算され前処理された理論XRD強度データのリスト。
- `labels_train`: list. 学習データ（x_train）の各列に対応するラベルのリスト。
- `x_train`: list. 学習データとなる背景多項式と理論XRD強度のリスト。
"""
module_input = inf["module_input"]
inf_input = inf["inf_input"]
inf_cif_list = inf["inf_cif_list"]
ncif = len(inf_cif_list)
nspectrum = inf_input["nspectrum"]
if ncif < 1:
app.terminate(f"\nError: CIF file is not found.\n", pause = True)
fwhm = cparams.fwhm_smear
Gf = cparams.Gfraction_smear
dQ2 = 12.0 * fwhm
nconv = int(dQ2 / cparams.xstep + 1.00001)
filter = [0.0] * (2 * nconv + 1)
if fwhm <= 0.0:
wGL = 1.0
else:
wGL = GaussLorentz(0.0, 0.0, fwhm / 2.0, C0 = 1.0, Gfraction = Gf, Gwratio = 1.0, A = None)
for i in range(-nconv, nconv + 1):
if fwhm <= 0.0:
filter[i] = 0.0
else:
filter[i] = GaussLorentz(0.0, i * cparams.xstep, fwhm / 2.0, C0 = 1.0, Gfraction = Gf, Gwratio = 1.0, A = None) / wGL
sample_infile = inf_input["sample_name"]
xQ2_infile = inf_input["data_list"][0]
nx_in = len(xQ2_infile)
dx = xQ2_infile[1] - xQ2_infile[0]
nskip = int(fwhm / dx / 5.0 + 1.0e-5)
if nskip == 0:
nskip = 1
minnx = 500
if int(nx_in / nskip) < minnx:
nskip = int(nx_in / minnx)
# nskip = 1
yobs_infile = inf_input["data_list"][1]
if nspectrum >= 3:
ysim_infile = inf_input["data_list"][2]
else:
ysim_infile = None
if cparams.yscale == 'log':
yobs_infile = log(yobs_infile)
if ysim_infile is not None:
ysim_infile = log(ysim_infile)
if fwhm > 0.0:
_xQ2_infile, yobs_infile = convolution(xQ2_infile, yobs_infile, filter, nskip = nskip)
if ysim_infile is not None:
_xQ2_infile, ysim_infile = convolution(xQ2_infile, ysim_infile, filter, nskip = nskip)
else:
_xQ2_infile = xQ2_infile
vmax = max(yobs_infile)
vmin = min(yobs_infile)
yobs_infile = normalize(yobs_infile, Amin = 0.0, Amax = 1.0, vmin = vmin, vmax = vmax)
if ysim_infile is not None:
ysim_infile = normalize(ysim_infile, Amin = 0.0, Amax = 1.0, vmin = vmin, vmax = vmax)
print(f"Read CIF data from {ncif} files")
xQ2 = None
ycif_list = []
sample_names = []
for i in range(1, ncif + 1):
inf_cif = inf_cif_list[i - 1]
xQ2, xrd_cal = inf_cif["conv_data"]
maxy = abs(max(xrd_cal))
if cparams.yscale == 'log':
xrd_cal = log(xrd_cal + maxy * 1.0e-7)
if fwhm > 0.0:
_xQ2, xrd_cal = convolution(xQ2, xrd_cal, filter, nskip = nskip)
else:
_xQ2 = xQ2
xrd_cal = normalize(xrd_cal, Amin = 0.0, Amax = 1.0, vmin = vmin, vmax = vmax)
xrd_cal = normalize(xrd_cal, Amin = 0.0, Amax = 1.0)
filename = inf_cif["filename"]
dirname, basename, filebody, ext = split_file_path(filename)
ycif_list.append(xrd_cal)
sample_names.append(filebody)
nx = len(xQ2)
xmin = min(xQ2)
xmax = max(xQ2)
xstep = (xmax - xmin) / (nx - 1)
x_list = []
Q2_list = []
for i in range(nx):
x_list.append(-1 + i * xstep)
Q2_list.append(xmin + i * xstep)
func1d = interp1d(_xQ2_infile, yobs_infile, bounds_error = False, fill_value = (0.0, 0.0))
# func1d = interp1d(xQ2_infile, yobs_infile, bounds_error = False, fill_value = (0.0, 0.0))
yobs_infile = func1d(Q2_list)
if ysim_infile is not None:
func1d = interp1d(_xQ2_infile, ysim_infile, bounds_error = False, fill_value = (0.0, 0.0))
# func1d = interp1d(xQ2_infile, ysim_infile, bounds_error = False, fill_value = (0.0, 0.0))
ysim_infile = func1d(Q2_list)
for i in range(1, ncif + 1):
func1d = interp1d(_xQ2, ycif_list[i-1], bounds_error = False, fill_value = (0.0, 0.0))
# func1d = interp1d(xQ2, ycif_list[i-1], bounds_error = False, fill_value = (0.0, 0.0))
ycif_list[i-1] = func1d(Q2_list)
bg_list = []
for order in range(cparams.BGorder + 1):
poly1d = legendre(order)
bg_list.append(poly1d(x_list))
if is_print:
print("# of cif data:", ncif)
print("BG order :", cparams.BGorder)
bg_names = [f"BG_order_{i}" for i in range(cparams.BGorder + 1)] # BG_order+1個の多項式基底
x_train = bg_list.copy()
labels_train = bg_names.copy()
x_train.extend(ycif_list)
labels_train.extend(sample_names) # CIFファイルの名前を追加
return Q2_list, yobs_infile, ysim_infile, bg_names, bg_list, sample_names, ycif_list, labels_train, x_train
[ドキュメント]
def check_overwrap(inf, app, cparams):
"""
異なるCIFファイル由来のXRDパターン間のオーバーラップ（重なり）をチェックし、可視化する。
詳細説明:
各CIFファイルの回折ピークが、他のCIFファイルから生成された回折パターンとどの程度重なっているかを計算し、
その結果をプロットする。これにより、複数の結晶相が混在する場合に、各相のピークが互いに干渉する度合いを評価できる。
:param inf: dict. 読み込まれたデータを含む辞書。
:param app: tkApplicationのインスタンス。
:param cparams: 設定パラメータを保持するオブジェクト。
:returns: なし
"""
print("")
print("#########################################")
print(" Check overwrap between cif XRD patterns")
print("#########################################")
print(f" smearing: fwhm={cparams.fwhm_smear} Gaussian fraction={cparams.Gfraction_smear}")
inf_cif_list = inf["inf_cif_list"]
Q2_list, yobs_infile, ysim_infile, bg_names, bg_list, sample_names, ycif_list, labels_train, x_train = collect_data(inf, app, cparams)
nx = len(Q2_list)
Qmin = Q2_list[0]
Qmax = Q2_list[nx - 1]
ncif = len(ycif_list)
nxtrain = len(x_train)
fwhm = cparams.fwhm + cparams.fwhm_smear
Gf = cparams.Gfraction_smear
dQ2 = fwhm # FWHMを考慮した畳み込み幅
nconv = int(dQ2 / cparams.xstep + 1.00001)
filter = [0.0] * (2 * nconv + 1)
if fwhm <= 0.0:
wGL = 1.0
else:
# FWHMの半分を半値幅として使用
wGL = GaussLorentz(0.0, 0.0, fwhm / 2.0, C0 = 1.0, Gfraction = Gf, Gwratio = 1.0, A = None)
for i in range(-nconv, nconv + 1):
if fwhm <= 0.0:
filter[i] = 0.0
else:
filter[i] = GaussLorentz(0.0, i * cparams.xstep, fwhm / 2.0, C0 = 1.0, Gfraction = Gf, Gwratio = 1.0, A = None) / wGL
# print("nconv=", nconv)
# print("filter=", filter)
# 構造のピークが他の構造のブロードニングされたパターンに与える影響を計算
corr = make_matrix3(ncif, ncif, nx)
for ic in range(ncif):
inf = inf_cif_list[ic]
Q2 = inf["diffractions"]["Q2"]
Int = inf["diffractions"]["intensity"]
nd = len(Q2)
ndiffractions = len(Q2)
for ic2 in range(ncif):
if ic == ic2:
continue
ycif = ycif_list[ic2] # 他のCIFパターン
for id in range(nd):
Q20 = Q2[id] # 現在のCIFのピーク位置
id0 = None
for i in range(nx):
if Q20 <= Q2_list[i]:
id0 = i
break
if id0 is None:
continue
# ピーク位置 Q20 を中心に、ブロードニングフィルターを適用
for idf in range(-nconv, nconv + 1):
i1 = id0 + idf
if 0 <= i1 < nx:
# 現在のピーク強度 * フィルタ値 * 他のCIFパターンの強度
# これはオーバーラップの度合いを示す
v = filter[idf] * Int[id] * ycif[i1]
if corr[ic][ic2][i1] is None:
corr[ic][ic2][i1] = v
else:
corr[ic][ic2][i1] += v
# print("corr=", corr)
fig, axes = plt.subplots(ncif, 1, figsize = figsize, sharex=True, gridspec_kw={'hspace': 0})
plot_event = tkPlotEvent(plt, distance = 'x')
ncolors = len(colors)
nmarkers = len(markers)
for ic in range(ncif):
inf = inf_cif_list[ic]
src = inf["diffractions"]["source"]
Q2 = inf["diffractions"]["Q2"]
dhkl = inf["diffractions"]["dhkl"]
hkl = inf["diffractions"]["hkl"]
mul = inf["diffractions"].get("mul", None)
if mul is None:
mul = [0 for i in range(len(src))]
Int = inf["diffractions"]["intensity"]
ndiffractions = len(Q2)
filename = inf["filename"]
dirname, basename, filebody, ext = split_file_path(filename)
phase = [filebody] * ndiffractions
ax = axes[ic]
ax.tick_params(labelsize = fontsize)
ax2 = ax.twinx()
ax2.tick_params(labelsize = 0)
ycif = ycif_list[ic]
ax.plot([], [], label = filebody, linestyle = '') # フェーズ名を表示するためのダミープロット
ax2.plot(Q2_list, ycif, picker = True, linestyle = 'dashed', color = 'black', linewidth = 0.5)
for ic2 in range(ncif):
if ic == ic2:
continue
inf_cif2 = inf_cif_list[ic2]
filename = inf_cif2["filename"]
dirname, basename, filebody, ext = split_file_path(filename)
color = colors[ic2 % ncolors]
marker = markers[ic2 % ncolors]
y = corr[ic][ic2]
vmax = max_none(y)
vmin = min_none(y)
y = normalize_none(y, Amin = 0.0, Amax = 1.0, vmin = vmin, vmax = vmax)
Int = normalize(Int, Amin = 0.0, Amax = 1.0) # ピーク強度を正規化
ax.plot(Q2_list, y, label = filebody, picker = True, color = color, linewidth = 1.0) # オーバーラッププロット
data = ax.plot(Q2, Int, linestyle = '', marker = marker, markerfacecolor = color, markeredgecolor = 'black', markersize = 2.0)
for j in range(ndiffractions):
ax.plot([Q2[j], Q2[j]], [0.0, Int[j]], linestyle = 'dashed', color = 'black', linewidth = 0.5)
hkl_str = [f"{hkl[i][0]} {hkl[i][1]} {hkl[i][2]}" for i in range(ndiffractions)]
plot_event.add_data({"label": "diffraction", "plot_type": "2D", "axis": ax, "data": data,
"x": Q2, 'y': Int, "xlabel": "2Theta", "ylabel": "Intensity",
"xlist": [src, phase, Q2, dhkl, hkl_str, mul, Int],
"xlabels": ['X-ray', 'phase', 'Q2', 'dhkl', 'hkl', 'multiplicity', 'Intensity']
})
ax.set_xlim([Qmin, Qmax])
ax.legend(fontsize = legend_fontsize)
if cparams.yscale == 'log':
ax.set_yscale('log')
ax.set_ylim([1.0e-4, ax.get_ylim()[1]])
ax2.set_yscale('linear') # ax2のスケールは常にリニア
ax.set_xlabel(r'2$\theta$', fontsize = fontsize)
plot_event.register_click(fig) # callback = lambda event: plot_event.onclick(event))
plot_event.register_pick(fig) # callback = lambda event: plot_event.onclick(event))
plt.tight_layout()
plt.pause(0.1)
input("Press ENTER to terminate>>")
[ドキュメント]
def CIF_correlation(inf, app, cparams):
"""
CIFファイルから生成されたXRDパターン間の相関係数を計算し、可視化する。
詳細説明:
すべてのCIFファイルから生成された理論XRDパターン間の相関係数を計算し、結果をテキストで出力する。
また、各CIFパターンと入力実験パターンとの相関係数を、スマアリングFWHMを変化させながらプロットし、
最も相関の高いスマアリング条件や、異なる相間の類似性を評価する。
:param inf: dict. 読み込まれたデータを含む辞書。
:param app: tkApplicationのインスタンス。
:param cparams: 設定パラメータを保持するオブジェクト。
:returns: なし
"""
print("")
print("#########################################")
print(" Check correlation between cif XRD patterns")
print("#########################################")
print(f" smearing: fwhm={cparams.fwhm_smear} Gaussian fraction={cparams.Gfraction_smear}")
inf_cif_list = inf["inf_cif_list"]
Q2_list, yobs_infile, ysim_infile, bg_names, bg_list, sample_names, ycif_list, labels_train, x_train = collect_data(inf, app, cparams)
xlist = [yobs_infile]
xlabels = ["obs"]
xlist.extend(ycif_list)
xlabels.extend(sample_names)
nx = len(xlist)
dx = Q2_list[1] - Q2_list[0]
#相関係数
print("")
print("Correlation coefficients:")
corr = np.zeros([nx, nx])
for i in range(nx):
for j in range(i, nx):
corr[i][j] = np.dot(xlist[i], xlist[j])
for i in range(nx):
for j in range(i + 1, nx):
corr[i][j] /= np.sqrt(corr[i][i] * corr[j][j])
if corr[i][j] >= 0.9:
print(f" identical? {corr[i][j]:8.4f}: ({i:2}-{j:2}): ({xlabels[i]:20})-({xlabels[j]:20})")
elif corr[i][j] >= 0.8:
print(f" similar {corr[i][j]:8.4f}: ({i:2}-{j:2}): ({xlabels[i]:20})-({xlabels[j]:20})")
elif corr[i][j] >= 0.5:
print(f" ***** {corr[i][j]:8.4f}: ({i:2}-{j:2}): ({xlabels[i]:20})-({xlabels[j]:20})")
elif corr[i][j] >= 0.5:
print(f" ***** {corr[i][j]:8.4f}: ({i:2}-{j:2}): ({xlabels[i]:20})-({xlabels[j]:20})")
elif corr[i][j] >= 0.4:
print(f" **** {corr[i][j]:8.4f}: ({i:2}-{j:2}): ({xlabels[i]:20})-({xlabels[j]:20})")
elif corr[i][j] >= 0.3:
print(f" *** {corr[i][j]:8.4f}: ({i:2}-{j:2}): ({xlabels[i]:20})-({xlabels[j]:20})")
else:
print(f" {corr[i][j]:8.4f}: ({i:2}-{j:2}): ({xlabels[i]:20})-({xlabels[j]:20})")
# input dataとの相関係数をsmearing fwhmを変えながらプロット
print("")
print("Correlation coefficients with input data with varied smearing fwhm:")
fwhm_smear_original = cparams.fwhm_smear
fwhm_list = np.arange(0.0, 1.001, 0.2)
corr_list = make_matrix2(nx, len(fwhm_list))
for i in range(len(fwhm_list)):
_fwhm = fwhm_list[i]
cparams.fwhm_smear = _fwhm
print(f" Calculating for fwhm={_fwhm}")
Q2_list0, yobs_infile0, ysim_infile0, bg_names0, bg_list0, sample_names0, ycif_list0, labels_train0, x_train0 \
= collect_data(inf, app, cparams, is_print = False)
nskip = int((_fwhm / 5.0) / dx + 1.0e-5)
# nskip = 1
xlist0 = [[v for i, v in enumerate(yobs_infile0) if i % nskip == 0]]
for y in ycif_list0:
xlist0.append([v for i, v in enumerate(y) if i % nskip == 0])
# xlist0 = [yobs_infile]
# xlist0.extend(ycif_list0)
nx_sampled = len(xlist0)
print(f" _fwhm={_fwhm} dx={dx} nskip={nskip} nx={len(xlist0[0])}")
# nskip = int((_fwhm / 4.0) / dx)
# print("nskip = ", nskip)
# print("")
norm = make_matrix1(nx_sampled)
for j in range(nx_sampled):
print(f" calculating {j}-th diagonal norm")
norm[j] = np.sqrt(np.dot(xlist0[j], xlist0[j]))
for j in range(nx_sampled):
print(f" calculating obs - {j}-th non-diagonal correlation")
corr_list[j][i] = np.dot(xlist0[0], xlist0[j]) / norm[0] / norm[j]
cparams.fwhm_smear = fwhm_smear_original
print("")
print("plot")
fig, ax = plt.subplots(1, 1, figsize = figsize_one)
plot_event = tkPlotEvent(plt, distance = 'x')
ncolors = len(colors)
for i in range(1, nx): # "obs"を除いたCIFパターンとobsの相関をプロット
ax.tick_params(labelsize = fontsize)
ax.plot(fwhm_list, corr_list[i], label = xlabels[i], picker = True, color = colors[i % ncolors], linewidth = 1.0)
ax.legend(fontsize = legend_fontsize)
ax.set_xlabel(r'fwhm ($\degree$)', fontsize = fontsize)
ax.set_ylabel(r'Correlation coefficient', fontsize = fontsize)
plot_event.register_pick(fig) # callback = lambda event: plot_event.onclick(event))
plt.tight_layout()
plt.pause(0.1)
print("")
input("Press ENTER to terminate>>")
[ドキュメント]
def correlation(inf, app, cparams):
"""
入力された実験XRDパターンとCIFファイルから生成された理論パターンとの相関係数を計算し、可視化する。
詳細説明:
実験XRDパターンと、読み込まれた各CIFファイルから計算された理論XRDパターン間の相関係数を計算し、
結果をテキストで出力する。さらに、これらの相関係数をスマアリングFWHMを変化させながらプロットし、
最も相関の高い相や適切なスマアリング条件を評価する。
:param inf: dict. 読み込まれたデータを含む辞書。
:param app: tkApplicationのインスタンス。
:param cparams: 設定パラメータを保持するオブジェクト。
:returns: なし
"""
print("")
print("#########################################")
print(" Check correlation between the input XRD and the cif XRD patterns")
print("#########################################")
print(f" smearing: fwhm={cparams.fwhm_smear} Gaussian fraction={cparams.Gfraction_smear}")
inf_cif_list = inf["inf_cif_list"]
Q2_list, yobs_infile, ysim_infile, bg_names, bg_list, sample_names, ycif_list, labels_train, x_train = collect_data(inf, app, cparams)
xlist = [yobs_infile]
xlabels = ["obs"]
xlist.extend(ycif_list)
xlabels.extend(sample_names)
nx = len(xlist)
dx = Q2_list[1] - Q2_list[0]
#相関係数
print("")
print("Correlation coefficients:")
corr = np.zeros([nx, nx])
for i in range(nx):
for j in range(i, nx):
corr[i][j] = np.dot(xlist[i], xlist[j])
for i in range(nx):
for j in range(i + 1, nx):
corr[i][j] /= np.sqrt(corr[i][i] * corr[j][j])
if corr[i][j] >= 0.9:
print(f" identical? {corr[i][j]:8.4f}: ({i:2}-{j:2}): ({xlabels[i]:20})-({xlabels[j]:20})")
elif corr[i][j] >= 0.8:
print(f" similar {corr[i][j]:8.4f}: ({i:2}-{j:2}): ({xlabels[i]:20})-({xlabels[j]:20})")
elif corr[i][j] >= 0.5:
print(f" ***** {corr[i][j]:8.4f}: ({i:2}-{j:2}): ({xlabels[i]:20})-({xlabels[j]:20})")
elif corr[i][j] >= 0.5:
print(f" ***** {corr[i][j]:8.4f}: ({i:2}-{j:2}): ({xlabels[i]:20})-({xlabels[j]:20})")
elif corr[i][j] >= 0.4:
print(f" **** {corr[i][j]:8.4f}: ({i:2}-{j:2}): ({xlabels[i]:20})-({xlabels[j]:20})")
elif corr[i][j] >= 0.3:
print(f" *** {corr[i][j]:8.4f}: ({i:2}-{j:2}): ({xlabels[i]:20})-({xlabels[j]:20})")
else:
print(f" {corr[i][j]:8.4f}: ({i:2}-{j:2}): ({xlabels[i]:20})-({xlabels[j]:20})")
# input dataとの相関係数をsmearing fwhmを変えながらプロット
print("")
print("Correlation coefficients with input data with varied smearing fwhm:")
fwhm_smear_original = cparams.fwhm_smear
fwhm_list = np.arange(0.0, 1.001, 0.2)
corr_list = make_matrix2(nx, len(fwhm_list))
for i in range(len(fwhm_list)):
_fwhm = fwhm_list[i]
cparams.fwhm_smear = _fwhm
print("")
print(f" Calculating for fwhm={_fwhm}")
Q2_list0, yobs_infile0, ysim_infile0, bg_names0, bg_list0, sample_names0, ycif_list0, labels_train0, x_train0 \
= collect_data(inf, app, cparams, is_print = False)
nskip = int((_fwhm / 5.0) / dx + 1.0e-5)
# nskip = 1
xlist0 = [[v for i, v in enumerate(yobs_infile0) if i % nskip == 0]]
for y in ycif_list0:
xlist0.append([v for i, v in enumerate(y) if i % nskip == 0])
# xlist0 = [yobs_infile]
# xlist0.extend(ycif_list0)
nx_sampled = len(xlist0)
print(f" _fwhm={_fwhm} dx={dx} nskip={nskip} nx={len(xlist0[0])}")
norm = make_matrix1(nx_sampled)
for j in range(nx_sampled):
print(f" calculating {j}-th diagonal norm")
norm[j] = np.sqrt(np.dot(xlist0[j], xlist0[j]))
for j in range(nx_sampled):
print(f" calculating obs - {j}-th non-diagonal correlation")
corr_list[j][i] = np.dot(xlist0[0], xlist0[j]) / norm[0] / norm[j]
cparams.fwhm_smear = fwhm_smear_original
print("")
print("plot")
fig, ax = plt.subplots(1, 1, figsize = figsize_one)
plot_event = tkPlotEvent(plt, distance = 'x')
ncolors = len(colors)
for i in range(1, nx): # "obs"を除いたCIFパターンとobsの相関をプロット
ax.tick_params(labelsize = fontsize)
ax.plot(fwhm_list, corr_list[i], label = xlabels[i], picker = True, color = colors[i % ncolors], linewidth = 1.0)
ax.legend(fontsize = legend_fontsize)
ax.set_xlabel(r'smearing FWHM ($\degree$)', fontsize = fontsize)
ax.set_ylabel(r'Correlation coefficient', fontsize = fontsize)
plot_event.register_pick(fig) # callback = lambda event: plot_event.onclick(event))
plt.tight_layout()
plt.pause(0.1)
print("")
input("Press ENTER to terminate>>")
[ドキュメント]
def fit(inf, app, cparams):
"""
LASSO回帰を用いて実験XRDパターンを背景とCIF由来の理論パターンでフィッティングする。
詳細説明:
実験XRDパターンを目的変数、背景多項式と各CIFファイルから計算された理論XRDパターンを説明変数として、
LASSO回帰を実行する。異なるalpha値での係数変化とRMSEの変化をプロットし、最終的に指定されたalpha値での
フィッティング結果（係数、RMSE、MAV）とフィッティングされたパターンをプロットで表示する。
:param inf: dict. 読み込まれたデータを含む辞書。
:param app: tkApplicationのインスタンス。
:param cparams: 設定パラメータを保持するオブジェクト。
:returns: なし
"""
print("")
print("#########################################")
print(" Fitting by LASSO")
print("#########################################")
print(f" alpha={cparams.alpha}")
print(f" smearing: FWHM={cparams.fwhm_smear} Gaussian fraction={cparams.Gfraction_smear}")
print("yscale:", cparams.yscale)
Q2_list, yobs_infile, ysim_infile, bg_names, bg_list, sample_names, ycif_list, labels_train, x_train \
= collect_data(inf, app, cparams)
nx = len(Q2_list)
ncif = len(ycif_list)
nxtrain = len(x_train)
print("")
print("LASSO alpha:", cparams.alpha)
x_train = np.array(x_train).T
alpha0 = 1.0e-8
ntry = 80
print("")
print(f"Lasso regression with varied alpha: start from {alpha0}")
# BG項のラベルも含むため、labels_trainの先頭から表示
print(f"{'':42}", labels_train)
_alpha = alpha0
alpha_list = []
c_list = []
RMSE_list = []
maxC = None
for i in range(ntry):
model = Lasso(alpha = _alpha, fit_intercept = False)
model.fit(x_train, yobs_infile)
yfit = model.predict(x_train)
MSE = mean_squared_error(yobs_infile, yfit)
RMSE = np.sqrt(MSE)
alpha_list.append(_alpha)
c_list.append(model.coef_)
RMSE_list.append(RMSE)
s = ["{:10.4g}".format(v) for v in model.coef_]
s = " ".join(s)
print(f" {i:2} alpha={_alpha:8.2g} RMSE={RMSE:12.4g} coeff={s}")
if maxC is None:
maxC = max(abs(max(model.coef_)), abs(min(model.coef_))) if len(model.coef_) > 0 else 0.0
else:
if len(model.coef_) > 0 and abs(max(model.coef_)) < 1.0e-5 * maxC and abs(min(model.coef_)) < 1.0e-5 * maxC:
break
_alpha *= 2.0
plot_event = tkPlotEvent(plt, distance = 'x')
print("")
print("plot LASSO analysis")
# coeffプロットはBG項も含めて全係数を表示するが、凡例はCIF名のみ
print("sample_names=", sample_names)
fig_lasso = plt.figure(figsize = figsize_coeff)
ax = fig_lasso.add_subplot(1, 1, 1)
ax2 = ax.twinx()
ax.tick_params(labelsize = fontsize)
ax2.tick_params(labelsize = fontsize)
ax.plot(alpha_list, RMSE_list, label = 'RMSE', picker = True, linestyle = 'dashed', color = 'black', linewidth = 1.5)
# BG項の係数もプロットする場合
# for j in range(len(bg_names)):
# ax2.plot(alpha_list, np.array(c_list).T[j], label = bg_names[j], picker = True, linewidth = 0.8, linestyle = 'dotted')
# CIF項の係数をプロット
for i in range(len(sample_names)):
# model.coef_ のうち、CIF成分に対応するインデックスは BGorder + 1 以降
ax2.plot(alpha_list, np.array(c_list).T[cparams.BGorder + 1 + i], label = sample_names[i], picker = True, linewidth = 1.0, color = colors[i % ncolors])
ax2.axhline(0.0, linestyle = 'dashed', color = 'red', linewidth = 0.5)
ax.set_xlabel('alpha', fontsize = fontsize)
ax.set_ylabel('RMSE', fontsize = fontsize)
ax2.set_ylabel('coefficient', fontsize = fontsize)
ax.set_xscale('log')
ax.set_yscale('log')
ax2.legend(fontsize = legend_fontsize_one)
plt.tight_layout()
plt.pause(0.1)
model = Lasso(alpha = cparams.alpha, fit_intercept = False)
model.fit(x_train, yobs_infile)
yfit = model.predict(x_train)
print("")
print(f"Lasso regression with alpha={cparams.alpha}")
print(f" Coefficients:")
# print(f" coefficients:", model.coef_)
yBG = np.zeros(nx)
for order in range(cparams.BGorder + 1):
print(f" {bg_names[order]}: {model.coef_[order]:10.6g}")
yBG += model.coef_[order] * bg_list[order] # BG多項式の寄与を計算
# CIF項の係数と、それらを乗じたパターンを更新
for i in range(len(sample_names)):
c = model.coef_[cparams.BGorder + 1 + i] # BG項の後にCIF項が続く
print(f" {sample_names[i]:20}: {c:10.6g}")
# collect_dataでycif_listは正規化された状態なので、ここで係数を乗じる
ycif_list[i] *= c
MAE = mean_absolute_error(yobs_infile, yfit)
MSE = mean_squared_error(yobs_infile, yfit)
RMSE = np.sqrt(MSE)
print(f" MAE : {MAE:12.4g}")
print(f" MSE : {MSE:12.4g}")
print(f" RMSE: {RMSE:12.4g}")
print("")
print("plot XRD profiles")
fig = plt.figure(figsize = figsize_one)
ncolors = len(colors)
ax = fig.add_subplot(1, 1, 1)
ax.tick_params(labelsize = fontsize)
ax.plot(Q2_list, yobs_infile, label = "obs", picker = True, linestyle = '',
marker = 'o', markersize = 2.0, markerfacecolor = colors[0], markeredgecolor = colors[0])
ax.plot(Q2_list, yBG, label = "background", picker = True, linestyle = '-', color = 'cyan', linewidth = 0.5)
# ax.plot(Q2_list, yobs_infile, label = "obs", color = colors[0], linewidth = 1.5)
ax.plot(Q2_list, yfit, label = "fit", picker = True, linestyle = 'dashed', color = colors[1], linewidth = 1.5)
# ax.plot(Q2_list, ysim_infile, label = "obs", color = colors[1], linewidth = 0.5)
# フィッティングされた各CIFパターンをプロット
for i in range(len(sample_names)):
ax.plot(Q2_list, ycif_list[i], label = f"{sample_names[i]}", picker = True, color = colors[(i+2) % ncolors ], linewidth = 1.2) # +2 は obs, fit, BGを避けるため
# for order in range(cparams.BGorder):
# ax.plot(Q2_list, bg_list[order], label = f"{order+1}-th order", color = colors[(ncif+2+order) % ncolors], linewidth = 0.5)
ax.set_xlabel('$2\\theta$', fontsize = fontsize)
if cparams.yscale == 'log':
ax.set_ylabel('log($y$)', fontsize = fontsize)
else:
ax.set_ylabel('$y$', fontsize = fontsize)
ax.legend(fontsize = legend_fontsize_one)
plot_event.register_pick(fig_lasso) # callback = lambda event: plot_event.onclick(event))
plot_event.register_pick(fig) # callback = lambda event: plot_event.onclick(event))
plt.tight_layout()
plt.pause(0.1)
input("Press ENTER to terminate>>")
[ドキュメント]
def main():
"""
メイン関数。XRD解析スクリプトのエントリポイント。
詳細説明:
アプリケーションの初期化、パラメータの読み込みと更新、プラグインモジュールのロード、
および指定されたモードに応じたXRD解析処理（プロット、相関分析、フィッティングなど）を実行する。
処理のログはファイルにも出力される。
:returns: なし
"""
global module_names, modules
#==================================================================
# Initialize parameters
#==================================================================
app = tkApplication(usage_str = usage_str, globals = globals(), locals = locals())
cparams = app.get_params()
logfile = app.replace_path(None, template = ["{dirname}", "{filebody}-out.txt"])
# logfile = app.replace_path(cparams.infile, template = ["{dirname}", "{filebody}-out.txt"])
print(f"Open logfile [{logfile}]")
app.redirect(targets = ["stdout", logfile], mode = 'w')
initialize(app, cparams)
update_vars(app, cparams)
cparams.outfile = replace_path(cparams.infile, template = os.path.join("{dirname}", "{filebody}-out.txt"))
print("")
print( "==========================================================================")
print(" Convert CIF file to powder XRD pattern")
print( "==========================================================================")
print(f"mode: {cparams.mode}")
print(f"Plug-in dir: {cparams.plugin_dir}")
print(f"Input file: {cparams.infile}")
print(f"Output file: {cparams.outfile}")
# Load modules
print("")
print(f"Load modules:")
module_names, modules = app.load_modules(cparams.plugin_dir, "*.py", target = "read_data", is_print = True)
for m in modules:
# if hasattr(m, "initialize"):
# print(f"initialize {m.name}")
# m.initialize(app, cparams)
input_type = m.get_input_type(app = app, cparams = cparams)
output_type = m.get_output_type(app = app, cparams = cparams)
print(f" {m.name}: input_type={input_type} output_type={output_type}")
if cparams.mode == 'plot':
inf = read_all_files(app, cparams, input_only = True)
plot_input(inf, app, cparams)
elif cparams.mode == 'overwrap':
inf = read_all_files(app, cparams)
check_overwrap(inf, app, cparams)
elif cparams.mode == 'CIFcorrelation':
inf = read_all_files(app, cparams)
CIF_correlation(inf, app, cparams)
elif cparams.mode == 'correlation':
inf = read_all_files(app, cparams)
correlation(inf, app, cparams)
elif cparams.mode == 'plot_all':
inf = read_all_files(app, cparams)
plot_all(inf, app, cparams)
elif cparams.mode == 'plot_one':
inf = read_all_files(app, cparams)
plot_one(inf, app, cparams)
elif cparams.mode == 'plot_one2':
inf = read_all_files(app, cparams)
plot_one2(inf, app, cparams)
elif cparams.mode == 'fit':
inf = read_all_files(app, cparams)
fit(inf, app, cparams)
app.terminate(usage = usage)
if __name__ == "__main__":
main()