import os
os.environ["KERAS_BACKEND"] = "jax" # JAXバックエンド
import keras
from keras import ops
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
plt.rcParams['font.family'] = 'Noto Sans CJK JP'
# ==========================
# データ読み込み
# ==========================
data_label = np.load("data_label.npz")
data = data_label["data"]
label = data_label["label"]
x_train, x_test, y_train, y_test = train_test_split(data, label, test_size=0.3, random_state=0)
# ==========================
# ラベルスケーリング
# ==========================
scaler_y = StandardScaler()
y_train_f = scaler_y.fit_transform(y_train)
y_test_f = scaler_y.transform(y_test)
# PINN内で逆変換するために平均と標準偏差をTensorとして保持
# float32型にキャストしておく
y_mean = keras.ops.convert_to_tensor(scaler_y.mean_, dtype="float32")
y_scale = keras.ops.convert_to_tensor(scaler_y.scale_, dtype="float32")
# ==========================
# 物理法則 (Circuit Theory)
# ==========================
def calculate_s_parameters_from_LC(freq_norm, lc_params):
"""
freq_norm: 正規化周波数 (batch, 200)
lc_params: 物理単位のLC値 (batch, 6) -> [L1, C1, L2, C2, L3, C3]
"""
f_stop = 20.0e9 # 20 GHz
z0 = 50.0 + 0j
# 周波数を物理単位(Hz)に戻して角周波数計算
freq_hz = freq_norm * f_stop
omega = 2.0 * np.pi * freq_hz
omega = ops.cast(omega, "complex64")
epsilon = 1e-9 # 0除算防止
# (batch, 1) に拡張してブロードキャスト
L1 = ops.expand_dims(lc_params[:, 0], axis=1) * 1e-9
C1 = ops.expand_dims(lc_params[:, 1], axis=1) * 1e-12
L2 = ops.expand_dims(lc_params[:, 2], axis=1) * 1e-9
C2 = ops.expand_dims(lc_params[:, 3], axis=1) * 1e-12
L3 = ops.expand_dims(lc_params[:, 4], axis=1) * 1e-9
C3 = ops.expand_dims(lc_params[:, 5], axis=1) * 1e-12
# Series (L1, C1)
Zs1 = 1j * omega * L1 + 1.0 / (1j * omega * C1 + epsilon)
# Shunt (L2, C2)
Yp2 = 1j * omega * C2 + 1.0 / (1j * omega * L2 + epsilon)
# Series (L3, C3)
Zs3 = 1j * omega * L3 + 1.0 / (1j * omega * C3 + epsilon)
# ABCD行列 (A, B, C, D) の要素計算
A = 1.0 + Zs1 * Yp2
B = (1.0 + Zs1 * Yp2) * Zs3 + Zs1
C = Yp2
D = Yp2 * Zs3 + 1.0
# Sパラメータ変換
delta = A + B/z0 + C*z0 + D
S11 = (A + B/z0 - C*z0 - D) / delta
S21 = 2.0 / delta
return S11, S21
# ==========================
# PINN モデル定義 (compute_lossのみオーバーライド)
# ==========================
class CircuitPINN(keras.Model):
def __init__(self, inputs, outputs, y_mean, y_scale, physics_weight=0.1, **kwargs):
super().__init__(inputs=inputs, outputs=outputs, **kwargs)
self.y_mean = y_mean
self.y_scale = y_scale
self.physics_weight = physics_weight
# Loss監視用メトリクス
self.data_loss_tracker = keras.metrics.Mean(name="data_loss")
self.phys_loss_tracker = keras.metrics.Mean(name="phys_loss")
# Keras 3 ではここをオーバーライドするのが推奨
def compute_loss(self, x=None, y=None, y_pred=None, sample_weight=None):
# 1. Data Loss (教師データとの誤差)
# y: 教師データ(正規化済み), y_pred: 予測値(正規化済み)
data_loss = ops.mean(ops.square(y - y_pred))
# 2. Physics Loss (物理法則との整合性)
# 正規化された予測値を物理単位(nH, pF)に戻す
y_pred_physical = y_pred * self.y_scale + self.y_mean
# L, Cは正の値であるべきなので絶対値をとる（またはReLU）
y_pred_physical = ops.abs(y_pred_physical)
# 入力 x から周波数と正解Sパラメータを抽出
# x shape: (batch, 200, 5) -> [freq, S11r, S11i, S21r, S21i]
freq_norm = x[:, :, 0]
s11_true_real = x[:, :, 1]
s11_true_imag = x[:, :, 2]
s21_true_real = x[:, :, 3]
s21_true_imag = x[:, :, 4]
# 予測されたLC値からSパラメータを計算
s11_pred, s21_pred = calculate_s_parameters_from_LC(freq_norm, y_pred_physical)
# Sパラメータの誤差を計算 (実部・虚部)
phys_loss = (
ops.mean(ops.square(s11_true_real - ops.real(s11_pred))) +
ops.mean(ops.square(s11_true_imag - ops.imag(s11_pred))) +
ops.mean(ops.square(s21_true_real - ops.real(s21_pred))) +
ops.mean(ops.square(s21_true_imag - ops.imag(s21_pred)))
)
total_loss = data_loss + self.physics_weight * phys_loss
# メトリクスの更新
self.data_loss_tracker.update_state(data_loss)
self.phys_loss_tracker.update_state(phys_loss)
return total_loss
@property
def metrics(self):
# fit() のプログレスバーに表示されるメトリクス一覧
return [self.data_loss_tracker, self.phys_loss_tracker]
# ==========================
# モデル構築
# ==========================
keras.utils.set_random_seed(1)
inputs = keras.Input(shape=(200, 5))
x = keras.layers.Flatten()(inputs)
x = keras.layers.Dense(256, activation="relu")(x)
x = keras.layers.BatchNormalization()(x)
x = keras.layers.Dense(256, activation="relu")(x)
x = keras.layers.BatchNormalization()(x)
#x = keras.layers.Dropout(0.1)(x)
outputs = keras.layers.Dense(6)(x)
# モデル初期化
model = CircuitPINN(
inputs=inputs,
outputs=outputs,
y_mean=y_mean,
y_scale=y_scale,
physics_weight=0.1 # 重みは適宜調整してください
)
# コンパイル (loss引数はcompute_lossで定義しているので不要)
model.compile(optimizer=keras.optimizers.Adam(learning_rate=0.001))
batch_size = 32
epochs = 300
# 学習
history = model.fit(
x_train,
y_train_f,
batch_size=batch_size,
epochs=epochs,
validation_split=0.15,
)
# ==========================
# 評価・可視化
# ==========================
# 推論
y_pred_f = model.predict(x_test)
y_pred = scaler_y.inverse_transform(y_pred_f)
# R2スコア
metric = keras.metrics.R2Score()
metric.update_state(y_test, y_pred)
print("\nR2 Score:", metric.result())
# 誤差率計算
eps = 1e-12
valid_mask = np.abs(y_test) > eps
pct_error = np.zeros_like(y_test)
pct_error[valid_mask] = np.abs(
(y_test[valid_mask] - y_pred[valid_mask]) / y_test[valid_mask] * 100
)
print("\n% Error per element (mean):", pct_error.mean(axis=0))
# プロット
row, column = 2, 3
legend = ["L1", "C1", "L2", "C2", "L3", "C3"]
fig, ax = plt.subplots(2, 3, figsize=(15,9))
for i in range(row):
for j in range(column):
count = column * i + j
maxvalue = max(y_pred[:, count].max(), y_test[:,count].max())
ax[i,j].scatter(y_pred[:, count], y_test[:,count], c="r", s=5)
ax[i,j].plot([0,maxvalue], [0,maxvalue], "--", c="black")
ax[i,j].set_xlabel("推定した値")
ax[i,j].set_ylabel("実際の値")
ax[i,j].set_xlim(0, maxvalue)
ax[i,j].set_ylim(0, maxvalue)
ax[i,j].grid()
ax[i,j].legend([legend[count] + f" 平均誤差{pct_error.mean(axis=0)[count]:.2f}%"])
fig.tight_layout()
plt.show()
|
最近のコメント