1. Overview

This Jupyter notebook implements a Physics-Informed Neural Network (PINN) prediction pipeline for Vapor-Liquid Equilibrium (VLE) data using COSMO-based molecular descriptors combined with the Margules activity coefficient model. The notebook loads pre-trained models and generates predictions for binary mixture thermodynamic properties (vapor composition y and system pressure P) across ~17,600 data points.

Domain Context

  • COSMO (COnductor-like Screening MOdel): A quantum chemistry method that produces sigma-profiles (charge density distributions) for molecules, used here as 51-dimensional molecular descriptors.
  • Margules Equation: A thermodynamic model for activity coefficients in binary mixtures, incorporated as a physics constraint layer in the neural network.
  • VLE Prediction: Given two compounds (A and B), temperature (T), and liquid composition (xA), the model predicts vapor composition (yA) and system pressure (P).

2. Notebook Structure (15 Code Cells)

Cell 1 — Initial Setup & Helper Functions

Purpose: Imports core libraries, defines file paths, and creates three helper functions.

Key Functions: - training_data(number): Reads a 51-point sigma-profile from a CSV file identified by compound ID. Returns zero vector if file is missing. - name_to_cosmo(COSMO_num): Looks up a compound ID in a mapping CSV (CID ↔ SMILES). Essentially an identity function that validates the compound exists. - train_datasets(fluids, T, x): Assembles a 110-dimensional feature vector by concatenating two sigma-profiles (51×2=102), temperature (1), mole fraction xA (1), and 6 placeholder zeros for P, VpA, VpB, gammaA, gammaB, and a spare slot.

Feature Vector Layout (110 dimensions):

Index Content
0–50 Sigma-profile of compound A
51–101 Sigma-profile of compound B
102 Temperature T (K)
103 Mole fraction xA
104 yA (target, placeholder)
105 ln(P) (pressure)
106 ln(VpA) (vapor pressure A)
107 ln(VpB) (vapor pressure B)
108 gammaA (activity coeff. A)
109 gammaB (activity coeff. B)

Cell 2 — Commented-Out Multi-Fold Prediction (MLP V9)

Purpose: An earlier, fully commented-out attempt to iterate over multiple cross-validation folds, load each model, and run predictions. Uses a simpler MLP (not PINN) with dual input structure [in_feat, in_phys].

Status: Entirely commented out. Not used.

Cell 3 — Single-Fold MLP V9 Prediction

Purpose: Runs row-by-row prediction using a simple MLP model (V9 version) without the Margules physics constraint layer. Writes results to CSV one row at a time. Uses MinMaxScaler for feature normalization and inverse-transforms predictions back to physical units.

Key Logic: - Pressure unit detection: P > 500 is assumed to be in Pa, otherwise kPa. - Inverse normalization: y_calc = pred * (max - min) + min for vapor composition; P = exp(pred * (max - min) + min) for pressure.

Cell 4 — Commented-Out Pathlib Version

Purpose: Another commented-out version attempting to use pathlib.Path for cleaner path handling, and a more explicit log transformation for vapor pressures. Not executed.

Cell 5 — Import Block for PINN Model

Purpose: Imports TensorFlow/Keras components needed for the PINN architecture: Model, Input, Dense, Concatenate, Lambda, and Keras backend K.

Cell 6 — Path Configuration for PINN (V8_7)

Purpose: Sets up all file paths for the PINN model version (V8, variant 7). Includes paths to: - Weight files (.weights.h5) - Hyperparameter pickle files (.pkl) - Preprocessing pipeline (.joblib) - Input/output CSV files

Cell 7 — Core PINN Architecture + First Prediction Attempt

Purpose: Defines the central components of the PINN model:

  1. hard_constraint_layer_v17(inputs, **kwargs) — The Margules physics constraint layer: - Receives MLP raw output [A12, A21, y_mlp, P_delta] and physical inputs [xA, lnPsatA, lnPsatB, T] - De-normalizes physical inputs from MinMax-scaled space - Computes Margules activity coefficients: (\ln \gamma_A = x_B^2 (A_{12} + 2(A_{21} - A_{12}) x_A)) - Calculates total pressure via modified Raoult's law: (P = x_A \gamma_A P_{sat,A} + x_B \gamma_B P_{sat,B}) - Returns [y_mlp, P_normalized + P_delta, gamma_A, gamma_B]

  2. build_and_load_model(params, fold_idx, scaler) — Builds a simpler model with uniform layer widths (all n_neurons) and tanh activation.

  3. First prediction loop — Loads model for fold 1, runs predictions, outputs results.

