From 6693f3916cb669b74b5e2f261d328eb65263e26c Mon Sep 17 00:00:00 2001 From: huangfu <3045324663@qq.com> Date: Sun, 19 Apr 2026 18:33:36 +0800 Subject: [PATCH] =?UTF-8?q?=E4=BF=AE=E6=94=B9=E8=84=9A=E6=9C=AC=E7=BB=93?= =?UTF-8?q?=E6=9E=84MOE?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- MZM_MoE_PINN_Model.ipynb | 949 +++++++++++++++++++++++++++++++++++++++ README.md | 203 ++++----- best_hyperparams.json | 11 + configs/default.yaml | 73 +-- src/config.py | 203 ++++----- src/evaluate.py | 83 ++-- src/infer.py | 23 +- src/main.py | 26 +- src/model.py | 129 ++++-- src/plots.py | 28 +- src/preprocess.py | 547 ++++++---------------- src/trainer.py | 262 ++++++----- tests/test_smoke.py | 234 ++++------ 13 files changed, 1675 insertions(+), 1096 deletions(-) create mode 100644 MZM_MoE_PINN_Model.ipynb create mode 100644 best_hyperparams.json diff --git a/MZM_MoE_PINN_Model.ipynb b/MZM_MoE_PINN_Model.ipynb new file mode 100644 index 0000000..760129e --- /dev/null +++ b/MZM_MoE_PINN_Model.ipynb @@ -0,0 +1,949 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "id": "E7c5HVjM1J5Q" + }, + "source": [ + "# MZM MoE PINN Model\n", + "\n", + "## Overview\n", + "This notebook implements a Physics-Informed Neural Network (PINN) architecture for modeling a Mach-Zehnder Modulator (MZM). The model is built using PyTorch and includes:\n", + "\n", + "- Data preprocessing\n", + "- Custom activation functions\n", + "- Flexible neural network design\n", + "- Training loop with learning rate scheduling\n", + "- Performance visualization\n", + "\n", + "The goal is to approximate nonlinear relationships between device parameters and output performance metrics.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 418 + }, + "id": "vazB3wCHBm1c", + "outputId": "19953917-22ca-48ed-fbb3-b3b976ce6819" + }, + "outputs": [], + "source": [ + "# ==============================\n", + "# Import Required Libraries\n", + "# ==============================\n", + "import numpy as np\n", + "import pandas as pd\n", + "import matplotlib.pyplot as plt\n", + "from io import StringIO\n", + "from sklearn.model_selection import train_test_split\n", + "from sklearn.preprocessing import StandardScaler\n", + "import re\n", + "import torch\n", + "import torch.nn as nn\n", + "from torch.optim.lr_scheduler import StepLR\n", + "from torch.utils.data import DataLoader, TensorDataset" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "id": "b459hcyW1J5S" + }, + "outputs": [], + "source": [ + "# random seed for selecting test set\n", + "random_state = 123\n", + "\n", + "# Custum functions are defined here\n", + "class GaussianActivation(nn.Module):\n", + " def forward(self, x):\n", + " return torch.exp(-x**2)\n", + "\n", + "def weights_init(layer_in):\n", + " # Initialize the parameters with He initialization\n", + " if isinstance(layer_in, nn.Linear):\n", + " nn.init.kaiming_uniform_(layer_in.weight)\n", + " # nn.init.xavier_uniform_(layer_in.weight)\n", + " layer_in.bias.data.fill_(0.0)\n", + "\n", + "def create_flexible_nn(input_dim, output_dim, hidden_dims, activation_fn=nn.ReLU(),positive_output=False):\n", + " \"\"\"\n", + " Create a flexible neural network with variable number of hidden layers\n", + "\n", + " Args:\n", + " input_dim (int): Number of input features\n", + " output_dim (int): Number of output targets\n", + " hidden_dims (list): List of integers specifying hidden layer dimensions\n", + " activation_fn: Activation function to use between layers\n", + " \"\"\"\n", + " layers = []\n", + "\n", + " # Input layer\n", + " layers.append(nn.Linear(input_dim, hidden_dims[0]))\n", + " layers.append(activation_fn)\n", + "\n", + " # Hidden layers\n", + " for i in range(len(hidden_dims) - 1):\n", + " layers.append(nn.Linear(hidden_dims[i], hidden_dims[i+1]))\n", + " layers.append(activation_fn)\n", + "\n", + " # Output layer\n", + " layers.append(nn.Linear(hidden_dims[-1], output_dim))\n", + "\n", + " # If positivity is requested, append a Softplus\n", + " if positive_output:\n", + " layers.append(nn.Softplus())\n", + "\n", + " return nn.Sequential(*layers)\n", + "\n", + "class ExpertNN(nn.Module):\n", + " def __init__(self, input_dim, output_dim, hidden_dims, activation_fn=nn.ReLU(),\n", + " dropout_rate=0.0, use_bn=False):\n", + " super().__init__()\n", + " layers = []\n", + "\n", + " prev_dim = input_dim\n", + " for h in hidden_dims:\n", + " layers.append(nn.Linear(prev_dim, h))\n", + " if use_bn:\n", + " layers.append(nn.BatchNorm1d(h))\n", + " layers.append(activation_fn)\n", + " if dropout_rate > 0:\n", + " layers.append(nn.Dropout(p=dropout_rate))\n", + " prev_dim = h\n", + " layers.append(nn.Linear(prev_dim, output_dim))\n", + "\n", + " self.net = nn.Sequential(*layers)\n", + "\n", + " def forward(self, x):\n", + " return self.net(x)\n", + "\n", + "\n", + "class MixtureOfExperts(nn.Module):\n", + " def __init__(self, input_dim, output_dim, hidden_dims, n_experts=3,\n", + " activation_fn=nn.ReLU(), gating_hidden=32,\n", + " dropout_rate=0.0, use_bn=False):\n", + " super().__init__()\n", + " self.experts = nn.ModuleList([\n", + " ExpertNN(input_dim, output_dim, hidden_dims,\n", + " activation_fn, dropout_rate, use_bn)\n", + " for _ in range(n_experts)\n", + " ])\n", + " self.gating = nn.Sequential(\n", + " nn.Linear(input_dim, gating_hidden),\n", + " nn.ReLU(),\n", + " nn.Linear(gating_hidden, n_experts),\n", + " nn.Softmax(dim=1)\n", + " )\n", + "\n", + " def forward(self, x):\n", + " # Gating weights\n", + " gate_weights = self.gating(x) # [batch, n_experts]\n", + " # Experts’ outputs\n", + " expert_outputs = torch.stack([expert(x) for expert in self.experts], dim=2) # [batch, output_dim, n_experts]\n", + " # Weighted sum\n", + " out = torch.bmm(expert_outputs, gate_weights.unsqueeze(2)).squeeze(2) # [batch, output_dim]\n", + " return out" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "SbevfhHm1J5T" + }, + "source": [ + "## Neural Network Architecture Design\n", + "\n", + "In this section, we define:\n", + "\n", + "- Input dimension (`D_i`)\n", + "- Output dimension (`D_o`)\n", + "- Hidden layer structure (`D_h`)\n", + "- Activation function\n", + "\n", + "The architecture is intentionally flexible so it can be tuned for bias-variance tradeoff optimization.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "id": "d0Hkww28B5Da" + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "MixtureOfExperts(\n", + " (experts): ModuleList(\n", + " (0-49): 50 x ExpertNN(\n", + " (net): Sequential(\n", + " (0): Linear(in_features=8, out_features=8, bias=True)\n", + " (1): BatchNorm1d(8, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True)\n", + " (2): ReLU()\n", + " (3): Linear(in_features=8, out_features=6, bias=True)\n", + " (4): BatchNorm1d(6, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True)\n", + " (5): ReLU()\n", + " (6): Linear(in_features=6, out_features=6, bias=True)\n", + " (7): BatchNorm1d(6, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True)\n", + " (8): ReLU()\n", + " (9): Linear(in_features=6, out_features=3, bias=True)\n", + " )\n", + " )\n", + " )\n", + " (gating): Sequential(\n", + " (0): Linear(in_features=8, out_features=8, bias=True)\n", + " (1): ReLU()\n", + " (2): Linear(in_features=8, out_features=50, bias=True)\n", + " (3): Softmax(dim=1)\n", + " )\n", + ")\n", + "\n", + "Number of parameters: 11972\n" + ] + } + ], + "source": [ + "###### Design the NN architecture\n", + "\n", + "acivation_fun = nn.ReLU()\n", + "# acivation_fun = nn.Tanh()\n", + "# acivation_fun = GaussianActivation()\n", + "\n", + "D_i = 8 # Input dimensions\n", + "D_o = 3 # Output dimensions\n", + "# Hidden layer dimensions\n", + "# D_h = [64, 128, 256, 128, 64]\n", + "# D_h = [200, 300, 350, 300, 200]\n", + "# D_h = [8,7,6,5,4]\n", + "# D_h = [8,6,6]\n", + "D_h = [64, 128, 64]\n", + "\n", + "#positive_output = False # True if we want a constraint to ensure +ve \"BW_3dB\", \"IL\", \"V_pi\" DOESN'T WORK WITH NORMALIZATION\n", + "#model = create_flexible_nn(D_i, D_o, D_h, acivation_fun,positive_output)\n", + "n_experts = 60 # number of experts in MoE model\n", + "gating_hidden_layer_width = 8 #width of the gating hidden layer\n", + "dropout_rate = 0.0 # dropout rate (regularization parameter)\n", + "use_bn = True # Batch normalization flag\n", + "\n", + "model_data = MixtureOfExperts(\n", + " D_i, D_o, D_h,\n", + " n_experts=n_experts,\n", + " activation_fn=acivation_fun,\n", + " gating_hidden=gating_hidden_layer_width,\n", + " dropout_rate=dropout_rate,\n", + " use_bn=use_bn\n", + ")\n", + "\n", + "print(model_data)\n", + "print(f\"\\nNumber of parameters: {sum(p.numel() for p in model_data.parameters())}\")\n", + "\n", + "# Design the NN training module\n", + "\n", + "# SGD parameters\n", + "batch_size = 128\n", + "learning_rate = 1e-3 #1e-3\n", + "weight_decay = 0.05 # Regularization parameter 0.05\n", + "momentum = 0.9 # used in momentum SGD optimizer\n", + "betas = (0.9, 0.999) # used in Adam optimizer beta1=0.9, beta2=0.999\n", + "\n", + "# Schedular parameters\n", + "LR_scheduler_gamma = 0.5 # decrease learning by 0.5 every N steps\n", + "LR_scheduler_step = 100 # step size parameter N of the LR scheduler\n", + "n_epoch = 100 # loop over the dataset n_epoch, e.g., 200 times\n", + "\n", + "# define MSE as loss function (regression problem)\n", + "loss_function = nn.MSELoss()\n", + "\n", + "# construct SGD optimizer\n", + "# optimizer = torch.optim.SGD(model.parameters(), lr = learning_rate, weight_decay=weight_decay, momentum=momentum)\n", + "optimizer = torch.optim.AdamW(model_data.parameters(), lr=learning_rate, weight_decay=weight_decay, betas=betas)\n", + "\n", + "# learning rate schedular by half every 10 epochs\n", + "# scheduler = StepLR(optimizer, step_size=LR_scheduler_step, gamma=LR_scheduler_gamma)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "id": "zYO1hW9L1J5U" + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " PN_offset Bias_V Core_width P+_width N+_width \\\n", + "0 -2.150000e-07 -10.0 4.500000e-07 2.062650e-07 1.083980e-07 \n", + "1 -2.150000e-07 -10.0 4.500000e-07 3.059360e-07 3.005950e-07 \n", + "2 -2.150000e-07 -10.0 5.081920e-07 2.085450e-07 1.936980e-07 \n", + "3 -2.150000e-07 -10.0 5.146500e-07 2.985000e-07 3.407150e-07 \n", + "4 -2.150000e-07 -10.0 5.180810e-07 1.589440e-07 2.348480e-07 \n", + "\n", + " P_width N_width Phase_length BW_3dB IL V_pi \n", + "0 1.000000e-06 7.551270e-07 0.002402 51.200000 2.58569 30.6728 \n", + "1 1.000000e-06 6.000000e-07 0.003196 36.179487 3.00849 39.7924 \n", + "2 8.262140e-07 7.634800e-07 0.002991 42.461539 2.82829 33.4881 \n", + "3 6.387450e-07 8.128130e-07 0.002953 38.692308 2.67160 40.2870 \n", + "4 6.000000e-07 6.000000e-07 0.001487 60.800000 1.39794 75.1647 \n" + ] + } + ], + "source": [ + "# Mount Drive\n", + "# from google.colab import drive\n", + "# drive.mount('/content/drive')\n", + "\n", + "# Load and scale data\n", + "# file_path = \"/content/drive/MyDrive/MZM Data/Sim_generated_dataset.txt\"\n", + "file_path = \"Sim_generated_dataset.txt\"\n", + "\n", + "with open(file_path) as f:\n", + " cleaned = [re.sub(r'[\\[\\]]', '', line.strip()) for line in f]\n", + "\n", + "df = pd.read_csv(StringIO(\"\\n\".join(cleaned)), header=None)\n", + "df.columns = [\n", + " \"PN_offset\", \"Bias_V\", \"Core_width\", \"P+_width\", \"N+_width\",\n", + " \"P_width\", \"N_width\", \"Phase_length\", \"BW_3dB\", \"IL\", \"V_pi\"\n", + "]\n", + "print(df.head()) # for debugging only to check the data" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "id": "JjcbkxaR1J5U" + }, + "outputs": [], + "source": [ + "# Remove the rows where Vpi is above a threshold (data cleaning)\n", + "df_cleaned = df[df['V_pi'] < 500].copy()\n", + "\n", + "# Split data into features/targets\n", + "feature_cols = df_cleaned.columns[:8] # first 8 columns\n", + "target_cols = df_cleaned.columns[8:] # last 3 columns\n", + "\n", + "x = df_cleaned[feature_cols].values\n", + "y = df_cleaned[target_cols].values\n", + "# print(y[:3]) # for debugging only to check the data\n", + "\n", + "# Train-test split\n", + "x_train_raw, x_test_raw, y_train_raw, y_test_raw = train_test_split(\n", + " x, y, test_size=0.1, random_state=random_state)\n", + "\n", + "# Fit scalers on TRAINING data only\n", + "scaler_x = StandardScaler()\n", + "\n", + "# Scale each target variable SEPARATELY\n", + "scaler_y1 = StandardScaler() # For BW_3dB\n", + "scaler_y2 = StandardScaler() # For IL\n", + "scaler_y3 = StandardScaler() # For V_pi\n", + "\n", + "# Scale features (x) - one scaler is fine\n", + "x_train_scaled = scaler_x.fit_transform(x_train_raw)\n", + "\n", + "# Scale each target dimension INDEPENDENTLY\n", + "y_train_scaled_1 = scaler_y1.fit_transform(y_train_raw[:, 0:1]) # BW_3dB\n", + "y_train_scaled_2 = scaler_y2.fit_transform(y_train_raw[:, 1:2]) # IL\n", + "y_train_scaled_3 = scaler_y3.fit_transform(y_train_raw[:, 2:3]) # V_pi\n", + "\n", + "# Combine targets back together\n", + "y_train_scaled = np.hstack([y_train_scaled_1, y_train_scaled_2, y_train_scaled_3])\n", + "\n", + "# Transform TEST data using training parameters\n", + "x_test_scaled = scaler_x.transform(x_test_raw)\n", + "\n", + "# Transform each test target separately\n", + "y_test_scaled_1 = scaler_y1.transform(y_test_raw[:, 0:1]) # BW_3dB\n", + "y_test_scaled_2 = scaler_y2.transform(y_test_raw[:, 1:2]) # IL\n", + "y_test_scaled_3 = scaler_y3.transform(y_test_raw[:, 2:3]) # V_pi\n", + "y_test_scaled = np.hstack([y_test_scaled_1, y_test_scaled_2, y_test_scaled_3])\n", + "\n", + "# Convert to tensors\n", + "x_train = torch.from_numpy(x_train_scaled).float()\n", + "y_train = torch.from_numpy(y_train_scaled).float()\n", + "x_test = torch.from_numpy(x_test_scaled).float()\n", + "y_test = torch.from_numpy(y_test_scaled).float()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "vZvTrgGptLKj" + }, + "source": [ + "The outputs are:\n", + "\n", + "- BW_3dB\n", + "\n", + "- IL\n", + "\n", + "- V_pi\n", + "\n", + "These are macroscopic performance metrics, not modeling the field profile directly, so solving full Maxwell is overkill. The most reasonable PINN constraint here is Phase Accumulation Relation enforced by gradients.\n", + "\n", + "## ✔ Constraint 1: Monotonic Bandwidth vs Length\n", + "\n", + "Physically:\n", + "\n", + "$$\n", + "BW \\propto \\frac{1}{L}\n", + "$$\n", + "\n", + "So:\n", + "\n", + "$$\n", + "\\frac{\\partial BW}{\\partial L} \\le 0\n", + "$$\n", + "\n", + "This must be enforced via autograd, not sorting.\n", + "\n", + "---\n", + "\n", + "## ✔ Constraint 2: Insertion Loss increases with Length\n", + "\n", + "$$\n", + "\\frac{\\partial IL}{\\partial L} \\ge 0\n", + "$$\n", + "\n", + "Longer waveguide → more absorption.\n", + "\n", + "---\n", + "\n", + "## ✔ Constraint 3: Vπ scales inversely with L\n", + "\n", + "$$\n", + "V_\\pi \\cdot L = constant\n", + "$$\n", + "\n", + "Instead of using dataset mean, enforce:\n", + "\n", + "$$\n", + "\\frac{\\partial (V_\\pi L)}{\\partial L} \\approx 0\n", + "$$\n", + "\n", + "---\n", + "\n", + "## ✔ Constraint 4: Smoothness Constraint\n", + "\n", + "Device metrics should be smooth functions of geometry.\n", + "\n", + "So penalize second derivatives:\n", + "\n", + "$$\n", + "\\left|\\frac{\\partial^2 y}{\\partial L^2}\\right|\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "id": "z03ciqUp1J5V" + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Using device: cuda\n" + ] + } + ], + "source": [ + "# ==============================\n", + "# GPU Setup\n", + "# ==============================\n", + "device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\n", + "print(f\"Using device: {device}\")\n", + "\n", + "# Move model and all data to GPU\n", + "model_data.to(device)\n", + "x_train = x_train.to(device)\n", + "y_train = y_train.to(device)\n", + "x_test = x_test.to(device)\n", + "y_test = y_test.to(device)\n", + "\n", + "# Precompute constant for Vpi*L on GPU\n", + "with torch.no_grad():\n", + " VpiL_const = (y_train[:, 2] * x_train[:, 7]).mean()\n", + "\n", + "# Store polynomial fit coefficients on GPU\n", + "L_train_np = x_train[:,7].cpu().numpy() # polyfit uses numpy, move to CPU\n", + "BW_train_np = y_train[:,0].cpu().numpy()\n", + "coeffs = np.polyfit(L_train_np, BW_train_np, deg=2) # [a2, a1, a0]\n", + "BW_coeffs = torch.tensor(coeffs, dtype=torch.float32, device=device)\n", + "\n", + "# Training errors storage\n", + "errors_train = np.zeros((n_epoch))\n", + "errors_test = np.zeros((n_epoch))" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "id": "Z5Fa_mZttHEc" + }, + "outputs": [ + { + "data": { + "text/plain": [ + "MixtureOfExperts(\n", + " (experts): ModuleList(\n", + " (0-49): 50 x ExpertNN(\n", + " (net): Sequential(\n", + " (0): Linear(in_features=8, out_features=8, bias=True)\n", + " (1): BatchNorm1d(8, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True)\n", + " (2): ReLU()\n", + " (3): Linear(in_features=8, out_features=6, bias=True)\n", + " (4): BatchNorm1d(6, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True)\n", + " (5): ReLU()\n", + " (6): Linear(in_features=6, out_features=6, bias=True)\n", + " (7): BatchNorm1d(6, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True)\n", + " (8): ReLU()\n", + " (9): Linear(in_features=6, out_features=3, bias=True)\n", + " )\n", + " )\n", + " )\n", + " (gating): Sequential(\n", + " (0): Linear(in_features=8, out_features=8, bias=True)\n", + " (1): ReLU()\n", + " (2): Linear(in_features=8, out_features=50, bias=True)\n", + " (3): Softmax(dim=1)\n", + " )\n", + ")" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# -- Precompute constant for VpiL\n", + "with torch.no_grad():\n", + " VpiL_const = (y_train[:, 2] * x_train[:, 7]).mean()\n", + "\n", + "# Define weights for each physics term\n", + "lambda_vpiL = 0.0\n", + "lambda_bw_mon = 0.0\n", + "lambda_IL_mon = 0.0\n", + "lambda_pos = 0.0\n", + "lambda_bw_poly = 0.2\n", + "\n", + "data_loader = DataLoader(\n", + " TensorDataset(x_train, y_train),\n", + " batch_size=batch_size,\n", + " shuffle=True,\n", + " worker_init_fn=np.random.seed(1)\n", + ")\n", + "\n", + "model_data.apply(weights_init)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "id": "zhkJ08bl1J5W" + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Epoch 0, train MSE 0.5561, test MSE 0.5682\n", + "Epoch 10, train MSE 0.0699, test MSE 0.0889\n", + "Epoch 20, train MSE 0.0549, test MSE 0.0779\n", + "Epoch 30, train MSE 0.0506, test MSE 0.0749\n", + "Epoch 40, train MSE 0.0478, test MSE 0.0732\n", + "Epoch 50, train MSE 0.0461, test MSE 0.0711\n", + "Epoch 60, train MSE 0.0444, test MSE 0.0705\n", + "Epoch 70, train MSE 0.0439, test MSE 0.0708\n", + "Epoch 80, train MSE 0.0425, test MSE 0.0699\n", + "Epoch 90, train MSE 0.0406, test MSE 0.0715\n" + ] + } + ], + "source": [ + "# ==============================\n", + "# Training Loop\n", + "# ==============================\n", + "for epoch in range(n_epoch):\n", + " model_data.train()\n", + " for x_batch, y_batch in data_loader:\n", + " # Move batch to GPU\n", + " x_batch = x_batch.to(device)\n", + " y_batch = y_batch.to(device)\n", + "\n", + " optimizer.zero_grad()\n", + " pred = model_data(x_batch)\n", + "\n", + " # Data loss\n", + " data_loss = loss_function(pred, y_batch)\n", + "\n", + " # Physics-informed losses\n", + " physics_loss_total = 0.0\n", + "\n", + " # 1) Vpi*L constant\n", + " if lambda_vpiL != 0:\n", + " Vpi_pred = pred[:,2]\n", + " L_batch = x_batch[:,7]\n", + " residual_vpi = Vpi_pred * L_batch - VpiL_const\n", + " physics_loss_vpiL = torch.mean(residual_vpi**2)\n", + " physics_loss_total += lambda_vpiL * physics_loss_vpiL\n", + "\n", + " # 2) BW monotonic decreasing w.r.t L\n", + " if lambda_bw_mon != 0:\n", + " L_sorted, idx = torch.sort(x_batch[:,7])\n", + " BW_sorted = pred[idx][:,0]\n", + " violation_bw = torch.relu(BW_sorted[1:] - BW_sorted[:-1])\n", + " physics_loss_bw = torch.mean(violation_bw**2)\n", + " physics_loss_total += lambda_bw_mon * physics_loss_bw\n", + "\n", + " # 3) IL monotonic increasing w.r.t L\n", + " if lambda_IL_mon != 0:\n", + " IL_sorted = pred[idx][:,1]\n", + " violation_il = torch.relu(IL_sorted[:-1] - IL_sorted[1:])\n", + " physics_loss_il = torch.mean(violation_il**2)\n", + " physics_loss_total += lambda_IL_mon * physics_loss_il\n", + "\n", + " # 4) Positivity penalty\n", + " if lambda_pos != 0:\n", + " neg_penalty = torch.relu(-pred)\n", + " physics_loss_pos = torch.mean(neg_penalty**2)\n", + " physics_loss_total += lambda_pos * physics_loss_pos\n", + "\n", + " # 5) BW polynomial fit consistency\n", + " if lambda_bw_poly != 0:\n", + " L_batch = x_batch[:,7]\n", + " BW_pred = pred[:,0]\n", + " BW_fit = BW_coeffs[0] * L_batch**2 + BW_coeffs[1] * L_batch + BW_coeffs[2]\n", + " residual_bw_poly = BW_pred - BW_fit\n", + " physics_loss_bw_poly = torch.mean(residual_bw_poly**2)\n", + " physics_loss_total += lambda_bw_poly * physics_loss_bw_poly\n", + "\n", + " # Total loss\n", + " loss = data_loss + physics_loss_total\n", + " loss.backward()\n", + " optimizer.step()\n", + "\n", + " # ==============================\n", + " # Evaluation (no_grad, GPU)\n", + " # ==============================\n", + " model_data.eval()\n", + " with torch.no_grad():\n", + " pred_train = model_data(x_train)\n", + " pred_test = model_data(x_test)\n", + "\n", + " train_loss = loss_function(pred_train, y_train).item()\n", + " test_loss = loss_function(pred_test, y_test).item()\n", + "\n", + " errors_train[epoch] = train_loss\n", + " errors_test[epoch] = test_loss\n", + "\n", + " # Print every 10 epochs + first epoch\n", + " if epoch % 10 == 0 or epoch == 0:\n", + " print(f\"Epoch {epoch:5d}, train MSE {train_loss:.4f}, test MSE {test_loss:.4f}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "P3fpdbKU1J5X" + }, + "source": [ + "## Training Performance Visualization\n", + "\n", + "This section plots the Mean Squared Error (MSE) for both training and testing sets across epochs.\n", + "\n", + "Monitoring both curves allows us to detect:\n", + "- Overfitting (train ↓, test ↑)\n", + "- Underfitting (both high)\n", + "- Proper convergence (both ↓ and stable)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "id": "Z35AHOJlCFRj" + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAHHCAYAAABDUnkqAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAXPhJREFUeJzt3Xl4E9X+BvA3SZuk6b5AS6HQslOQIhTKDkoVgauCIMtFgYqoP0HlVryKXEFxAZeLuIKoiDsCgriCWNnkIvu+ylKKQDeW7k3b5Pz+OCRp6EJLk0ybvp/nmSfNZDI5mQby9nvOnFEJIQSIiIiI3IRa6QYQERERORLDDREREbkVhhsiIiJyKww3RERE5FYYboiIiMitMNwQERGRW2G4ISIiIrfCcENERERuheGGiIiI3ArDDbmFCRMmIDIyUulmEBFRLcBwQ06lUqmqtGzYsEHpptrZsGEDVCoVVqxYoXRTrstoNOLpp59GeHg4vLy8EBcXh3Xr1lX5+efOncPIkSMREBAAPz8/3H333Th16pTdNgUFBZg4cSI6dOgAf39/+Pj4ICYmBm+99RaKi4vL7HPdunXo3bs3DAYDAgMDMWLECCQnJ9ttc/HiRbz++uvo27cvGjRogICAAHTv3h3ffPNNmf0dOnQI9957L5o3bw6DwYCQkBD07dsXP/zwQ5XfZ2n9+/ev0ufy+eefv6H9X+v999/HkiVLqry9SqXClClTHPLazvbxxx+jXbt20Ov1aNWqFd55550qP7cqn93k5ORKf0eTJk2q9j4BoLi4GC+88AKaN28OnU6H5s2b46WXXkJJSUmlbX755ZehUqnQoUOHKr9Pcj0PpRtA7u3zzz+3u//ZZ59h3bp1Zda3a9euRq/z4Ycfwmw212gfddWECROwYsUKTJ06Fa1atcKSJUswePBgrF+/Hr179670ubm5ubjllluQlZWFZ599Fp6ennjzzTfRr18/7N27F8HBwQBkuDl06BAGDx6MyMhIqNVq/O9//8O//vUvbNu2DV999ZV1nz/++CPuvvtudO7cGXPnzkV2djbeeust9O7dG3v27EGDBg0AAFu3bsWMGTMwePBg/Oc//4GHhwe+/fZbjB49GocPH8YLL7xg3eeZM2eQk5OD8ePHIzw8HPn5+fj2229x11134YMPPsBDDz1UrWM2Y8YMPPjgg9b7O3bswNtvv41nn33W7rPYsWPHau23Iu+//z5CQkIwYcIEh+yvtvjggw/wyCOPYPjw4UhMTMTmzZvx+OOPIz8/H08//fR1n1+Vz26DBg3K/H8BAGvWrMGXX36J22+/vdr7BID77rsPy5cvxwMPPIDY2Fj8+eefeO6555CSkoJFixaV296///4br7zyCry9vatzmEgJgsiFJk+eLKryscvLy3NBayq2fv16AUAsX75c0XZcz7Zt2wQA8frrr1vXFRQUiBYtWogePXpc9/mvvvqqACC2b99uXXfkyBGh0WjE9OnTr/v8KVOmCADiwoUL1nXR0dGiZcuWwmg0Wtft3btXqNVqkZiYaF136tQpkZycbLc/s9ksbr31VqHT6URubm6lr11SUiJiYmJEmzZtrtvO61m+fLkAINavX1/jfZWnffv2ol+/flXeHoCYPHmyU9riKPn5+SI4OFgMGTLEbv3YsWOFt7e3uHTpUqXPr+lnd8CAAcLPz08UFBRUe5/bt28XAMRzzz1nt88nn3xSqFQqsW/fvnJfc9SoUeLWW28V/fr1E+3bt79uG0k57JYixfXv3x8dOnTArl270LdvXxgMBjz77LMAgNWrV2PIkCEIDw+HTqdDixYt8OKLL8JkMtnt49oxN5ZS9htvvIFFixahRYsW0Ol06Nq1K3bs2OGwtp86dQr33nsvgoKCYDAY0L17d/z0009ltnvnnXfQvn17azdNbGysXbUjJycHU6dORWRkJHQ6HRo2bIjbbrsNu3fvrvT1V6xYAY1GY1e50Ov1mDhxIrZu3YqzZ89e9/ldu3ZF165drevatm2LAQMGYNmyZdd9/5ZjfuXKFQDApUuXcPjwYQwbNgxarda6XUxMDNq1a4elS5da10VFRaFZs2Z2+1OpVBg6dCiMRmOZrrFraTQaREREWF/bGX755Rf06dMH3t7e8PX1xZAhQ3Do0CG7bVJTU5GQkIAmTZpAp9OhUaNGuPvuu63dcJGRkTh06BA2btxo7Urp379/jduWl5eHJ598EhEREdDpdGjTpg3eeOMNCCHstrN0EQYEBMDHxwdt2rSx/vuyuN7nszzr16/HxYsX8eijj9qtnzx5MvLy8sr9d1BaTT67Fy5cwPr163HPPfdAr9dXe5+bN28GAIwePdpuv6NHj4YQotyu0U2bNmHFihWYP39+pe+Lagd2S1GtcPHiRQwaNAijR4/Gfffdh9DQUADAkiVL4OPjg8TERPj4+OD333/HzJkzkZ2djddff/26+/3qq6+Qk5ODhx9+GCqVCq+99hruuecenDp1Cp6enjVqc1paGnr27In8/Hw8/vjjCA4Oxqeffoq77roLK1aswLBhwwDILrPHH38cI0aMwBNPPIHCwkLs378f27Ztwz//+U8AwCOPPIIVK1ZgypQpiI6OxsWLF/HHH3/gyJEj6Ny5c4Vt2LNnD1q3bg0/Pz+79d26dQMA7N27FxEREeU+12w2Y//+/XjggQfKPNatWzf8+uuvyMnJga+vr3V9UVERsrOzUVBQgJ07d+KNN95As2bN0LJlSwByvAMAeHl5ldmnwWDAoUOHkJqairCwsArfU2pqKgAgJCSkzGN5eXkoKChAVlYWvv/+e/zyyy8YNWpUhfuqic8//xzjx4/HwIED8eqrryI/Px8LFiywdq9Zgt3w4cNx6NAhPPbYY4iMjER6ejrWrVuHlJQUREZGYv78+Xjsscfg4+ODGTNmAID1832jhBC46667sH79ekycOBGdOnXC2rVr8dRTT+HcuXN48803AcixSv/4xz/QsWNHzJ49GzqdDidOnMCWLVus+6rK57M8e/bsAQDExsbare/SpQvUajX27NmD++67r9Ln3+hnd+nSpTCbzRg7duwN7bOiz6nBYAAA7Nq1y269yWTCY489hgcffBA33XRThe+JahGFK0dUz5TXLdWvXz8BQCxcuLDM9vn5+WXWPfzww8JgMIjCwkLruvHjx4tmzZpZ758+fVoAEMHBwXbl8dWrVwsA4ocffqi0nVXplpo6daoAIDZv3mxdl5OTI6KiokRkZKQwmUxCCCHuvvvu65aw/f39b6gbon379uLWW28ts/7QoUMVHlOLjIwMAUDMnj27zGPvvfeeACCOHj1qt/7rr78WAKxLbGys2L9/v/Vxk8kkAgICxIABA+yel5mZKby9vQUAsXPnzgrbdPHiRdGwYUPRp0+fch9/+OGHra+tVqvFiBEjrtv9URXXdkvl5OSIgIAAMWnSJLvtUlNThb+/v3X95cuXy3SDlMfR3VLfffedACBeeuklu/UjRowQKpVKnDhxQgghxJtvvikAiIyMjAr3VZXPZ3kmT54sNBpNuY81aNBAjB49utLn1+Sz26VLF9GoUSPrv7Hq7vPbb78VAMTnn39ut93ChQsFANGhQwe79e+++67w9/cX6enpQgjBbqk6gN1SVCvodDokJCSUWV/6L6ucnBxkZmaiT58+yM/Px9GjR6+731GjRiEwMNB6v0+fPgBw3S6Pqvj555/RrVs3u0GKPj4+eOihh5CcnIzDhw8DAAICAvD3339X2h0WEBCAbdu24fz589VqQ0FBAXQ6XZn1llJ9QUFBpc8FUK3n33LLLVi3bh2WL1+ORx55BJ6ensjLy7M+rlar8fDDDyMpKQnTp0/HX3/9hV27dmHkyJEoKiqqtE2Wv8SvXLlS4Rk3U6dOxbp16/Dpp59i0KBBMJlM1v060rp163DlyhWMGTMGmZmZ1kWj0SAuLg7r168HID+fWq0WGzZswOXLlx3ejor8/PPP0Gg0ePzxx+3WP/nkkxBC4JdffgEgP1eA7N6taMB9VT6f5SkoKLDreixNr9dX+tmzPP9GPrvHjx/Hrl27MHr0aKjV9l9hVd3n4MGD0axZM0ybNg0rV67EmTNnsGzZMsyYMQMeHh52r33x4kXMnDkTzz33nHUwPNV+DDdUKzRu3Ljc/ygPHTqEYcOGwd/fH35+fmjQoIG11J2VlXXd/TZt2tTuviXoOOKL6MyZM2jTpk2Z9Zazbc6cOQMAePrpp+Hj44Nu3bqhVatWmDx5sl23AAC89tprOHjwICIiItCtWzc8//zzVQpgXl5e1hJ7aYWFhdbHK3sugGo9PzQ0FPHx8RgxYgQWLFiAf/zjH7jtttusXUkAMHv2bEycOBGvvfYaWrdujdjYWHh4eGDixIkAZAAsz2OPPYY1a9bgo48+QkxMTLnbtG3bFvHx8Rg3bhx+/PFH5Obm4s477ywzzqSm/vrrLwDArbfeigYNGtgtv/76K9LT0wHIYPjqq6/il19+QWhoKPr27YvXXnvN7ng4w5kzZxAeHm7XZQiU/eyNGjUKvXr1woMPPojQ0FCMHj0ay5Ytsws6Vfl8lsfLy6vCYFlYWFjpZ8/y/Bv57H755ZcAUKZLqjr71Ov1+OmnnxAcHIzhw4cjMjIS48aNw8yZMxEUFGT3Gf3Pf/6DoKAgPPbYY5W+H6pdGG6oVijvP7IrV66gX79+2LdvH2bPno0ffvgB69atw6uvvgoAVTr1W6PRlLve0V+GlWnXrh2OHTuGpUuXonfv3vj222/Ru3dvzJo1y7rNyJEjcerUKbzzzjsIDw/H66+/jvbt21v/Aq9Io0aNcOHChTLrLevCw8MrfG5QUBB0Ot0NPx8ARowYgdzcXKxevdq6TqvV4qOPPsL58+exadMmHDt2DGvXrkVWVhbUarV1fE5pL7zwAt5//33MnTsX999/f6Wvee3r79ixA8ePH6/yc6rC8tn6/PPPsW7dujJL6fc7depUHD9+HHPmzIFer8dzzz2Hdu3aWcekKMnLywubNm3Cb7/9hvvvvx/79+/HqFGjcNttt1kH5Vfl81meRo0awWQyWYOeRVFRES5evHjdz86Nfna/+uortGnTBl26dKnRPtu3b4+DBw/i4MGD2Lx5M86fP49JkyYhMzMTrVu3BiBD7qJFi/D444/j/PnzSE5ORnJyMgoLC1FcXIzk5GRcunSp0vdJClG4W4zqmYrG3JTXf71q1SoBQGzcuNFu/aJFi8qctlvRmJvyxkIAELNmzaq0nVUZc9O6dWvRrVu3Muvnzp0rAIgDBw6U+zyj0SiGDBkiNBqN3WmspaWlpYnGjRuLXr16VdrOadOmCY1GI7KysuzWv/zyywKASElJqfT5sbGxomvXrmXW33bbbaJ58+aVPlcIeYo3APHqq69Wul1JSYlo1KhRuaf4vvvuuwKAmDp16nVf71rz588XAMS2bduq/dzSrh1zs2zZMgFArF27ttr7On78uDAYDGLs2LHWdR06dHDomJuHHnpIaDQakZ2dbbf+zz//FADEO++8U+FzLZ+NdevWlft4VT6fQgjx448/CgDip59+slu/ZcsWAUB89tlnFT5XiBv77FreX3njxG50n6X99NNPAoD44IMPhBC2/wcqW5544olK90nKYOWGai1L1UWUqrIUFRXh/fffV6pJdgYPHozt27dj69at1nV5eXlYtGgRIiMjER0dDUD22Zem1WoRHR0NIQSKi4thMpnKdLE1bNgQ4eHh5ZbYSxsxYgRMJpPdpGNGoxGffPIJ4uLi7M42SUlJKTNOyVL52Llzp3XdsWPH8Pvvv+Pee++1rsvMzCy32vXRRx8BKHvGzLXeeOMNXLhwAU8++aTd+m+++QaPP/44xo4di3nz5lX4/GurA4CcYfazzz6Dl5eX9Vg7ysCBA+Hn54dXXnml3BmYMzIyAAD5+fnWLg+LFi1awNfX1+535+3t7dBT1gcPHgyTyYR3333Xbv2bb74JlUqFQYMGAUC5VYVOnToBsHVHXu/zWZFbb70VQUFBWLBggd36BQsWwGAwYMiQIdZ1mZmZOHr0KPLz863rqvPZtbCcnl7RWVw3sk+LgoICPPfcc2jUqBHGjBkDAOjQoQNWrVpVZmnfvj2aNm2KVatWWbtbqXbhqeBUa/Xs2ROBgYEYP348Hn/8cahUKnz++ecu7VL69ttvyx24PH78eDzzzDP4+uuvMWjQIDz++OMICgrCp59+itOnT+Pbb7+1Dna8/fbbERYWhl69eiE0NBRHjhzBu+++iyFDhsDX1xdXrlxBkyZNMGLECMTExMDHxwe//fYbduzYgf/+97+Vti8uLg733nsvpk+fjvT0dLRs2RKffvopkpOT8fHHH9ttO27cOGzcuNHu+D366KP48MMPMWTIEEybNg2enp6YN28eQkND7YLIF198gYULF2Lo0KFo3rw5cnJysHbtWqxbtw533nknbr31Vrttv/32W/Tt29f6XpYtW4YHH3wQw4cPt263fft2jBs3DsHBwRgwYIB1LIVFz5490bx5cwDAww8/jOzsbPTt2xeNGzdGamoqvvzySxw9ehT//e9/7cZILFmyBAkJCfjkk09ueEZgPz8/LFiwAPfffz86d+6M0aNHo0GDBkhJScFPP/2EXr164d1338Xx48cxYMAAjBw5EtHR0fDw8MCqVauQlpZmN4dKly5dsGDBArz00kto2bIlGjZsaHfMyrNz50689NJLZdb3798fd955J2655RbMmDEDycnJiImJwa+//orVq1dj6tSpaNGiBQA5/mnTpk0YMmQImjVrhvT0dLz//vto0qSJdSD89T6fFfHy8sKLL76IyZMn495778XAgQOxefNmfPHFF3j55ZcRFBRk3fbdd9/FCy+8gPXr11vn+KnOZxeQp2N/88036N69u/X9Xas6+xw5ciTCw8MRHR2N7OxsLF68GKdOncJPP/1kfd8hISEYOnRomdexzHVT3mNUSyhZNqL6pzrdUkLIEnf37t2Fl5eXCA8PF//+97/F2rVrXdYtVdFiOf375MmTYsSIESIgIEDo9XrRrVs38eOPP9rt64MPPhB9+/YVwcHBQqfTiRYtWoinnnrKWjo3Go3iqaeeEjExMcLX11d4e3uLmJgY8f7771faRouCggIxbdo0ERYWJnQ6nejatatYs2ZNme0sp9xf6+zZs2LEiBHCz89P+Pj4iH/84x/ir7/+sttmx44d4t577xVNmzYVOp1OeHt7i86dO4t58+aJ4uJiu223bdsm+vbtKwIDA4VerxcxMTFi4cKFwmw22233ySefVHqMP/nkE+u2X3/9tYiPjxehoaHCw8NDBAYGivj4eLF69eoy7+edd94RAMo9BhWpaIbi9evXi4EDBwp/f3+h1+tFixYtxIQJE6yns2dmZorJkyeLtm3bCm9vb+Hv7y/i4uLEsmXL7PaTmpoqhgwZInx9fQWA63ZRVXZcXnzxRSGEPF39X//6lwgPDxeenp6iVatW4vXXX7c7zklJSeLuu+8W4eHhQqvVivDwcDFmzBhx/Phx6zbX+3xez6JFi0SbNm2EVqsVLVq0EG+++WaZ3/WsWbPKPb5V/ewKIcSaNWsEAPH2229X2p6q7vPVV18Vbdu2FXq9XgQGBoq77rpL7Nmzp0rvmaeC134qIVz4ZzARkZONHDkSycnJ2L59u9JNISKFsFuKiNyGEAIbNmzAF198oXRTiEhBrNwQERGRW+HZUkRERORWGG6IiIjIrTDcEBERkVthuCEiIiK3Uu/OljKbzTh//jx8fX2hUqmUbg4RERFVgRACOTk5CA8PL3NF+GvVu3Bz/vz5SqfgJiIiotrr7NmzaNKkSaXb1LtwY5lW++zZs/Dz81O4NURERFQV2dnZiIiIqPSyIBb1LtxYuqL8/PwYboiIiOqYqgwp4YBiIiIicisMN0RERORWGG6IiIjIrdS7MTdERETOZDKZUFxcrHQz6iStVnvd07yrguGGiIjIAYQQSE1NxZUrV5RuSp2lVqsRFRUFrVZbo/0w3BARETmAJdg0bNgQBoOBE8VWk2WS3QsXLqBp06Y1On4MN0RERDVkMpmswSY4OFjp5tRZDRo0wPnz51FSUgJPT88b3g8HFBMREdWQZYyNwWBQuCV1m6U7ymQy1Wg/DDdEREQOwq6omnHU8WO4ISIiIrfCcENEREQOERkZifnz5yvdDA4oJiIiqs/69++PTp06OSSU7NixA97e3jVvVA0x3DhIURGQlgYIATRtqnRriIiIHEMIAZPJBA+P60eGBg0auKBF18duKQfZvl2Gmvh4pVtCRERUNRMmTMDGjRvx1ltvQaVSQaVSYcmSJVCpVPjll1/QpUsX6HQ6/PHHHzh58iTuvvtuhIaGwsfHB127dsVvv/1mt79ru6VUKhU++ugjDBs2DAaDAa1atcL333/v9PfFcOMgOp28NRqVbQcREdUCQgB5ecosQlS5mW+99RZ69OiBSZMm4cKFC7hw4QIiIiIAAM888wzmzp2LI0eOoGPHjsjNzcXgwYORlJSEPXv24I477sCdd96JlJSUSl/jhRdewMiRI7F//34MHjwYY8eOxaVLl2p0eK+H3VIOotfL28JCZdtBRES1QH4+4OOjzGvn5gJVHPfi7+8PrVYLg8GAsLAwAMDRo0cBALNnz8Ztt91m3TYoKAgxMTHW+y+++CJWrVqF77//HlOmTKnwNSZMmIAxY8YAAF555RW8/fbb2L59O+64445qv7WqYuXGQVi5ISIidxIbG2t3Pzc3F9OmTUO7du0QEBAAHx8fHDly5LqVm44dO1p/9vb2hp+fH9LT053SZgtWbhyElRsiIrIyGGQFRanXdoBrz3qaNm0a1q1bhzfeeAMtW7aEl5cXRowYgaKiokr3c+1lFFQqFcxms0PaWBGGGwcpXbkRAuAklURE9ZhKVeWuIaVptdoqXe5gy5YtmDBhAoYNGwZAVnKSk5Od3Lobw24pB7FUbgB5WjgREVFdEBkZiW3btiE5ORmZmZkVVlVatWqFlStXYu/evdi3bx/++c9/Or0Cc6MYbhxEl3nO+jPH3RARUV0xbdo0aDQaREdHo0GDBhWOoZk3bx4CAwPRs2dP3HnnnRg4cCA6d+7s4tZWjUqIapwz5gays7Ph7++PrKws+Pn5OWy/Ysv/oO7dE4CczK9hQ4ftmoiIarnCwkKcPn0aUVFR0Jcu5VO1VHYcq/P9zcqNg6j0OmghSzas3BARESmH4cZR9HroIU+V4hlTREREymG4cRSdDjpWboiIiBTHcOMorNwQERHVCgw3jlK6cpN//fkCiIiIyDkYbhyldOUmp1jhxhAREdVfDDeOUrpyk8dwQ0REpBSGG0fx9CxVuSlRuDFERET1F8ONo6hU0KllxcaYx3BDRESkFIYbB9JfDTeFuQw3RERESmG4cSCd5mrlhmdLERFRHdG/f39MnTrVYfubMGEChg4d6rD93QiGGwfSa2TFpjCP4YaIiEgpDDcOpPOQ4cZYwHBDRES134QJE7Bx40a89dZbUKlUUKlUSE5OxsGDBzFo0CD4+PggNDQU999/PzIzM63PW7FiBW666SZ4eXkhODgY8fHxyMvLw/PPP49PP/0Uq1evtu5vw4YNLn9fHi5/RTem95ChpjC/Xl1onYiIriEEkJ+vzGsbDIBKVbVt33rrLRw/fhwdOnTA7NmzAQCenp7o1q0bHnzwQbz55psoKCjA008/jZEjR+L333/HhQsXMGbMGLz22msYNmwYcnJysHnzZgghMG3aNBw5cgTZ2dn45JNPAABBQUHOeqsVYrhxIJ2nDDfGArPCLSEiIiXl5wM+Psq8dm4u4O1dtW39/f2h1WphMBgQFhYGAHjppZdw880345VXXrFut3jxYkREROD48ePIzc1FSUkJ7rnnHjRr1gwAcNNNN1m39fLygtFotO5PCQw3DqS/Gm4KC1i5ISKiumnfvn1Yv349fMpJZydPnsTtt9+OAQMG4KabbsLAgQNx++23Y8SIEQgMDFSgteVjuHEgnacMNcZChhsiovrMYJAVFKVeuyZyc3Nx55134tVXXy3zWKNGjaDRaLBu3Tr873//w6+//op33nkHM2bMwLZt2xAVFVWzF3cQhhsH0mtld1ShUeGGEBGRolSqqncNKU2r1cJksp0I07lzZ3z77beIjIyEh0f5MUGlUqFXr17o1asXZs6ciWbNmmHVqlVITEwssz8l8GwpB9Jpr1ZuGG6IiKiOiIyMxLZt25CcnIzMzExMnjwZly5dwpgxY7Bjxw6cPHkSa9euRUJCAkwmE7Zt24ZXXnkFO3fuREpKClauXImMjAy0a9fOur/9+/fj2LFjyMzMRHGx66+3yHDjQHqdDDeFxioOUyciIlLYtGnToNFoEB0djQYNGqCoqAhbtmyByWTC7bffjptuuglTp05FQEAA1Go1/Pz8sGnTJgwePBitW7fGf/7zH/z3v//FoEGDAACTJk1CmzZtEBsbiwYNGmDLli0uf0/slnIgnV7eGosYboiIqG5o3bo1tm7dWmb9ypUry92+Xbt2WLNmTYX7a9CgAX799VeHte9GsHLjQHqdvC008rASEREphd/CDqTTy4qNsZiHlYiISCn8FnYgvZcMN4UMN0RERIrht7ADWSs3JRqFW0JERFR/Mdw4kN4gD2dhMcMNEVF9JAQnca0JRx0/hhsH0nnJw2k08SQ0IqL6xNPTEwCQr9TVMt1EUVERAECjqVmRgN/CDqT3lr+MwhIeViKi+kSj0SAgIADp6ekAAIPBAFVVL81NAACz2YyMjAwYDIYKZ0auKn4LO5DuargxmjwVbgkREbma5SrYloBD1adWq9G0adMaB0OGGwfSe8vDWWhmuCEiqm9UKhUaNWqEhg0bKnLJAXeg1WqhVtd8xAzDjQPproYbo9kTQsgLpxERUf2i0WhqPGaEaqZWDCh+7733EBkZCb1ej7i4OGzfvr3CbZcsWQKVSmW36PV6F7a2YnofGW4E1GBoJyIiUobi4eabb75BYmIiZs2ahd27dyMmJgYDBw6stM/Sz88PFy5csC5nzpxxYYsrpvPVWn/mlcGJiIiUoXi4mTdvHiZNmoSEhARER0dj4cKFMBgMWLx4cYXPUalUCAsLsy6hoaEubHHFdD62sTaFhQo2hIiIqB5TNNwUFRVh165diI+Pt65Tq9WIj48v9wqlFrm5uWjWrBkiIiJw991349ChQxVuazQakZ2dbbc4i9qghyeKrr6u016GiIiIKqFouMnMzITJZCpTeQkNDUVqamq5z2nTpg0WL16M1atX44svvoDZbEbPnj3x999/l7v9nDlz4O/vb10iIiIc/j6sdDroIUs2rNwQEREpQ/Fuqerq0aMHxo0bh06dOqFfv35YuXIlGjRogA8++KDc7adPn46srCzrcvbsWec1Tq+HDrJkw8oNERGRMhQ9FTwkJAQajQZpaWl269PS0qyTIV2Pp6cnbr75Zpw4caLcx3U6HXQ6XY3bWiWs3BARESlO0cqNVqtFly5dkJSUZF1nNpuRlJSEHj16VGkfJpMJBw4cQKNGjZzVzKorXbkp5MXTiIiIlKD4JH6JiYkYP348YmNj0a1bN8yfPx95eXlISEgAAIwbNw6NGzfGnDlzAACzZ89G9+7d0bJlS1y5cgWvv/46zpw5gwcffFDJtyHpdNAjAwBQmFMMQFv59kRERORwioebUaNGISMjAzNnzkRqaio6deqENWvWWAcZp6Sk2E3FfPnyZUyaNAmpqakIDAxEly5d8L///Q/R0dFKvQWb0pWbvBIw3BAREbmeSghRr/pPsrOz4e/vj6ysLPj5+Tl252Yzemm24n/ohZWfXMGwCQGO3T8REVE9VZ3v7zp3tlStplZDp7o6z01eicKNISIiqp8YbhxMr5YXlSrMZbghIiJSAsONg+k0MtwY80wKt4SIiKh+YrhxML1GVmwKGW6IiIgUwXDjYDoPGW6MBWaFW0JERFQ/Mdw4mN5DVmwK8xluiIiIlMBw42C6q+GGlRsiIiJlMNw4mN7zauWmoF5NH0RERFRrMNw4mM5TVmx4bSkiIiJlMNw4mF4rww2vCk5ERKQMhhsH02llxabQqHBDiIiI6imGGwfT62S4MRpVCreEiIiofmK4cTCdTt4WMtwQEREpguHGwfR6eWssZrghIiJSAsONg+n0MtQUFvHQEhERKYHfwA5mq9zw0BIRESmB38AOpvOSh7SwWKNwS4iIiOonhhsH03vJbiljCcMNERGREhhuHExnkKGmsMRD4ZYQERHVTww3DqY3yENqNDHcEBERKYHhxsGslRuTp8ItISIiqp8YbhxM7yMrNqzcEBERKYPhxsEs4abQrIXghcGJiIhcjuHGwXQ+sjtKQI2SEoUbQ0REVA8x3DiYpXIDAIWFCjaEiIionmK4cTCdr9b6s9GoYEOIiIjqKYYbB1N76eCJIgCs3BARESmB4cbR9HroIEs2rNwQERG5HsONo+l00EOWbFi5ISIicj2GG0dj5YaIiEhRDDeOxsoNERGRohhuHE2ns1Vu8k0KN4aIiKj+YbhxNL3eVrnJLlK4MURERPUPw42jla7c5BUr3BgiIqL6h+HG0Tw8oL8abgpzef0FIiIiV2O4cQKdWlZsjAw3RERELsdw4wR6jQw3hXkcUExERORqDDdOoNPIig3PliIiInI9hhsn0Htcrdww3BAREbkcw40T6DxkqDHmmxVuCRERUf3DcOME+qvhprCA4YaIiMjVGG6cQOcpQ42xQCjcEiIiovqH4cYJ9J5XKzeFDDdERESuxnDjBDqtDDVGXjiTiIjI5RhunECvld1SvCo4ERGR6zHcOIFOd7Vyw+tmEhERuRzDjRPor3ZLFRpVCreEiIio/mG4cQKdXoYaYxHDDRERkasx3DiBXi9vC4t4eImIiFyN375OYK3cFPPwEhERuRq/fZ1A7yXDTWGxRuGWEBER1T8MN06g85KH1VjCcENERORqtSLcvPfee4iMjIRer0dcXBy2b99epectXboUKpUKQ4cOdW4Dq0lvkIe1sMRD4ZYQERHVP4qHm2+++QaJiYmYNWsWdu/ejZiYGAwcOBDp6emVPi85ORnTpk1Dnz59XNTSqmPlhoiISDmKh5t58+Zh0qRJSEhIQHR0NBYuXAiDwYDFixdX+ByTyYSxY8fihRdeQPPmzV3Y2qrR+8iKTaHJU+GWEBER1T+KhpuioiLs2rUL8fHx1nVqtRrx8fHYunVrhc+bPXs2GjZsiIkTJ7qimdWmM8iKjdHMcENERORqig4KyczMhMlkQmhoqN360NBQHD16tNzn/PHHH/j444+xd+/eKr2G0WiE0Wi03s/Ozr7h9lYVKzdERETKUbxbqjpycnJw//3348MPP0RISEiVnjNnzhz4+/tbl4iICCe3EtB5y3BjhgYlJU5/OSIiIipF0cpNSEgINBoN0tLS7NanpaUhLCyszPYnT55EcnIy7rzzTus6s1legdvDwwPHjh1DixYt7J4zffp0JCYmWu9nZ2c7PeDofW0Vm8JCwMfHqS9HREREpShaudFqtejSpQuSkpKs68xmM5KSktCjR48y27dt2xYHDhzA3r17rctdd92FW265BXv37i03tOh0Ovj5+dktzqbzsYWbUj1iRERE5AKKT8SSmJiI8ePHIzY2Ft26dcP8+fORl5eHhIQEAMC4cePQuHFjzJkzB3q9Hh06dLB7fkBAAACUWa8kjbceHihGCTxRWKh0a4iIiOoXxcPNqFGjkJGRgZkzZyI1NRWdOnXCmjVrrIOMU1JSoFbXqaFBgE4HHYwogScrN0RERC6mEkIIpRvhStnZ2fD390dWVpbzuqgOHkTITWG4iBAcOgRERzvnZYiIiOqL6nx/17GSSB2h10MHWbJh5YaIiMi1GG6cQaeDHnKwDcfcEBERuRbDjTOUrtwUmBVuDBERUf3CcOMMpSs3OcUKN4aIiKh+YbhxhtKVmzyGGyIiIldiuHEGT89SlRtef4GIiMiVGG6cQaWCTi0rNsY8hhsiIiJXYrhxEv3VcFOYy3BDRETkSgw3TqLTyFBjzDcp3BIiIqL6heHGSfQeVys3+TwVnIiIyJUYbpxEp5EVG1ZuiIiIXIvhxkn0nrJbqrCgXl26i4iISHEMN06i85DdUZyhmIiIyLUYbpxEr5XdUYWFrNwQERG5EsONk+g8ZagxMtwQERG5FMONk+i1sjuqsFClcEuIiIjqF4YbJ9HprlZujAo3hIiIqJ5huHESvU7eFhpZuSEiInIlhhsn0V0NN8YihhsiIiJXYrhxEv3VbqnCIh5iIiIiV+I3r5PovOShNRbzEBMREbkSv3mdRK+Xt4UMN0RERC7Fb14n0Rk0AABjiUbhlhAREdUvDDdOoveSA4kLiz0UbgkREVH9wnDjJJbKTaGJ4YaIiMiVGG6cRO99tVuK4YaIiMilGG6cROctQ02hyVPhlhAREdUvDDdOYq3cmBluiIiIXInhxkkslRuT0KCkROHGEBER1SMMN06i97VVbHjxTCIiItdhuHESnY8t3BQWKtgQIiKieobhxkk8vHXQQPZHsXJDRETkOgw3zqLTQQeZali5ISIich2GG2fR66GHTDWs3BAREbkOw42zsHJDRESkCIYbZ2HlhoiISBEMN85SunJTIBRuDBERUf3BcOMspSs3eZzFj4iIyFUYbpyldOUml+GGiIjIVRhunEWns1VucosVbgwREVH9wXDjLGo1dKoiAKzcEBERuRLDjRPp1bJiwzE3RERErsNw40QGjRxzk5djUrglRERE9QfDjRMFeuYAAC5f4qngRERErlKtcPPaa6+hoKDAen/Lli0wlpqhLicnB48++qjjWlfHBWtluLl4UaVwS4iIiOqPaoWb6dOnIycnx3p/0KBBOHfunPV+fn4+PvjgA8e1ro4L0uYBAC5dYYGMiIjIVar1rSuEqPQ+2QvS5wMALmUx3BAREbkKv3WdyBJuLmZ5KNwSIiKi+oPhxomCveUkfpeyPRVuCRERUf1R7ZLCRx99BB8fHwBASUkJlixZgpCQEACwG49DQJDharjJZbghIiJylWqFm6ZNm+LDDz+03g8LC8Pnn39eZhuSgnzkDMVX8nUwmQCNRuEGERER1QPVCjfJyclOaoZ7Cgq0Dbi+fBm4WuAiIiIiJ+KYGyfyCAmAH7IAAJcuKdwYIiKieqJa4Wbr1q348ccf7dZ99tlniIqKQsOGDfHQQw/ZTepXVe+99x4iIyOh1+sRFxeH7du3V7jtypUrERsbi4CAAHh7e6NTp05lusZqjaAgBEGmGoYbIiIi16hWuJk9ezYOHTpkvX/gwAFMnDgR8fHxeOaZZ/DDDz9gzpw51WrAN998g8TERMyaNQu7d+9GTEwMBg4ciPT09HK3DwoKwowZM7B161bs378fCQkJSEhIwNq1a6v1ui5RKtxcvKhwW4iIiOqJaoWbvXv3YsCAAdb7S5cuRVxcHD788EMkJibi7bffxrJly6rVgHnz5mHSpElISEhAdHQ0Fi5cCIPBgMWLF5e7ff/+/TFs2DC0a9cOLVq0wBNPPIGOHTvijz/+qNbrukRQEIIhUw0rN0RERK5RrXBz+fJlhIaGWu9v3LgRgwYNst7v2rUrzp49W+X9FRUVYdeuXYiPj7c1SK1GfHw8tm7det3nCyGQlJSEY8eOoW/fvuVuYzQakZ2dbbe4DLuliIiIXK5a4SY0NBSnT58GIIPJ7t270b17d+vjOTk58PSs+pwumZmZMJlMdoHJ8jqpqakVPi8rKws+Pj7QarUYMmQI3nnnHdx2223lbjtnzhz4+/tbl4iIiCq3r8aCgxluiIiIXKxa4Wbw4MF45plnsHnzZkyfPh0GgwF9+vSxPr5//360aNHC4Y28lq+vL/bu3YsdO3bg5ZdfRmJiIjZs2FDuttOnT0dWVpZ1qU5lqcZKdUtxzA0REZFrVGuemxdffBH33HMP+vXrBx8fHyxZsgRardb6+OLFi3H77bdXeX8hISHQaDRIS0uzW5+WloawsLAKn6dWq9GyZUsAQKdOnXDkyBHMmTMH/fv3L7OtTqeDTqercpscqnS3VKYZPPOeiIjI+aoVbkJCQrBp0yZrt5Dmmil3ly9fDl9f3yrvT6vVokuXLkhKSsLQoUMBAGazGUlJSZgyZUqV92M2m2/oFHSnCwiwhZv0YgAKhSwiIqJ6pFrh5oEHHqjSdhWd6VSexMREjB8/HrGxsejWrRvmz5+PvLw8JCQkAADGjRuHxo0bW08xnzNnDmJjY9GiRQsYjUb8/PPP+Pzzz7FgwYLqvBXX0GgQ5F0E5AEXM8xKt4aIiKheqFa4WbJkCZo1a4abb74ZQojrP6EKRo0ahYyMDMycOROpqano1KkT1qxZYx1knJKSArXa1p2Tl5eHRx99FH///Te8vLzQtm1bfPHFFxg1apRD2uNowf4lQB5w6ZJK6aYQERHVCypRjZQyefJkfP3112jWrBkSEhJw3333ISgoyJntc7js7Gz4+/sjKysLfn5+Tn+9ox1GoN2hFQjwLsLlXO31n0BERERlVOf7u1ojXN977z1cuHAB//73v/HDDz8gIiICI0eOxNq1ax1WyXE3QSHyEF/J08JkUrgxRERE9UC1T9/R6XQYM2YM1q1bh8OHD6N9+/Z49NFHERkZidzcXGe0sU4LCrNVay5fVrAhRERE9USNzk1Wq9VQqVQQQsDEskS5eGVwIiIi16p2uDEajfj6669x2223oXXr1jhw4ADeffddpKSkwMfHxxltrNt4CQYiIiKXqtbZUo8++iiWLl2KiIgIPPDAA/j6668REhLirLa5h6vhJhlRnKWYiIjIBaoVbhYuXIimTZuiefPm2LhxIzZu3FjuditXrnRI49wCrwxORETkUtUKN+PGjYNKxflaqoXdUkRERC5V7Un8qJqCghCEFAAMN0RERK7AKzk6G68MTkRE5FIMN84WHGzrlrrI60sRERE5G8ONswUG2sJNWonCjSEiInJ/DDfO5uGBIK9CALwyOBERkSsw3LhAsL+s2Fzi5ReIiIicjuHGBSwXTr90RaNsQ4iIiOoBhhsXsF0Z3JNXBiciInIyhhsXCGpom06IVwYnIiJyLoYbF+CVwYmIiFyH4cYVSs91w3BDRETkVAw3rlDq+lKcpZiIiMi5GG5cgVcGJyIichmGG1fglcGJiIhchuHGFRhuiIiIXIbhxhV4ZXAiIiKXYbhxhdKVm4tC4cYQERG5N4YbVygdbtJ5ZXAiIiJnYrhxBa0WQfoCAMDFDF5/gYiIyJkYblzEemVwDigmIiJyKoYbFwkKlGNtLmXxyuBERETOxHDjItYrg+fyyuBERETOxHDjIrwyOBERkWsw3LgIrwxORETkGgw3rsJZiomIiFyC4cZVeGVwIiIil2C4cZXgYF4ZnIiIyAUYblyF3VJEREQuwXDjKgw3RERELsFw4yq8MjgREZFLMNy4Cq8MTkRE5BIMN64SGGgLNxm8MjgREZGzMNy4ipcXgrR5AICL6WaFG0NEROS+GG5cKNivGAAHFBMRETkTw40LBQXJ24tXeNiJiIichd+yLhTaQHZHZeV5Ii9P4cYQERG5KYYbFwoM1cIfVwAAycmKNoWIiMhtMdy4UlAQmuMUAOD0aYXbQkRE5KYYblwpKAhRkKnm1CmF20JEROSmGG5ciZUbIiIip2O4cSVWboiIiJyO4caVWLkhIiJyOoYbV7qmciN4iSkiIiKHY7hxpeBgNMMZqGBGXh6Qmal0g4iIiNwPw40rBQdDDyPCcR4Ax90QERE5A8ONK4WFAZ6eHHdDRETkRLUi3Lz33nuIjIyEXq9HXFwctm/fXuG2H374Ifr06YPAwEAEBgYiPj6+0u1rFY0GaNaMZ0wRERE5keLh5ptvvkFiYiJmzZqF3bt3IyYmBgMHDkR6enq522/YsAFjxozB+vXrsXXrVkREROD222/HuXPnXNzyGxQZycoNERGREykebubNm4dJkyYhISEB0dHRWLhwIQwGAxYvXlzu9l9++SUeffRRdOrUCW3btsVHH30Es9mMpKQkF7f8BkVFsXJDRETkRIqGm6KiIuzatQvx8fHWdWq1GvHx8di6dWuV9pGfn4/i4mIEBQWV+7jRaER2drbdoqioKFZuiIiInEjRcJOZmQmTyYTQ0FC79aGhoUhNTa3SPp5++mmEh4fbBaTS5syZA39/f+sSERFR43bXSKnKTUoKUFKibHOIiIjcjeLdUjUxd+5cLF26FKtWrYJery93m+nTpyMrK8u6nD171sWtvEZUFBrhAnQohMkEKN0cIiIid+Oh5IuHhIRAo9EgLS3Nbn1aWhrCwsIqfe4bb7yBuXPn4rfffkPHjh0r3E6n00Gn0zmkvQ4RFQU1BCKRjGNoi1OngKgopRtFRETkPhSt3Gi1WnTp0sVuMLBlcHCPHj0qfN5rr72GF198EWvWrEFsbKwrmuo4DRoABgPH3RARETmJopUbAEhMTMT48eMRGxuLbt26Yf78+cjLy0NCQgIAYNy4cWjcuDHmzJkDAHj11Vcxc+ZMfPXVV4iMjLSOzfHx8YGPj49i76PKVCogMhJRh3nGFBERkTMoHm5GjRqFjIwMzJw5E6mpqejUqRPWrFljHWSckpICtdpWYFqwYAGKioowYsQIu/3MmjULzz//vCubfuOiotD8MCs3REREzqB4uAGAKVOmYMqUKeU+tmHDBrv7ycnJzm+Qs5U6Y4rhhoiIyLHq9NlSdRYn8iMiInIahhsllJrILyMDyM1VuD1ERERuhOFGCVFR8Ec2AlWXAbBrioiIyJEYbpQQGQkAaC5OAmC4ISIiciSGGyUEBAABARx3Q0RE5AQMN0rhBTSJiIicguFGKTxjioiIyCkYbpTCyg0REZFTMNwo5ZqJ/IRQuD1ERERuguFGKVFRaIYzUMGM/HwgPV3pBhEREbkHhhulREVBi2I0UZ0DwHE3REREjsJwo5RmzQBwrhsiIiJHY7hRisEAhIZax92cPKlwe4iIiNwEw42SoqIQjcMAgL17lW0KERGRu2C4UVJUFOKwDQCwbZvCbSEiInITDDdKiopCF+yCWmXGuXPAuXNKN4iIiKjuY7hRUlQUvJGPDj7JAIDt25VtDhERkTtguFFSVBQAIE69EwC7poiIiByB4UZJkZEAgLi83wEw3BARETkCw42SmjYF1Gp0K9kCANi5EzCZFG4TERFRHcdwoyRPT6BJE0TjMHy8TMjNBY4cUbpRREREdRvDjdKioqCBGbGRGQDYNUVERFRTDDdKuzqouFuwvLgUz5giIiKqGYYbpbVrBwCIM8lxN6zcEBER1QzDjdK6dAEAxP29EgBw4ACQl6dkg4iIiOo2hhulde4MAGh89k+Eh5lhNgO7dyvcJiIiojqM4UZpgYFAixYAgLgWHFRMRERUUww3tUFsLAAgzldeIZzhhoiI6MYx3NQGV8fddMvbAIBnTBEREdUEw01tcDXcxJ75FioVkJICpKYq3CYiIqI6iuGmNrg6qNg35RDaty0BwK4pIiKiG8VwUxsEBAAtWwIAujWVJRt2TREREd0YhpvawjKo2OsAAFZuiIiIbhTDTW1hmcwvex0AYMcOwGxWskFERER1E8NNbXG1ctP+xGr4+gLZ2cDWrQq3iYiIqA5iuKktbr4ZAOCRcgpDBxUCAJYuVbJBREREdRPDTW3h7w+0bg0AGN3hEABg2TKgpETJRhEREdU9DDe1ydVxN/HmXxEUBKSnAxs3KtwmIiKiOobhpja5Gm60+3Zg+HC5il1TRERE1cNwU5tcHVSMXbswerT8ceVKoKhIuSYRERHVNQw3tcnVQcVISUG/6AyEhgKXLgG//aZss4iIiOoShpvaxM8PaNMGAKDZuwsjR8rV7JoiIiKqOoab2ubquBvs3GntmvruO6CgQLEWERER1SkMN7VNqXE33bsDTZsCOTnAL78o2ywiIqK6guGmtilVuVGrgVGj5F12TREREVUNw01tc/PNgEoF/P03kJxs7Zr68UdZwSEiIqLKMdzUNr6+QP/+8uelS3HzzUCrVnLMzfffK9oyIiKiOoHhpjYaO1befvklVCpYqzf//jdw+rRyzSIiIqoLGG5qo+HDAa0WOHgQ2L8f//oX0L49cP48EB8PXLigdAOJiIhqL4ab2iggABgyRP781VcIDAR+/RVo3hw4dQq4/XY5uR8RERGVxXBTW/3zn/L2q68Asxnh4XKm4kaNZEFn0CAOMCYiIioPw01tNWSInLH47FlgyxYAQFQUsG4dEBQEbN8ODB0KFBYq20wiIqLaRvFw89577yEyMhJ6vR5xcXHYvn17hdseOnQIw4cPR2RkJFQqFebPn++6hrqalxdwzz3y5y+/tK5u3x5Yswbw8QF+/x0YMwYoKVGojURERLWQouHmm2++QWJiImbNmoXdu3cjJiYGAwcORHp6ernb5+fno3nz5pg7dy7CwsJc3FoFWM6aWr7c7tLgXbvK08J1OnlphkmTALNZmSYSERHVNoqGm3nz5mHSpElISEhAdHQ0Fi5cCIPBgMWLF5e7fdeuXfH6669j9OjR0Ol0Lm6tAm65BQgLk6OH164t89DSpYBaDSxZAjz1FCCEMs0kIiKqTRQLN0VFRdi1axfi4+NtjVGrER8fj61btzrsdYxGI7Kzs+2WOkOjsU1y89VXZR4eOhT4+GP587x5wJw5rmsaERFRbaVYuMnMzITJZEJoaKjd+tDQUKSmpjrsdebMmQN/f3/rEhER4bB9u4TlrKnVq8s9PWrCBOC//5U/z5ghf2YFh4iI6jPFBxQ72/Tp05GVlWVdzp49q3STqic21nb9he++K3eTxETg2Wflz9OmyWIPTxMnIqL6SrFwExISAo1Gg7S0NLv1aWlpDh0srNPp4OfnZ7fUKSqVbWDxv/8NJCeXu9lLLwFvvgl4eADLlslBxwcPuq6ZREREtYVi4Uar1aJLly5ISkqyrjObzUhKSkKPHj2Ualbt9K9/ATfdBKSmAnfcAVy8WGYTlQqYOhXYtAlo0gQ4dgzo1g348EMgLY1dVUREVH8o2i2VmJiIDz/8EJ9++imOHDmC//u//0NeXh4SEhIAAOPGjcP06dOt2xcVFWHv3r3Yu3cvioqKcO7cOezduxcnTpxQ6i24hp8f8MsvQESETC133SW7qcrRowewe7e8RENBAfDQQ/KEq5AQoE8fef+aM8uJiIjcikoIZf+mf/fdd/H6668jNTUVnTp1wttvv424uDgAQP/+/REZGYklS5YAAJKTkxEVFVVmH/369cOGDRuq9HrZ2dnw9/dHVlZW3euiOnQI6N0buHIFGDZMphSNptxNTSbgtdfk2VSnTpWt3ISGyvlxHnpIZiYiIqLarDrf34qHG1er0+EGADZulGWZoiJg8mTgnXdkn1QlCgpkwefIEVnV+eIL2cMFyHly7rwTePhhudsKshIREZGiGG4qUefDDSBHDI8eLcsxr74qBxpXQ3GxPPFqwQJg/Xrb+qZNgYkTgQcekON2iIiIaguGm0q4RbgBgPnz5UBjQF57yjIfTjUdOSIHHX/6qZwIGZDVnNtuk5e2uvtu2YVlYTIBe/cCGzYAWi0wahTQsGFN3ggREdH1MdxUwm3CDQA8+aScmtjTUw44HjDghndVWAisXAksWiR7vixUKjnMp39/YP9++diVK7bHPTxkt9bEicDAgfI+ERGRozHcVMKtwo3ZLC8LvmyZPKNq82agY8ca7/avv4AVK2TY2bmz7ON+fkDfvkB6OlD6Iu7h4fLlw8JsS7NmQJcuspvrOkODiIiIKsRwUwm3CjeALLkMHCgnuAkPB7ZulYNnHCQlBVi1CtixQ061c+utwM032yo0Bw/KM7I++8zWrVWehg3lZMtdusihQunpcv6d9HSgpASIiZGPdekiA1J2NrBnjxwAvXu3nNpn7FjgvvtkdxgREdUvDDeVcLtwAwCXL8tJbA4dkgNkliyRk/25kNEo89Xff8szsSzL8ePAgQNyrE5VaTQVb9+kieyNmzQJ8PZ2TNuJiKj2Y7iphFuGGwA4exYYNEgGHAB44glg7lxAr1e2XZCnou/bJ7u49u2TlZfQUFnNCQ2VvWu7dwO7dsltLBMwt2oFdO4sK0Vmszzr/cIF+VhwsJywMDtbjgG6cgXIzweCguQ+LfuPiADatJFLy5bytS9elAWu//1PLpcuAVFRQPPmQIsW8jY4GPDxsS3e3oBOx641IiKlMNxUwm3DDSBTxNNPyxQAAB06AF99JfuT6gghgHPn5Liea389hYWy++u114CTJ6u/b7Vahh5LQLqR53t7AwaDvG3USAawli3lEhUFNGggA5afX8VBSAjg/HmZQw8flvsbPJin3xMRVYbhphJuHW4sfv4ZSEiQA1o0Gnle96hRwNChQECA0q2rMZMJ+OknOWYnIADw95eLwSCrMJaxPGlpwOnTcgLDY8fsr5Teti3Qs6dcwsPl9UhPnpSzOZ8+LStBublyKSysfhs9PGTI8fGRFR/LYjbLrrqsrLLP6dxZnnk2aJB8/qVLssp08aKsOMXGypzKM9KIqD5iuKlEvQg3gPx2nzQJ+P572zqtVg4+HjdOTmDj6alc+1xMCFmxSUmRVZaQkKo/t6QEyMuT3V6W29xc2RN44oQ8u+zECRmQLl2q8LJfdjQa2Y7oaBnCtm6t2sVNvbzkoOuuXeXPBQW2xWwGfH3l4ucnb4uKZNddVpa8BWTWHTKkbGUMkKFu715ZWUpLs42dysqS+zIa5WI2y4HfffrI5UbGsJeUyHmWvLyAyEiGNiKqHMNNJepNuLE4fhz45hu5WMbjALJc8cgjMgCFhSnXPjdUUGCruOTn2wKBJRS0bCm7s3Q623PS02U16ocf5MBsvV6O+wkKkrdZWfKMtfIqPjdCq7VN1OjhAWzZIpfSH5HqaNoU6N5dvrcWLeRt8+YyPxcXy2BUXCwHnG/ZAvzxhwx0ubny+R4ecvtWreRtUBAQGCgrcwEB8hhYxmhd2+VnMsljnpMjw1lWlry1VOo0GtmlqNHIY24JgL6+stp3/rwMpydOyOqd0SinOrjtNvvQZpnAcv16OXbs0iU5lv/yZfl6Pj7247aaNZOv4eVlW8xmW0UwN9c2TswydULDhtf/m8Nslm0sKbEtJpP8nRoM9mPDhJDHwdLGwEA5Dq0qY8fy8+XvKyNDhvDAwPK3KymR/82YTLYKpVYr72dkyM92ejqQmSl/dxER8rhGRFR+UkBJifxjxPIHxPHjtj8kAgPl9YOHDQPatSv7fgoK5DHy96/5OLmSEnksSv8hodHYxuJ5e1ctmJtMwJkz8n0cOyZv8/Jsn+vQUNmt7eEhf2+WpbDQvop7+bL8XIWHy6VRI/k+L1+Wj2dmyu1btpR/xHh5lW3LyZPyK6GwEIiLk0tV/uArLJR/9Pj6yt/Btce2uFgOK0hJke+jZ8+qHeOqYripRL0LN6UdPCjH4Hz8sfzfBpD/k44YAUyfXqfG5tRHli6tbdvkafJms/0Xp1otvzBzcmSVJjtbftH4+cn//Pz85Bfcd9/J/1wrYvliDguT/+GGhcmAUbp7raREznG0aZMcDF6ds+FK8/WV/yFWp+tPq5X/sRYV2cKjs7RqJac/uHBBvtfSE1g6i6+vfI+enrZbo1G+1/z86x8ry9gwDw8Z9Mxm+8d9fGRYiY6WQSw/3xbSLl+WlbqzZ+XPFiqVrNT17Qv06yfbaAnE27bZQmp1BQTIL1VLiA8Kkq974oTsHi4uvv4+WreWJ4fm5Mgv7ZMn5RcsII9dgwa2xdtb/uGg09nfWn7WauV/jcnJtiUz8/pt0OnkZzIw0BbMhZDvxRKCL12q2vtxJB8fORrhn/8EunWT//aXLJF/XFyrRQtZFTYY5GfGsmRny5D799/2x8LT0xbKPTzkZ+b8edvnrX9/+8v7OALDTSXqdbixMBqBb78F3n1X/vkMyP+9Ro4Enn9eDkght3b4sJyk8ccf5V+hvXrZxiBV93IaubnyY7R/v/xisVRAzpyxVRQsX9QBAfIst9695Wu2by8/eufOyb/K//pLPu/KFVu14fJl+Z9qerr9uKlrqdUyxFnGYVn+eZvNsh0mkwwGOTm2xWiUX6qlK05mM5CUJMPbtaHNMoFlz54y+Fm+0Pz9ZZA4dcq2pKTIv8xL/8WvVstg4OMjb3U62zixtDQZGm+EWl02xJSm08njUt0vWB8f+by//77+dgaDfZVSrZaBwlKZCAmRxyglRf6OK/tdlm63paLXqpUMMi1byuO7ahXw228y5LqKl5cMQiaT/N1WN9TrdLL9ljM4fX1ldcsyTjAjQ+5TpbItWq38jFqWwEAZOM6fty3Z2baAGBJim9P1zJny26FWA/HxsuqzbRtw9GjV34NWW/kx12plVa5HD+Dzz6t3fK6H4aYSDDfX2L1bnjK+fLm8r1bLmfKefVb+6yO6QWaz7T9oRykokF8Cly/LLxmDwbbo9fLjWx0lJRV3KWRlyWuobd4sv6BvucV+AktHM5tl+Lhyxb4rr7hYfil6edm/V09P2RaNRh5jS/dJXp5ciott3XqWroniYhk+Dx2Sy5kzti4GSzegZQqFiAhbF2BqqjwOGzfKClZ+vuyG7NXLFlI1Gtt7sXyrVPa7z8qSofbiRVu3y6VLsj2WMNO4ceW/0+xsYM0aWYlo0EAGVMvi4yPDgqVrLCPDvpu4sFAupe8bjTIcREbKpVkz2QZLxaf0+xFC/o7y8mxdf5bl0iW5belqTlCQ7EYqfZycSQj5R8fXX8tJ7NPTZRfe+PHyv/jGjW3bXr4su70tc5Kp1XJRqeR7b9JEfh6aNJGfEaNR7s8yJq+oSHY1Nm0qPz/V/XdYVQw3lWC4qcC+fcCsWcDq1bZ13brJaYFHjbK/eiYREdUZJSUyjDRqVLfn6mK4qQTDzXXs2AG8+KI8ndxSc9VoZAdqVJTtvGs/P/mnTZ8+8k8SIiIiJ2K4qQTDTRWlpcnh9F9+aX91zGupVPLCUP37y6Vr17r/5wEREdU6DDeVYLi5AX/9BaxbJzuSs7Jsy4EDcqKSa4WEyMATEyPPwLJM49uwIUMPERHdEIabSjDcOFhamhx1aRl5efRoxacQ+PnZJkCJjJTdXJGRcuBy8+YMPkREVCGGm0ow3DhZQYE8DWPfPjnj2ZEjsvJz9mzlU/C2bAn84x/y+gN9+sj9bNokz8lNSpLPv+ce4LHHgE6dXPVuiIiolmC4qQTDjUIKC+XkFH/9JWfGOn3aNkvW4cP2k2/4+MhwU1EFqHdvGXL69pXbWWY3KymR525yzA8RkdthuKkEw00tlJMjx/T8+KO8BoFl9uSWLeX0sAMGyEksFi0CVqy4/kxnvr5ytq/WreV4n9JdYI0by1nnSl9ds6BABqrSS1iYfE2GJCKiWoHhphIMN7Wc2SwHKgcEyCrMtc6fBz74QAadtDTbrGZeXjKI/P33jV8L4Fo6nW32qtBQ24xpHh5yadFCzmAWG2t/oSiLrCzb/O5ERFQjDDeVYLhxExVNgVpUJOf+t1yd7uRJWzfYmTO2ecP9/GRgadhQhiPLtK65uXLa08zMql2mG5ABJjYW6NBBTtdp6W7LypJTdTZvLqdwjY6W1aO0NPn4mTPyVquVU9927iyXmBg5/WxBgW0aVb1eBi1nTf1JRFTLMdxUguGmHjOb5Rzvvr7Xr6YUFckq0dmzcrFc9MVkkt1iRqO8mNKWLbZuNGczGOSZZe3ayet/abX2c8ir1bLK1KyZbSnv0r01IYR8v2fPym7DgADH7ZuIqBIMN5VguCGHEkJerGfLFnnbuLHtojTNmsnxRIcPyzPIDh+WVw0MC7PfJj9fXuPLsiQn2/Zv6dbKz7/xSwpbLkyk18vxRJarGTZsKH82mWwX4cnIkJWr0hccCgyU7+PYMblkZdn23aYNEBcnLzQUHW27gqJlsVxG3NPTVnWyXHUwL0++L61Wvp6vr+suvENEdQ7DTSUYbqjWy8+XQUCrtQWC4mJ5ttmRI3IuoePHZbDS620BqKjIdsnlM2dkUHEGlUpO1Fjd/Xt6yudWdklhLy95eQ/LpaQti4eHLQzl5clqlUolw5Bl0WptV5f08rIter3tVq22VbosVa9r/wv08pID0Js3l0vpq05mZMiux/R02XWZkSFvMzPlMeneXYa98i6tbrlsN7sWiW5Idb6/nXR9WyK6YQZD2XWenrJKUp0rtRcUyIqLZexOQYEcU2S5TLJl8fCQFZyQEHnr71/2Msc6ne31W7aUX/gZGfLSHNu2yeX0afkalhBiNNq359rKk0ol92O59LWlzQUFMkDUFo0ayfZlZlb9OVFRQMeO8jhajrNlHFdgIBAcbFtKXzo6MFD+Pq5csS1ZWfL3V1RkWwDbJb8tz7dUzby9bZUzHx9ZEbOcBajTyXCl0dhuPTzsK2slJfJ3m5Yml8uX5eWsW7Ysf5qFkhLZTj8/GTBLEwK4cEF24R4+LKuWPXrIyiXPRCQnYuWGiJzDZLJ9GRuN8tZstn356vW2LzijUQaBnBz5RWn5YrUsludZFr1efnFaxkFZXssy75Hl1lKdsYQ7s1k+11LxsnzZl5aTI6tkJ0/KLrrS1Gr7Lj1LKAwOllWzbdvkl3hdpNHIkGM0VjyY3mCQZwn6+dkqV5cv27YPCpIBJixMrjtwoPxQGBoqQ05srPy5dNjTam1j20wm+Xu0/D5OnJA/GwxAz57ybMXu3WWAA+TvNyNDjpfLz7eFOstnzmSyH6em0chxauWNwRPC9h4tAdByK4T8vFmCuRDyfQQFVR7ahJCfj99+k8vly/J4tmxpu23dmmPZKsBuqUow3BBRlQghr6d2+rQMQZYv4euNC8rKAnbskN2HgYG2MGS5ttrFi/ZL6QrZ5cvySz0wUFbQAgLkrZeXbINWKxezWYbA0s+zdNlZJrW0nP2Xk2O7reo0CWq1DG6hobINf/8tx4JZutaqQ62WFb/oaDkQfc+eGx8/Vtn+c3Nlleh682CVxzIOrmlTGWgtXbv5+dXbj14vq1zh4fJ3qNXKMGQZ/L9pU9WqkmFh8qSBtm1l4GncWC7h4bJ6ZjbL32d2tq066+lpey2tVn7eAgOrfyxqMYabSjDcEFG9JYT8YjSb7c/+K12F8PIqP8QVFckv/BMnZHAqXbkKCpJftBcuyC/v1FS53w4dZKixjFsC5Bfx7t3A1q3AwYP2QS8zUz7PMqeURiNDXWSkrbLRvLkMnVu2yKX0AHxABsjQUFmtsQS8vDz76SMsYdFovH6ACQyUx6u4WLatuFgGqtJhwmyWAbMq9Hp5iZnbbpNVo1On5DG1LBcuVG0/VREQII9XixYyMF2+bKuGpqfLwGWZt8vDQ76XBg1kgLIs3t62wGxZLJUvy9g1T0/ZFduihVyaN5f3S//eHYDhphIMN0REbuT8eTmmJzBQVjdCQ+WXbWlC2H+RW7qOhJChylKpSUmRXVzNmslAFRFR/gSd5TEaZVvOnwfOnZNhr/Q4KSGArl1ld1plU1FkZ8uzEo8elcupU/b7LSiQ26nVtrMMvbzsQ2pRkazsKaltW3kChAMx3FSC4YaIiOokIWQ3lEYjxxBVNr4nL09WtU6elAEpNVVW5EqfhWgJRZalsFBWdC5csC0FBfaD1C1nIVrO0rRMVXHqlG1s1MmT8hqAP/3k0LfPs6WIiIjcjUolB3NXhbe3nBm9fXvntqk8lkqZgjjhAhERETmOZVyTghhuiIiIyK0w3BAREZFbYbghIiIit8JwQ0RERG6F4YaIiIjcCsMNERERuRWGGyIiInIrDDdERETkVhhuiIiIyK0w3BAREZFbYbghIiIit8JwQ0RERG6F4YaIiIjciofSDXA1IQQAIDs7W+GWEBERUVVZvrct3+OVqXfhJicnBwAQERGhcEuIiIiounJycuDv71/pNipRlQjkRsxmM86fPw9fX1+oVCqH7js7OxsRERE4e/Ys/Pz8HLpvssdj7To81q7DY+06PNau46hjLYRATk4OwsPDoVZXPqqm3lVu1Go1mjRp4tTX8PPz4z8WF+Gxdh0ea9fhsXYdHmvXccSxvl7FxoIDiomIiMitMNwQERGRW2G4cSCdTodZs2ZBp9Mp3RS3x2PtOjzWrsNj7To81q6jxLGudwOKiYiIyL2xckNERERuheGGiIiI3ArDDREREbkVhhsiIiJyKww3DvLee+8hMjISer0ecXFx2L59u9JNqvPmzJmDrl27wtfXFw0bNsTQoUNx7Ngxu20KCwsxefJkBAcHw8fHB8OHD0daWppCLXYfc+fOhUqlwtSpU63reKwd59y5c7jvvvsQHBwMLy8v3HTTTdi5c6f1cSEEZs6ciUaNGsHLywvx8fH466+/FGxx3WQymfDcc88hKioKXl5eaNGiBV588UW7axPxWN+4TZs24c4770R4eDhUKhW+++47u8ercmwvXbqEsWPHws/PDwEBAZg4cSJyc3Nr3jhBNbZ06VKh1WrF4sWLxaFDh8SkSZNEQECASEtLU7ppddrAgQPFJ598Ig4ePCj27t0rBg8eLJo2bSpyc3Ot2zzyyCMiIiJCJCUliZ07d4ru3buLnj17Ktjqum/79u0iMjJSdOzYUTzxxBPW9TzWjnHp0iXRrFkzMWHCBLFt2zZx6tQpsXbtWnHixAnrNnPnzhX+/v7iu+++E/v27RN33XWXiIqKEgUFBQq2vO55+eWXRXBwsPjxxx/F6dOnxfLly4WPj4946623rNvwWN+4n3/+WcyYMUOsXLlSABCrVq2ye7wqx/aOO+4QMTEx4s8//xSbN28WLVu2FGPGjKlx2xhuHKBbt25i8uTJ1vsmk0mEh4eLOXPmKNgq95Oeni4AiI0bNwohhLhy5Yrw9PQUy5cvt25z5MgRAUBs3bpVqWbWaTk5OaJVq1Zi3bp1ol+/ftZww2PtOE8//bTo3bt3hY+bzWYRFhYmXn/9deu6K1euCJ1OJ77++mtXNNFtDBkyRDzwwAN26+655x4xduxYIQSPtSNdG26qcmwPHz4sAIgdO3ZYt/nll1+ESqUS586dq1F72C1VQ0VFRdi1axfi4+Ot69RqNeLj47F161YFW+Z+srKyAABBQUEAgF27dqG4uNju2Ldt2xZNmzblsb9BkydPxpAhQ+yOKcBj7Ujff/89YmNjce+996Jhw4a4+eab8eGHH1ofP336NFJTU+2Otb+/P+Li4nisq6lnz55ISkrC8ePHAQD79u3DH3/8gUGDBgHgsXamqhzbrVu3IiAgALGxsdZt4uPjoVarsW3bthq9fr27cKajZWZmwmQyITQ01G59aGgojh49qlCr3I/ZbMbUqVPRq1cvdOjQAQCQmpoKrVaLgIAAu21DQ0ORmpqqQCvrtqVLl2L37t3YsWNHmcd4rB3n1KlTWLBgARITE/Hss89ix44dePzxx6HVajF+/Hjr8Szv/xQe6+p55plnkJ2djbZt20Kj0cBkMuHll1/G2LFjAYDH2omqcmxTU1PRsGFDu8c9PDwQFBRU4+PPcEN1wuTJk3Hw4EH88ccfSjfFLZ09exZPPPEE1q1bB71er3Rz3JrZbEZsbCxeeeUVAMDNN9+MgwcPYuHChRg/frzCrXMvy5Ytw5dffomvvvoK7du3x969ezF16lSEh4fzWLs5dkvVUEhICDQaTZmzRtLS0hAWFqZQq9zLlClT8OOPP2L9+vVo0qSJdX1YWBiKiopw5coVu+157Ktv165dSE9PR+fOneHh4QEPDw9s3LgRb7/9Njw8PBAaGspj7SCNGjVCdHS03bp27dohJSUFAKzHk/+n1NxTTz2FZ555BqNHj8ZNN92E+++/H//6178wZ84cADzWzlSVYxsWFob09HS7x0tKSnDp0qUaH3+GmxrSarXo0qULkpKSrOvMZjOSkpLQo0cPBVtW9wkhMGXKFKxatQq///47oqKi7B7v0qULPD097Y79sWPHkJKSwmNfTQMGDMCBAwewd+9e6xIbG4uxY8daf+axdoxevXqVmdLg+PHjaNasGQAgKioKYWFhdsc6Ozsb27Zt47Gupvz8fKjV9l9zGo0GZrMZAI+1M1Xl2Pbo0QNXrlzBrl27rNv8/vvvMJvNiIuLq1kDajQcmYQQ8lRwnU4nlixZIg4fPiweeughERAQIFJTU5VuWp32f//3f8Lf319s2LBBXLhwwbrk5+dbt3nkkUdE06ZNxe+//y527twpevToIXr06KFgq91H6bOlhOCxdpTt27cLDw8P8fLLL4u//vpLfPnll8JgMIgvvvjCus3cuXNFQECAWL16tdi/f7+4++67eXryDRg/frxo3Lix9VTwlStXipCQEPHvf//bug2P9Y3LyckRe/bsEXv27BEAxLx588SePXvEmTNnhBBVO7Z33HGHuPnmm8W2bdvEH3/8IVq1asVTwWuTd955RzRt2lRotVrRrVs38eeffyrdpDoPQLnLJ598Yt2moKBAPProoyIwMFAYDAYxbNgwceHCBeUa7UauDTc81o7zww8/iA4dOgidTifatm0rFi1aZPe42WwWzz33nAgNDRU6nU4MGDBAHDt2TKHW1l3Z2dniiSeeEE2bNhV6vV40b95czJgxQxiNRus2PNY3bv369eX+Hz1+/HghRNWO7cWLF8WYMWOEj4+P8PPzEwkJCSInJ6fGbVMJUWqqRiIiIqI6jmNuiIiIyK0w3BAREZFbYbghIiIit8JwQ0RERG6F4YaIiIjcCsMNERERuRWGGyIiInIrDDdEVO+pVCp89913SjeDiByE4YaIFDVhwgSoVKoyyx133KF004iojvJQugFERHfccQc++eQTu3U6nU6h1hBRXcfKDREpTqfTISwszG4JDAwEILuMFixYgEGDBsHLywvNmzfHihUr7J5/4MAB3HrrrfDy8kJwcDAeeugh5Obm2m2zePFitG/fHjqdDo0aNcKUKVPsHs/MzMSwYcNgMBjQqlUrfP/9985900TkNAw3RFTrPffccxg+fDj27duHsWPHYvTo0Thy5AgAIC8vDwMHDkRgYCB27NiB5cuX47fffrMLLwsWLMDkyZPx0EMP4cCBA/j+++/RsmVLu9d44YUXMHLkSOzfvx+DBw/G2LFjcenSJZe+TyJykBpfepOIqAbGjx8vNBqN8Pb2tltefvllIYS8Ovwjjzxi95y4uDjxf//3f0IIIRYtWiQCAwNFbm6u9fGffvpJqNVqkZqaKoQQIjw8XMyYMaPCNgAQ//nPf6z3c3NzBQDxyy+/OOx9EpHrcMwNESnulltuwYIFC+zWBQUFWX/u0aOH3WM9evTA3r17AQBHjhxBTEwMvL29rY/36tULZrMZx44dg0qlwvnz5zFgwIBK29CxY0frz97e3vDz80N6evqNviUiUhDDDREpztvbu0w3kaN4eXlVaTtPT0+7+yqVCmaz2RlNIiIn45gbIqr1/vzzzzL327VrBwBo164d9u3bh7y8POvjW7ZsgVqtRps2beDr64vIyEgkJSW5tM1EpBxWbohIcUajEampqXbrPDw8EBISAgBYvnw5YmNj0bt3b3z55ZfYvn07Pv74YwDA2LFjMWvWLIwfPx7PP/88MjIy8Nhjj+H+++9HaGgoAOD555/HI488goYNG2LQoEHIycnBli1b8Nhjj7n2jRKRSzDcEJHi1qxZg0aNGtmta9OmDY4ePQpAnsm0dOlSPProo2jUqBG+/vprREdHAwAMBgPWrl2LJ554Al27doXBYMDw4cMxb948677Gjx+PwsJCvPnmm5g2bRpCQkIwYsQI171BInIplRBCKN0IIqKKqFQqrFq1CkOHDlW6KURUR3DMDREREbkVhhsiIiJyKxxzQ0S1GnvOiai6WLkhIiIit8JwQ0RERG6F4YaIiIjcCsMNERERuRWGGyIiInIrDDdERETkVhhuiIiIyK0w3BAREZFbYbghIiIit/L/sLHtRbQ6Q/4AAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Plot the results\n", + "fig, ax = plt.subplots()\n", + "ax.plot(errors_train,'r-',label='train')\n", + "ax.plot(errors_test,'b-',label='test')\n", + "#ax.set_ylim(0,100); ax.set_xlim(0,n_epoch)\n", + "ax.set_xlabel('Epoch'); ax.set_ylabel('MSE')\n", + "ax.set_title('Train Loss %3.5f, Test Loss %3.5f'%(errors_train[-1],errors_test[-1]))\n", + "ax.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [], + "source": [ + "# Physics Informed Model\n", + "model_pinn = MixtureOfExperts(\n", + " D_i, D_o, D_h,\n", + " n_experts=n_experts,\n", + " activation_fn=acivation_fun,\n", + " gating_hidden=gating_hidden_layer_width,\n", + " dropout_rate=dropout_rate,\n", + " use_bn=use_bn\n", + ")\n", + "model_pinn.apply(weights_init)\n", + "model_pinn.to(device)\n", + "\n", + "lambda_bw_mon = 0.1\n", + "lambda_IL_mon = 0.1\n", + "lambda_vpiL = 0.05\n", + "lambda_smooth = 0.01\n", + "\n", + "optimizer = torch.optim.AdamW(model_pinn.parameters(), lr=learning_rate, weight_decay=weight_decay, betas=betas)" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Epoch 0 | Train 0.577587 | Test 0.602078\n", + "Epoch 10 | Train 0.078204 | Test 0.101025\n", + "Epoch 20 | Train 0.053335 | Test 0.077052\n", + "Epoch 30 | Train 0.046076 | Test 0.070806\n", + "Epoch 40 | Train 0.043260 | Test 0.068207\n", + "Epoch 50 | Train 0.040504 | Test 0.065854\n", + "Epoch 60 | Train 0.038483 | Test 0.064187\n", + "Epoch 70 | Train 0.037348 | Test 0.064485\n", + "Epoch 80 | Train 0.037371 | Test 0.066420\n", + "Epoch 90 | Train 0.034940 | Test 0.064381\n" + ] + } + ], + "source": [ + "# ==============================\n", + "# Training Loop\n", + "# ==============================\n", + "\n", + "for epoch in range(n_epoch):\n", + "\n", + " model_pinn.train()\n", + "\n", + " for x_batch, y_batch in data_loader:\n", + "\n", + " x_batch = x_batch.to(device)\n", + " y_batch = y_batch.to(device)\n", + "\n", + " x_batch.requires_grad_(True)\n", + "\n", + " optimizer.zero_grad()\n", + "\n", + " pred = model_pinn(x_batch)\n", + "\n", + " data_loss = loss_function(pred, y_batch)\n", + "\n", + " BW_pred = pred[:, 0]\n", + " IL_pred = pred[:, 1]\n", + " Vpi_pred = pred[:, 2]\n", + " L_batch = x_batch[:, 7]\n", + "\n", + " physics_loss_total = 0.0\n", + "\n", + " # =========================================\n", + " # 1) dBW/dL <= 0\n", + " # =========================================\n", + " if lambda_bw_mon != 0:\n", + "\n", + " grads = torch.autograd.grad(\n", + " BW_pred,\n", + " x_batch,\n", + " grad_outputs=torch.ones_like(BW_pred),\n", + " create_graph=True\n", + " )[0]\n", + "\n", + " dBW_dL = grads[:, 7]\n", + "\n", + " physics_bw = torch.mean(torch.relu(dBW_dL) ** 2)\n", + " physics_loss_total += lambda_bw_mon * physics_bw\n", + "\n", + "\n", + " # =========================================\n", + " # 2) dIL/dL >= 0\n", + " # =========================================\n", + " if lambda_IL_mon != 0:\n", + "\n", + " grads = torch.autograd.grad(\n", + " IL_pred,\n", + " x_batch,\n", + " grad_outputs=torch.ones_like(IL_pred),\n", + " create_graph=True\n", + " )[0]\n", + "\n", + " dIL_dL = grads[:, 7]\n", + "\n", + " physics_il = torch.mean(torch.relu(-dIL_dL) ** 2)\n", + " physics_loss_total += lambda_IL_mon * physics_il\n", + "\n", + "\n", + " # =========================================\n", + " # 3) d(Vpi*L)/dL ≈ 0\n", + " # =========================================\n", + " if lambda_vpiL != 0:\n", + "\n", + " VpiL = Vpi_pred * L_batch\n", + "\n", + " grads = torch.autograd.grad(\n", + " VpiL,\n", + " x_batch,\n", + " grad_outputs=torch.ones_like(VpiL),\n", + " create_graph=True\n", + " )[0]\n", + "\n", + " dVpiL_dL = grads[:, 7]\n", + "\n", + " physics_vpiL = torch.mean(dVpiL_dL ** 2)\n", + " physics_loss_total += lambda_vpiL * physics_vpiL\n", + "\n", + "\n", + " # =========================================\n", + " # 4) Smoothness (second derivative of BW)\n", + " # =========================================\n", + " if lambda_smooth != 0:\n", + "\n", + " grads1 = torch.autograd.grad(\n", + " BW_pred,\n", + " x_batch,\n", + " grad_outputs=torch.ones_like(BW_pred),\n", + " create_graph=True\n", + " )[0]\n", + "\n", + " dBW_dL = grads1[:, 7]\n", + "\n", + " grads2 = torch.autograd.grad(\n", + " dBW_dL,\n", + " x_batch,\n", + " grad_outputs=torch.ones_like(dBW_dL),\n", + " create_graph=True\n", + " )[0]\n", + "\n", + " d2BW_dL2 = grads2[:, 7]\n", + "\n", + " physics_smooth = torch.mean(d2BW_dL2 ** 2)\n", + " physics_loss_total += lambda_smooth * physics_smooth\n", + "\n", + "\n", + " loss = data_loss + physics_loss_total\n", + "\n", + " loss.backward()\n", + " optimizer.step()\n", + "\n", + "\n", + " # =========================\n", + " # Evaluation\n", + " # =========================\n", + " model_pinn.eval()\n", + " with torch.no_grad():\n", + "\n", + " pred_train = model_pinn(x_train)\n", + " pred_test = model_pinn(x_test)\n", + "\n", + " train_loss = loss_function(pred_train, y_train).item()\n", + " test_loss = loss_function(pred_test, y_test).item()\n", + "\n", + " errors_train[epoch] = train_loss\n", + " errors_test[epoch] = test_loss\n", + "\n", + " if epoch % 10 == 0 or epoch == 0:\n", + " print(f\"Epoch {epoch:5d} | Train {train_loss:.6f} | Test {test_loss:.6f}\")" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAHHCAYAAABDUnkqAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAXVtJREFUeJzt3XlcVOXiP/DPDDAzLLLDIAqCW+6YomiaWuJuqWVaWS7tqaU/rrf0etOyb2nZ7Wrq1ZtldkuzNG1PM1xKM3fcd0VIZXNhh4GZ5/fH4wyMLILMzIHh8369zmvgzJkzzxwG5sOzqoQQAkREREROQq10AYiIiIhsieGGiIiInArDDRERETkVhhsiIiJyKgw3RERE5FQYboiIiMipMNwQERGRU2G4ISIiIqfCcENEREROheGGnML48eMRERGhdDGIiKgWYLghu1KpVFXatm3bpnRRrWzbtg0qlQrr1q1Tuii3VVhYiFdffRWhoaFwd3dHTEwMNm/eXOXHX7p0CaNGjYKvry+8vb0xbNgwnD9/3uqY/Px8PP3002jXrh18fHzg5eWFqKgoLFy4EEVFRZWe/9lnn4VKpcLQoUPLvf+7775Dp06doNPpEB4ejtmzZ6O4uLjMcZs3b0bPnj3h4eEBPz8/jBw5EomJiVV+naU58n2Zl5eH119/vcrn4nvvfLnHpqam4vnnn0ejRo2g0+kQERGBp59+2uqY119/vdyfo06nszquOu/nPn36VPj+cHNzq/JrJcdyVboA5Nw+++wzq+//97//YfPmzWX2t27dukbPs3z5cphMphqdo64aP3481q1bh6lTp6JFixZYuXIlBg8ejK1bt6Jnz56VPjYnJwf33XcfMjMz8Y9//ANubm7497//jd69eyMhIQEBAQEA5IfBsWPHMHjwYERERECtVuOPP/7A//t//w+7d+/G6tWryz3/vn37sHLlyjIfLmY///wzhg8fjj59+mDRokU4cuQI/u///g9paWlYunSp5bgffvgBw4YNQ6dOnTBv3jxkZWVh4cKF6NmzJw4ePIigoKBqXTNHvS8BGW7eeOMNAPKD0pk44r0HAMnJyejRowcA4IUXXkCjRo1w+fJl7Nmzp9xzL126FF5eXpbvXVxcrO6vzvt55syZeOaZZ6wen5ubixdeeAH9+/ev2oUixxNEDjRp0iRRlbddbm6uA0pTsa1btwoAYu3atYqW43Z2794tAIj58+db9uXn54tmzZqJ7t273/bx77zzjgAg9uzZY9l34sQJ4eLiImbMmHHbx0+ePFkAEFeuXClzn8lkEt27dxdPPfWUaNKkiRgyZEiZY9q0aSOioqJEUVGRZd/MmTOFSqUSJ06csDquefPmorCw0LIvISFBqNVqERcXd9ty3k5V35d3Ij09XQAQs2fPrtLxfO+Vfe8NGjRIREZGioyMjErPOXv2bAFApKenV/PVSJW9n0v77LPPBACxatWqO3oesj82S5Hi+vTpg3bt2mH//v3o1asXPDw88I9//AMA8O2332LIkCEIDQ2FVqtFs2bN8Oabb8JoNFqd49Y+N4mJiVCpVHjvvffw4YcfolmzZtBqtejSpQv27t1rs7KfP38ejzzyCPz9/eHh4YFu3brhxx9/LHPcokWL0LZtW0uTSnR0tNV/h9nZ2Zg6dSoiIiKg1WoRHByMfv364cCBA5U+/7p16+Di4oLnnnvOsk+n0+Hpp5/Grl27kJycfNvHd+nSBV26dLHsa9WqFfr27Yuvvvrqtq/ffM1v3LhR5r7PPvsMR48exVtvvVXuY48fP47jx4/jueeeg6trSSXyxIkTIYSwNMtcu3YNx48fx4gRI6DRaCzHRUVFoXXr1lizZs1ty3knTCYTFixYgLZt20Kn00Gv1+P555/H9evXrY7bt28fBgwYgMDAQLi7uyMyMhJPPfUUAPk+NNcqvfHGG5bmjNdff73G5asv772TJ0/i559/xt///ncEBASgoKDgtk2hQghkZWVBCFHpcbeq7P1c2urVq+Hp6Ylhw4ZV6/zkOGyWolrh6tWrGDRoEB599FE88cQT0Ov1AICVK1fCy8sLcXFx8PLywpYtWzBr1ixkZWVh/vz5tz3v6tWrkZ2djeeffx4qlQrvvvsuHnroIZw/f77G7eWpqam45557kJeXh5dffhkBAQH49NNP8eCDD2LdunUYMWIEANlk9vLLL2PkyJGYMmUKCgoKcPjwYezevRuPP/44AFnVvm7dOkyePBlt2rTB1atXsWPHDpw4cQKdOnWqsAwHDx5Ey5Yt4e3tbbW/a9euAICEhASEhYWV+1iTyYTDhw9bPohvffwvv/yC7OxsNGjQwLLfYDAgKysL+fn52LdvH9577z00adIEzZs3t3p8dnY2Xn31VfzjH/9ASEhIhWUHgOjoaKv9oaGhaNy4seX+wsJCAIC7u3uZc3h4eODYsWNISUmp8Hnu1PPPP4+VK1diwoQJePnll3HhwgUsXrwYBw8exM6dO+Hm5oa0tDT0798fQUFBmD59Onx9fZGYmIj169cDAIKCgrB06VK8+OKLGDFiBB566CEAQIcOHWpUtvr03vv1118BAHq9Hn379sWWLVvg4uKCfv36YenSpeUOJGjatClycnLg6emJ4cOH41//+pflb0ppVX0/l5aeno7Nmzdj9OjR8PT0rPA4UpjCNUdUz5RX/d+7d28BQCxbtqzM8Xl5eWX2Pf/888LDw0MUFBRY9o0bN040adLE8v2FCxcEABEQECCuXbtm2f/tt98KAOL777+vtJxVaRqYOnWqACB+//13y77s7GwRGRkpIiIihNFoFEIIMWzYMNG2bdtKn8/Hx0dMmjSp0mPK07ZtW3H//feX2X/s2LEKr6mZublkzpw5Ze5bsmSJACBOnjxptf+LL74QACxbdHS0OHz4cJnHT5s2TURGRlp+RuU1S82fP18AEElJSWUe36VLF9GtWzchhBBGo1H4+vqKvn37Wh2TkZEhPD09BQCxb9++Cl9nVdz6vvz999/LbXbYuHGj1f4NGzYIAGLv3r0VntsezVL16b338ssvW36XBw4cKL788ksxf/584eXlJZo1a2bVhL1gwQIxefJksWrVKrFu3ToxZcoU4erqKlq0aCEyMzPLPFdV38+lLVq0SAAQP/30U6XHkbLYLEW1glarxYQJE8rsL/3fenZ2NjIyMnDvvfciLy8PJ0+evO15R48eDT8/P8v39957LwBUOCKjOn766Sd07drVquOkl5cXnnvuOSQmJuL48eMAAF9fX/z111+VNof5+vpi9+7duHz5crXKkJ+fD61WW2a/uQNvfn5+pY8FUK3H33fffdi8eTPWrl2LF154AW5ubsjNzbU65vTp01i4cCHmz59f7rmr+vzm+9VqNZ5//nnEx8djxowZOHPmDPbv349Ro0bBYDDc9nXeibVr18LHxwf9+vVDRkaGZevcuTO8vLywdetWAPLnBsgOz7drKrGl+vTey8nJAQCEhITgxx9/xKhRozBt2jQsX74c586ds2pimzJlChYtWoTHH38cDz/8MBYsWIBPP/0UZ86cwX/+858yz1WV9/OtVq9ejaCgIPTr16/S40hZDDdUKzRq1MiqP4XZsWPHMGLECPj4+MDb2xtBQUF44oknAACZmZm3PW94eLjV9+agc2u/iTtx8eJF3HXXXWX2m0fYXLx4EQDw6quvwsvLC127dkWLFi0wadIk7Ny50+ox7777Lo4ePYqwsDB07doVr7/+epUCmLu7u6XZprSCggLL/ZU9FkC1Hq/X6xEbG4uRI0di6dKlGDp0KPr164eUlBTLMVOmTME999yDhx9++LZlr+z5Sz/3nDlz8PTTT+Pdd99Fy5YtER0dDVdXV8tQ4NIjY2zhzJkzyMzMRHBwMIKCgqy2nJwcpKWlAQB69+6Nhx9+GG+88QYCAwMxbNgwfPLJJ+W+JluqT+898+2oUaOgVpd8ZD3yyCNwdXXFH3/8UWk5H3/8cYSEhFiat0qryvu5tPPnz2PXrl0YPXq0VT8xqn0YbqhWKO8P4Y0bN9C7d28cOnQIc+bMwffff4/NmzfjnXfeAYAqDf2+dQiomahmR8OaaN26NU6dOoU1a9agZ8+e+Prrr9GzZ0/Mnj3bcsyoUaNw/vx5LFq0CKGhoZg/fz7atm2Ln3/+udJzN2zYEFeuXCmz37wvNDS0wsf6+/tDq9Xe8eMBYOTIkcjJycG3334LANiyZQs2btyIKVOmIDEx0bIVFxcjPz8fiYmJyMrKspS99HPd+vyln1uj0eCjjz7C5cuX8dtvv+HUqVPYtGkTMjMzoVarK+0jcSdMJhOCg4OxefPmcrc5c+YAgGU+ml27dmHy5Mm4dOkSnnrqKXTu3NlS46AkZ3jvmW9v7TPj4uKCgICAKv2jEhYWhmvXrt32uFvfz7cy1xKNGTPmtuciZTHcUK21bds2XL16FStXrsSUKVMwdOhQxMbGWjUzKalJkyY4depUmf3m5rImTZpY9nl6emL06NH45JNPkJSUhCFDhuCtt96y/JcKyA+LiRMn4ptvvsGFCxcQEBBQ4Ugjs44dO+L06dOWwGC2e/duy/0VUavVaN++Pfbt21fmvt27d6Np06ZWnYnLY246MNeiJSUlAQAeeughREZGWrZLly5hy5YtiIyMxIoVK6zKduvzX758GX/99Ve5Zdfr9bj33nvRsmVLGI1GbNu2DTExMTavuWnWrBmuXr2KHj16IDY2tswWFRVldXy3bt3w1ltvYd++fVi1ahWOHTtmGcWlUqlsWjagfr33OnfuDEBO+FeawWBARkbGbec4EkJYjVqrzK3v51utXr0azZo1Q7du3W57LlIWww3VWuZal9K1LAaDody2cyUMHjwYe/bswa5duyz7cnNz8eGHHyIiIgJt2rQBIEeClabRaNCmTRsIIVBUVASj0Vjmj2lwcDBCQ0Nv27wxcuRIGI1GfPjhh5Z9hYWF+OSTTxATE2M1WiUpKalMP6WRI0di7969Vh8yp06dwpYtW/DII49Y9mVkZJRb2/XRRx8BKBnxdP/992PDhg1ltqCgIERHR2PDhg144IEHAABt27ZFq1at8OGHH1oN7V+6dClUKhVGjhxZ6Wt/7733cOXKFfztb3+r9Lg7MWrUKBiNRrz55ptl7isuLrYMFb5+/XqZ62L+UDf/7Dw8PADcfnhxddSn916fPn0QHByMVatWWQWylStXwmg0WvV9SU9PL1POpUuXIj09HQMHDrTsq+r7ubSDBw/ixIkTllFmVLux0ZBqrXvuuQd+fn4YN24cXn75ZahUKnz22WcObVL6+uuvy+24PG7cOEyfPh1ffPEFBg0ahJdffhn+/v749NNPceHCBXz99deW/gH9+/dHSEgIevToAb1ejxMnTmDx4sUYMmQIGjRogBs3bqBx48YYOXIkoqKi4OXlhV9//RV79+7Fv/71r0rLFxMTg0ceeQQzZsxAWloamjdvjk8//RSJiYn4+OOPrY4dO3Ystm/fbnX9Jk6ciOXLl2PIkCGYNm0a3Nzc8P7770Ov11uFhs8//xzLli3D8OHD0bRpU2RnZ2PTpk3YvHkzHnjgAdx///0AZB+nW/s5AcDUqVOh1+sxfPhwq/3z58/Hgw8+iP79++PRRx/F0aNHsXjxYjzzzDNWswN//vnn+Prrr9GrVy/L9fnqq6/wzDPPlOnbM378eMvP4U7XG+vduzeef/55zJ07FwkJCejfvz/c3Nxw5swZrF27FgsXLsTIkSPx6aef4j//+Q9GjBiBZs2aITs7G8uXL4e3tzcGDx4MQDa5tmnTBl9++SVatmwJf39/tGvXDu3atau0DHzvSVqtFvPnz8e4cePQq1cvPPnkk0hKSsLChQtx7733WobXA7LGavTo0Wjfvj10Oh127NiBNWvWoGPHjnj++ectx1X1/VzaqlWrALBJqs5QaJQW1VMVDQWvaLjqzp07Rbdu3YS7u7sIDQ0Vr7zyiti0aZMAILZu3Wo5rqKh4KVnTzVDFYblmofjVrSZh+CeO3dOjBw5Uvj6+gqdTie6du0qfvjhB6tz/fe//xW9evUSAQEBQqvVimbNmom///3vlqGphYWF4u9//7uIiooSDRo0EJ6eniIqKkr85z//qbSMZvn5+WLatGkiJCREaLVa0aVLF7Fx48Yyx5mH3N8qOTlZjBw5Unh7ewsvLy8xdOhQcebMGatj9u7dKx555BERHh4utFqt8PT0FJ06dRLvv/++1ezCFalohmIh5HDqjh07Cq1WKxo3biz++c9/CoPBYHXM7t27Ra9evYSfn5/Q6XQiKipKLFu2TJhMpjLne/jhh4W7u7u4fv36bctlVtEMxR9++KHo3LmzcHd3Fw0aNBDt27cXr7zyirh8+bIQQogDBw6Ixx57zHJdgoODxdChQ8sMTf/jjz9E586dhUajue37j++9M2WOE0IO246KihJarVbo9XoxefJkkZWVZXXMM888I9q0aSMaNGgg3NzcRPPmzcWrr75a5rjqvp+NRqNo1KiR6NSpU5WuCylPJYQD/w0mIrIzvV6PsWPHVmmSRyJyTgw3ROQ0jh07hu7du+P8+fMIDAxUujhEpBCGGyIiInIqHC1FREREToXhhoiIiJwKww0RERE5FYYbIiIicir1bhI/k8mEy5cvo0GDBnaZFp2IiIhsTwiB7OxshIaGWi2iWp56F24uX75sNS04ERER1R3Jyclo3LhxpcfUu3BjXowtOTkZ3t7eCpeGiIiIqiIrKwthYWG3XdAXqIfhxtwU5e3tzXBDRERUx1SlSwk7FBMREZFTYbghIiIip8JwQ0RERE6l3vW5ISIisiej0YiioiKli1EnaTSa2w7zrgqGGyIiIhsQQiAlJQU3btxQuih1llqtRmRkJDQaTY3Ow3BDRERkA+ZgExwcDA8PD04UW03mSXavXLmC8PDwGl0/hhsiIqIaMhqNlmATEBCgdHHqrKCgIFy+fBnFxcVwc3O74/OwQzEREVENmfvYeHh4KFySus3cHGU0Gmt0HoYbIiIiG2FTVM3Y6vox3BAREZFTYbghIiIim4iIiMCCBQuULoby4WbJkiWIiIiATqdDTEwM9uzZU+nxN27cwKRJk9CwYUNotVq0bNkSP/30k4NKS0RE5Fz69OmDqVOn2uRce/fuxXPPPWeTc9WEoqOlvvzyS8TFxWHZsmWIiYnBggULMGDAAJw6dQrBwcFljjcYDOjXrx+Cg4Oxbt06NGrUCBcvXoSvr6/jC1+mbEBqKiAEEB6udGmIiIhsQwgBo9EIV9fbR4agoCAHlOj2FK25ef/99/Hss89iwoQJaNOmDZYtWwYPDw+sWLGi3ONXrFiBa9eu4ZtvvkGPHj0QERGB3r17IyoqysElL2vPHhlqYmOVLgkREVHVjB8/Htu3b8fChQuhUqmgUqmwcuVKqFQq/Pzzz+jcuTO0Wi127NiBc+fOYdiwYdDr9fDy8kKXLl3w66+/Wp3v1mYplUqFjz76CCNGjICHhwdatGiB7777zu6vS7FwYzAYsH//fsSWSgNqtRqxsbHYtWtXuY/57rvv0L17d0yaNAl6vR7t2rXD22+/XemQscLCQmRlZVlt9qDVmp/PLqcnIqK6RAggN1eZTYgqF3PhwoXo3r07nn32WVy5cgVXrlxBWFgYAGD69OmYN28eTpw4gQ4dOiAnJweDBw9GfHw8Dh48iIEDB+KBBx5AUlJSpc/xxhtvYNSoUTh8+DAGDx6MMWPG4Nq1azW6vLejWLNURkYGjEYj9Hq91X69Xo+TJ0+W+5jz589jy5YtGDNmDH766SecPXsWEydORFFREWbPnl3uY+bOnYs33njD5uW/lU4nbwsK7P5URERU2+XlAV5eyjx3Tg7g6VmlQ318fKDRaODh4YGQkBAAsHwGz5kzB/369bMc6+/vb9VS8uabb2LDhg347rvvMHny5AqfY/z48XjssccAAG+//TY++OAD7NmzBwMHDqz2S6sqxTsUV4fJZEJwcDA+/PBDdO7cGaNHj8bMmTOxbNmyCh8zY8YMZGZmWrbk5GS7lI01N0RE5Eyio6Otvs/JycG0adPQunVr+Pr6wsvLCydOnLhtzU2HDh0sX3t6esLb2xtpaWl2KbOZYjU3gYGBcHFxQWpqqtX+1NRUS3q8VcOGDeHm5gYXFxfLvtatWyMlJQUGg6Hchba0Wi205uRhR6y5ISIiCw8PWYOi1HPbgOcttT/Tpk3D5s2b8d5776F58+Zwd3fHyJEjYTAYKj3PrcsoqFQqmEwmm5SxIoqFG41Gg86dOyM+Ph7Dhw8HIGtm4uPjK6ze6tGjB1avXg2TyWRZEv306dNo2LBhjVcQranSNTdCAJykkoioHlOpqtw0pDSNRlOl5Q527tyJ8ePHY8SIEQBkTU5iYqKdS3dnFG2WiouLw/Lly/Hpp5/ixIkTePHFF5Gbm4sJEyYAAMaOHYsZM2ZYjn/xxRdx7do1TJkyBadPn8aPP/6It99+G5MmTVLqJViYa24AOSyciIioLoiIiMDu3buRmJiIjIyMCmtVWrRogfXr1yMhIQGHDh3C448/bvcamDulaLgZPXo03nvvPcyaNQsdO3ZEQkICNm7caOlknJSUhCtXrliODwsLw6ZNm7B371506NABL7/8MqZMmYLp06cr9RIsSocbNk0REVFdMW3aNLi4uKBNmzYICgqqsA/N+++/Dz8/P9xzzz144IEHMGDAAHTq1MnBpa0alRDVGDPmBLKysuDj44PMzEx4e3vb7LxCADdbypCaCpQzByERETmpgoICXLhwAZGRkdCV/m+XqqWy61idz+86NVqqNlOpSvrdsOaGiIhIOQw3NsTh4ERERMpjuLEhDgcnIiJSHsONDbHmhoiISHkMNzbEmhsiIiLlMdzYyqlT0KbLpR1Yc0NERKQchhtbuXYNuhspAFhzQ0REpCSGG1vR6aCDTDUMN0RERMphuLEVrRZayPYoNksREREph+HGVlhzQ0REVCsw3NgKa26IiKgO6tOnD6ZOnWqz840fPx7Dhw+32fnuBMONrWi1JTU3ubdfOp6IiIjsg+HGVnS6kpqbPIYbIiKq/caPH4/t27dj4cKFUKlUUKlUSExMxNGjRzFo0CB4eXlBr9fjySefREZGhuVx69atQ/v27eHu7o6AgADExsYiNzcXr7/+Oj799FN8++23lvNt27bN4a/L1eHP6KxYc0NERDcJAeTlKfPcHh5yMeeqWLhwIU6fPo127dphzpw5AAA3Nzd07doVzzzzDP79738jPz8fr776KkaNGoUtW7bgypUreOyxx/Duu+9ixIgRyM7Oxu+//w4hBKZNm4YTJ04gKysLn3zyCQDA39/fXi+1Qgw3tuLqCi0MAIDCvGKFC0NERErKywO8vJR57pwcwNOzasf6+PhAo9HAw8MDISEhAID/+7//w9133423337bctyKFSsQFhaG06dPIycnB8XFxXjooYfQpEkTAED79u0tx7q7u6OwsNByPiUw3NiKSgWdaxFQDBTkmpQuDRER0R05dOgQtm7dCq9y0tm5c+fQv39/9O3bF+3bt8eAAQPQv39/jBw5En5+fgqUtnwMNzakcy2W4SaP4YaIqD7z8JA1KEo9d03k5OTggQcewDvvvFPmvoYNG8LFxQWbN2/GH3/8gV9++QWLFi3CzJkzsXv3bkRGRtbsyW2E4caGtK4y1BQWMNwQEdVnKlXVm4aUptFoYDSW9BXt1KkTvv76a0RERMDVtfyYoFKp0KNHD/To0QOzZs1CkyZNsGHDBsTFxZU5nxI4WsqGdG7yh1mQLxQuCRERUdVERERg9+7dSExMREZGBiZNmoRr167hsccew969e3Hu3Dls2rQJEyZMgNFoxO7du/H2229j3759SEpKwvr165Geno7WrVtbznf48GGcOnUKGRkZKCoqcvhrYrixIa2bueaG4YaIiOqGadOmwcXFBW3atEFQUBAMBgN27twJo9GI/v37o3379pg6dSp8fX2hVqvh7e2N3377DYMHD0bLli3xz3/+E//6178waNAgAMCzzz6Lu+66C9HR0QgKCsLOnTsd/prYLGVDOo0MN1x+gYiI6oqWLVti165dZfavX7++3ONbt26NjRs3Vni+oKAg/PLLLzYr351gzY0NaTWyxobLLxARESmH4caGdFoZbgoKqjh7EhEREdkcw40NaW+Gm0IDww0REZFSGG5sSKeTtwUMN0RERIphuLEhnVbeFhh4WYmI6iMhOFq2Jmx1/fgpbENanayxKSziZSUiqk/c3NwAAHlKrZbpJAwGuUaji4tLjc7DoeA2pPOQoaagqGY/FCIiqltcXFzg6+uLtLQ0AICHhwdUVV2amwAAJpMJ6enp8PDwqHBm5KpiuLEhS81NMcMNEVF9Y14F2xxwqPrUajXCw8NrHAwZbmxI5ylDTUExLysRUX2jUqnQsGFDBAcHK7LkgDPQaDRQq2vetYOfwjakdZc/EIPRFULIhdOIiKh+cXFxqXGfEaoZ9ny1IXPNDcBZiomIiJTCcGNDWs+SijCGGyIiImUw3NiQxqMk3HDxTCIiImUw3NiQyl0HHfIBMNwQEREpheHGlrRaaCHbo9gsRUREpAyGG1vS6aCDrLJhzQ0REZEyGG5siTU3REREimO4sSWtljU3RERECmO4sSWdjjU3RERECmO4sSXW3BARESmO4caWWHNDRESkOIYbW2LNDRERkeIYbmyJ4YaIiEhxDDe2xGYpIiIixTHc2FLpmpt8oXBhiIiI6ieGG1sqXXOTb1S4MERERPUTw40tla65yWG4ISIiUkKtCDdLlixBREQEdDodYmJisGfPngqPXblyJVQqldWm0+kcWNpKlF5+IY/hhoiISAmKh5svv/wScXFxmD17Ng4cOICoqCgMGDAAaWlpFT7G29sbV65csWwXL150YIkroVZDpy4CABTkMtwQEREpQfFw8/777+PZZ5/FhAkT0KZNGyxbtgweHh5YsWJFhY9RqVQICQmxbHq93oElrpzWtRgAUJhvUrgkRERE9ZOi4cZgMGD//v2IjY217FOr1YiNjcWuXbsqfFxOTg6aNGmCsLAwDBs2DMeOHavw2MLCQmRlZVlt9qRzlTU2BQw3REREilA03GRkZMBoNJapedHr9UhJSSn3MXfddRdWrFiBb7/9Fp9//jlMJhPuuece/PXXX+UeP3fuXPj4+Fi2sLAwm7+O0rRuN8NNHoeCExERKUHxZqnq6t69O8aOHYuOHTuid+/eWL9+PYKCgvDf//633ONnzJiBzMxMy5acnGzX8uluhpvCAoYbIiIiJbgq+eSBgYFwcXFBamqq1f7U1FSEhIRU6Rxubm64++67cfbs2XLv12q10Gq1NS5rVek0MtQUMNwQEREpQtGaG41Gg86dOyM+Pt6yz2QyIT4+Ht27d6/SOYxGI44cOYKGDRvaq5jVonWTfW24/AIREZEyFK25AYC4uDiMGzcO0dHR6Nq1KxYsWIDc3FxMmDABADB27Fg0atQIc+fOBQDMmTMH3bp1Q/PmzXHjxg3Mnz8fFy9exDPPPKPky7DQaW/W3BSqFC4JERFR/aR4uBk9ejTS09Mxa9YspKSkoGPHjti4caOlk3FSUhLU6pIKpuvXr+PZZ59FSkoK/Pz80LlzZ/zxxx9o06aNUi/BirkFrJDhhoiISBEqIUS96hySlZUFHx8fZGZmwtvb2+bn39LlVfTd9w7aNrqOo3/52fz8RERE9VF1Pr/r3Gip2k6rkzU2hQZeWiIiIiXwE9jGzMtcFRTx0hIRESmBn8A2pnWXl7SgyEXhkhAREdVPDDc2pnO/2SxVzHBDRESkBIYbG9N5sOaGiIhISQw3Nqb1kKGmyOQKE9fOJCIicjiGGxsz19wAnKWYiIhICQw3Nqb1LJkXkeGGiIjI8RhubMzNww0qyPaoggKFC0NERFQPMdzYmEqnhRayyoY1N0RERI7HcGNrOh10kFU2rLkhIiJyPIYbW9OW1Nww3BARETkew42tlaq5YbMUERGR4zHc2JpWy2YpIiIiBTHc2JqWHYqJiIiUxHBja+xQTEREpCiGG1tjzQ0REZGiGG5sjTU3REREimK4sTXW3BARESmK4cbWWHNDRESkKIYbW+MkfkRERIpiuLG1UvPcFBaYFC4MERFR/cNwY2ulm6VyjAoXhoiIqP5huLG10h2K8xhuiIiIHI3hxtbc3KAz97nJZbghIiJyNIYbW1OpoHUtBgAU5rPPDRERkaMx3NiB7ma4KchjuCEiInI0hhs70LrJUFNYIBQuCRERUf3DcGMHOlfZ16aAzVJEREQOx3BjB1qNrLHhJH5ERESOx3BjBzqNuVlK4YIQERHVQww3dmAONwVcOJOIiMjhGG7sQKuVt4UGlbIFISIiqocYbuxAp73Z56aQ4YaIiMjRGG7sQKuToaawiJeXiIjI0fjpawc6nbwtMPDyEhERORo/fe1A686aGyIiIqXw09cOdO7yshYUuShcEiIiovqH4cYOtOZwU+yqcEmIiIjqH4YbO9B5yhqbYpMLjEaFC0NERFTPMNzYgc6j5LIWciI/IiIih2K4sQOtR0lfG4YbIiIix2K4sQNXDw3UuLkyONeXIiIiciiGGztQ6bTQQlbZsOaGiIjIsRhu7EGngw6yyoY1N0RERI7FcGMPWtbcEBERKYXhxh60WtbcEBERKYThxh50OkvNDcMNERGRY9WKcLNkyRJERERAp9MhJiYGe/bsqdLj1qxZA5VKheHDh9u3gNVVquaGzVJERESOpXi4+fLLLxEXF4fZs2fjwIEDiIqKwoABA5CWllbp4xITEzFt2jTce++9DippNbDmhoiISDGKh5v3338fzz77LCZMmIA2bdpg2bJl8PDwwIoVKyp8jNFoxJgxY/DGG2+gadOmDixtFbHmhoiISDGKhhuDwYD9+/cjNjbWsk+tViM2Nha7du2q8HFz5sxBcHAwnn766ds+R2FhIbKysqw2u2OHYiIiIsUoGm4yMjJgNBqh1+ut9uv1eqSkpJT7mB07duDjjz/G8uXLq/Qcc+fOhY+Pj2ULCwurcblvq1SzFGtuiIiIHEvxZqnqyM7OxpNPPonly5cjMDCwSo+ZMWMGMjMzLVtycrKdSwnW3BARESnIVcknDwwMhIuLC1JTU632p6amIiQkpMzx586dQ2JiIh544AHLPpPJBABwdXXFqVOn0KxZM6vHaLVaaLVaO5S+Eqy5ISIiUoyiNTcajQadO3dGfHy8ZZ/JZEJ8fDy6d+9e5vhWrVrhyJEjSEhIsGwPPvgg7rvvPiQkJDimyakqStfc5AuFC0NERFS/KFpzAwBxcXEYN24coqOj0bVrVyxYsAC5ubmYMGECAGDs2LFo1KgR5s6dC51Oh3bt2lk93tfXFwDK7FdUqeUXCvJMAFyULQ8REVE9oni4GT16NNLT0zFr1iykpKSgY8eO2Lhxo6WTcVJSEtTqOtU1yGrhzMI8IxhuiIiIHEfxcAMAkydPxuTJk8u9b9u2bZU+duXKlbYvUE2VrrnJNSpcGCIiovqljlWJ1BEuLtCpDQCAwnyTwoUhIiKqXxhu7ETnUgzA3OeGiIiIHIXhxk60bjLUFBYw3BARETkSw42d6NxkXxsOBSciInIshhs7Kam5YbghIiJyJIYbO9FpZbjh8gtERESOxXBjJ1o3WWNTUKBSuCRERET1C8ONnei0MtwUGhQuCBERUT3DcGMn5rU6CwpZc0NERORIDDd2otPJ20IDLzEREZEj8ZPXTnTu8ragiJeYiIjIkfjJaydanWyOKmS4ISIicih+8tqJzl1e2oIirghORETkSAw3dqK9GW4KixluiIiIHInhxk50HvLSFptcUFyscGGIiIjqEYYbO9F6lNTYFBYqWBAiIqJ6huHGTnSeDDdERERKYLixE1cPDdS4uTI415ciIiJyGIYbe9FqoYNMNay5ISIichyGG3spFW5Yc0NEROQ4DDf2otNBC1llw5obIiIix2G4sRfW3BARESmC4cZeWHNDRESkCIYbe2HNDRERkSIYbuxFq7XU3DDcEBEROQ7Djb3odBwKTkREpACGG3thzQ0REZEiGG7shTU3REREimC4sRetFh7IAwDk5ChcFiIionqE4cZetFoEIgMAkJGhcFmIiIjqEYYbe9HpGG6IiIgUUK1w8+677yI/P9/y/c6dO1FYqkNJdnY2Jk6caLvS1WWla27ShcKFISIiqj+qFW5mzJiB7Oxsy/eDBg3CpUuXLN/n5eXhv//9r+1KV5fpdAhCOgCGGyIiIkeqVrgRQlT6PZVSquYmneGGiIjIYdjnxl40mlJ9blQKF4aIiKj+YLixF5UKgRrZhHf1ugomk8LlISIiqidcq/uAjz76CF5eXgCA4uJirFy5EoGBgQBg1R+HgEBtNmAATCYVbtwA/P2VLhEREZHzU4lqdJyJiIiASnX7JpYLFy7UqFD2lJWVBR8fH2RmZsLb29u+T6bXwyftNLLgg5Mngbvusu/TEREROavqfH5Xq+YmMTGxJuWqf252Ks6CDzIyGG6IiIgcgX1u7EmrLRkOzon8iIiIHKJa4WbXrl344YcfrPb973//Q2RkJIKDg/Hcc89ZTepX73GWYiIiIoerVriZM2cOjh07Zvn+yJEjePrppxEbG4vp06fj+++/x9y5c21eyDrLaq4bhctCRERUT1Qr3CQkJKBv376W79esWYOYmBgsX74ccXFx+OCDD/DVV1/ZvJB1FmtuiIiIHK5a4eb69evQ6/WW77dv345BgwZZvu/SpQuSk5NtV7q6jn1uiIiIHK5a4Uav11uGeRsMBhw4cADdunWz3J+dnQ03NzfblrAuK714JsMNERGRQ1Qr3AwePBjTp0/H77//jhkzZsDDwwP33nuv5f7Dhw+jWbNmNi9knVWqWYp9boiIiByjWvPcvPnmm3jooYfQu3dveHl5YeXKldBoNJb7V6xYgf79+9u8kHWWtzcCcQYAa26IiIgcpVrhJjAwEL/99hsyMzPh5eUFFxcXq/vXrl2LBg0a2LSAdVpQEILwBwCGGyIiIkepVrh56qmnqnTcihUrqlWIJUuWYP78+UhJSUFUVBQWLVqErl27lnvs+vXr8fbbb+Ps2bMoKipCixYt8Le//Q1PPvlktZ7TIQIDLc1SWVmAwQCUqugiIiIiO6hWuFm5ciWaNGmCu+++G9VYkqpSX375JeLi4rBs2TLExMRgwYIFGDBgAE6dOoXg4OAyx/v7+2PmzJlo1aoVNBoNfvjhB0yYMAHBwcEYMGCATcpkM0FB8MUNqGGECS7IyABCQ5UuFBERkXOr1sKZkyZNwhdffIEmTZpgwoQJeOKJJ+Bfw6WuY2Ji0KVLFyxevBgAYDKZEBYWhpdeegnTp0+v0jk6deqEIUOG4M0337ztsQ5dOPO774BhwxDsehXpxf44dAjo0MG+T0lEROSMqvP5Xa3RUkuWLMGVK1fwyiuv4Pvvv0dYWBhGjRqFTZs23VFNjsFgwP79+xEbG1tSILUasbGx2LVr120fL4RAfHw8Tp06hV69epV7TGFhIbKysqw2hwkKkjcqDgcnIiJylGovnKnVavHYY49h8+bNOH78ONq2bYuJEyciIiICOTk51TpXRkYGjEaj1cSAgJxPJyUlpcLHmTs0azQaDBkyBIsWLUK/fv3KPXbu3Lnw8fGxbGFhYdUqY40EBsobYxoAhhsiIiJHqNGq4Gq1GiqVCkIIGI1GW5Xptho0aICEhATs3bsXb731FuLi4rBt27Zyj50xYwYyMzMtm0NnUL5ZcxNoSgXAuW6IiIgcoVodigHZzLN+/XqsWLECO3bswNChQ7F48WIMHDgQanX1slJgYCBcXFyQmppqtT81NRUhISEVPk6tVqN58+YAgI4dO+LEiROYO3cu+vTpU+ZYrVYLrVZbrXLZjI8P4OqKwGI2SxERETlKtdLIxIkT0bBhQ8ybNw9Dhw5FcnIy1q5di8GDB1c72ACARqNB586dER8fb9lnMpkQHx+P7t27V/k8JpMJhYWF1X5+u1OpgMBAri9FRETkQNWquVm2bBnCw8PRtGlTbN++Hdu3by/3uPXr11f5nHFxcRg3bhyio6PRtWtXLFiwALm5uZgwYQIAYOzYsWjUqBHmzp0LQPahiY6ORrNmzVBYWIiffvoJn332GZYuXVqdl+I4gYEITGHNDRERkaNUK9yMHTsWKpXKpgUYPXo00tPTMWvWLKSkpKBjx47YuHGjpZNxUlKSVa1Qbm4uJk6ciL/++gvu7u5o1aoVPv/8c4wePdqm5bKZUhP5sc8NERGR/VVrnhtn4NB5bgBg1ChsXJuFQdiIqCggIcH+T0lERORs7DbPDd0B9rkhIiJyKIYbewsKsmqWql/1ZERERI7HcGNvpfrcGAxANec5JCIiompiuLG3oCB4IA86tRyqzqYpIiIi+2K4sbfAQKgABLlcA8BwQ0REZG8MN/ZmXoJByE7FHA5ORERkXww39mZePLNYLjHBmhsiIiL7YrixN3O44XBwIiIih2C4sTetFmjQgHPdEBEROQjDjSNwCQYiIiKHYbhxhFIT+bHmhoiIyL4YbhyhVM0Nww0REZF9Mdw4QlAQ+9wQERE5CMONI7DPDRERkcMw3DhCqT43164BRqPC5SEiInJiDDeOEBiIAFwFIFcFv35d4fIQERE5MYYbRwgKghuK4euSBYD9boiIiOyJ4cYRzLMUq+Timex3Q0REZD8MN45gXjzTlAaANTdERET2xHDjCOaaGxMXzyQiIrI3hhtH8PUFXFw41w0REZEDMNw4gkrFuW6IiIgchOHGUbgEAxERkUMw3DgKl2AgIiJyCIYbR2HNDRERkUMw3DhKqSUY2OeGiIjIfhhuHOWWDsVCKFweIiIiJ8Vw4yhBQWiESwCA3FyuL0VERGQvDDeOEhgID+QjRCMX0LxwQeHyEBEROSmGG0e5uQRDU9ckAMD580oWhoiIyHkx3DjKzSUYmprOAWC4ISIisheGG0cx19wUngDAcENERGQvDDeOEhAAAIgUrLkhIiKyJ4YbR9HpAC8vNIVMNQw3RERE9sFw40hBQZZwc/EiUFyscHmIiIicEMONIwUGIhSXoXE1wmgE/vpL6QIRERE5H4YbRwoKghoCkYHZANg0RUREZA8MN45kHg7uIyfyY7ghIiKyPYYbRzIPB3dPAcBwQ0REZA8MN450s+Ym0jUZAMMNERGRPTDcOJK55oZz3RAREdkNw40jmfvccJZiIiIiu2G4caSbNTeROUcAAFevAllZShaIiIjI+TDcONLNmhvvqxfMX+LCBQXLQ0RE5IQYbhypUSN5m52NpuFyemI2TREREdkWw40jeXoCDRsCAJoGZgJguCEiIrI1hhtHa94cAOe6ISIisheGG0dr0QIAEMnVwYmIiOyiVoSbJUuWICIiAjqdDjExMdizZ0+Fxy5fvhz33nsv/Pz84Ofnh9jY2EqPr3XMNTe5RwEw3BAREdma4uHmyy+/RFxcHGbPno0DBw4gKioKAwYMQFpaWrnHb9u2DY899hi2bt2KXbt2ISwsDP3798elS5ccXPI7ZA43V/cCABITAZNJwfIQERE5GZUQQihZgJiYGHTp0gWLFy8GAJhMJoSFheGll17C9OnTb/t4o9EIPz8/LF68GGPHjr3t8VlZWfDx8UFmZia8vb1rXP5qS0gA7r4bxYEhcL9xBcXFQHIy0Lix44tCRERUV1Tn81vRmhuDwYD9+/cjNjbWsk+tViM2Nha7du2q0jny8vJQVFQEf39/exXTtpo1AwC4ZqSgSZgRAJumiIiIbEnRcJORkQGj0Qi9Xm+1X6/XIyUlpUrnePXVVxEaGmoVkEorLCxEVlaW1aaoBg2AkBAAQNPgXAAMN0RERLakeJ+bmpg3bx7WrFmDDRs2QKfTlXvM3Llz4ePjY9nCwsIcXMpymPvdNJD9ihhuiIiIbEfRcBMYGAgXFxekpqZa7U9NTUXIzdqNirz33nuYN28efvnlF3To0KHC42bMmIHMzEzLlpycbJOy18jNcBOpugiA4YaIiMiWFA03Go0GnTt3Rnx8vGWfyWRCfHw8unfvXuHj3n33Xbz55pvYuHEjoqOjK30OrVYLb29vq01xN+e6aWo4CYDhhoiIyJZclS5AXFwcxo0bh+joaHTt2hULFixAbm4uJkyYAAAYO3YsGjVqhLlz5wIA3nnnHcyaNQurV69GRESEpW+Ol5cXvLy8FHsd1WJulrpxAAAXzyQiIrIlxcPN6NGjkZ6ejlmzZiElJQUdO3bExo0bLZ2Mk5KSoFaXVDAtXboUBoMBI0eOtDrP7Nmz8frrrzuy6HfOHG4u/Q4ASEkB8vIADw8lC0VEROQcFJ/nxtEUn+dGFgLw8QEA+PmYcCNThaNHgbZtlSkOERFRbVdn5rmpt7y9geBgAEDThvkA2O+GiIjIVhhulGJumvK5CoDhhoiIyFYYbpRyc8RUS10SAODwYSULQ0RE5DwYbpRys+amh5tc0fz335UsDBERkfNguFGKOdzkbIJKBZw5A1y5onCZiIiInADDjVJuNkv5JB5Cx45yF2tviIiIao7hRik3a26QkoJe3QwAgN9+U7A8REREToLhRik+PkBQEACgV7NLABhuiIiIbIHhRkk3a2/u9T0CADhyBLh2TckCERER1X0MN0q6GW6C0o6hdWu5a8cOBctDRETkBBhulHSzUzHOnkWvXvJLNk0RERHVDMONksydihluiIiIbIbhRknmcHPmDO69V3554ACQna1ckYiIiOo6hhslmcPNlSsI889FZCRgNAK7dilbLCIiorqM4UZJfn5AQID8+tw59O4tv9y+XbkiERER1XUMN0or1TTFfjdEREQ1x3CjNPOIqdOnLeFmzx4gP1+5IhEREdVlDDdKi4qStzt3omlTIDQUMBhkwCEiIqLqY7hRWr9+8nbrVqgMhWyaIiIiqiGGG6V16ACEhAB5ecDOnQw3RERENcRwozSVCujfX369aZMl3Pzxh2yeIiIiouphuKkNBgyQt7/8gtatSypyNm5UtlhERER1EcNNbWDud5OQAHV6Kh5/XH77v/8pVyQiIqK6iuGmNggKAjp1kl//8gvGjZNffv89cO2acsUiIiKqixhuagtz09SmTejQQY4QNxiAr75StlhERER1DcNNbWEON5s3AyYTxo6V37JpioiIqHoYbmqL7t0BLy8gLQ04dAiPPw6o1XIRzTNnlC4cERFR3cFwU1toNMB998mvN21CSEhJZc5nnylXLCIiorqG4aY2KdXvBoClaeqzzwCTSaEyERER1TEMN7WJOdzs3Ank5GDYMMDbG0hMBH7/XdGSERER1RkMN7VJ8+ZA06ZAURGwbRvc3YFHHpF3sWMxERFR1TDc1DYVNE2tXStnLSYiIqLKMdzUNqXWmQKAnj2BiAggOxv49lvlikVERFRXMNzUNvffD7i6yvHfR49CrS6pvXn7bdliRURERBVjuKltvL2BBx+UX3/0EQDg5ZeBwEDg6FHgX/9SsGxERER1AMNNbfTMM/L2s8+AggIEBADvvy93vfEGcP68ckUjIiKq7RhuaqP+/YHwcLlq5oYNAIAnngD69gUKCoCJEwEhFC4jERFRLcVwUxu5uABPPSW/Xr4cAKBSAUuXAlqt7Gu8Zo2C5SMiIqrFGG5qqwkTZKLZuhU4exYA0KIF8M9/yrunTgWuX1eueERERLUVw01tFR4ODBwov/74Y8vuV14BWreW62tOn65Q2YiIiGoxhpva7Nln5e0nn1jGgGs0wH//K3d/+CEX1SQiIroVw01tNnQooNcDqanADz9Ydt97LzBtmvx6/Hjg66+VKR4REVFtxHBTm7m5yfQCWDoWm73zjuyWYzIBjz0G/Pyz44tHRERUGzHc1HbmOW82bgSSky271WqZd0aPli1WDz0EbN+uUBmJiIhqEYab2q55c+C+++TENubONje5uMg+N0OHyvlvhg4Fdu9WqJxERES1BMNNXTBpkrx9/30gKcnqLjc3uWL4/fcDOTly/r8//1SgjERERLUEw01d8NBDQK9eQH4+EBdX5m6dTq4Y3qsXkJUlA87OnQqUk4iIqBZguKkLVCpg8WLZDvX118DmzWUO8fICfvpJtmBlZwMDBgC//aZAWYmIiBSmeLhZsmQJIiIioNPpEBMTgz179lR47LFjx/Dwww8jIiICKpUKCxYscFxBlda+PfDSS/Lrl14CDIYyh3h6yhHj/foBubnAoEFygmMiIqL6RNFw8+WXXyIuLg6zZ8/GgQMHEBUVhQEDBiAtLa3c4/Py8tC0aVPMmzcPISEhDi5tLfD663Lem1OngH//u9xDPDxkE9XAgUBenlxsMyxMzo3z5JPAa6/JGp5yshEREZFTUAmh3PrSMTEx6NKlCxYvXgwAMJlMCAsLw0svvYTpt1lbICIiAlOnTsXUqVOr9ZxZWVnw8fFBZmYmvL2977ToyvnsM2DsWJliTp6UyaUcBQUyzKxbV/5p/PyAESPkUPL77wdcXe1YZiIiohqqzue3YjU3BoMB+/fvR2xsbElh1GrExsZi165dShWr9nviCaBnT1kt87e/VXiYTidHUaWmytFTa9YA8+YBTz8tK3+uXwdWrJB9cxo1Aj79VI42JyIiqusU+389IyMDRqMRer3ear9er8fJkydt9jyFhYUoLCy0fJ+VlWWzcyvC3Lm4UyeZXj78EHjuuQoPDw6WW0xMyT6jUXY2/uorWbOTllayjMN//ws0bGj/l0FERGQvincotre5c+fCx8fHsoVV0IxTp0RFAbNny69ffBH4/vtqPdzFRY6qWroUuHwZmDtXLsj5/fdA27bAF1+wFoeIiOouxcJNYGAgXFxckJqaarU/NTXVpp2FZ8yYgczMTMuWXGoJgzrttddkG5PJJDvO3OHMfW5uwPTpwP79sjLo+nXg8ceBhx+WwYeIiKiuUSzcaDQadO7cGfHx8ZZ9JpMJ8fHx6N69u82eR6vVwtvb22pzCioVsGwZMHiwnNxv6FA5iuoOtWsn89Ebb8jOxRs2AK1by6cwmWxYbiIiIjtTtFkqLi4Oy5cvx6effooTJ07gxRdfRG5uLiZMmAAAGDt2LGbMmGE53mAwICEhAQkJCTAYDLh06RISEhJw9uxZpV6CslxdZceZrl2Bq1fl+O8rV+74dG5uwKxZshana1c52/GLL8qZj48ft2G5iYiI7EjRoeAAsHjxYsyfPx8pKSno2LEjPvjgA8Tc7P3ap08fREREYOXKlQCAxMREREZGljlH7969sW3btio9X50fCl6e9HSgRw/gzBk52d+2bYC/f41OaTQCS5YA//iHnBDQzQ145BHgqadkfx210/fWIiKi2qQ6n9+KhxtHc8pwAwDnz8sh4leuyGqXX38FGjSo8WmTkuS6nT/8ULKvSRNgwgTZN6d5c9lCRkREZE8MN5Vw2nADAMeOAb17yyaqXr2An3+Wk/3VkBDA3r3AJ58Aq1fL5iqzRo3kU5q3li0ZdoiIyPYYbirh1OEGkB1m7r9fJpCBA4FvvgG0WpudPj9fdjb+5BNg+3agqMj6/tBQueSDeWvcWC71kJ4u59PJzATuvhvw8bFZkYiIqB5guKmE04cbANixQ049nJcHPPggMGcO0KGDzatU8vLkCKvt2+X2559AqfkSAciWsexs633+/nKZrBdekH15iIiIbofhphL1ItwAss/NkCElK2Tq9XK58P795fDxgACbP2V+PvDHH0B8vHz6/ftLhpG7uABBQTJfmQd03XUX8N57sphsyiIiosow3FSi3oQbQI6amj9f3ubllezXaOQsfc89JzvK2ClZXL8u17YKCpILdarVQHEx8PHHcg7C9HR5XLdusmKpcWPZh6dxY9kn2tfXLsUiIqI6iOGmEvUq3JgVFgK7dgG//AL8+CNw+HDJfS1bAs8+K9uIvLwcVqTMTLnsw7//XVK5VJq3NzBtGjB1qk0GfRERUR3HcFOJehlubrV/P7B8ObBqFZCTI/eFhsplw8eMcegkNhcvAps2AZcuye2vv+REy4mJ8v7AQLk8xMSJchDYwYPAgQMyn0VEAHFxsraHiIicG8NNJRhuSsnJkatkzpsn58kBZHvQwoWyrUghJpOceHn2bOD0abnPza3syCxADgR74QUZgMxLkhUWArt3A1u3yhqixo2BsDC5NWnCVc+JiOoihptKMNyUo6AAWLAAeOutkpqcsWPlPj8/xYpVXAz8739yvaukJNkpuU0bOZS8XTvg22+BnTvlse7ustLp4kU5WCw/v+LzxsTIGp+HHpIrWBARUe3HcFMJhptKXLkCzJwJrFwpZ+5r1Ej2/h0wQNFiGQyyYikiAtDpSvYLIUdlvfaarKkpTa+X0/2EhsqmruRkuV26VDKCKzwceOklueZocjJw7pzcLl6Ug8latJBdklq0kOe7ehXIyCjZcnJkiCookLdFRTIsubmVbN7e8rHBwfI2JIR9iIiI7gTDTSUYbqpg1y5g3Di5VhUAPP+8HLPtwA7H1SEE8NNPsq9069Zy8sDWrcsfBJaaCvznP8DSpSWjtRwtKMg6OLVoATRtCjRrducjxISQA+IKCuRWWCg3tVoGQq1W3np6cm4hIqqbGG4qwXBTRXl5siPLokXy+6ZNZa3OqFG1NuRUR0GB7E/9wQeyX09kpAwXzZrJGqKMDLn/9GmZ8fLy5EoWQUGyk3NAgKyV0elkk5hOJ0OD0ShrcMzbjRtyZubUVLmZW/0q4u8vyxIQIGdx9vaWm7+/7CsUGiq3wEDg5ElZY2XeUlNv/7rd3OTCpyNGAMOGVd7/KDMTSEiQnbjz8uSxpTd/f2Wa9QwGOc2Avz+DGtV+QsgBEL/+Kmfh8POT/8T4+cm1+fR6pUtYdzDcVILhppq2bJGrZCYlye+9vIBHHwWeflp2XnGC2feEqPxlCCHDkLt7zZ8rOxs4e7YkNJ0+Lb8/f75q4aQqXFxk2NJoZBNcYaEs/61UKtlvvF07eZzJJMNZTg5w6JBsorudBg3kH2nz5u0t95lvi4tlk11eXkk/qODgkma64GC5kH3TpmV/BkajnBRywwYZ5C5fllvpGrfAwJKw1aiRbGps0kTehofL83t7lz13QYFshU1JkWHJ/FdQCPm85povg0FuXl4y2Jo3X1/52sz3Gwxyny3eI7aUlwf89huwebO8DQwERo+W4fbWJVCuXpV92PLySq5jw4ZlB0+ar1V1f/Vv93t2q8REOU1XfDzQqhVwzz1y69y54uuckSHD+MGD8vliY2UfvVtfw/XrcvqvpCT5HgkJkVtwsAz1ycnyvuRkuZJNeLj8pyMyUv7zk5cnf3/PnpW3aWnyH6O2beXWuLFsDl+9Gvj8c+Do0fLL6+YmZ+KYOVP+01JdQsjRpVeuAJ063dmyNkLIf44+/VQGME9P+R4PDpa3LVvKpvvw8Oqf29YYbirBcHMHsrJkO87HH5c0VQHyt7l3b+Dee+VW3icUVVlODnDhgvyjfv26vOxZWfKPbUaG/AN2+bK8TU+XHz4xMXKAW0yMDAkeHuXXpggha5IuXJDLjW3YULafUnnCw+UfTV9f+bzmzdZNemFhQJ8+cmvUCPj+e+Drr2X4qClXV1nLExgo356XL8vra2sqlbxeLVvK2bebNpXPbTLJ628yyQDYtq3sGH9rBajJJH/OV6/K90JubsmtORyag6KLi/yQbdpUbkFB8tiTJ4Hjx+W2d68MK+XNI6XVyg+sfv1krcJvv5X/AezmJj90TSbrkOruLl+n+bU2aybLfvZsyZaaKgOg0VjSzy0yEujRQ4aUHj3ktXBxsX7O48flAM7Vq+VjyytTRIT8EDZvarVcNzg5uezxer3sNnjfffKfiVtnT7eHW5ed0Wjk5PDu7vK9d+OGvF7mKS90OmDyZODVV+X7tDzFxfK9e/68/N3duVOG/6tX5f1qtfxd7dNH/llu0ED2MTT/U3DjRsk/AyEh8rqYQ82pU7d/TR07ytreBx6QNcvm8F9YCFy7VtJn0by1by9rx22J4aYSDDc1IATw++8y5KxdW3ZIUsOG8reqTx/5l6RFC4YdOzGZaj4d0aVLsp9Sero8l3nTaktGpVW0Soe5ye369ZLtxg35B90cynJy5Ie7u3vJJoT8L9e8Xb4sa4nKG+YPyP9Ehw+X2blRo5JmOT8/+Qe1dOD66y/53/bFiyX/defmVvz6tVr5R16nk29T8+biIu8zb25u8nWVXvz1Vmp19T8smzSRvyJZWfI6pKTID7A74e5e8QjB8HAZYvr2lR+Mq1YBJ06Uf2zr1vID8OJF+f4oL1zYkkYjf8YNGsjN1VUGD7N+/eRE6hcvyg/ynTtvX8PZvLn8kC8okLU+Fb0H7rpL1lpmZMhzpqTI97BOVzJ1RHi4LFdSkvzH4Pz5kqblRo3kc7VoIcPlmTMyYJ0+XXLdevcGnnhCTghf3sDTbdtkrc0ff8jvvbxkYNRq5bXRamV4uHhRvp/L+3nodLKWxVy5fifc3WUZH31U/gzM7/XUVNkFc+fO6r+/27UDjhy58zKVh+GmEgw3NpKVJcdc//abDDx795b9hGrYEIiOLmmv8PGRVQADB8p3PhFkbcAff8g/9OamgthY4JFH5AeyRnPn5y4okP/Zmke6mUwlfZd8fe8sexsM8u2v0cjNzU2Gm4wM+R+webt4UR5vDo0qlfwAPXas4hoplUqWy8tLbuaaCQ8PuZlDosFQ8mF76VJJU1FIiAymbdrIX7H775cfwKVfp7kPyOrVwJ49cumTXr1kgAwOLjmuuLjkP383N+uQmp1d8jpPn5blCAyUz9W8uazJadRIPs7FRW7m5zXXOPz5Z/l90FQq2Ww2Y4b881GaELK249Il65otg0GGgo4d5Z+a0j+rHTuAn3+Wz9msmXxP9e0rm47K+9m6uVX8vhBChmpz5/zyGAwy6Pj6Vm2CUSFk+f75T9mcVhlXVxm6OnaUNV89e8p/QjQaGe63b5e/Qzt2yPOa/xkwv9/NNcDmrVEjGb5GjrS+brdKT5f/CH37reypUFxc8v7XauVjzYMizFuLFvK9YEsMN5VguLGT/PySJcK3bZNxv7z6cED+5Rg1Si4N3qqVI0tJVCtcvSqbX86dk9k/NFSGLr2++p2kCwvlf/X+/nKrK4qLZblzcmRYysmRW/v28oOxvjGZ5P+I165Z9/lydbXuA3VrM159wnBTCYYbB8nPlw26p06VdBzJypJ/zX/6SR6jVsuZ92bOlHXEREREFWC4qQTDTS1w6JBcW+Hbb0v26fWyfrVjR3nbuzfHSBIRkQXDTSUYbmqR/ftlyPnpp5JOA2Zqtex4MWaMbIDntL5ERPUaw00lGG5qodxc2a3ePGPc3r3WPevc3WXQ0emsZ8iLiACGDJG9Ayvq3UdERE6B4aYSDDd1xNmzcjjHqlUlS4NXRKuVQ8/79ZMBqHRvPJNJ9sBzdS2Z3e7++zlai4iojmG4qQTDTR0jhGy++uMP2VRlXpHSxUXu/+EHOSa2ujp3BsaPBx57zHoyF/OMVEFBXDKciKgWYbipBMONkxFCTsv6ww9y+Ll5FjrzLFhqtZz5yjxValqanIvePFuamxsQFSVnoUtPlyO6AFnD06GDnA2sUyc5iYZKVTLVrMkkm9Nu3CjZCgpkJ2jzbHPmWy6ARERUYww3lWC4IaSnA198AaxceftZs2rKPE++eXazyEg545V5ZjZPTzk5SePGsraoptMOExE5KYabSjDckJWjR2X/nsDAkpXivL1lU9eBAyVbYmLJNLPmWy8vOe2nedNo5NSzly+XTO1a0USG5dFoZG1P48Yy8Pj4lMzq7Osrm8/Mm7+/LEdxsdyKimTNUem1EMxzyZuPDwiQNUtNm9bvmcCIqE5iuKkEww05jMkkw07p1QQTE2VzVuktI0Me56hfRQ8POZ9Qp05yTiHzKoTm+f7VajnZonluffNEjKUXXNLpSlajNG8NG8raqdBQrilGRDbHcFMJhhuqlYqKZE1PcrKs9blxQ87qnJlZUgtjXiTp6lXZ6RmQnZ7Nm1Yr5/L39ZW3Pj6yg3Tpx12+XPEKi7bi7l6yoqC/v/zevDiSRmM9nL+4WAYj8/LWTZvK8ufmlqzel54ur4l5RcyLF+U+Hx9ZG1U6YAUFlWyBgfIYb29ZhjsNXEVF8vpW9njzn1GGOiK7qc7nN4eDENUGbm5y8ZgmTez7PEajrI0xN7cdPCgDhHkFQvMqhBERckkM8xYQUDLEvrBQBqTr12Wtk3n76y/ZnJefL+ctutMlgd3cKl4m/E6p1TLkeHjI85tXvDSvegmUBJPiYrnYkXmJc4NB1lSV7iTu4yOXTL5ypWRJb51O9qkyb2FhJc1/5vDj6yvXU2vVqvxlos3Pf+WKDLrJyfLc5tUSIyOrHqAyMuTP1/xzPnRIvv5OneRowU6d5EJO5mXRiZwIa26IyJoQNavlSEyUyyKfOyfDQV5eyWYwWAcLFxcZEs6fl8enppacS6cr6QfVqJFcPdC86fXy3Oblvs1benrJlpEhj6mtf+KCg2UNl8lUsmpkTo4MjUZj+Y/x9paj+5o1k6+rdJ+rzEz5ms3XpKCg6mW5tS9ZcHDJFhAgn8tgkFtRkfzZBQfLn4NeL39GJpMMtnl58tZgkOczn1utLlkl1LyVN/O4EDIo794tly0/fVqGydK1cuHhMiA2alT2vWoyycCuVlsvc16a0ShHWd64IUdF1mQGdINBBtAGDWQ5Sw8KMK9qmpQkA6uLi6xF1OnkbViYDKxUJWyWqgTDDVEtZu6DFBAg+//UtEbBZJIftllZsiYmL09+OJs/pIuK5IepeQPkB1CDBnIzj2zLzJQ1NObO4pmZ8kPdvJx3w4ay7BculGyXLpUERfPrSEuTH6p//VV5uV1dSwKdXi/Pd+RI9TqoA7Jp0Ny3qmNHGZ4OHJBzRO3fL6+1ktzd5TVu0KCkz9f58zIIVIWnpww5EREy1CUlyWtrvk7BwTK8dOggjzt3Tgamffvk+wGQP5tWrYAuXYDoaFkOcxg312RqNNbb5cvAiRPyZ3nuXEkYValkwPHzkwEvJeX2r+Guu4AHHpDbPfdUfX6t/Hx5na5ckWVt2VKGpTsZcSmEfL+ePSu/d3WVAdbVFQgJke/FWoDhphIMN0SkuOxsWSNx/rz8sPTyKtn8/GSguXVEW1GR7Nx98KD8ADfPvG3ezH2QzP2QgoMrX5ZECFlLVFxcMn+T0SgDUFqa3FJTZWhwcSmpcdNoZI1EamrJlp4uy2DuX+XuLo8zP495bihzv6/Ll0vmlCqPi4tsMouJkbc5OSU1cmlpMuydPVtxDZc5TFb28ebpKa/17YJmVbi6lsyddSt3d9ncHBoqy5OfL2vV8vNlMCr9OF9fOVrSPE2Ep6c8d+mwlZsrr8GNG+W/ptat5abTlfTby8yUz+fnZ91H7do1GZqPHpXvhYqEhwM9ewI9egDduskAdf26fPz167Js5jBkrpXV6+WyOTbEcFMJhhsiolrAHKLM/ZtycuRtw4aytsnDo/LHGwwyHJ48KTuZBwaWNFuGhsoweOwYcPiw7G906pQMGTExQNeuQJs2Jc2i+/aVrGlnNJZ0gC/dCd68pEthoXwuc4ho1ark+Up/4Gu18vkCAiqugczMBH75Bfj+e7mA8NWr1buG5r5gWq0MezXpq+biIjv0u7qWdPYvKpK1TxWFyMp06yYnVrUhhptKMNwQEVGtYzTKEHb9uqydMdfUFBeXTPpp3oKCZAj08SkJTkVFMuwdPy6bzIxG67mydLqSQQDmPmmenrJmrH17GdJ0urLlysmR/Z927JDb/v0y8Pn7y5ogf39ZJnMYMt+2bg188IFNLxHDTSUYboiIiOqe6nx+c653IiIicioMN0RERORUGG6IiIjIqTDcEBERkVNhuCEiIiKnwnBDREREToXhhoiIiJwKww0RERE5FYYbIiIicioMN0RERORUGG6IiIjIqTDcEBERkVNhuCEiIiKnwnBDRERETsVV6QI4mhACgFw6nYiIiOoG8+e2+XO8MvUu3GRnZwMAwsLCFC4JERERVVd2djZ8fHwqPUYlqhKBnIjJZMLly5fRoEEDqFQqm547KysLYWFhSE5Ohre3t03PTdZ4rR2H19pxeK0dh9facWx1rYUQyM7ORmhoKNTqynvV1LuaG7VajcaNG9v1Oby9vfnL4iC81o7Da+04vNaOw2vtOLa41rersTFjh2IiIiJyKgw3RERE5FQYbmxIq9Vi9uzZ0Gq1ShfF6fFaOw6vtePwWjsOr7XjKHGt612HYiIiInJurLkhIiIip8JwQ0RERE6F4YaIiIicCsMNERERORWGGxtZsmQJIiIioNPpEBMTgz179ihdpDpv7ty56NKlCxo0aIDg4GAMHz4cp06dsjqmoKAAkyZNQkBAALy8vPDwww8jNTVVoRI7j3nz5kGlUmHq1KmWfbzWtnPp0iU88cQTCAgIgLu7O9q3b499+/ZZ7hdCYNasWWjYsCHc3d0RGxuLM2fOKFjiusloNOK1115DZGQk3N3d0axZM7z55ptWaxPxWt+53377DQ888ABCQ0OhUqnwzTffWN1flWt77do1jBkzBt7e3vD19cXTTz+NnJycmhdOUI2tWbNGaDQasWLFCnHs2DHx7LPPCl9fX5Gamqp00eq0AQMGiE8++UQcPXpUJCQkiMGDB4vw8HCRk5NjOeaFF14QYWFhIj4+Xuzbt09069ZN3HPPPQqWuu7bs2ePiIiIEB06dBBTpkyx7Oe1to1r166JJk2aiPHjx4vdu3eL8+fPi02bNomzZ89ajpk3b57w8fER33zzjTh06JB48MEHRWRkpMjPz1ew5HXPW2+9JQICAsQPP/wgLly4INauXSu8vLzEwoULLcfwWt+5n376ScycOVOsX79eABAbNmywur8q13bgwIEiKipK/Pnnn+L3338XzZs3F4899liNy8ZwYwNdu3YVkyZNsnxvNBpFaGiomDt3roKlcj5paWkCgNi+fbsQQogbN24INzc3sXbtWssxJ06cEADErl27lCpmnZadnS1atGghNm/eLHr37m0JN7zWtvPqq6+Knj17Vni/yWQSISEhYv78+ZZ9N27cEFqtVnzxxReOKKLTGDJkiHjqqaes9j300ENizJgxQghea1u6NdxU5doeP35cABB79+61HPPzzz8LlUolLl26VKPysFmqhgwGA/bv34/Y2FjLPrVajdjYWOzatUvBkjmfzMxMAIC/vz8AYP/+/SgqKrK69q1atUJ4eDiv/R2aNGkShgwZYnVNAV5rW/ruu+8QHR2NRx55BMHBwbj77ruxfPlyy/0XLlxASkqK1bX28fFBTEwMr3U13XPPPYiPj8fp06cBAIcOHcKOHTswaNAgALzW9lSVa7tr1y74+voiOjrackxsbCzUajV2795do+evdwtn2lpGRgaMRiP0er3Vfr1ej5MnTypUKudjMpkwdepU9OjRA+3atQMApKSkQKPRwNfX1+pYvV6PlJQUBUpZt61ZswYHDhzA3r17y9zHa20758+fx9KlSxEXF4d//OMf2Lt3L15++WVoNBqMGzfOcj3L+5vCa10906dPR1ZWFlq1agUXFxcYjUa89dZbGDNmDADwWttRVa5tSkoKgoODre53dXWFv79/ja8/ww3VCZMmTcLRo0exY8cOpYvilJKTkzFlyhRs3rwZOp1O6eI4NZPJhOjoaLz99tsAgLvvvhtHjx7FsmXLMG7cOIVL51y++uorrFq1CqtXr0bbtm2RkJCAqVOnIjQ0lNfaybFZqoYCAwPh4uJSZtRIamoqQkJCFCqVc5k8eTJ++OEHbN26FY0bN7bsDwkJgcFgwI0bN6yO57Wvvv379yMtLQ2dOnWCq6srXF1dsX37dnzwwQdwdXWFXq/ntbaRhg0bok2bNlb7WrdujaSkJACwXE/+Tam5v//975g+fToeffRRtG/fHk8++ST+3//7f5g7dy4AXmt7qsq1DQkJQVpamtX9xcXFuHbtWo2vP8NNDWk0GnTu3Bnx8fGWfSaTCfHx8ejevbuCJav7hBCYPHkyNmzYgC1btiAyMtLq/s6dO8PNzc3q2p86dQpJSUm89tXUt29fHDlyBAkJCZYtOjoaY8aMsXzNa20bPXr0KDOlwenTp9GkSRMAQGRkJEJCQqyudVZWFnbv3s1rXU15eXlQq60/5lxcXGAymQDwWttTVa5t9+7dcePGDezfv99yzJYtW2AymRATE1OzAtSoOzIJIeRQcK1WK1auXCmOHz8unnvuOeHr6ytSUlKULlqd9uKLLwofHx+xbds2ceXKFcuWl5dnOeaFF14Q4eHhYsuWLWLfvn2ie/fuonv37gqW2nmUHi0lBK+1rezZs0e4urqKt956S5w5c0asWrVKeHh4iM8//9xyzLx584Svr6/49ttvxeHDh8WwYcM4PPkOjBs3TjRq1MgyFHz9+vUiMDBQvPLKK5ZjeK3vXHZ2tjh48KA4ePCgACDef/99cfDgQXHx4kUhRNWu7cCBA8Xdd98tdu/eLXbs2CFatGjBoeC1yaJFi0R4eLjQaDSia9eu4s8//1S6SHUegHK3Tz75xHJMfn6+mDhxovDz8xMeHh5ixIgR4sqVK8oV2oncGm54rW3n+++/F+3atRNarVa0atVKfPjhh1b3m0wm8dprrwm9Xi+0Wq3o27evOHXqlEKlrbuysrLElClTRHh4uNDpdKJp06Zi5syZorCw0HIMr/Wd27p1a7l/o8eNGyeEqNq1vXr1qnjssceEl5eX8Pb2FhMmTBDZ2dk1LptKiFJTNRIRERHVcexzQ0RERE6F4YaIiIicCsMNERERORWGGyIiInIqDDdERETkVBhuiIiIyKkw3BAREZFTYbghonpPpVLhm2++UboYRGQjDDdEpKjx48dDpVKV2QYOHKh00YiojnJVugBERAMHDsQnn3xitU+r1SpUGiKq61hzQ0SK02q1CAkJsdr8/PwAyCajpUuXYtCgQXB3d0fTpk2xbt06q8cfOXIE999/P9zd3REQEIDnnnsOOTk5VsesWLECbdu2hVarRcOGDTF58mSr+zMyMjBixAh4eHigRYsW+O677+z7oonIbhhuiKjWe+211/Dwww/j0KFDGDNmDB599FGcOHECAJCbm4sBAwbAz88Pe/fuxdq1a/Hrr79ahZelS5di0qRJeO6553DkyBF89913aN68udVzvPHGGxg1ahQOHz6MwYMHY8yYMbh27ZpDXycR2UiNl94kIqqBcePGCRcXF+Hp6Wm1vfXWW0IIuTr8Cy+8YPWYmJgY8eKLLwohhPjwww+Fn5+fyMnJsdz/448/CrVaLVJSUoQQQoSGhoqZM2dWWAYA4p///Kfl+5ycHAFA/PzzzzZ7nUTkOOxzQ0SKu++++7B06VKrff7+/pavu3fvbnVf9+7dkZCQAAA4ceIEoqKi4Onpabm/R48eMJlMOHXqFFQqFS5fvoy+fftWWoYOHTpYvvb09IS3tzfS0tLu9CURkYIYbohIcZ6enmWaiWzF3d29Sse5ublZfa9SqWAymexRJCKyM/a5IaJa788//yzzfevWrQEArVu3xqFDh5Cbm2u5f+fOnVCr1bjrrrvQoEEDREREID4+3qFlJiLlsOaGiBRXWFiIlJQUq32urq4IDAwEAKxduxbR0dHo2bMnVq1ahT179uDjjz8GAIwZMwazZ8/GuHHj8PrrryM9PR0vvfQSnnzySej1egDA66+/jhdeeAHBwcEYNGgQsrOzsXPnTrz00kuOfaFE5BAMN0SkuI0bN6Jhw4ZW++666y6cPHkSgBzJtGbNGkycOBENGzbEF198gTZt2gAAPDw8sGnTJkyZMgVdunSBh4cHHn74Ybz//vuWc40bNw4FBQX497//jWnTpiEwMBAjR4503AskIodSCSGE0oUgIqqISqXChg0bMHz4cKWLQkR1BPvcEBERkVNhuCEiIiKnwj43RFSrseWciKqLNTdERETkVBhuiIiIyKkw3BAREZFTYbghIiIip8JwQ0RERE6F4YaIiIicCsMNERERORWGGyIiInIqDDdERETkVP4/g9UU6raNFKgAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Plot the results\n", + "fig, ax = plt.subplots()\n", + "ax.plot(errors_train,'r-',label='train')\n", + "ax.plot(errors_test,'b-',label='test')\n", + "#ax.set_ylim(0,100); ax.set_xlim(0,n_epoch)\n", + "ax.set_xlabel('Epoch'); ax.set_ylabel('MSE')\n", + "ax.set_title('Train Loss %3.5f, Test Loss %3.5f'%(errors_train[-1],errors_test[-1]))\n", + "ax.legend()\n", + "plt.show()" + ] + } + ], + "metadata": { + "accelerator": "GPU", + "colab": { + "gpuType": "T4", + "provenance": [] + }, + "jupytext": { + "main_language": "python" + }, + "kernelspec": { + "display_name": "Python (mzmvenv)", + "language": "python", + "name": "mzmvenv" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.14.3" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/README.md b/README.md index cfc9aaa..eaf3130 100644 --- a/README.md +++ b/README.md @@ -1,163 +1,146 @@ -# MZM 器件性能 MLP 回归基线(PyTorch) +# MZM MoE PINN(PyTorch) -本项目实现一个**多输出回归**基线模型:用 **8 个器件/偏置参数**预测 **3 个射频性能指标**。代码面向科研复现:可配置 YAML、固定随机种子、训练集拟合标准化器、完整日志与可视化产物。 +本仓库已从原先的 MLP 基线迁移为 **Mixture-of-Experts + Physics-Informed Neural Network** 训练流程,目标是尽量对齐 `MZM_MoE_PINN_Model.ipynb` 的训练行为,同时保留仓库化的命令行入口、结果目录和复现产物。 -## 项目简介 +## 当前训练管线 -- **任务类型**:监督学习,多输出回归(非分类)。 -- **输入(8 维)**:工艺与偏置相关参数。 -- **输出(3 维)**:`BW_3dB`、`IL`、`V_pi`。 -- **模型**:原生 PyTorch MLP,可选 BatchNorm / Dropout / 残差(同维时相加)。 -- **损失**:默认在**标准化后的输出空间**使用加权 `SmoothL1Loss`(Huber);可选加权 MSE。 -- **v1 目标**:先把数据清洗、划分、训练、评估、日志与可视化流程跑通;**不引入 physics loss**。 +- **输入**:8 个器件/偏置参数 +- **输出**:`BW_3dB`、`IL`、`V_pi` +- **模型**:MoE,包含多个专家网络与一个 gating 网络 +- **数据清洗**:保留 `V_pi < 500` +- **数据划分**:`train_test_split(test_size=0.1, random_state=123)` +- **标准化**: + - `X` 使用一个 `StandardScaler` + - `Y` 的三个目标分别使用独立的 `StandardScaler` +- **损失**: + - 数据项:标准化空间 `MSE` + - 物理项:`dBW/dL <= 0`、`dIL/dL >= 0`、`d(V_pi*L)/dL ~= 0`、`d2BW/dL2` 平滑项 +- **优化器**:`AdamW(lr=1e-3, weight_decay=0.05, betas=(0.9, 0.999))` +- **训练方式**:固定 `100` epoch,无早停;每轮计算全量 train/test MSE + +默认物理约束权重来自 `best_hyperparams.json`: + +```json +{ + "lambda_bw_mon": 0.0, + "lambda_IL_mon": 0.3, + "lambda_vpiL": 0.005, + "lambda_smooth": 0.1 +} +``` ## 数据格式 -数据为 **txt 或 csv**,每行 **11 个逗号分隔的浮点数**,无表头(txt)或表头与下列字段一致(csv)。 +数据文件为 11 列逗号分隔浮点数,列含义如下: -| 顺序 | 列名 | 含义 | 作为 | -| --- | --- | --- | --- | -| 1 | `PN_offset` | PN 偏移 | 输入 | -| 2 | `Bias_V` | 偏置电压 | 输入 | -| 3 | `Core_width` | 芯区宽度 | 输入 | -| 4 | `P+_width` | P+ 区宽度 | 输入 | -| 5 | `N+_width` | N+ 区宽度 | 输入 | -| 6 | `P_width` | P 区宽度 | 输入 | -| 7 | `N_width` | N 区宽度 | 输入 | -| 8 | `Phase_length` | 相位区长度 | 输入 | -| 9 | `BW_3dB` | 3 dB 带宽 | 目标 | -| 10 | `IL` | 插入损耗 | 目标 | -| 11 | `V_pi` | 半波电压 | 目标 | +| 顺序 | 列名 | 作为 | +| --- | --- | --- | +| 1 | `PN_offset` | 输入 | +| 2 | `Bias_V` | 输入 | +| 3 | `Core_width` | 输入 | +| 4 | `P+_width` | 输入 | +| 5 | `N+_width` | 输入 | +| 6 | `P_width` | 输入 | +| 7 | `N_width` | 输入 | +| 8 | `Phase_length` | 输入 | +| 9 | `BW_3dB` | 输出 | +| 10 | `IL` | 输出 | +| 11 | `V_pi` | 输出 | -- 自动忽略空行与行首行尾空格。 -- 每行必须恰好 **11 列**;也支持仿真导出的 **整行方括号** 写法,例如 `[a, b, ..., k]`(与无括号的 `a, b, ..., k` 等价)。 -- 否则整文件解析失败并给出错误行号提示。 +支持两种文本格式: -## TXT 数据清洗流程(以 V_pi 为准) +```text +a,b,c,...,k +[a, b, c, ..., k] +``` -本仓库约定:**txt 每行从左到右第 11 个逗号分隔浮点数**即半波电压 **`V_pi`**(与表头列名一致)。清洗时以该列为**物理可信区间**的主门控,避免异常仿真/标注污染训练。 - -建议按以下顺序理解流水线(与 `src/preprocess.py` 中 `clean_dataframe` 实现一致): - -1. **解析与建表**:读取 txt → 校验每行 11 列 → 转为 `float` → 构建 `DataFrame`(最后一列为 `V_pi`)。 -2. **(可选)去重**:`remove_duplicate_rows: true` 时删除 11 列完全相同的重复行。 -3. **V_pi 区间门控(主清洗)**:默认启用 `filter_v_pi_range: true`,仅保留 - `v_pi_min <= V_pi <= v_pi_max`(默认 **`[0, 500]`**)。**区间之外整行剔除**。 - 该步骤专门针对「以最后一列 `V_pi` 为正常范围」的需求。 -4. **(可选)严格正电压**:`remove_nonpositive_vpi: true` 时,在区间过滤之后再删除 `V_pi <= 0`(若需保留 `V_pi = 0` 且仍在 `[0,500]` 内,请保持为 `false`)。 -5. **后续步骤**:默认采用**按 8 个输入字段分组**的 train/val/test 切分,避免「同输入异输出」同时落入不同集合;再按 `split_stratify_target`(默认 `V_pi`)做组级近似分层;之后才做(可选)训练集离群策略与仅在训练集上拟合 `StandardScaler`。 - -清洗前会在日志与 `data_report.md` 中报告:给定 `[v_pi_min, v_pi_max]` 下 **`V_pi` 越界行数**、重复样本、同输入异输出等统计,便于核对。 - -## 环境要求 - -- Python **3.10+**(已在 3.13 下通过冒烟测试)。 -- 推荐使用虚拟环境。 - -### 安装依赖 +## 安装 ```bash cd /path/to/photonAI python -m venv .venv -source .venv/bin/activate # Windows 使用 .venv\Scripts\activate +source .venv/bin/activate pip install -U pip pip install -r requirements.txt ``` -## 放置数据 - -1. 将原始 txt(例如仓库根目录下的 `Sim_MZM_dataset.txt`)复制或软链接到 `data/dataset.txt`。 -2. 或在 `configs/default.yaml` 中修改 `data_path` 为绝对路径或相对项目根目录的路径。 - -若路径不存在,程序会给出明确报错,不会静默失败。 - ## 训练 ```bash python -m src.main train --config configs/default.yaml ``` -或使用脚本: +训练完成后会在 `results/run_时间戳/` 下生成: -```bash -bash scripts/train.sh -``` +- `config_snapshot.yaml` +- `data_report.md` +- `data_stats.csv` +- `cleaning_meta.json` +- `split_indices.json` +- `x_scaler.pkl` +- `y_scalers.pkl` +- `train_log.csv` +- `checkpoints/best.pt` +- `checkpoints/last.pt` +- `metrics.csv` +- `summary.json` +- `summary.md` +- `test_predictions.csv` +- `figures/*.png` -训练会在 `results/run_时间戳/` 下生成: +说明: -- `config_snapshot.yaml`:本次运行配置快照。 -- `split_indices.json`:对**清洗后**样本行的 train/val/test 索引,便于 `eval` 完全复现划分。 -- `x_scaler.pkl` / `y_scaler.pkl`:`StandardScaler`,推理阶段用于反标准化。 -- `data_report.md` / `data_stats.csv`:数据统计与清洗说明。 -- `cleaning_meta.json`:清洗与划分元信息。 -- `train_log.csv`:逐 epoch 的 train/val loss 与学习率。 -- `checkpoints/best.pt`、`checkpoints/last.pt`:最优与最后一轮权重。 -- 训练结束后:`metrics.csv`、`summary.json`、`summary.md`、`test_predictions.csv`、`figures/*.png`。 +- `train_log.csv` 记录每轮的全量 `train_loss` / `test_loss` +- `summary.*` 与 `metrics.csv` 中的 `loss` 为**标准化空间 MSE** +- 物理空间指标仍输出 `MAE / RMSE / R²` -**说明(损失列)**:`metrics.csv` / `summary.*` 中的 `loss` 与 `*_loss` 均在**标准化输出空间**按训练准则(Huber / 加权 MSE)计算;物理量空间以 **MAE / RMSE / R²** 为主指标。 +## 评估 -## 评估(复现划分与 scaler) - -在**同一数据文件**与 `config_snapshot.yaml` 前提下,可仅运行评估: +按训练时保存的切分索引与 scaler 重算 train/test 指标: ```bash python -m src.main eval --config configs/default.yaml --run-dir results/run_YYYYMMDD_HHMMSS ``` -若不指定 `--run-dir`,将在 `configs/default.yaml` 的 `output_dir`(默认 `results`)下自动选择**最近修改时间**的 `run_*` 目录。 - -```bash -bash scripts/eval.sh --run-dir results/run_某次训练 -``` - ## 推理 -输入文件需包含上述 **8 个输入列**(csv 带表头,或 8 列无表头 txt)。 +输入文件需包含 8 个输入列(csv 带表头,或 8 列 txt): ```bash python -m src.main infer --config configs/default.yaml --input path/to/inputs.csv --output path/to/preds.csv ``` -脚本封装: +输出列为原始 8 个输入 + `pred_BW_3dB`、`pred_IL`、`pred_V_pi`。 -```bash -bash scripts/infer.sh path/to/inputs.csv --run-dir results/run_某次训练 --output preds.csv -``` +## 默认配置 -输出列为 8 个输入 + `pred_BW_3dB`、`pred_IL`、`pred_V_pi`(**物理量空间**,已反标准化)。 +`configs/default.yaml` 目前对应 notebook 风格的默认 MoE PINN 参数: + +- `data.test_size: 0.1` +- `data.random_state: 123` +- `data.filter_v_pi_max: 500.0` +- `model.hidden_dims: [64, 128, 64]` +- `model.n_experts: 60` +- `model.gating_hidden: 8` +- `model.dropout_rate: 0.0` +- `model.use_bn: true` +- `optimizer.lr: 0.001` +- `optimizer.weight_decay: 0.05` +- `training.batch_size: 128` +- `training.epochs: 100` +- `physics.*` 默认由 `best_hyperparams.json` 提供,再由 YAML 显式值覆盖 ## 测试 ```bash -pip install pytest -pytest -q tests/test_smoke.py +PYTHONPATH=. PYTEST_DISABLE_PLUGIN_AUTOLOAD=1 pytest -q tests/test_smoke.py ``` -## 配置说明(`configs/default.yaml`) +## 说明 -主要字段: - -- **数据与清洗**:`data_path`、`remove_duplicate_rows`、**`filter_v_pi_range` / `v_pi_min` / `v_pi_max`**(默认按 **`V_pi ∈ [0, 500]`** 剔除越界行,对应 txt **第 11 列**)、`remove_nonpositive_vpi`、`outlier_strategy`(`none` / `iqr` / `zscore` / `quantile_clip`)及 `outlier_apply_to`(`targets` / `all`)。 -- **划分**:`split_ratios`、`random_seed`、`split_mode`、`split_stratify_target`、`split_stratify_bins`。默认 `grouped_stratified`:先按 8 维输入分组,再按指定目标(默认 `V_pi`)做组级近似分层;也可切回 `random`。**仅在训练子集**上拟合标准化器;离群阈值(若启用)也在训练子集上统计。 -- **模型**:`hidden_dims`、`batchnorm`、`dropout`、`residual`。 -- **训练**:`AdamW`、`lr`、`weight_decay`、`batch_size`、`epochs`、早停 `early_stopping_patience`。 -- **调度器**:`cosine`(默认)或 `plateau`。 -- **损失**:`huber`(默认)或 `weighted_mse`,`target_weights` 长度须为 3。 - -默认策略刻意**不删除**仅因统计极端的样本(`outlier_strategy: none`),但在报告中给出极端值计数;**默认以 `V_pi` 物理区间 `[0,500]` 删除越界行**;`remove_nonpositive_vpi` 默认为 `false`,以便与「0 属于合法下界」一致,需要时可改为 `true`。 - -## 项目结构 - -```text -. -├── README.md -├── requirements.txt -├── .gitignore -├── configs -│ └── default.yaml -├── data -├── reports -├── results +- 现在的主流程优先保证与 notebook 的 **数据切分、标准化、模型结构、物理损失和训练循环** 一致。 +- 为了适配仓库化使用,仍保留了 `train / eval / infer` CLI 与 `run_*` 结果目录结构。 +- 旧的 MLP baseline 文档与配置已不再是当前默认路径。 ├── scripts │ ├── train.sh │ ├── eval.sh diff --git a/best_hyperparams.json b/best_hyperparams.json new file mode 100644 index 0000000..13856cd --- /dev/null +++ b/best_hyperparams.json @@ -0,0 +1,11 @@ +{ + "best_config": { + "lambda_bw_mon": 0.0, + "lambda_IL_mon": 0.3, + "lambda_vpiL": 0.005, + "lambda_smooth": 0.1 + }, + "final_train_loss": 0.011141298338770866, + "final_test_loss": 0.013557782396674156, + "n_params": 619503 +} \ No newline at end of file diff --git a/configs/default.yaml b/configs/default.yaml index 6c645f4..7f687a4 100644 --- a/configs/default.yaml +++ b/configs/default.yaml @@ -1,70 +1,39 @@ -# 默认配置:MZM MLP 多输出回归基线 -# 将数据 txt 放到 data/ 下并修改 data_path,或保持路径指向你的文件 +# 默认配置:MZM MoE PINN(按 notebook 迁移) data_path: data/dataset.txt +best_hyperparams_path: best_hyperparams.json -split_ratios: [0.7, 0.15, 0.15] # train, val, test;可改为 [0.8, 0.1, 0.1] -random_seed: 42 -# 切分策略:按 8 维输入分组,避免“同输入异输出”跨集合泄漏;再按目标分桶近似分层 -split_mode: grouped_stratified # grouped_stratified | random -split_stratify_target: V_pi -split_stratify_bins: 10 - -remove_duplicate_rows: true - -# 异常值处理策略:none | iqr | zscore | quantile_clip -# 默认仅报告极端值,不删除;物理上不可信的 V_pi 由下方区间门控剔除 -outlier_strategy: none -# 启用非 none 策略时,在训练子集上拟合阈值;iqr/zscore 仅删训练集离群行;quantile_clip 按训练分位数 winsorize -outlier_apply_to: targets # targets | all -outlier_config: - iqr_k: 1.5 - zscore_threshold: 4.0 - quantile_lower: 0.001 - quantile_upper: 0.999 - -# 以 txt 第 11 列(列名 V_pi)为物理门控:仅保留闭区间 [v_pi_min, v_pi_max] 内样本 -filter_v_pi_range: true -v_pi_min: 0.0 -v_pi_max: 500.0 - -# 在区间过滤之后,是否再剔除 V_pi<=0;若需保留 V_pi=0(仍在 [0,500] 内),请设为 false -remove_nonpositive_vpi: false +data: + test_size: 0.1 + random_state: 123 + filter_v_pi_max: 500.0 model: input_dim: 8 - hidden_dims: [200, 300, 350, 300, 200] output_dim: 3 - batchnorm: false - # 温和 dropout,实测略优于全 0(见 results/run_20260419_163305) - dropout: 0.05 - residual: false + hidden_dims: [64, 128, 64] + n_experts: 60 + gating_hidden: 8 + dropout_rate: 0.0 + use_bn: true + activation: relu optimizer: - name: adamw lr: 0.001 - weight_decay: 0.0001 - -scheduler: - type: cosine # cosine | plateau - plateau_factor: 0.5 - plateau_patience: 10 - plateau_min_lr: 1.0e-6 + weight_decay: 0.05 + betas: [0.9, 0.999] training: batch_size: 128 - epochs: 300 - early_stopping_patience: 30 + epochs: 100 num_workers: 0 -loss: - type: huber # huber | weighted_mse - huber_delta: 1.0 - # BW_3dB, IL, V_pi;略加重 V_pi 以小幅提升其测试 R² - target_weights: [1.0, 1.0, 1.2] +# 默认会先从 best_hyperparams.json 读取这些系数,再用此处显式值覆盖 +physics: + lambda_bw_mon: 0.0 + lambda_IL_mon: 0.3 + lambda_vpiL: 0.005 + lambda_smooth: 0.1 -# 总输出目录;每次训练会在其下创建 run_时间戳/ output_dir: results - -# 评估/推理时若未指定 run_dir,可填最近一次 run 的路径(可选) last_run_dir: null diff --git a/src/config.py b/src/config.py index 770a7a0..9af2dc7 100644 --- a/src/config.py +++ b/src/config.py @@ -2,6 +2,7 @@ from __future__ import annotations +import json from dataclasses import dataclass, field from pathlib import Path from typing import Any, List, Optional @@ -9,137 +10,118 @@ from typing import Any, List, Optional import yaml +@dataclass +class DataConfig: + test_size: float = 0.1 + random_state: int = 123 + filter_v_pi_max: float = 500.0 + + @dataclass class ModelConfig: input_dim: int = 8 - hidden_dims: List[int] = field(default_factory=lambda: [200, 300, 350, 300, 200]) output_dim: int = 3 - batchnorm: bool = False - dropout: float = 0.0 - residual: bool = False + hidden_dims: List[int] = field(default_factory=lambda: [64, 128, 64]) + n_experts: int = 60 + gating_hidden: int = 8 + dropout_rate: float = 0.0 + use_bn: bool = True + activation: str = "relu" @dataclass class OptimizerConfig: - name: str = "adamw" lr: float = 1e-3 - weight_decay: float = 1e-4 - - -@dataclass -class SchedulerConfig: - type: str = "cosine" # cosine | plateau - plateau_factor: float = 0.5 - plateau_patience: int = 10 - plateau_min_lr: float = 1e-6 + weight_decay: float = 0.05 + betas: List[float] = field(default_factory=lambda: [0.9, 0.999]) @dataclass class TrainingConfig: batch_size: int = 128 - epochs: int = 300 - early_stopping_patience: int = 30 + epochs: int = 100 num_workers: int = 0 @dataclass -class LossConfig: - type: str = "huber" # huber | weighted_mse - huber_delta: float = 1.0 - target_weights: List[float] = field(default_factory=lambda: [1.0, 1.0, 1.0]) - - -@dataclass -class OutlierConfig: - iqr_k: float = 1.5 - zscore_threshold: float = 4.0 - quantile_lower: float = 0.001 - quantile_upper: float = 0.999 +class PhysicsConfig: + lambda_bw_mon: float = 0.0 + lambda_IL_mon: float = 0.3 + lambda_vpiL: float = 0.005 + lambda_smooth: float = 0.1 @dataclass class AppConfig: data_path: str - split_ratios: List[float] - random_seed: int - split_mode: str - split_stratify_target: str - split_stratify_bins: int - remove_duplicate_rows: bool - outlier_strategy: str - outlier_config: OutlierConfig - outlier_apply_to: str # targets | all - remove_nonpositive_vpi: bool - filter_v_pi_range: bool - v_pi_min: float - v_pi_max: float + data: DataConfig model: ModelConfig optimizer: OptimizerConfig - scheduler: SchedulerConfig training: TrainingConfig - loss: LossConfig + physics: PhysicsConfig output_dir: str + best_hyperparams_path: Optional[str] = None last_run_dir: Optional[str] = None @staticmethod - def from_dict(raw: dict[str, Any]) -> "AppConfig": - m = raw.get("model", {}) - o = raw.get("optimizer", {}) - s = raw.get("scheduler", {}) - t = raw.get("training", {}) - l = raw.get("loss", {}) - oc = raw.get("outlier_config", {}) + def from_dict(raw: dict[str, Any], cfg_dir: Path) -> "AppConfig": + data_raw = raw.get("data", {}) + model_raw = raw.get("model", {}) + optimizer_raw = raw.get("optimizer", {}) + training_raw = raw.get("training", {}) + physics_raw = raw.get("physics", {}) + best_hyperparams_path = raw.get("best_hyperparams_path") + + if best_hyperparams_path: + hp_path = Path(best_hyperparams_path) + if not hp_path.is_absolute(): + cand = (cfg_dir / hp_path).resolve() + if cand.is_file(): + hp_path = cand + else: + hp_path = (cfg_dir.parent / hp_path).resolve() + with hp_path.open("r", encoding="utf-8") as f: + hp_raw = json.load(f) + physics_raw = {**hp_raw.get("best_config", {}), **physics_raw} + best_hyperparams_path = str(hp_path) + return AppConfig( data_path=str(raw["data_path"]), - split_ratios=list(raw["split_ratios"]), - random_seed=int(raw["random_seed"]), - split_mode=str(raw.get("split_mode", "grouped_stratified")), - split_stratify_target=str(raw.get("split_stratify_target", "V_pi")), - split_stratify_bins=int(raw.get("split_stratify_bins", 10)), - remove_duplicate_rows=bool(raw["remove_duplicate_rows"]), - outlier_strategy=str(raw.get("outlier_strategy", "none")), - outlier_config=OutlierConfig( - iqr_k=float(oc.get("iqr_k", 1.5)), - zscore_threshold=float(oc.get("zscore_threshold", 4.0)), - quantile_lower=float(oc.get("quantile_lower", 0.001)), - quantile_upper=float(oc.get("quantile_upper", 0.999)), + data=DataConfig( + test_size=float(data_raw.get("test_size", raw.get("test_size", 0.1))), + random_state=int(data_raw.get("random_state", raw.get("random_seed", 123))), + filter_v_pi_max=float( + data_raw.get("filter_v_pi_max", raw.get("v_pi_max", 500.0)) + ), ), - outlier_apply_to=str(raw.get("outlier_apply_to", "targets")), - remove_nonpositive_vpi=bool(raw.get("remove_nonpositive_vpi", False)), - filter_v_pi_range=bool(raw.get("filter_v_pi_range", True)), - v_pi_min=float(raw.get("v_pi_min", 0.0)), - v_pi_max=float(raw.get("v_pi_max", 500.0)), model=ModelConfig( - input_dim=int(m.get("input_dim", 8)), - hidden_dims=list(m.get("hidden_dims", [200, 300, 350, 300, 200])), - output_dim=int(m.get("output_dim", 3)), - batchnorm=bool(m.get("batchnorm", False)), - dropout=float(m.get("dropout", 0.0)), - residual=bool(m.get("residual", False)), + input_dim=int(model_raw.get("input_dim", 8)), + output_dim=int(model_raw.get("output_dim", 3)), + hidden_dims=list(model_raw.get("hidden_dims", [64, 128, 64])), + n_experts=int(model_raw.get("n_experts", 60)), + gating_hidden=int(model_raw.get("gating_hidden", 8)), + dropout_rate=float(model_raw.get("dropout_rate", 0.0)), + use_bn=bool(model_raw.get("use_bn", True)), + activation=str(model_raw.get("activation", "relu")), ), optimizer=OptimizerConfig( - name=str(o.get("name", "adamw")), - lr=float(o.get("lr", 1e-3)), - weight_decay=float(o.get("weight_decay", 1e-4)), - ), - scheduler=SchedulerConfig( - type=str(s.get("type", "cosine")), - plateau_factor=float(s.get("plateau_factor", 0.5)), - plateau_patience=int(s.get("plateau_patience", 10)), - plateau_min_lr=float(s.get("plateau_min_lr", 1e-6)), + lr=float(optimizer_raw.get("lr", 1e-3)), + weight_decay=float(optimizer_raw.get("weight_decay", 0.05)), + betas=[float(x) for x in optimizer_raw.get("betas", [0.9, 0.999])], ), training=TrainingConfig( - batch_size=int(t.get("batch_size", 128)), - epochs=int(t.get("epochs", 300)), - early_stopping_patience=int(t.get("early_stopping_patience", 30)), - num_workers=int(t.get("num_workers", 0)), + batch_size=int(training_raw.get("batch_size", 128)), + epochs=int(training_raw.get("epochs", 100)), + num_workers=int(training_raw.get("num_workers", 0)), ), - loss=LossConfig( - type=str(l.get("type", "huber")), - huber_delta=float(l.get("huber_delta", 1.0)), - target_weights=[float(x) for x in l.get("target_weights", [1.0, 1.0, 1.0])], + physics=PhysicsConfig( + lambda_bw_mon=float(physics_raw.get("lambda_bw_mon", 0.0)), + lambda_IL_mon=float(physics_raw.get("lambda_IL_mon", 0.3)), + lambda_vpiL=float(physics_raw.get("lambda_vpiL", 0.005)), + lambda_smooth=float(physics_raw.get("lambda_smooth", 0.1)), ), output_dir=str(raw.get("output_dir", "results")), + best_hyperparams_path=best_hyperparams_path, last_run_dir=raw.get("last_run_dir"), ) @@ -153,24 +135,25 @@ def load_config(path: str | Path) -> AppConfig: raw = yaml.safe_load(f) if not isinstance(raw, dict): raise ValueError("YAML 根节点必须是字典") - cfg = AppConfig.from_dict(raw) - sr = cfg.split_ratios - if len(sr) != 3: - raise ValueError("split_ratios 必须为长度为 3 的列表 [train, val, test]") - if abs(sum(sr) - 1.0) > 1e-6: - raise ValueError(f"split_ratios 之和必须为 1,当前为 {sum(sr)}") - if cfg.split_mode not in ("random", "grouped_stratified"): - raise ValueError("split_mode 必须为 random 或 grouped_stratified") - if cfg.split_stratify_target not in ("BW_3dB", "IL", "V_pi"): - raise ValueError("split_stratify_target 必须为 BW_3dB、IL 或 V_pi") - if cfg.split_stratify_bins < 2: - raise ValueError("split_stratify_bins 必须 >= 2") - if cfg.outlier_strategy not in ("none", "iqr", "zscore", "quantile_clip"): - raise ValueError(f"未知 outlier_strategy: {cfg.outlier_strategy}") - if cfg.outlier_apply_to not in ("targets", "all"): - raise ValueError("outlier_apply_to 必须为 targets 或 all") - if len(cfg.loss.target_weights) != 3: - raise ValueError("loss.target_weights 长度必须为 3") - if cfg.filter_v_pi_range and cfg.v_pi_min >= cfg.v_pi_max: - raise ValueError("启用 filter_v_pi_range 时须满足 v_pi_min < v_pi_max") + cfg = AppConfig.from_dict(raw, path.parent.resolve()) + if not 0.0 < cfg.data.test_size < 1.0: + raise ValueError("data.test_size 必须在 (0, 1) 之间") + if cfg.data.filter_v_pi_max <= 0: + raise ValueError("data.filter_v_pi_max 必须 > 0") + if cfg.model.input_dim != 8: + raise ValueError("model.input_dim 必须为 8") + if cfg.model.output_dim != 3: + raise ValueError("model.output_dim 必须为 3") + if not cfg.model.hidden_dims: + raise ValueError("model.hidden_dims 不能为空") + if cfg.model.n_experts < 1: + raise ValueError("model.n_experts 必须 >= 1") + if cfg.model.gating_hidden < 1: + raise ValueError("model.gating_hidden 必须 >= 1") + if cfg.model.activation not in ("relu", "gaussian"): + raise ValueError("model.activation 必须为 relu 或 gaussian") + if len(cfg.optimizer.betas) != 2: + raise ValueError("optimizer.betas 长度必须为 2") + if cfg.training.batch_size < 1 or cfg.training.epochs < 1: + raise ValueError("training.batch_size 与 training.epochs 必须 >= 1") return cfg diff --git a/src/evaluate.py b/src/evaluate.py index b1f8494..9db8186 100644 --- a/src/evaluate.py +++ b/src/evaluate.py @@ -1,12 +1,12 @@ -"""加载最优模型并在各划分上评估,导出 CSV / JSON。""" +"""加载模型并在 train/test 上评估,导出 CSV / JSON。""" from __future__ import annotations import csv import json import logging -from pathlib import Path from dataclasses import asdict +from pathlib import Path from typing import Dict, Tuple import numpy as np @@ -16,15 +16,10 @@ import torch.nn as nn from src.config import AppConfig from src.data import INPUT_COLUMNS, TARGET_COLUMNS -from src.losses import build_loss -from src.metrics import ( - FullMetricsReport, - compute_full_report, - report_to_flat_dict, -) -from src.model import MLPRegressor -from src.preprocess import ProcessedDataBundle -from src.trainer import evaluate_loss_loader, load_weights +from src.metrics import FullMetricsReport, compute_full_report, report_to_flat_dict +from src.model import create_model_from_config +from src.preprocess import ProcessedDataBundle, inverse_transform_targets +from src.trainer import load_weights logger = logging.getLogger(__name__) @@ -40,25 +35,32 @@ def predict_all( for xb, yb in loader: xb = xb.to(device) pr = model(xb).detach().cpu().numpy() - yt = yb.numpy() preds.append(pr) - trues.append(yt) + trues.append(yb.numpy()) return np.concatenate(preds, axis=0), np.concatenate(trues, axis=0) +def _select_checkpoint(run_dir: Path) -> Path: + ckpt_last = run_dir / "checkpoints" / "last.pt" + if ckpt_last.is_file(): + return ckpt_last + ckpt_best = run_dir / "checkpoints" / "best.pt" + if ckpt_best.is_file(): + return ckpt_best + raise FileNotFoundError(f"未找到 {run_dir}/checkpoints/last.pt 或 best.pt") + + def evaluate_split( model: nn.Module, - criterion: nn.Module, loader: torch.utils.data.DataLoader, device: torch.device, - y_scaler, + y_scalers, split_name: str, ) -> Tuple[FullMetricsReport, FullMetricsReport, float]: - """返回 (标准化空间报告, 物理空间报告, 平均损失)。""" - loss = evaluate_loss_loader(model, loader, criterion, device) pred_n, true_n = predict_all(model, loader, device) - pred_p = y_scaler.inverse_transform(pred_n) - true_p = y_scaler.inverse_transform(true_n) + loss = float(nn.functional.mse_loss(torch.from_numpy(pred_n), torch.from_numpy(true_n)).item()) + pred_p = inverse_transform_targets(pred_n, y_scalers) + true_p = inverse_transform_targets(true_n, y_scalers) rep_n = compute_full_report(split_name, loss, true_n, pred_n, TARGET_COLUMNS) rep_p = compute_full_report(split_name, loss, true_p, pred_p, TARGET_COLUMNS) return rep_n, rep_p, loss @@ -70,31 +72,18 @@ def run_full_evaluation( run_dir: Path, device: torch.device, ) -> Tuple[nn.Module, Dict]: - """载入 best.pt,在 train/val/test 上评估并写 metrics.csv 与 summary.json;返回模型与摘要。""" - model = MLPRegressor( - input_dim=cfg.model.input_dim, - hidden_dims=cfg.model.hidden_dims, - output_dim=cfg.model.output_dim, - batchnorm=cfg.model.batchnorm, - dropout=cfg.model.dropout, - residual=cfg.model.residual, - ).to(device) - ckpt_best = run_dir / "checkpoints" / "best.pt" - load_weights(model, ckpt_best, device) - - criterion = build_loss(cfg.loss).to(device) - y_scaler = bundle.y_scaler + model = create_model_from_config(cfg).to(device) + load_weights(model, _select_checkpoint(run_dir), device) rows = [] summary: Dict = {"splits": {}} - - for name, loader in ( - ("train", bundle.train_loader), - ("val", bundle.val_loader), - ("test", bundle.test_loader), - ): + for name, loader in (("train", bundle.train_loader), ("test", bundle.test_loader)): rep_n, rep_p, loss = evaluate_split( - model, criterion, loader, device, y_scaler, name + model, + loader, + device, + bundle.y_scalers, + name, ) summary["splits"][name] = { "loss": loss, @@ -106,11 +95,10 @@ def run_full_evaluation( rows.append(row) metrics_path = run_dir / "metrics.csv" - if rows: - with metrics_path.open("w", newline="", encoding="utf-8") as f: - writer = csv.DictWriter(f, fieldnames=list(rows[0].keys())) - writer.writeheader() - writer.writerows(rows) + with metrics_path.open("w", newline="", encoding="utf-8") as f: + writer = csv.DictWriter(f, fieldnames=list(rows[0].keys())) + writer.writeheader() + writer.writerows(rows) (run_dir / "summary.json").write_text( json.dumps(summary, indent=2, ensure_ascii=False, default=str), encoding="utf-8", @@ -125,11 +113,10 @@ def export_test_predictions_csv( device: torch.device, path: Path, ) -> None: - """导出测试集物理空间真值、预测与误差。""" model.eval() pred_n, true_n = predict_all(model, bundle.test_loader, device) - pred_p = bundle.y_scaler.inverse_transform(pred_n) - true_p = bundle.y_scaler.inverse_transform(true_n) + pred_p = inverse_transform_targets(pred_n, bundle.y_scalers) + true_p = inverse_transform_targets(true_n, bundle.y_scalers) err = pred_p - true_p cols: Dict[str, np.ndarray] = {} for j, name in enumerate(INPUT_COLUMNS): diff --git a/src/infer.py b/src/infer.py index 365e60c..b26f9f8 100644 --- a/src/infer.py +++ b/src/infer.py @@ -6,13 +6,14 @@ import logging from pathlib import Path from typing import Optional +import numpy as np import pandas as pd import torch from src.config import AppConfig, load_config from src.data import INPUT_COLUMNS, TARGET_COLUMNS, _strip_optional_list_brackets -from src.model import MLPRegressor -from src.preprocess import load_scalers +from src.model import create_model_from_config +from src.preprocess import inverse_transform_targets, load_scalers from src.trainer import load_weights logger = logging.getLogger(__name__) @@ -68,23 +69,19 @@ def run_inference( 载入 best 模型与 scaler,对输入表进行批量推理并写出 CSV(物理量空间)。 """ X = _read_inputs_table(Path(input_path)).to_numpy(dtype=np.float32) - X_scaler, y_scaler = load_scalers(run_dir) + X_scaler, y_scalers = load_scalers(run_dir) Xn = X_scaler.transform(X) - model = MLPRegressor( - input_dim=cfg.model.input_dim, - hidden_dims=cfg.model.hidden_dims, - output_dim=cfg.model.output_dim, - batchnorm=cfg.model.batchnorm, - dropout=cfg.model.dropout, - residual=cfg.model.residual, - ).to(device) - load_weights(model, run_dir / "checkpoints" / "best.pt", device) + model = create_model_from_config(cfg).to(device) + ckpt = run_dir / "checkpoints" / "last.pt" + if not ckpt.is_file(): + ckpt = run_dir / "checkpoints" / "best.pt" + load_weights(model, ckpt, device) model.eval() with torch.no_grad(): pred_n = model(torch.from_numpy(Xn).float().to(device)).cpu().numpy() - pred_p = y_scaler.inverse_transform(pred_n) + pred_p = inverse_transform_targets(pred_n, y_scalers) out = pd.DataFrame(X, columns=INPUT_COLUMNS) for j, name in enumerate(TARGET_COLUMNS): diff --git a/src/main.py b/src/main.py index 87a38b6..42d516e 100644 --- a/src/main.py +++ b/src/main.py @@ -13,7 +13,7 @@ import torch from src.config import load_config from src.data import load_raw_txt, quality_report_before_clean, summarize_for_console from src.evaluate import export_test_predictions_csv, run_full_evaluation -from src.model import MLPRegressor +from src.model import create_model_from_config, weights_init from src.plots import generate_all_figures from src.preprocess import prepare_training_data, rebuild_bundle_for_eval from src.trainer import fit @@ -61,7 +61,7 @@ def _write_summary_md(run_dir: Path, summary: dict) -> None: lines.append(f"## {split}") lines.append("") lines.append( - f"- **损失(标准化输出空间 Huber/MSE 准则)**: {block['loss']:.6f}" + f"- **损失(标准化输出空间 MSE)**: {block['loss']:.6f}" ) for space, label in ("normalized", "标准化空间"), ("physical", "物理量空间"): sub = block[space] @@ -78,7 +78,7 @@ def _write_summary_md(run_dir: Path, summary: dict) -> None: def cmd_train(args: argparse.Namespace) -> None: cfg_path = _resolve_cfg_path(args.config) cfg = load_config(cfg_path) - set_global_seed(cfg.random_seed) + set_global_seed(cfg.data.random_state) run_dir = make_run_dir(resolve_path(cfg.output_dir, _project_root())) shutil.copy2(cfg_path, run_dir / "config_snapshot.yaml") @@ -86,21 +86,15 @@ def cmd_train(args: argparse.Namespace) -> None: data_path = resolve_path(cfg.data_path, _project_root()) df = load_raw_txt(data_path) - q = quality_report_before_clean(df, cfg.v_pi_min, cfg.v_pi_max) + q = quality_report_before_clean(df, 0.0, cfg.data.filter_v_pi_max) logger.info("数据质量(清洗前): %s", summarize_for_console(df, q)) bundle = prepare_training_data(df, cfg, run_dir) device = torch.device("cuda" if torch.cuda.is_available() else "cpu") - model = MLPRegressor( - input_dim=cfg.model.input_dim, - hidden_dims=cfg.model.hidden_dims, - output_dim=cfg.model.output_dim, - batchnorm=cfg.model.batchnorm, - dropout=cfg.model.dropout, - residual=cfg.model.residual, - ).to(device) + model = create_model_from_config(cfg).to(device) + model.apply(weights_init) - history = fit(model, cfg, bundle.train_loader, bundle.val_loader, run_dir, device) + history = fit(model, cfg, bundle, run_dir, device) model_eval, summary = run_full_evaluation(cfg, bundle, run_dir, device) export_test_predictions_csv( bundle, model_eval, device, run_dir / "test_predictions.csv" @@ -115,7 +109,7 @@ def cmd_eval(args: argparse.Namespace) -> None: run_dir = _resolve_run_dir(cfg_path, args.run_dir, args.output_dir) snap = run_dir / "config_snapshot.yaml" cfg = load_config(snap if snap.is_file() else cfg_path) - set_global_seed(cfg.random_seed) + set_global_seed(cfg.data.random_state) setup_logging(run_dir / "eval.log") data_path = resolve_path(cfg.data_path, _project_root()) @@ -135,7 +129,7 @@ def cmd_infer(args: argparse.Namespace) -> None: run_dir = _resolve_run_dir(cfg_path, args.run_dir, args.output_dir) snap = run_dir / "config_snapshot.yaml" cfg = load_config(snap if snap.is_file() else cfg_path) - set_global_seed(cfg.random_seed) + set_global_seed(cfg.data.random_state) setup_logging(None) from src.infer import run_inference @@ -147,7 +141,7 @@ def cmd_infer(args: argparse.Namespace) -> None: def build_parser() -> argparse.ArgumentParser: - p = argparse.ArgumentParser(description="MZM MLP 训练 / 评估 / 推理") + p = argparse.ArgumentParser(description="MZM MoE PINN 训练 / 评估 / 推理") sub = p.add_subparsers(dest="command", required=True) pt = sub.add_parser("train", help="训练模型") diff --git a/src/model.py b/src/model.py index c2c2fc0..7645ca7 100644 --- a/src/model.py +++ b/src/model.py @@ -1,4 +1,4 @@ -"""可配置 MLP 回归模型。""" +"""MoE PINN 模型定义。""" from __future__ import annotations @@ -8,54 +8,103 @@ import torch import torch.nn as nn -def kaiming_init_module(m: nn.Module) -> None: - """对 Linear 使用 Kaiming uniform(ReLU),偏置置零。""" - if isinstance(m, nn.Linear): - nn.init.kaiming_uniform_(m.weight, nonlinearity="relu") - if m.bias is not None: - nn.init.zeros_(m.bias) - elif isinstance(m, nn.BatchNorm1d): - nn.init.ones_(m.weight) - nn.init.zeros_(m.bias) +class GaussianActivation(nn.Module): + def forward(self, x: torch.Tensor) -> torch.Tensor: + return torch.exp(-(x**2)) -class MLPRegressor(nn.Module): - """ - 多层感知机回归:输入 8 维,输出 3 维。 +def weights_init(layer_in: nn.Module) -> None: + """与 notebook 保持一致的 Kaiming 初始化。""" + if isinstance(layer_in, nn.Linear): + nn.init.kaiming_uniform_(layer_in.weight) + if layer_in.bias is not None: + layer_in.bias.data.fill_(0.0) - 可选 BatchNorm1d、Dropout、以及在相邻层维度相等时的残差相加。 - """ +def build_activation(name: str) -> nn.Module: + if name == "relu": + return nn.ReLU() + if name == "gaussian": + return GaussianActivation() + raise ValueError(f"未知激活函数: {name}") + + +class ExpertNN(nn.Module): def __init__( self, input_dim: int, - hidden_dims: List[int], output_dim: int, - batchnorm: bool = False, - dropout: float = 0.0, - residual: bool = False, + hidden_dims: List[int], + activation_fn: nn.Module, + dropout_rate: float = 0.0, + use_bn: bool = False, ) -> None: super().__init__() - self.residual = residual - dims = [input_dim] + list(hidden_dims) + [output_dim] - self._hidden_blocks = nn.ModuleList() - for i in range(len(dims) - 2): - in_d, out_d = dims[i], dims[i + 1] - seq_layers: list[nn.Module] = [nn.Linear(in_d, out_d)] - if batchnorm: - seq_layers.append(nn.BatchNorm1d(out_d)) - seq_layers.append(nn.ReLU(inplace=True)) - if dropout and dropout > 0: - seq_layers.append(nn.Dropout(p=dropout)) - self._hidden_blocks.append(nn.Sequential(*seq_layers)) - self._head = nn.Linear(dims[-2], dims[-1]) - self.apply(kaiming_init_module) + layers: list[nn.Module] = [] + prev_dim = input_dim + for h in hidden_dims: + layers.append(nn.Linear(prev_dim, h)) + if use_bn: + layers.append(nn.BatchNorm1d(h)) + layers.append(type(activation_fn)() if isinstance(activation_fn, nn.ReLU) else activation_fn.__class__()) + if dropout_rate > 0: + layers.append(nn.Dropout(p=dropout_rate)) + prev_dim = h + layers.append(nn.Linear(prev_dim, output_dim)) + self.net = nn.Sequential(*layers) def forward(self, x: torch.Tensor) -> torch.Tensor: - h = x - for block in self._hidden_blocks: - inp = h - h = block(inp) - if self.residual and inp.shape[-1] == h.shape[-1]: - h = h + inp - return self._head(h) + return self.net(x) + + +class MixtureOfExperts(nn.Module): + def __init__( + self, + input_dim: int, + output_dim: int, + hidden_dims: List[int], + n_experts: int = 3, + activation_fn: nn.Module | None = None, + gating_hidden: int = 32, + dropout_rate: float = 0.0, + use_bn: bool = False, + ) -> None: + super().__init__() + act = activation_fn if activation_fn is not None else nn.ReLU() + self.experts = nn.ModuleList( + [ + ExpertNN( + input_dim, + output_dim, + hidden_dims, + act, + dropout_rate=dropout_rate, + use_bn=use_bn, + ) + for _ in range(n_experts) + ] + ) + self.gating = nn.Sequential( + nn.Linear(input_dim, gating_hidden), + nn.ReLU(), + nn.Linear(gating_hidden, n_experts), + nn.Softmax(dim=1), + ) + + def forward(self, x: torch.Tensor) -> torch.Tensor: + gate_weights = self.gating(x) + expert_outputs = torch.stack([expert(x) for expert in self.experts], dim=2) + return torch.bmm(expert_outputs, gate_weights.unsqueeze(2)).squeeze(2) + + +def create_model_from_config(cfg) -> MixtureOfExperts: + return MixtureOfExperts( + input_dim=cfg.model.input_dim, + output_dim=cfg.model.output_dim, + hidden_dims=cfg.model.hidden_dims, + n_experts=cfg.model.n_experts, + activation_fn=build_activation(cfg.model.activation), + gating_hidden=cfg.model.gating_hidden, + dropout_rate=cfg.model.dropout_rate, + use_bn=cfg.model.use_bn, + ) diff --git a/src/plots.py b/src/plots.py index a26fa63..523e880 100644 --- a/src/plots.py +++ b/src/plots.py @@ -13,22 +13,22 @@ import torch from src.config import AppConfig from src.data import TARGET_COLUMNS from src.evaluate import predict_all -from src.model import MLPRegressor -from src.preprocess import ProcessedDataBundle +from src.model import create_model_from_config +from src.preprocess import ProcessedDataBundle, inverse_transform_targets from src.trainer import TrainHistory, load_weights logger = logging.getLogger(__name__) def plot_loss_curves(history: TrainHistory, out_path: Path) -> None: - """绘制 train/val loss 曲线。""" + """绘制 train/test loss 曲线。""" sns.set_theme(style="whitegrid", context="talk") fig, ax = plt.subplots(figsize=(8, 5)) ax.plot(history.epoch, history.train_loss, label="Train loss", linewidth=2) - ax.plot(history.epoch, history.val_loss, label="Val loss", linewidth=2) + ax.plot(history.epoch, history.test_loss, label="Test loss", linewidth=2) ax.set_xlabel("Epoch") ax.set_ylabel("Loss (normalized target space)") - ax.set_title("Training / Validation Loss") + ax.set_title("Training / Test Loss") ax.legend() fig.tight_layout() out_path.parent.mkdir(parents=True, exist_ok=True) @@ -104,19 +104,15 @@ def generate_all_figures( fig_dir.mkdir(parents=True, exist_ok=True) plot_loss_curves(history, fig_dir / "loss_curve.png") - model = MLPRegressor( - input_dim=cfg.model.input_dim, - hidden_dims=cfg.model.hidden_dims, - output_dim=cfg.model.output_dim, - batchnorm=cfg.model.batchnorm, - dropout=cfg.model.dropout, - residual=cfg.model.residual, - ).to(device) - load_weights(model, run_dir / "checkpoints" / "best.pt", device) + model = create_model_from_config(cfg).to(device) + ckpt = run_dir / "checkpoints" / "last.pt" + if not ckpt.is_file(): + ckpt = run_dir / "checkpoints" / "best.pt" + load_weights(model, ckpt, device) pred_n, true_n = predict_all(model, bundle.test_loader, device) - pred_p = bundle.y_scaler.inverse_transform(pred_n) - true_p = bundle.y_scaler.inverse_transform(true_n) + pred_p = inverse_transform_targets(pred_n, bundle.y_scalers) + true_p = inverse_transform_targets(true_n, bundle.y_scalers) for i, name in enumerate(TARGET_COLUMNS): plot_scatter_true_pred( diff --git a/src/preprocess.py b/src/preprocess.py index a6a64af..cb96ba3 100644 --- a/src/preprocess.py +++ b/src/preprocess.py @@ -7,7 +7,7 @@ import logging import pickle from dataclasses import dataclass from pathlib import Path -from typing import Dict, List, Optional, Tuple +from typing import List, Sequence, Tuple import numpy as np import pandas as pd @@ -17,331 +17,55 @@ from sklearn.preprocessing import StandardScaler from torch.utils.data import DataLoader, TensorDataset from src.config import AppConfig -from src.data import ALL_COLUMNS, INPUT_COLUMNS, TARGET_COLUMNS, quality_report_before_clean +from src.data import INPUT_COLUMNS, TARGET_COLUMNS, quality_report_before_clean logger = logging.getLogger(__name__) @dataclass class ProcessedDataBundle: - """训练用张量与 DataLoader,以及划分后的 numpy(含测试集原始物理量用于导出)。""" + """训练用张量与 DataLoader,以及 notebook 风格的 train/test 划分数据。""" train_loader: DataLoader - val_loader: DataLoader test_loader: DataLoader X_train: np.ndarray - X_val: np.ndarray X_test: np.ndarray y_train: np.ndarray - y_val: np.ndarray y_test: np.ndarray + X_train_raw: np.ndarray X_test_raw: np.ndarray + y_train_raw: np.ndarray y_test_raw: np.ndarray X_scaler: StandardScaler - y_scaler: StandardScaler + y_scalers: List[StandardScaler] feature_names: List[str] target_names: List[str] -def _mask_outliers_iqr( - values: np.ndarray, col_names: List[str], k: float -) -> np.ndarray: - """返回 True 表示该行在任一选定列上超出训练集 IQR 范围(基于传入的 values 统计)。""" - mask = np.zeros(len(values), dtype=bool) - for j, _ in enumerate(col_names): - col = values[:, j] - q1, q3 = np.percentile(col, [25, 75]) - iqr = q3 - q1 - lo, hi = q1 - k * iqr, q3 + k * iqr - mask |= (col < lo) | (col > hi) - return mask - - -def _mask_outliers_zscore(values: np.ndarray, threshold: float) -> np.ndarray: - mask = np.zeros(len(values), dtype=bool) - for j in range(values.shape[1]): - col = values[:, j] - mu, sig = col.mean(), col.std(ddof=0) - if sig < 1e-12: - continue - z = np.abs((col - mu) / sig) - mask |= z > threshold - return mask - - -def _winsorize_train_apply_all( - train: np.ndarray, - val: np.ndarray, - test: np.ndarray, - ql: float, - qu: float, -) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: - """按训练集分位数对 train/val/test 同步裁剪(列方向)。""" - lo = np.quantile(train, ql, axis=0) - hi = np.quantile(train, qu, axis=0) - def clip_arr(a: np.ndarray) -> np.ndarray: - return np.clip(a, lo, hi) - return clip_arr(train), clip_arr(val), clip_arr(test) - - def clean_dataframe( df: pd.DataFrame, cfg: AppConfig, report_lines: List[str], ) -> pd.DataFrame: - """ - 清洗流程(顺序固定,便于复现与审计): - - 1. 可选:完全重复行去重。 - 2. 可选:以 **最后一列对应字段 V_pi**(txt 第 11 个逗号分隔字段)为门控,仅保留 - ``v_pi_min <= V_pi <= v_pi_max``(默认 [0, 500])。 - 3. 可选:再移除 ``V_pi <= 0``(与区间门控独立,由配置控制)。 - """ + """按 notebook 逻辑清洗:仅保留 V_pi < 阈值。""" out = df.copy() n0 = len(out) - if cfg.remove_duplicate_rows: - out = out.drop_duplicates() - report_lines.append(f"去完全重复行: {n0} -> {len(out)}") - if cfg.filter_v_pi_range: - n1 = len(out) - lo, hi = float(cfg.v_pi_min), float(cfg.v_pi_max) - mask = (out["V_pi"] >= lo) & (out["V_pi"] <= hi) - out = out[mask].reset_index(drop=True) - report_lines.append( - f"V_pi 物理区间过滤 [{lo}, {hi}](txt 第 11 列 / 列名 V_pi): {n1} -> {len(out)}" - ) - if cfg.remove_nonpositive_vpi: - n2 = len(out) - out = out[out["V_pi"] > 0].reset_index(drop=True) - report_lines.append(f"移除 V_pi<=0: {n2} -> {len(out)}") + vmax = float(cfg.data.filter_v_pi_max) + out = out[out["V_pi"] < vmax].reset_index(drop=True) + report_lines.append(f"V_pi 阈值过滤 (< {vmax}):{n0} -> {len(out)}") if len(out) == 0: - raise ValueError( - "清洗后样本数为 0:请检查 V_pi 区间配置、数据源或是否过度去重。" - ) + raise ValueError("清洗后样本数为 0,请检查数据源或 V_pi 阈值设置。") return out -def _random_split_indices( - n: int, - ratios: List[float], - seed: int, -) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: - """返回 train/val/test 的整数索引(先 shuffle 再按比例切分)。""" - rng = np.random.default_rng(seed) - idx = np.arange(n) - rng.shuffle(idx) - tr, va, te = ratios - n_test = int(round(n * te)) - n_val = int(round(n * va)) - n_train = n - n_val - n_test - if n_train <= 0 or n_val <= 0 or n_test <= 0: - raise ValueError( - f"划分后样本过少: train={n_train}, val={n_val}, test={n_test},请调整比例或数据量" - ) - i_train = idx[:n_train] - i_val = idx[n_train : n_train + n_val] - i_test = idx[n_train + n_val :] - return i_train, i_val, i_test - - -def _quantile_bin_labels(values: np.ndarray, n_bins: int) -> np.ndarray | None: - """ - 基于秩做近似等频分桶,避免重复值导致的 qcut 退化。 - 返回每个样本所属桶标签;若样本过少则返回 None。 - """ - if len(values) < 2: - return None - q = min(int(n_bins), len(values)) - if q < 2: - return None - ranks = pd.Series(values).rank(method="first") - labels = pd.qcut(ranks, q=q, labels=False, duplicates="drop") - if labels is None: - return None - arr = np.asarray(labels, dtype=int) - if len(np.unique(arr)) < 2: - return None - return arr - - -def _grouped_split_indices( - df: pd.DataFrame, - ratios: List[float], - seed: int, - stratify_target: str, - stratify_bins: int, - report_lines: List[str], -) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: - """ - 先按输入 8 维分组,再在组级别按目标统计量近似分层切分。 - 这样可避免“同输入异输出”同时落在 train/val/test,提升评估稳定性。 - 若分层条件不足,则退化为组级随机切分。 - """ - if len(df) < 3: - raise ValueError("样本数过少,无法做 train/val/test 切分。") - - group_ids = df.groupby(INPUT_COLUMNS, sort=False, dropna=False).ngroup().to_numpy() - n_groups = int(group_ids.max()) + 1 - group_df = df.copy() - group_df["_group_id"] = group_ids - group_stat = ( - group_df.groupby("_group_id", sort=True) - .agg(group_size=("V_pi", "size"), strat_value=(stratify_target, "median")) - .reset_index() - ) - group_id_arr = group_stat["_group_id"].to_numpy(dtype=int) - labels = _quantile_bin_labels( - group_stat["strat_value"].to_numpy(dtype=np.float64), - stratify_bins, - ) - - tr, va, te = ratios - holdout_ratio = va + te - val_ratio_in_holdout = va / holdout_ratio - - def _split_groups(use_stratify: bool) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: - strat = labels if use_stratify and labels is not None else None - train_groups, holdout_groups = train_test_split( - group_id_arr, - train_size=tr, - test_size=holdout_ratio, - random_state=seed, - shuffle=True, - stratify=strat, - ) - holdout_strat = None - if strat is not None: - label_map = dict(zip(group_id_arr.tolist(), labels.tolist())) - holdout_labels = np.asarray( - [label_map[int(g)] for g in holdout_groups], - dtype=int, - ) - if len(np.unique(holdout_labels)) >= 2: - holdout_strat = holdout_labels - val_groups, test_groups = train_test_split( - holdout_groups, - train_size=val_ratio_in_holdout, - test_size=1.0 - val_ratio_in_holdout, - random_state=seed + 1, - shuffle=True, - stratify=holdout_strat, - ) - return ( - np.asarray(train_groups, dtype=int), - np.asarray(val_groups, dtype=int), - np.asarray(test_groups, dtype=int), - ) - - split_note = ( - f"按输入分组切分,共 {n_groups} 个唯一输入组;" - f"组级按 {stratify_target} 中位数分 {min(stratify_bins, n_groups)} 桶近似分层。" - ) - try: - train_groups, val_groups, test_groups = _split_groups(use_stratify=True) - report_lines.append(split_note) - except ValueError as e: - train_groups, val_groups, test_groups = _split_groups(use_stratify=False) - report_lines.append(f"{split_note} 但分层条件不足,退化为组级随机切分:{e}") - - i_train = np.flatnonzero(np.isin(group_ids, train_groups)) - i_val = np.flatnonzero(np.isin(group_ids, val_groups)) - i_test = np.flatnonzero(np.isin(group_ids, test_groups)) - return i_train, i_val, i_test - - -def build_split_indices( - df: pd.DataFrame, - cfg: AppConfig, - report_lines: List[str], -) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: - """根据配置生成 train/val/test 行索引。""" - if cfg.split_mode == "random": - report_lines.append("切分策略:随机打乱后按比例切分。") - return _random_split_indices(len(df), cfg.split_ratios, cfg.random_seed) - return _grouped_split_indices( - df=df, - ratios=cfg.split_ratios, - seed=cfg.random_seed, - stratify_target=cfg.split_stratify_target, - stratify_bins=cfg.split_stratify_bins, - report_lines=report_lines, - ) - - -def apply_train_only_outliers( - X_train: np.ndarray, - y_train: np.ndarray, - X_val: np.ndarray, - y_val: np.ndarray, - X_test: np.ndarray, - y_test: np.ndarray, - cfg: AppConfig, - report_lines: List[str], -) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray]: - """ - 仅在训练集上估计阈值: - - iqr/zscore: 从训练集删除离群行(val/test 不动) - - quantile_clip: 对 train/val/test 同步 winsorize(阈值来自 train) - """ - strat = cfg.outlier_strategy - if strat == "none": - report_lines.append("outlier_strategy=none:不对数值做裁剪/删除(除配置项外)。") - return X_train, y_train, X_val, y_val, X_test, y_test - - cols = cfg.outlier_apply_to - if cols == "all": - train_mat = np.hstack([X_train, y_train]) - val_mat = np.hstack([X_val, y_val]) - test_mat = np.hstack([X_test, y_test]) - names = INPUT_COLUMNS + TARGET_COLUMNS - else: - train_mat = y_train.copy() - val_mat = y_val.copy() - test_mat = y_test.copy() - names = TARGET_COLUMNS - - if strat == "quantile_clip": - ql = cfg.outlier_config.quantile_lower - qu = cfg.outlier_config.quantile_upper - tr2, va2, te2 = _winsorize_train_apply_all(train_mat, val_mat, test_mat, ql, qu) - report_lines.append( - f"quantile_clip: 按训练集分位数 [{ql}, {qu}] 对 {cols} 列 winsorize。" - ) - if cols == "all": - d = len(INPUT_COLUMNS) - X_train, y_train = tr2[:, :d], tr2[:, d:] - X_val, y_val = va2[:, :d], va2[:, d:] - X_test, y_test = te2[:, :d], te2[:, d:] - else: - y_train, y_val, y_test = tr2, va2, te2 - return X_train, y_train, X_val, y_val, X_test, y_test - - if strat == "iqr": - mask = _mask_outliers_iqr(train_mat, names, cfg.outlier_config.iqr_k) - elif strat == "zscore": - mask = _mask_outliers_zscore(train_mat, cfg.outlier_config.zscore_threshold) - else: - raise ValueError(f"未知 outlier_strategy: {strat}") - - removed = int(mask.sum()) - kept = ~mask - X_train, y_train = X_train[kept], y_train[kept] - report_lines.append( - f"{strat}: 在训练子集上检测 {cols} 离群,删除训练行 {removed},保留 {len(X_train)}。" - ) - return X_train, y_train, X_val, y_val, X_test, y_test - - def build_dataloaders( X_train: np.ndarray, y_train: np.ndarray, - X_val: np.ndarray, - y_val: np.ndarray, X_test: np.ndarray, y_test: np.ndarray, batch_size: int, num_workers: int, -) -> Tuple[DataLoader, DataLoader, DataLoader]: +) -> Tuple[DataLoader, DataLoader]: def to_loader(X: np.ndarray, y: np.ndarray, shuffle: bool) -> DataLoader: ds = TensorDataset( torch.from_numpy(X).float(), @@ -357,127 +81,77 @@ def build_dataloaders( return ( to_loader(X_train, y_train, shuffle=True), - to_loader(X_val, y_val, shuffle=False), to_loader(X_test, y_test, shuffle=False), ) +def _fit_target_scalers(y_train: np.ndarray) -> Tuple[List[StandardScaler], np.ndarray]: + scalers: List[StandardScaler] = [] + scaled_cols = [] + for i in range(y_train.shape[1]): + scaler = StandardScaler() + scaled_cols.append(scaler.fit_transform(y_train[:, i : i + 1])) + scalers.append(scaler) + return scalers, np.hstack(scaled_cols) + + +def transform_targets(y: np.ndarray, y_scalers: Sequence[StandardScaler]) -> np.ndarray: + cols = [scaler.transform(y[:, i : i + 1]) for i, scaler in enumerate(y_scalers)] + return np.hstack(cols) + + +def inverse_transform_targets( + y_scaled: np.ndarray, y_scalers: Sequence[StandardScaler] +) -> np.ndarray: + cols = [scaler.inverse_transform(y_scaled[:, i : i + 1]) for i, scaler in enumerate(y_scalers)] + return np.hstack(cols) + + def save_scalers( X_scaler: StandardScaler, - y_scaler: StandardScaler, + y_scalers: Sequence[StandardScaler], run_dir: Path, ) -> None: with (run_dir / "x_scaler.pkl").open("wb") as f: pickle.dump(X_scaler, f) - with (run_dir / "y_scaler.pkl").open("wb") as f: - pickle.dump(y_scaler, f) + with (run_dir / "y_scalers.pkl").open("wb") as f: + pickle.dump(list(y_scalers), f) -def load_scalers(run_dir: Path) -> Tuple[StandardScaler, StandardScaler]: +def load_scalers(run_dir: Path) -> Tuple[StandardScaler, List[StandardScaler]]: with (run_dir / "x_scaler.pkl").open("rb") as f: X_scaler = pickle.load(f) - with (run_dir / "y_scaler.pkl").open("rb") as f: - y_scaler = pickle.load(f) - return X_scaler, y_scaler + with (run_dir / "y_scalers.pkl").open("rb") as f: + y_scalers = pickle.load(f) + return X_scaler, list(y_scalers) def save_split_indices( run_dir: Path, i_train: np.ndarray, - i_val: np.ndarray, i_test: np.ndarray, ) -> None: - """保存对清洗后矩阵行的划分索引,便于 eval 阶段完全复现。""" payload = { "train": i_train.astype(int).tolist(), - "val": i_val.astype(int).tolist(), "test": i_test.astype(int).tolist(), } (run_dir / "split_indices.json").write_text( - json.dumps(payload, indent=2), encoding="utf-8" + json.dumps(payload, indent=2), + encoding="utf-8", ) -def load_split_indices(run_dir: Path) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: +def load_split_indices(run_dir: Path) -> Tuple[np.ndarray, np.ndarray]: path = run_dir / "split_indices.json" if not path.is_file(): - raise FileNotFoundError( - f"未找到 {path}。请使用本仓库训练产生的 run 目录,或先完成一次训练。" - ) + raise FileNotFoundError(f"未找到 {path}") data = json.loads(path.read_text(encoding="utf-8")) return ( np.asarray(data["train"], dtype=int), - np.asarray(data["val"], dtype=int), np.asarray(data["test"], dtype=int), ) -def rebuild_bundle_for_eval( - df: pd.DataFrame, - cfg: AppConfig, - run_dir: Path, -) -> ProcessedDataBundle: - """ - 与训练阶段相同的清洗、划分与离群处理,但使用已保存的 StandardScaler 仅做 transform。 - 用于独立 eval / infer 流程,避免重新拟合 scaler 造成分布偏移。 - """ - report_lines: List[str] = [] - cleaned = clean_dataframe(df, cfg, report_lines) - X = cleaned[INPUT_COLUMNS].to_numpy(dtype=np.float64) - y = cleaned[TARGET_COLUMNS].to_numpy(dtype=np.float64) - i_tr, i_va, i_te = load_split_indices(run_dir) - for name, idx in ("train", i_tr), ("val", i_va), ("test", i_te): - if len(idx) == 0 or int(idx.max()) >= len(X) or int(idx.min()) < 0: - raise ValueError( - f"split_indices.json 与当前数据不兼容({name} 索引越界或为空)。" - f"请确认 data_path 指向与训练相同的清洗后样本空间。" - ) - - X_train, y_train = X[i_tr], y[i_tr] - X_val, y_val = X[i_va], y[i_va] - X_test, y_test = X[i_te], y[i_te] - X_train, y_train, X_val, y_val, X_test, y_test = apply_train_only_outliers( - X_train, y_train, X_val, y_val, X_test, y_test, cfg, report_lines - ) - - X_scaler, y_scaler = load_scalers(run_dir) - X_train_s = X_scaler.transform(X_train) - y_train_s = y_scaler.transform(y_train) - X_val_s = X_scaler.transform(X_val) - y_val_s = y_scaler.transform(y_val) - X_test_s = X_scaler.transform(X_test) - y_test_s = y_scaler.transform(y_test) - - train_loader, val_loader, test_loader = build_dataloaders( - X_train_s, - y_train_s, - X_val_s, - y_val_s, - X_test_s, - y_test_s, - cfg.training.batch_size, - cfg.training.num_workers, - ) - - return ProcessedDataBundle( - train_loader=train_loader, - val_loader=val_loader, - test_loader=test_loader, - X_train=X_train_s, - X_val=X_val_s, - X_test=X_test_s, - y_train=y_train_s, - y_val=y_val_s, - y_test=y_test_s, - X_test_raw=X_test, - y_test_raw=y_test, - X_scaler=X_scaler, - y_scaler=y_scaler, - feature_names=list(INPUT_COLUMNS), - target_names=list(TARGET_COLUMNS), - ) - - def write_data_report_md( path: Path, raw_quality: dict, @@ -515,12 +189,9 @@ def prepare_training_data( cfg: AppConfig, run_dir: Path, ) -> ProcessedDataBundle: - """ - 完整预处理流水线:质量报告 -> 清洗 -> 划分 -> 训练集离群处理 -> 标准化 -> DataLoader。 - 将 data_report.md 与 cleaning 元数据写入 run_dir。 - """ + """按 notebook 一致逻辑准备 train/test、scaler 与 DataLoader。""" report_lines: List[str] = [] - raw_q = quality_report_before_clean(df, cfg.v_pi_min, cfg.v_pi_max) + raw_q = quality_report_before_clean(df, 0.0, cfg.data.filter_v_pi_max) stats_before = df.describe().T cleaned = clean_dataframe(df, cfg, report_lines) @@ -528,39 +199,35 @@ def prepare_training_data( X = cleaned[INPUT_COLUMNS].to_numpy(dtype=np.float64) y = cleaned[TARGET_COLUMNS].to_numpy(dtype=np.float64) + all_idx = np.arange(len(X)) - i_tr, i_va, i_te = build_split_indices(cleaned, cfg, report_lines) - save_split_indices(run_dir, i_tr, i_va, i_te) - X_train, y_train = X[i_tr], y[i_tr] - X_val, y_val = X[i_va], y[i_va] - X_test, y_test = X[i_te], y[i_te] + X_train_raw, X_test_raw, y_train_raw, y_test_raw, i_train, i_test = train_test_split( + X, + y, + all_idx, + test_size=cfg.data.test_size, + random_state=cfg.data.random_state, + ) report_lines.append( - f"划分 train/val/test = {cfg.split_ratios},样本数 " - f"{len(X_train)}/{len(X_val)}/{len(X_test)}" - ) - - X_train, y_train, X_val, y_val, X_test, y_test = apply_train_only_outliers( - X_train, y_train, X_val, y_val, X_test, y_test, cfg, report_lines + f"train_test_split(test_size={cfg.data.test_size}, random_state={cfg.data.random_state}) " + f"-> {len(X_train_raw)}/{len(X_test_raw)}" ) + save_split_indices(run_dir, i_train, i_test) X_scaler = StandardScaler() - y_scaler = StandardScaler() - X_train_s = X_scaler.fit_transform(X_train) - y_train_s = y_scaler.fit_transform(y_train) - X_val_s = X_scaler.transform(X_val) - y_val_s = y_scaler.transform(y_val) - X_test_s = X_scaler.transform(X_test) - y_test_s = y_scaler.transform(y_test) - - save_scalers(X_scaler, y_scaler, run_dir) + X_train = X_scaler.fit_transform(X_train_raw) + X_test = X_scaler.transform(X_test_raw) + y_scalers, y_train = _fit_target_scalers(y_train_raw) + y_test = transform_targets(y_test_raw, y_scalers) + save_scalers(X_scaler, y_scalers, run_dir) meta = { "raw_quality": raw_q, "cleaning_steps": report_lines, - "split_ratios": cfg.split_ratios, - "n_train": int(len(X_train_s)), - "n_val": int(len(X_val_s)), - "n_test": int(len(X_test_s)), + "test_size": cfg.data.test_size, + "random_state": cfg.data.random_state, + "n_train": int(len(X_train)), + "n_test": int(len(X_test)), } (run_dir / "cleaning_meta.json").write_text( json.dumps(meta, indent=2, ensure_ascii=False, default=str), @@ -576,31 +243,77 @@ def prepare_training_data( stats_after.to_csv(run_dir / "data_stats.csv", encoding="utf-8") logger.info("预处理完成:%s", run_dir / "data_report.md") - train_loader, val_loader, test_loader = build_dataloaders( - X_train_s, - y_train_s, - X_val_s, - y_val_s, - X_test_s, - y_test_s, + train_loader, test_loader = build_dataloaders( + X_train, + y_train, + X_test, + y_test, cfg.training.batch_size, cfg.training.num_workers, ) return ProcessedDataBundle( train_loader=train_loader, - val_loader=val_loader, test_loader=test_loader, - X_train=X_train_s, - X_val=X_val_s, - X_test=X_test_s, - y_train=y_train_s, - y_val=y_val_s, - y_test=y_test_s, - X_test_raw=X_test, - y_test_raw=y_test, + X_train=X_train, + X_test=X_test, + y_train=y_train, + y_test=y_test, + X_train_raw=X_train_raw, + X_test_raw=X_test_raw, + y_train_raw=y_train_raw, + y_test_raw=y_test_raw, X_scaler=X_scaler, - y_scaler=y_scaler, + y_scalers=y_scalers, + feature_names=list(INPUT_COLUMNS), + target_names=list(TARGET_COLUMNS), + ) + + +def rebuild_bundle_for_eval( + df: pd.DataFrame, + cfg: AppConfig, + run_dir: Path, +) -> ProcessedDataBundle: + """使用训练时保存的切分索引与 scaler 重新构造 train/test 数据。""" + report_lines: List[str] = [] + cleaned = clean_dataframe(df, cfg, report_lines) + X = cleaned[INPUT_COLUMNS].to_numpy(dtype=np.float64) + y = cleaned[TARGET_COLUMNS].to_numpy(dtype=np.float64) + i_train, i_test = load_split_indices(run_dir) + for name, idx in (("train", i_train), ("test", i_test)): + if len(idx) == 0 or int(idx.max()) >= len(X) or int(idx.min()) < 0: + raise ValueError(f"split_indices.json 与当前数据不兼容({name} 索引越界或为空)") + + X_train_raw, X_test_raw = X[i_train], X[i_test] + y_train_raw, y_test_raw = y[i_train], y[i_test] + X_scaler, y_scalers = load_scalers(run_dir) + X_train = X_scaler.transform(X_train_raw) + X_test = X_scaler.transform(X_test_raw) + y_train = transform_targets(y_train_raw, y_scalers) + y_test = transform_targets(y_test_raw, y_scalers) + train_loader, test_loader = build_dataloaders( + X_train, + y_train, + X_test, + y_test, + cfg.training.batch_size, + cfg.training.num_workers, + ) + + return ProcessedDataBundle( + train_loader=train_loader, + test_loader=test_loader, + X_train=X_train, + X_test=X_test, + y_train=y_train, + y_test=y_test, + X_train_raw=X_train_raw, + X_test_raw=X_test_raw, + y_train_raw=y_train_raw, + y_test_raw=y_test_raw, + X_scaler=X_scaler, + y_scalers=y_scalers, feature_names=list(INPUT_COLUMNS), target_names=list(TARGET_COLUMNS), ) diff --git a/src/trainer.py b/src/trainer.py index 5d32660..80931e6 100644 --- a/src/trainer.py +++ b/src/trainer.py @@ -1,4 +1,4 @@ -"""训练循环、早停、调度器与 checkpoint。""" +"""Notebook 风格的 MoE + autograd PINN 训练循环。""" from __future__ import annotations @@ -6,17 +6,15 @@ import csv import logging from dataclasses import dataclass from pathlib import Path -from typing import Dict, List, Optional, Tuple +from typing import Dict, List, Tuple +import numpy as np import torch import torch.nn as nn from torch.optim import AdamW -from torch.optim.lr_scheduler import CosineAnnealingLR, ReduceLROnPlateau -from tqdm import tqdm from src.config import AppConfig -from src.losses import build_loss -from src.model import MLPRegressor +from src.preprocess import ProcessedDataBundle logger = logging.getLogger(__name__) @@ -25,8 +23,7 @@ logger = logging.getLogger(__name__) class TrainHistory: epoch: List[int] train_loss: List[float] - val_loss: List[float] - lr: List[float] + test_loss: List[float] def _move_batch( @@ -36,145 +33,172 @@ def _move_batch( return x.to(device), y.to(device) -def train_one_epoch( - model: nn.Module, - loader: torch.utils.data.DataLoader, - criterion: nn.Module, - optimizer: torch.optim.Optimizer, - device: torch.device, -) -> float: - model.train() - total, n = 0.0, 0 - for batch in loader: - xb, yb = _move_batch(batch, device) - optimizer.zero_grad(set_to_none=True) - pred = model(xb) - loss = criterion(pred, yb) - loss.backward() - optimizer.step() - total += float(loss.detach().cpu()) * xb.size(0) - n += xb.size(0) - return total / max(n, 1) - - -@torch.no_grad() -def evaluate_loss_loader( - model: nn.Module, - loader: torch.utils.data.DataLoader, - criterion: nn.Module, - device: torch.device, -) -> float: - model.eval() - total, n = 0.0, 0 - for batch in loader: - xb, yb = _move_batch(batch, device) - pred = model(xb) - loss = criterion(pred, yb) - total += float(loss.detach().cpu()) * xb.size(0) - n += xb.size(0) - return total / max(n, 1) - - -def build_optimizer_and_scheduler( - model: nn.Module, cfg: AppConfig -) -> Tuple[AdamW, object]: - opt = AdamW( +def build_optimizer(model: nn.Module, cfg: AppConfig) -> AdamW: + beta1, beta2 = cfg.optimizer.betas + return AdamW( model.parameters(), lr=cfg.optimizer.lr, weight_decay=cfg.optimizer.weight_decay, + betas=(beta1, beta2), ) - if cfg.scheduler.type == "cosine": - sched: torch.optim.lr_scheduler._LRScheduler = CosineAnnealingLR( - opt, T_max=cfg.training.epochs, eta_min=cfg.scheduler.plateau_min_lr - ) - elif cfg.scheduler.type == "plateau": - sched = ReduceLROnPlateau( - opt, - mode="min", - factor=cfg.scheduler.plateau_factor, - patience=cfg.scheduler.plateau_patience, - min_lr=cfg.scheduler.plateau_min_lr, - ) + + +def compute_pinn_loss( + model: nn.Module, + x_batch: torch.Tensor, + y_batch: torch.Tensor, + criterion: nn.Module, + cfg: AppConfig, +) -> Tuple[torch.Tensor, torch.Tensor, Dict[str, torch.Tensor]]: + x_in = x_batch.detach().clone().requires_grad_(True) + pred = model(x_in) + data_loss = criterion(pred, y_batch) + + bw_pred = pred[:, 0] + il_pred = pred[:, 1] + vpi_pred = pred[:, 2] + length_batch = x_in[:, 7] + + losses: Dict[str, torch.Tensor] = {} + total = data_loss + + if cfg.physics.lambda_bw_mon != 0: + grads = torch.autograd.grad( + bw_pred, + x_in, + grad_outputs=torch.ones_like(bw_pred), + create_graph=True, + )[0] + d_bw_d_l = grads[:, 7] + losses["bw_mon"] = torch.mean(torch.relu(d_bw_d_l) ** 2) + total = total + cfg.physics.lambda_bw_mon * losses["bw_mon"] + + if cfg.physics.lambda_IL_mon != 0: + grads = torch.autograd.grad( + il_pred, + x_in, + grad_outputs=torch.ones_like(il_pred), + create_graph=True, + )[0] + d_il_d_l = grads[:, 7] + losses["IL_mon"] = torch.mean(torch.relu(-d_il_d_l) ** 2) + total = total + cfg.physics.lambda_IL_mon * losses["IL_mon"] + + if cfg.physics.lambda_vpiL != 0: + vpi_l = vpi_pred * length_batch + grads = torch.autograd.grad( + vpi_l, + x_in, + grad_outputs=torch.ones_like(vpi_l), + create_graph=True, + )[0] + d_vpi_l_d_l = grads[:, 7] + losses["vpiL"] = torch.mean(d_vpi_l_d_l**2) + total = total + cfg.physics.lambda_vpiL * losses["vpiL"] + + if cfg.physics.lambda_smooth != 0: + grads1 = torch.autograd.grad( + bw_pred, + x_in, + grad_outputs=torch.ones_like(bw_pred), + create_graph=True, + )[0] + d_bw_d_l = grads1[:, 7] + grads2 = torch.autograd.grad( + d_bw_d_l, + x_in, + grad_outputs=torch.ones_like(d_bw_d_l), + create_graph=True, + )[0] + d2_bw_d_l2 = grads2[:, 7] + losses["smooth"] = torch.mean(d2_bw_d_l2**2) + total = total + cfg.physics.lambda_smooth * losses["smooth"] + + return total, data_loss, losses + + +@torch.no_grad() +def evaluate_full_batch_mse( + model: nn.Module, + x: np.ndarray | torch.Tensor, + y: np.ndarray | torch.Tensor, + device: torch.device, +) -> float: + model.eval() + if isinstance(x, torch.Tensor): + x_t = x.to(device) else: - raise ValueError(f"未知 scheduler.type: {cfg.scheduler.type}") - return opt, sched + x_t = torch.from_numpy(x).float().to(device) + if isinstance(y, torch.Tensor): + y_t = y.to(device) + else: + y_t = torch.from_numpy(y).float().to(device) + pred = model(x_t) + return float(nn.functional.mse_loss(pred, y_t).item()) def fit( model: nn.Module, cfg: AppConfig, - train_loader: torch.utils.data.DataLoader, - val_loader: torch.utils.data.DataLoader, + bundle: ProcessedDataBundle, run_dir: Path, device: torch.device, ) -> TrainHistory: - """ - 训练模型:早停依据验证集损失;保存 best / last 权重到 run_dir/checkpoints。 - 同步写入 train_log.csv。 - """ - criterion = build_loss(cfg.loss).to(device) - optimizer, scheduler = build_optimizer_and_scheduler(model, cfg) + """按 notebook 风格训练,并记录每轮全量 train/test MSE。""" + criterion = nn.MSELoss().to(device) + optimizer = build_optimizer(model, cfg) ckpt_dir = run_dir / "checkpoints" ckpt_dir.mkdir(parents=True, exist_ok=True) log_path = run_dir / "train_log.csv" - best_val = float("inf") - best_epoch = -1 - patience_left = cfg.training.early_stopping_patience + x_train_full = torch.from_numpy(bundle.X_train).float().to(device) + y_train_full = torch.from_numpy(bundle.y_train).float().to(device) + x_test_full = torch.from_numpy(bundle.X_test).float().to(device) + y_test_full = torch.from_numpy(bundle.y_test).float().to(device) - hist = TrainHistory(epoch=[], train_loss=[], val_loss=[], lr=[]) + hist = TrainHistory(epoch=[], train_loss=[], test_loss=[]) + best_test = float("inf") with log_path.open("w", newline="", encoding="utf-8") as fcsv: writer = csv.writer(fcsv) - writer.writerow(["epoch", "train_loss", "val_loss", "lr", "best_val"]) + writer.writerow(["epoch", "train_loss", "test_loss"]) - for epoch in range(1, cfg.training.epochs + 1): - tr_loss = train_one_epoch(model, train_loader, criterion, optimizer, device) - va_loss = evaluate_loss_loader(model, val_loader, criterion, device) + for epoch in range(cfg.training.epochs): + model.train() + for x_batch, y_batch in bundle.train_loader: + x_batch, y_batch = _move_batch((x_batch, y_batch), device) + optimizer.zero_grad() + loss, _, _ = compute_pinn_loss(model, x_batch, y_batch, criterion, cfg) + loss.backward() + optimizer.step() - if cfg.scheduler.type == "cosine": - scheduler.step() - elif cfg.scheduler.type == "plateau": - scheduler.step(va_loss) + train_loss = evaluate_full_batch_mse(model, x_train_full, y_train_full, device) + test_loss = evaluate_full_batch_mse(model, x_test_full, y_test_full, device) - lr_now = float(optimizer.param_groups[0]["lr"]) hist.epoch.append(epoch) - hist.train_loss.append(tr_loss) - hist.val_loss.append(va_loss) - hist.lr.append(lr_now) - - improved = va_loss + 1e-12 < best_val - if improved: - best_val = va_loss - best_epoch = epoch - patience_left = cfg.training.early_stopping_patience - torch.save( - {"epoch": epoch, "model_state": model.state_dict(), "val_loss": va_loss}, - ckpt_dir / "best.pt", - ) - else: - patience_left -= 1 - - writer.writerow([epoch, tr_loss, va_loss, lr_now, best_val]) + hist.train_loss.append(train_loss) + hist.test_loss.append(test_loss) + writer.writerow([epoch, train_loss, test_loss]) fcsv.flush() - logger.info( - "Epoch %d | train_loss=%.6f val_loss=%.6f | best_val=%.6f @%d", - epoch, - tr_loss, - va_loss, - best_val, - best_epoch, - ) + if epoch % 10 == 0 or epoch == 0: + logger.info( + "Epoch %5d | Train %.6f | Test %.6f", + epoch, + train_loss, + test_loss, + ) - torch.save( - {"epoch": epoch, "model_state": model.state_dict(), "val_loss": va_loss}, - ckpt_dir / "last.pt", - ) - - if patience_left <= 0: - logger.info("早停触发于 epoch %d,最佳 epoch=%d", epoch, best_epoch) - break + payload = { + "epoch": epoch, + "model_state": model.state_dict(), + "train_loss": train_loss, + "test_loss": test_loss, + } + torch.save(payload, ckpt_dir / "last.pt") + if test_loss < best_test: + best_test = test_loss + torch.save(payload, ckpt_dir / "best.pt") return hist diff --git a/tests/test_smoke.py b/tests/test_smoke.py index 7ac23c8..253999f 100644 --- a/tests/test_smoke.py +++ b/tests/test_smoke.py @@ -1,4 +1,4 @@ -"""最小冒烟测试:模型 shape 与单 epoch 训练不报错。""" +"""最小冒烟测试:MoE + PINN 流程可运行。""" from __future__ import annotations @@ -10,82 +10,93 @@ import torch import yaml from src.config import load_config -from src.data import INPUT_COLUMNS from src.data import load_raw_txt -from src.model import MLPRegressor -from src.preprocess import prepare_training_data -from src.trainer import fit +from src.model import create_model_from_config +from src.preprocess import inverse_transform_targets, prepare_training_data +from src.trainer import compute_pinn_loss, fit -def _write_synthetic_txt(path: Path, n: int = 64) -> None: +def _write_synthetic_txt(path: Path, n: int = 96) -> None: rng = np.random.default_rng(0) x = rng.normal(size=(n, 8)) y = np.zeros((n, 3)) - y[:, 0] = rng.normal(size=n) - y[:, 1] = rng.normal(size=n) - # 第 11 列 V_pi 落在默认物理门控 [0, 500] 内 - y[:, 2] = rng.uniform(1.0, 400.0, size=n) + length = np.abs(x[:, 7]) + 0.5 + y[:, 0] = 2.0 - 0.2 * length + rng.normal(scale=0.05, size=n) + y[:, 1] = 0.3 + 0.1 * length + rng.normal(scale=0.03, size=n) + y[:, 2] = 20.0 / length + rng.normal(scale=0.2, size=n) mat = np.hstack([x, y]) - lines = [",".join(str(v) for v in row) for row in mat] - path.write_text("\n".join(lines), encoding="utf-8") + path.write_text("\n".join(",".join(str(v) for v in row) for row in mat), encoding="utf-8") -def test_mlp_forward_shape() -> None: - m = MLPRegressor(8, [16, 16], 3, batchnorm=False, dropout=0.0, residual=False) +def _make_cfg(tmp_path: Path, data_txt: Path, epochs: int = 1) -> Path: + cfg_dict = { + "data_path": str(data_txt), + "data": { + "test_size": 0.1, + "random_state": 123, + "filter_v_pi_max": 500.0, + }, + "model": { + "input_dim": 8, + "output_dim": 3, + "hidden_dims": [16, 16], + "n_experts": 4, + "gating_hidden": 4, + "dropout_rate": 0.0, + "use_bn": False, + "activation": "relu", + }, + "optimizer": { + "lr": 0.001, + "weight_decay": 0.01, + "betas": [0.9, 0.999], + }, + "training": { + "batch_size": 16, + "epochs": epochs, + "num_workers": 0, + }, + "physics": { + "lambda_bw_mon": 0.1, + "lambda_IL_mon": 0.1, + "lambda_vpiL": 0.05, + "lambda_smooth": 0.01, + }, + "output_dir": str(tmp_path / "results"), + } + cfg_path = tmp_path / "cfg.yaml" + cfg_path.write_text(yaml.safe_dump(cfg_dict), encoding="utf-8") + return cfg_path + + +def test_moe_forward_shape(tmp_path: Path) -> None: + data_txt = tmp_path / "data.txt" + _write_synthetic_txt(data_txt, n=32) + cfg = load_config(_make_cfg(tmp_path, data_txt)) + model = create_model_from_config(cfg) x = torch.randn(5, 8) - y = m(x) + y = model(x) assert y.shape == (5, 3) +def test_pinn_loss_backpropagates(tmp_path: Path) -> None: + data_txt = tmp_path / "data.txt" + _write_synthetic_txt(data_txt, n=40) + cfg = load_config(_make_cfg(tmp_path, data_txt)) + model = create_model_from_config(cfg) + x = torch.randn(8, 8) + y = torch.randn(8, 3) + loss, data_loss, terms = compute_pinn_loss(model, x, y, torch.nn.MSELoss(), cfg) + loss.backward() + assert float(loss.detach()) >= float(data_loss.detach()) + assert "IL_mon" in terms + assert any(p.grad is not None for p in model.parameters()) + + def test_one_epoch_training_pipeline(tmp_path: Path) -> None: data_txt = tmp_path / "data.txt" _write_synthetic_txt(data_txt, n=80) - - cfg_dict = { - "data_path": str(data_txt), - "split_ratios": [0.7, 0.15, 0.15], - "random_seed": 1, - "remove_duplicate_rows": False, - "outlier_strategy": "none", - "outlier_apply_to": "targets", - "outlier_config": { - "iqr_k": 1.5, - "zscore_threshold": 4.0, - "quantile_lower": 0.001, - "quantile_upper": 0.999, - }, - "remove_nonpositive_vpi": False, - "filter_v_pi_range": True, - "v_pi_min": 0.0, - "v_pi_max": 500.0, - "model": { - "input_dim": 8, - "hidden_dims": [32, 32], - "output_dim": 3, - "batchnorm": False, - "dropout": 0.0, - "residual": False, - }, - "optimizer": {"name": "adamw", "lr": 0.01, "weight_decay": 0.0}, - "scheduler": { - "type": "cosine", - "plateau_factor": 0.5, - "plateau_patience": 10, - "plateau_min_lr": 1e-6, - }, - "training": { - "batch_size": 16, - "epochs": 1, - "early_stopping_patience": 1, - "num_workers": 0, - }, - "loss": {"type": "huber", "huber_delta": 1.0, "target_weights": [1.0, 1.0, 1.0]}, - "output_dir": str(tmp_path / "results"), - } - cfg_path = tmp_path / "cfg.yaml" - cfg_path.write_text(yaml.safe_dump(cfg_dict), encoding="utf-8") - - cfg = load_config(cfg_path) + cfg = load_config(_make_cfg(tmp_path, data_txt, epochs=1)) df = load_raw_txt(data_txt) run_dir = tmp_path / "run0" @@ -95,99 +106,12 @@ def test_one_epoch_training_pipeline(tmp_path: Path) -> None: bundle = prepare_training_data(df, cfg, run_dir) device = torch.device("cpu") - model = MLPRegressor( - input_dim=cfg.model.input_dim, - hidden_dims=cfg.model.hidden_dims, - output_dim=cfg.model.output_dim, - batchnorm=cfg.model.batchnorm, - dropout=cfg.model.dropout, - residual=cfg.model.residual, - ).to(device) - fit(model, cfg, bundle.train_loader, bundle.val_loader, run_dir, device) - assert (run_dir / "checkpoints" / "best.pt").is_file() + model = create_model_from_config(cfg).to(device) + history = fit(model, cfg, bundle, run_dir, device) + assert (run_dir / "checkpoints" / "last.pt").is_file() + assert len(history.train_loss) == 1 meta = json.loads((run_dir / "cleaning_meta.json").read_text(encoding="utf-8")) assert meta["n_train"] > 0 - - -def test_grouped_split_keeps_same_inputs_together(tmp_path: Path) -> None: - data_txt = tmp_path / "grouped_data.txt" - rng = np.random.default_rng(7) - rows = [] - base_inputs = rng.normal(size=(24, 8)) - for x in base_inputs: - for _ in range(3): - y0 = float(x[0] * 2.0 + rng.normal(scale=0.01)) - y1 = float(x[1] * -1.5 + rng.normal(scale=0.01)) - y2 = float(abs(x[2]) * 20.0 + 10.0 + rng.normal(scale=0.1)) - rows.append(np.concatenate([x, [y0, y1, y2]])) - mat = np.asarray(rows, dtype=float) - data_txt.write_text( - "\n".join(",".join(str(v) for v in row) for row in mat), - encoding="utf-8", - ) - - cfg_dict = { - "data_path": str(data_txt), - "split_ratios": [0.7, 0.15, 0.15], - "random_seed": 3, - "split_mode": "grouped_stratified", - "split_stratify_target": "V_pi", - "split_stratify_bins": 6, - "remove_duplicate_rows": False, - "outlier_strategy": "none", - "outlier_apply_to": "targets", - "outlier_config": { - "iqr_k": 1.5, - "zscore_threshold": 4.0, - "quantile_lower": 0.001, - "quantile_upper": 0.999, - }, - "remove_nonpositive_vpi": False, - "filter_v_pi_range": True, - "v_pi_min": 0.0, - "v_pi_max": 500.0, - "model": { - "input_dim": 8, - "hidden_dims": [16, 16], - "output_dim": 3, - "batchnorm": False, - "dropout": 0.0, - "residual": False, - }, - "optimizer": {"name": "adamw", "lr": 0.01, "weight_decay": 0.0}, - "scheduler": { - "type": "cosine", - "plateau_factor": 0.5, - "plateau_patience": 10, - "plateau_min_lr": 1e-6, - }, - "training": { - "batch_size": 16, - "epochs": 1, - "early_stopping_patience": 1, - "num_workers": 0, - }, - "loss": {"type": "huber", "huber_delta": 1.0, "target_weights": [1.0, 1.0, 1.0]}, - "output_dir": str(tmp_path / "results"), - } - cfg_path = tmp_path / "cfg_grouped.yaml" - cfg_path.write_text(yaml.safe_dump(cfg_dict), encoding="utf-8") - - cfg = load_config(cfg_path) - df = load_raw_txt(data_txt) - run_dir = tmp_path / "run_grouped" - run_dir.mkdir() - bundle = prepare_training_data(df, cfg, run_dir) - - split_data = json.loads((run_dir / "split_indices.json").read_text(encoding="utf-8")) - split_name_by_row = {} - for split_name, indices in split_data.items(): - for idx in indices: - split_name_by_row[int(idx)] = split_name - - cleaned = df.reset_index(drop=True) - for _, sub in cleaned.groupby(INPUT_COLUMNS, dropna=False): - assigned = {split_name_by_row[int(i)] for i in sub.index.to_list()} - assert len(assigned) == 1 - - assert len(bundle.X_train) > 0 + assert meta["n_test"] > 0 + restored = inverse_transform_targets(bundle.y_test[:3], bundle.y_scalers) + assert restored.shape == (3, 3)