Cell 8 — PKL-Based Model Loading with Fallback

Purpose: Loads hyperparameters from a pickle file, handles key name mismatches (n_neurons vs neurons vs units), and loads weights with a fallback path attempt.

Cell 9 — Complete PINN Pipeline (V8_7, Fold 1)

Purpose: Self-contained version with: - build_pinn_model_from_params() using per-layer neuron counts (u0, u1, ...) and swish activation - Output branches: tanh(A12, A21) * 3.0, clip(yA, 0, 1), P_delta - Diagnostic print statements for normalization verification - Row-by-row prediction loop

Cell 10 — Near-Identical Variant

Purpose: Same as Cell 9 with minor differences: - Different output filename (_Fix.csv) - Uses np.asarray() to fix a potential InvalidIndexError - No diagnostic prints

Cell 11 — Synchronized Prediction with Real-Time Output

Purpose: Same prediction logic with: - Version/fold changed to fold_idx = 5 - Real-time console output with formatted table - Cumulative error tracking (total_ad_y, total_ard_p) - Final summary statistics (AAD-y, AARD-P)

Note: Has a bug in p_args — 'max_gammaB': sc.data_min_[-1] uses data_min_ instead of data_max_.

Cell 12 — Final Executed Prediction (17,656 rows)

Purpose: The only cell that was actually executed to completion. Produces the main results file. Uses fold 5 with gamma placeholders set to 1.0 (which differs from other cells that use 0). This cell generated 17,656 predictions.

Final Results: AAD-y = 15.98%, AARD-P = 159.43%

Cell 13 — Batch Prediction (Vectorized)

Purpose: Attempts batch prediction using model.predict() on the entire dataset at once (batch_size=256) rather than row-by-row. Achieves the same AARD-P of 159.43%.

Cell 14 — Scaler Boundary Check

Purpose: Reloads the scaler and prints boundary values to verify feature alignment. Uses grouping_alldata_V3.csv (different version from other cells).

Cell 15 — Gamma Training Memory Check

Purpose: Prints min/max values for gamma columns (indices 108, 109) from the scaler, revealing the training data range: gammaA ∈ [0.0001, 45354.2], gammaB ∈ [0.025, 1083.9].

3. Data Flow Diagram

Input CSV (grouping_alldata_V2.csv)
    │
    ├── comp.A, comp.B → Sigma-profile CSVs (51-dim each)
    ├── T, xA → Physical conditions
    ├── P_exp, VpA, VpB → Experimental values
    │
    ▼
Feature Vector Assembly (110 dim)
    │  [sigma_A(51) | sigma_B(51) | T | xA | yA | lnP | lnVpA | lnVpB | gammaA | gammaB]
    │
    ▼
MinMaxScaler (from training pipeline)
    │
    ├── X_f = scaled[:, :104]          → MLP Input (sigma-profiles + T + xA)
    └── X_p = scaled[:, [103,106,107,102]] → Physics Input [xA, lnVpA, lnVpB, T]
         │
         ▼
    PINN Model
    ┌─────────────────────────────────────────┐
    │  MLP Backbone (swish, variable layers)  │
    │       ↓                                 │
    │  Raw Output: [A12, A21, yA, P_delta]    │
    │       ↓                                 │
    │  hard_constraint_layer_v17              │
    │  (Margules equation + Raoult's law)     │
    │       ↓                                 │
    │  Output: [y_pred, P_pred, γA, γB]       │
    └─────────────────────────────────────────┘
         │
         ▼
    Inverse Normalization
         │
         ├── y_cal = pred[0] * (max_y - min_y) + min_y
         └── P_cal = exp(pred[1] * (max_p - min_p) + min_p)
         │
         ▼
    Error Metrics (AARD-P, AAD-y) → Output CSV

4. Key Issues & Code Smells

4.1 Massive Code Duplication

The most critical issue. The same logic is copy-pasted across 10+ cells with minor variations: - hard_constraint_layer_v17 is defined 6 times (Cells 7, 9, 10, 11, 12, 13) - build_pinn_model_from_params is defined 5 times - get_desc / training_data is defined 4 times - Path configuration is repeated in 8 cells - The prediction loop is written 7 times

This makes the notebook extremely difficult to maintain and debug. A single bug fix must be applied in multiple places.

4.2 Hardcoded Absolute Windows Paths

All paths use r"C:\Users\Hung\Desktop\2025_COSMO". This makes the notebook: - Impossible to run on another machine without editing - Impossible to run on macOS/Linux - Fragile to directory structure changes

4.3 Magic Numbers Everywhere

Column indices like 103, 105, 106, 107, [-6], [-7] appear throughout without named constants. The 110-dimensional feature vector layout is undocumented and relies on implicit indexing conventions.

4.4 Fragile Pressure Unit Detection

The heuristic P > 500 → Pa, else kPa is used to auto-detect pressure units. This fails for any system with pressure between 500–1000 kPa (which is physically plausible) and creates silent data corruption.

4.5 Row-by-Row Prediction Loop

Most cells call model.predict() inside a Python for loop for each of the 17,600+ rows. This is extremely slow compared to batch prediction. Cell 13 attempted batch prediction but was added as an afterthought.

4.6 Bug in Cell 10/11 — gammaB Normalization

'max_gammaB': sc.data_min_[-1]  # BUG: should be sc.data_max_[-1]

This causes incorrect gamma_B denormalization.

4.7 Inconsistent Gamma Placeholder Values

  • Cells 1–11 use 0 as placeholder for gammaA, gammaB in the feature vector
  • Cell 12 uses 1.0 as placeholder
  • The correct value depends on what the scaler expects, and this inconsistency can cause different predictions

4.8 Repeated CSV File Reading in Loops

training_data() and name_to_cosmo() open and read CSV files on every call. For 17,600 iterations, this means ~35,200 file open/close operations for sigma-profiles alone.

4.9 No Progress Bar or Time Estimation

Uses print(f"已完成 {i+1}...") every 500 rows. A proper progress bar (e.g., tqdm) would be more informative.

4.10 Commented-Out Code Not Cleaned Up

Cells 2 and 4 are entirely commented-out blocks (hundreds of lines), left in the notebook as "historical reference" but adding noise and confusion.

5. Improvement Suggestions

5.1 Restructure into a Clean Python Module

Extract all reusable logic into a Python module (e.g., cosmo_vle_predictor.py):

# cosmo_vle_predictor.py

import os
import numpy as np
import pandas as pd
import joblib
import csv
import pickle
import tensorflow as tf
from tensorflow.keras.models import Model
from tensorflow.keras.layers import Input, Dense, Lambda, Concatenate
import tensorflow.keras.backend as K
from pathlib import Path
from dataclasses import dataclass

# --- Constants ---
SIGMA_PROFILE_DIM = 51
FEATURE_DIM = 110
MLP_INPUT_DIM = 104
PHYS_INPUT_DIM = 4

# Feature vector column indices
IDX_T = 102
IDX_XA = 103
IDX_YA = 104
IDX_LNP = 105
IDX_LNVPA = 106
IDX_LNVPB = 107
IDX_GAMMA_A = 108
IDX_GAMMA_B = 109


@dataclass
class ModelConfig:
    root: Path
    results_version: int
    version: int
    fold_idx: int

    @property
    def model_train_dir(self) -> Path:
        return self.root / "model_train" / f"results_15_V{self.results_version}_{self.version}_Validation"

    @property
    def pipeline_path(self) -> Path:
        return self.model_train_dir / "model" / "pipelineMLP_COSMO_merged.joblib"

    @property
    def weight_path(self) -> Path:
        return self.model_train_dir / "model" / "h5" / f"PINN_V{self.results_version}_Fold_{self.fold_idx}.weights.h5"

    @property
    def pkl_path(self) -> Path:
        return self.model_train_dir / "model" / "pkl" / f"PINN_V{self.results_version}_Fold_{self.fold_idx}.pkl"

5.2 Cache Sigma-Profiles in Memory

Load all sigma-profiles once at startup instead of reading files in every iteration:

class SigmaProfileCache:
    def __init__(self, folder_path: Path):
        self.folder_path = folder_path
        self._cache: dict[str, list[float]] = {}
        self._zero_profile = np.zeros(SIGMA_PROFILE_DIM).tolist()

    def get(self, compound_id) -> list[float]:
        key = str(compound_id)
        if key not in self._cache:
            file_path = self.folder_path / f"{key}-opt-b3lyp.csv"
            if file_path.exists():
                with open(file_path, "r", encoding="utf-8") as f:
                    self._cache[key] = [float(row[1]) for row in csv.reader(f)][:SIGMA_PROFILE_DIM]
            else:
                self._cache[key] = self._zero_profile
        return self._cache[key]

5.3 Use Batch Prediction Instead of Row-by-Row

Assemble all feature vectors first, then predict in one call:

def predict_batch(model, scaler, feature_matrix: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    scaled = np.asarray(scaler.transform(feature_matrix))
    X_f = scaled[:, :MLP_INPUT_DIM]
    X_p = scaled[:, [IDX_XA, IDX_LNVPA, IDX_LNVPB, IDX_T]]

    preds = model.predict([X_f, X_p], batch_size=256, verbose=0)

    y_cal = preds[:, 0] * (scaler.data_max_[-6] - scaler.data_min_[-6]) + scaler.data_min_[-6]
    ln_p_cal = preds[:, 1] * (scaler.data_max_[-5] - scaler.data_min_[-5]) + scaler.data_min_[-5]
    p_cal = np.exp(ln_p_cal)

    return y_cal, p_cal

5.4 Use Configuration Files Instead of Hardcoded Paths

Replace hardcoded paths with a YAML/JSON configuration file:

# config.yaml
root: ./2025_COSMO
results_version: 8
version: 7
fold_idx: 5
input_file: grouping/data/grouping_alldata_V2.csv
sigma_profile_dir: grouping/data/s-profile-area

5.5 Define the Physics Layer as a Proper Keras Layer

Replace the Lambda with a custom tf.keras.layers.Layer for better serialization and clarity:

class MarguleConstraintLayer(tf.keras.layers.Layer):
    def __init__(self, scaler_params: dict, **kwargs):
        super().__init__(**kwargs)
        self.p = scaler_params

    def call(self, inputs):
        mlp_out, phys_in = inputs
        A12, A21 = mlp_out[:, 0:1], mlp_out[:, 1:2]
        y_mlp, P_delta = mlp_out[:, 2:3], mlp_out[:, 3:4]

        xA = phys_in[:, 0:1] * (self.p["max_xA"] - self.p["min_xA"]) + self.p["min_xA"]
        xB = 1.0 - xA
        lnPsatA = phys_in[:, 1:2] * (self.p["max_VpA"] - self.p["min_VpA"]) + self.p["min_VpA"]
        lnPsatB = phys_in[:, 2:3] * (self.p["max_VpB"] - self.p["min_VpB"]) + self.p["min_VpB"]

        ln_gammaA = tf.square(xB) * (A12 + 2.0 * (A21 - A12) * xA)
        ln_gammaB = tf.square(xA) * (A21 + 2.0 * (A12 - A21) * xB)

        P_phys = xA * tf.exp(ln_gammaA + lnPsatA) + xB * tf.exp(ln_gammaB + lnPsatB)
        P_log = tf.math.log(tf.maximum(P_phys, 1e-7))
        P_norm = (P_log - self.p["min_p"]) / (self.p["max_p"] - self.p["min_p"] + 1e-8)

        return tf.concat([y_mlp, P_norm + P_delta, tf.exp(ln_gammaA), tf.exp(ln_gammaB)], axis=-1)

    def get_config(self):
        config = super().get_config()
        config["scaler_params"] = self.p
        return config

5.6 Explicit Pressure Unit Handling

Replace the fragile P > 500 heuristic with explicit metadata:

def convert_pressure(value: float, unit: str) -> float:
    """Convert pressure to kPa with explicit unit specification."""
    conversion = {"Pa": 1e-3, "kPa": 1.0, "bar": 100.0, "atm": 101.325, "mmHg": 0.133322}
    if unit not in conversion:
        raise ValueError(f"Unknown pressure unit: {unit}")
    return value * conversion[unit]

Alternatively, add a P_unit column to the input CSV.

5.7 Add Progress Bars with tqdm

from tqdm import tqdm

for i in tqdm(range(len(df_input)), desc="Predicting"):
    ...

5.8 Use Type Hints and Docstrings

All functions should have clear type hints and docstrings explaining: - What the function does - Expected input/output shapes and units - Physical meaning of parameters

5.9 Remove Dead Code

Delete all commented-out cells (Cells 2, 4) and consolidate the prediction logic into a single, well-tested cell or function.

5.10 Add Validation and Sanity Checks

def validate_feature_vector(vec: np.ndarray, scaler) -> None:
    """Check that feature values are within expected scaler bounds."""
    scaled = scaler.transform(vec.reshape(1, -1))
    out_of_range = (scaled < -0.5) | (scaled > 1.5)
    if out_of_range.any():
        bad_cols = np.where(out_of_range[0])[0]
        raise ValueError(f"Features out of scaler range at indices: {bad_cols}")

6. Suggested Refactored Notebook Structure

A clean version of this notebook should have 5–6 cells total:

Cell Purpose
1 Configuration & imports (load config from YAML, define constants)
2 Model architecture definition (single hard_constraint_layer, single build_model)
3 Data loading & feature assembly (batch, with caching)
4 Batch prediction & inverse normalization
5 Results output (CSV + summary statistics)
6 (Optional) Visualization / error analysis

Or Better Yet: Convert to a Python Script

For a production prediction pipeline, convert the notebook into a command-line script:

python predict_vle.py --config config.yaml --input data.csv --output results.csv

This provides: - Version control friendliness (no JSON notebook diffs) - Reproducibility via configuration files - Easy integration into automated workflows

7. Performance Considerations

Aspect Current Suggested
Prediction Row-by-row (17,600 calls to model.predict) Batch prediction (1 call, batch_size=256)
File I/O Read sigma-profiles on every iteration Cache in memory at startup
CSV writing Row-by-row append with mode='a' Collect results in list, write once at end
Model loading Rebuilt from scratch in every cell Load once, reuse
Expected speedup — 50–100x for prediction alone

8. Identified Bugs

# Location Bug Impact
1 Cell 10/11 p_args 'max_gammaB': sc.data_min_[-1] instead of sc.data_max_[-1] Incorrect gamma_B denormalization
2 Cell 12 Gamma placeholders set to 1.0 vs 0 in other cells Inconsistent feature scaling
3 Cell 1 name_to_cosmo Hardcoded range range(0, 1960) — will silently fail for datasets with >1960 compounds Missing compound IDs
4 Cell 3 No log transform on P_exp before scaling — differs from later cells Feature misalignment with training
5 Cell 12 p_args Missing 'min_gammaA', 'max_gammaA', 'min_gammaB', 'max_gammaB' keys hard_constraint_layer_v17 will use kwargs.get() which returns None, causing runtime error or silent incorrect output

9. Summary

This notebook is a working but highly experimental research codebase that has evolved through many iterations of trial-and-error. It successfully implements a PINN-based VLE prediction system combining COSMO molecular descriptors with Margules activity coefficient physics constraints. However, it suffers from severe code duplication (~10 near-identical copies of the core logic), hardcoded paths, magic numbers, and inconsistent handling of physical units and placeholder values.

The most impactful improvements would be: 1. Consolidate all definitions into a single module (eliminate ~90% of duplicated code) 2. Use batch prediction instead of row-by-row (~50-100x speedup) 3. Cache sigma-profiles in memory (~2x I/O speedup) 4. Fix the gamma normalization bug in the physics constraint layer 5. Replace hardcoded paths with relative paths or configuration files