{ "nbformat": 4, "nbformat_minor": 0, "metadata": { "colab": { "provenance": [] }, "kernelspec": { "name": "python3", "display_name": "Python 3" }, "language_info": { "name": "python" } }, "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "gLbOycylHye0", "outputId": "383d9be6-53ea-47d7-a9e1-876d6f6ad2ea" }, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "DATA AUDIT\n", "==============================================================================\n", "\n", "Rows as loaded: 94\n", "Sources: 3 -> ['Abdulwahab+et+al', 'Diego+et+al', 'Wang+et+al']\n", "SCM types: ['CPA', 'Control', 'FA10', 'FA30', 'GBS10', 'GBS30', 'NS1', 'NS3', 'RHA', 'SF10', 'SF5']\n", "\n", "-- Missing values --\n", " none\n", "\n", "-- Exact duplicate rows --\n", " 6 duplicate row(s) <-- investigate, these inflate R^2\n", "\n", "-- SCM_% vs dosage implied by SCM_Type label --\n", " consistent\n", "\n", "-- Control rows with non-zero sulphate concentration --\n", " 4 of 4 control rows <-- check whether these are unexposed controls\n", "\n", "-- Ranges --\n", " exposure n=94 min=0 max=500 mean=96.37 sd=142.2\n", " sulphate n=94 min=0.352 max=0.704 mean=0.447 sd=0.1282\n", " strength n=94 min=6.12 max=66 mean=36.39 sd=22.65\n", "\n", "-- Strength distribution (the bimodality problem) --\n", " below 30 MPa: 40 at/above 30 MPa: 54\n", " by source:\n", " Abdulwahab+et+al n=40 range 6.12-21.50 MPa\n", " Diego+et+al n=18 range 43.00-64.00 MPa\n", " Wang+et+al n=36 range 47.00-66.00 MPa\n", " If sources occupy separate strength bands, a high R^2 mostly\n", " reflects the model identifying the source, not the chemistry.\n", " The leave-one-source-out result is the honest test.\n", "\n", "\n", "Modelling on 94 rows, 12 predictors. Reference category: Control\n", "\n", "Multicollinearity (VIF > 10 is the usual flag):\n", " feature VIF\n", " exposure 1.27\n", " sulphate 6.72\n", " SCM_CPA 7.02\n", " SCM_FA10 1.91\n", " SCM_FA30 1.91\n", "SCM_GBS10 1.91\n", "SCM_GBS30 1.91\n", " SCM_NS1 1.91\n", " SCM_NS3 1.91\n", " SCM_RHA 5.85\n", " SCM_SF10 5.58\n", " SCM_SF5 1.91 \n", "\n", "Selected alpha over grid 1e-4..1e4: 0.07149\n", " NOTE: alpha this small means ridge behaves almost exactly\n", " like OLS. Say so in the paper rather than claiming ridge\n", " overcame multicollinearity.\n", "\n", "TABLE 2 -- model comparison\n", " Model R2_train R2_test MAE_test RMSE_test CV_R2_mean CV_R2_sd\n", "OLS (baseline) 0.9751 0.9681 3.030 3.841 0.9636 0.0128\n", " Ridge 0.9750 0.9676 3.040 3.867 0.9639 0.0123\n", " Lasso 0.9748 0.9662 3.133 3.948 0.9606 0.0086\n", " Random Forest 0.9913 0.9828 2.188 2.819 0.9704 0.0132 \n", "\n", "LEAVE-ONE-SOURCE-OUT (report this prominently)\n", " held_out_source n_test R2 MAE RMSE\n", "Abdulwahab+et+al 40 -157.9154 42.3141 42.4472\n", " Diego+et+al 18 -1.8350 7.9881 9.3279\n", " Wang+et+al 36 -14.5548 17.0399 19.3843\n", "\n", "TABLE 3 -- coefficients\n", " Variable Coef_standardised Coef_raw_units Type\n", " exposure -0.7352 -0.0050 continuous\n", " sulphate -2.9277 -23.0379 continuous\n", " SCM_CPA -41.9981 -41.9981 indicator\n", " SCM_FA10 2.2617 2.2617 indicator\n", " SCM_FA30 -2.6219 -2.6219 indicator\n", "SCM_GBS10 7.1005 7.1005 indicator\n", "SCM_GBS30 0.7657 0.7657 indicator\n", " SCM_NS1 2.8680 2.8680 indicator\n", " SCM_NS3 3.5640 3.5640 indicator\n", " SCM_RHA 10.2800 10.2800 indicator\n", " SCM_SF10 11.0715 11.0715 indicator\n", " SCM_SF5 9.3795 9.3795 indicator\n", "Intercept 50.4380 61.1340 intercept \n", "\n", "Equation in RAW units (directly applicable, reproduces Table 4):\n", "\n", "CS = 61.1340\n", " - 0.0050 x exposure\n", " - 23.0379 x sulphate\n", " - 41.9981 x SCM_CPA\n", " + 2.2617 x SCM_FA10\n", " - 2.6219 x SCM_FA30\n", " + 7.1005 x SCM_GBS10\n", " + 0.7657 x SCM_GBS30\n", " + 2.8680 x SCM_NS1\n", " + 3.5640 x SCM_NS3\n", " + 10.2800 x SCM_RHA\n", " + 11.0715 x SCM_SF10\n", " + 9.3795 x SCM_SF5\n", "\n", "Reference category: Control (coefficient = 0 by construction)\n", "Selected alpha: 0.0714943\n", "\n", "Standardisation applied to continuous predictors before fitting:\n", " exposure: mean=93.8000, sd=147.5027\n", " sulphate: mean=0.4440, sd=0.1271\n", "\n", "State this standardisation explicitly in Section 2.4.\n", "\n", "TABLE 4 -- CPA at 0.367 mol/L\n", " Exposure_days Actual_MPa n_obs Predicted_MPa In_training_range\n", " 7 6.86 4 10.65 yes\n", " 14 8.61 4 10.61 yes\n", " 28 11.69 4 10.54 yes\n", " 56 13.29 5 10.40 yes\n", " 90 9.81 5 10.23 yes\n", " 180 NaN 0 9.78 NO -- extrapolation \n", "\n", "Done. Everything written to /content/run_output\n", "Read audit_report.txt before you touch the manuscript.\n" ] } ], "source": [ "#!/usr/bin/env python3\n", "\"\"\"\n", "run_analysis.py\n", "================================================================================\n", "Frozen, single-execution analysis pipeline for:\n", "\n", " \"Predicting compressive strength of SCM-blended mortar under sulphate\n", " exposure using ridge regression\"\n", "\n", "Run this ONCE. Every table and figure in the manuscript comes from this one\n", "execution. Do not hand-edit any number afterwards.\n", "\n", " python run_analysis.py --data final_cleaned_sulfate_data.csv\n", "\n", "Outputs land in ./run_output/ :\n", " audit_report.txt data quality flags -- READ THIS FIRST\n", " table2_model_comparison.csv\n", " table3_coefficients.csv\n", " table4_cpa_predictions.csv\n", " loso_validation.csv leave-one-source-out results\n", " equation.txt raw-unit equation you can paste into Section 2.4\n", " figure3_actual_vs_pred.png\n", " figure4_cpa_timeseries.png\n", " figure5_strength_by_scm.png\n", " figure6_coefficients.png\n", " manifest.txt seed, versions, alpha, row counts -- for Methods\n", "================================================================================\n", "\"\"\"\n", "\n", "import argparse\n", "import json\n", "import platform\n", "import sys\n", "from datetime import datetime\n", "from pathlib import Path\n", "\n", "import numpy as np\n", "import pandas as pd\n", "import matplotlib\n", "matplotlib.use(\"Agg\")\n", "import matplotlib.pyplot as plt\n", "\n", "import sklearn\n", "from sklearn.compose import ColumnTransformer\n", "from sklearn.ensemble import RandomForestRegressor\n", "from sklearn.linear_model import LinearRegression, Ridge, RidgeCV, LassoCV\n", "from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score\n", "from sklearn.model_selection import (KFold, LeaveOneGroupOut, cross_val_score,\n", " train_test_split)\n", "from sklearn.pipeline import Pipeline\n", "from sklearn.preprocessing import StandardScaler\n", "\n", "# ==============================================================================\n", "# CONFIGURATION -- edit this block to match your CSV, then leave it alone\n", "# ==============================================================================\n", "\n", "SEED = 42\n", "TEST_SIZE = 0.20\n", "CV_FOLDS = 5\n", "ALPHA_GRID = np.logspace(-4, 4, 200) # wide grid: this is the scaling fix\n", "\n", "# Map your CSV's column headers to the names the script uses internally.\n", "# Left = internal name, right = whatever your CSV actually calls it.\n", "COLUMN_MAP = {\n", " \"source\": \"Source Title\",\n", " \"scm_type\": \"SCM_Type\",\n", " \"scm_pct\": \"SCM_%\",\n", " \"exposure\": \"Exposure_Days\",\n", " \"sulphate\": \"Sulfate_Conc_mol/L\",\n", " \"strength\": \"Comp_Strength_MPa\",\n", "}\n", "\n", "# Which SCM_Type value is the no-replacement control. All coefficients are\n", "# interpreted relative to this. Blank/NaN cells are treated as control.\n", "REFERENCE_SCM = \"Control\"\n", "\n", "# Continuous predictors -- these get standardised.\n", "CONTINUOUS = [\"exposure\", \"sulphate\"]\n", "\n", "# Table 4 scenario\n", "TABLE4_SCM = \"CPA\"\n", "TABLE4_SULPHATE = 0.367\n", "TABLE4_DAYS = [7, 14, 28, 56, 90]\n", "TABLE4_EXTRAPOLATE = [180] # reported separately as out-of-range\n", "\n", "# Figure 5 projection range. Keep the ceiling near your data max unless you\n", "# have the long-duration rows back in.\n", "FIG5_MAX_DAYS = 90\n", "FIG5_SULPHATE = 0.367\n", "\n", "OUTDIR = Path(\"run_output\")\n", "\n", "# Expected dosage implied by each SCM_Type label, used only for the audit.\n", "IMPLIED_DOSAGE = {\n", " \"SF5\": 5, \"SF10\": 10, \"NS1\": 1, \"NS3\": 3,\n", " \"FA10\": 10, \"FA30\": 30, \"GBS10\": 10, \"GBS30\": 30,\n", "}\n", "\n", "plt.rcParams.update({\n", " \"figure.dpi\": 200, \"savefig.dpi\": 300, \"savefig.bbox\": \"tight\",\n", " \"font.size\": 9, \"axes.grid\": True, \"grid.alpha\": 0.3,\n", "})\n", "\n", "\n", "# ==============================================================================\n", "# 1. LOAD AND NORMALISE\n", "# ==============================================================================\n", "\n", "def load(path):\n", " df = pd.read_csv(path) if str(path).lower().endswith(\".csv\") \\\n", " else pd.read_excel(path)\n", " df.columns = [str(c).strip() for c in df.columns]\n", "\n", " missing = [src for src in COLUMN_MAP.values() if src not in df.columns]\n", " if missing:\n", " sys.exit(\n", " f\"\\nERROR: these columns are not in your file: {missing}\\n\"\n", " f\"Your file has: {list(df.columns)}\\n\"\n", " f\"Fix COLUMN_MAP at the top of this script.\\n\"\n", " )\n", "\n", " out = pd.DataFrame({k: df[v] for k, v in COLUMN_MAP.items()})\n", " out[\"scm_type\"] = (out[\"scm_type\"].astype(str).str.strip()\n", " .replace({\"nan\": REFERENCE_SCM, \"\": REFERENCE_SCM,\n", " \"None\": REFERENCE_SCM}))\n", " for c in CONTINUOUS + [\"strength\"]:\n", " out[c] = pd.to_numeric(out[c], errors=\"coerce\")\n", " out[\"source\"] = out[\"source\"].astype(str).str.strip()\n", " return out, df\n", "\n", "\n", "# ==============================================================================\n", "# 2. DATA AUDIT -- catches the inconsistencies a reviewer would find\n", "# ==============================================================================\n", "\n", "def audit(df, raw):\n", " L = [\"DATA AUDIT\", \"=\" * 78, \"\"]\n", " L.append(f\"Rows as loaded: {len(df)}\")\n", " L.append(f\"Sources: {df['source'].nunique()} -> \"\n", " f\"{sorted(df['source'].unique())}\")\n", " L.append(f\"SCM types: {sorted(df['scm_type'].unique())}\")\n", " L.append(\"\")\n", "\n", " L.append(\"-- Missing values --\")\n", " nulls = df.isna().sum()\n", " L += [f\" {k}: {v}\" for k, v in nulls.items() if v] or [\" none\"]\n", " L.append(\"\")\n", "\n", " L.append(\"-- Exact duplicate rows --\")\n", " d = df.duplicated().sum()\n", " L.append(f\" {d} duplicate row(s)\"\n", " + (\" <-- investigate, these inflate R^2\" if d else \"\"))\n", " L.append(\"\")\n", "\n", " L.append(\"-- SCM_% vs dosage implied by SCM_Type label --\")\n", " bad = []\n", " for lbl, dose in IMPLIED_DOSAGE.items():\n", " sub = df[df[\"scm_type\"] == lbl]\n", " if len(sub) and not np.allclose(sub[\"scm_pct\"].dropna(), dose):\n", " vals = sorted(sub[\"scm_pct\"].dropna().unique().tolist())\n", " bad.append(f\" {lbl}: label implies {dose}%, CSV has {vals}\")\n", " L += bad or [\" consistent\"]\n", " if bad:\n", " L.append(\" <-- FIX THIS. SF5 rows recorded at 10% is a real error and\")\n", " L.append(\" it also means SCM_Type and SCM_% encode the same thing\")\n", " L.append(\" twice, which is a collinearity problem in itself.\")\n", " L.append(\"\")\n", "\n", " L.append(\"-- Control rows with non-zero sulphate concentration --\")\n", " ctrl = df[df[\"scm_type\"] == REFERENCE_SCM]\n", " odd = ctrl[ctrl[\"sulphate\"] > 0]\n", " L.append(f\" {len(odd)} of {len(ctrl)} control rows\"\n", " + (\" <-- check whether these are unexposed controls\" if len(odd)\n", " else \"\"))\n", " L.append(\"\")\n", "\n", " L.append(\"-- Ranges --\")\n", " for c in CONTINUOUS + [\"strength\"]:\n", " s = df[c].dropna()\n", " L.append(f\" {c:<10} n={len(s):<4} min={s.min():<10.4g} \"\n", " f\"max={s.max():<10.4g} mean={s.mean():<10.4g} sd={s.std():.4g}\")\n", " L.append(\"\")\n", "\n", " L.append(\"-- Strength distribution (the bimodality problem) --\")\n", " s = df[\"strength\"].dropna()\n", " lo, hi = (s < 30).sum(), (s >= 30).sum()\n", " L.append(f\" below 30 MPa: {lo} at/above 30 MPa: {hi}\")\n", " L.append(\" by source:\")\n", " for src, g in df.groupby(\"source\")[\"strength\"]:\n", " L.append(f\" {src:<28} n={len(g):<4} \"\n", " f\"range {g.min():.2f}-{g.max():.2f} MPa\")\n", " L.append(\" If sources occupy separate strength bands, a high R^2 mostly\")\n", " L.append(\" reflects the model identifying the source, not the chemistry.\")\n", " L.append(\" The leave-one-source-out result is the honest test.\")\n", " L.append(\"\")\n", " return \"\\n\".join(L)\n", "\n", "\n", "# ==============================================================================\n", "# 3. DESIGN MATRIX\n", "# ==============================================================================\n", "\n", "def build_design(df):\n", " df = df.dropna(subset=CONTINUOUS + [\"strength\"]).reset_index(drop=True)\n", "\n", " dummies = pd.get_dummies(df[\"scm_type\"], prefix=\"SCM\", dtype=float)\n", " ref_col = f\"SCM_{REFERENCE_SCM}\"\n", " if ref_col in dummies.columns:\n", " dummies = dummies.drop(columns=[ref_col])\n", " ref_used = REFERENCE_SCM\n", " else:\n", " ref_used = dummies.columns[0].replace(\"SCM_\", \"\")\n", " dummies = dummies.drop(columns=[dummies.columns[0]])\n", "\n", " X = pd.concat([df[CONTINUOUS], dummies], axis=1)\n", " y = df[\"strength\"].astype(float)\n", " return X, y, df[\"source\"], ref_used, list(dummies.columns)\n", "\n", "\n", "def vif_table(X):\n", " \"\"\"Multicollinearity check -- justifies ridge, or shows it isn't needed.\"\"\"\n", " rows = []\n", " Xv = X.astype(float).values\n", " for i, name in enumerate(X.columns):\n", " others = np.delete(Xv, i, axis=1)\n", " if np.ptp(Xv[:, i]) == 0:\n", " rows.append({\"feature\": name, \"VIF\": np.nan}); continue\n", " r2 = LinearRegression().fit(others, Xv[:, i]).score(others, Xv[:, i])\n", " rows.append({\"feature\": name,\n", " \"VIF\": np.inf if r2 >= 1 - 1e-12 else 1.0 / (1.0 - r2)})\n", " return pd.DataFrame(rows)\n", "\n", "\n", "def make_pipe(model, X):\n", " \"\"\"Standardise continuous predictors only; leave 0/1 dummies alone.\"\"\"\n", " cont = [c for c in CONTINUOUS if c in X.columns]\n", " dums = [c for c in X.columns if c not in cont]\n", " pre = ColumnTransformer(\n", " [(\"scale\", StandardScaler(), cont), (\"pass\", \"passthrough\", dums)],\n", " remainder=\"drop\")\n", " return Pipeline([(\"pre\", pre), (\"model\", model)]), cont + dums\n", "\n", "\n", "# ==============================================================================\n", "# 4. FIT, COMPARE, VALIDATE\n", "# ==============================================================================\n", "\n", "def metrics(y, p):\n", " return {\"R2\": r2_score(y, p),\n", " \"MAE\": mean_absolute_error(y, p),\n", " \"RMSE\": float(np.sqrt(mean_squared_error(y, p)))}\n", "\n", "\n", "def main():\n", " ap = argparse.ArgumentParser()\n", " ap.add_argument(\"--data\", required=True, help=\"path to your CSV or XLSX\")\n", " ap.add_argument(\"--outdir\", default=str(OUTDIR))\n", " # Simulate command-line arguments for Colab environment\n", " # Replace 'final_cleaned_sulfate_data.csv' with your actual data file name\n", " args = ap.parse_args(['--data', 'final_cleaned_sulfate_data_v2.csv']) # Use v2 data\n", "\n", " out = Path(args.outdir); out.mkdir(exist_ok=True)\n", " np.random.seed(SEED)\n", "\n", " df, raw = load(args.data)\n", " (out / \"audit_report.txt\").write_text(audit(df, raw))\n", " print(audit(df, raw))\n", "\n", " X, y, groups, ref_used, dummy_cols = build_design(df)\n", " print(f\"\\nModelling on {len(X)} rows, {X.shape[1]} predictors. \"\n", " f\"Reference category: {ref_used}\\n\")\n", "\n", " vif = vif_table(X)\n", " vif.to_csv(out / \"vif.csv\", index=False)\n", " print(\"Multicollinearity (VIF > 10 is the usual flag):\")\n", " print(vif.round(2).to_string(index=False), \"\\n\")\n", "\n", " X_tr, X_te, y_tr, y_te = train_test_split(\n", " X, y, test_size=TEST_SIZE, random_state=SEED)\n", "\n", " # ---- alpha selection over a wide grid ------------------------------------\n", " ridge_cv, feat_order = make_pipe(\n", " RidgeCV(alphas=ALPHA_GRID, cv=KFold(CV_FOLDS, shuffle=True,\n", " random_state=SEED)), X)\n", " ridge_cv.fit(X_tr, y_tr)\n", " alpha = float(ridge_cv.named_steps[\"model\"].alpha_)\n", " print(f\"Selected alpha over grid 1e-4..1e4: {alpha:.4g}\")\n", " if alpha < 0.1:\n", " print(\" NOTE: alpha this small means ridge behaves almost exactly\")\n", " print(\" like OLS. Say so in the paper rather than claiming ridge\")\n", " print(\" overcame multicollinearity.\\n\")\n", "\n", " # ---- Table 2: model comparison ------------------------------------------\n", " candidates = {\n", " \"OLS (baseline)\": LinearRegression(),\n", " \"Ridge\": Ridge(alpha=alpha, random_state=None),\n", " \"Lasso\": LassoCV(cv=CV_FOLDS, random_state=SEED,\n", " max_iter=100000),\n", " \"Random Forest\": RandomForestRegressor(n_estimators=500,\n", " random_state=SEED),\n", " }\n", " rows, fitted = [], {}\n", " cv = KFold(CV_FOLDS, shuffle=True, random_state=SEED)\n", " for name, est in candidates.items():\n", " pipe, _ = make_pipe(est, X)\n", " pipe.fit(X_tr, y_tr)\n", " fitted[name] = pipe\n", " tr = metrics(y_tr, pipe.predict(X_tr))\n", " te = metrics(y_te, pipe.predict(X_te))\n", " cvs = cross_val_score(pipe, X_tr, y_tr, cv=cv, scoring=\"r2\")\n", " rows.append({\n", " \"Model\": name,\n", " \"R2_train\": round(tr[\"R2\"], 4),\n", " \"R2_test\": round(te[\"R2\"], 4),\n", " \"MAE_test\": round(te[\"MAE\"], 3),\n", " \"RMSE_test\": round(te[\"RMSE\"], 3),\n", " \"CV_R2_mean\": round(float(cvs.mean()), 4),\n", " \"CV_R2_sd\": round(float(cvs.std()), 4),\n", " })\n", " table2 = pd.DataFrame(rows)\n", " table2.to_csv(out / \"table2_model_comparison.csv\", index=False)\n", " print(\"TABLE 2 -- model comparison\\n\" + table2.to_string(index=False), \"\\n\")\n", "\n", " ridge = fitted[\"Ridge\"]\n", "\n", " # ---- Leave-one-source-out: the validation that actually matters ----------\n", " loso_rows = []\n", " if groups.nunique() > 1:\n", " for tr_i, te_i in LeaveOneGroupOut().split(X, y, groups):\n", " held = groups.iloc[te_i].unique()[0]\n", " p, _ = make_pipe(Ridge(alpha=alpha), X)\n", " p.fit(X.iloc[tr_i], y.iloc[tr_i])\n", " m = metrics(y.iloc[te_i], p.predict(X.iloc[te_i]))\n", " loso_rows.append({\"held_out_source\": held, \"n_test\": len(te_i),\n", " **{k: round(v, 4) for k, v in m.items()}})\n", " loso = pd.DataFrame(loso_rows)\n", " loso.to_csv(out / \"loso_validation.csv\", index=False)\n", " print(\"LEAVE-ONE-SOURCE-OUT (report this prominently)\\n\"\n", " + (loso.to_string(index=False) if len(loso) else \" only one source\")\n", " + \"\\n\")\n", "\n", " # ---- Table 3 + raw-unit equation ----------------------------------------\n", " scaler = ridge.named_steps[\"pre\"].named_transformers_[\"scale\"]\n", " coef = ridge.named_steps[\"model\"].coef_\n", " intercept = float(ridge.named_steps[\"model\"].intercept_)\n", " cont = [c for c in CONTINUOUS if c in X.columns]\n", " n_cont = len(cont)\n", "\n", " raw_coef, raw_intercept = {}, intercept\n", " for i, name in enumerate(cont):\n", " mu, sd = scaler.mean_[i], scaler.scale_[i]\n", " raw_coef[name] = coef[i] / sd\n", " raw_intercept -= coef[i] * mu / sd\n", " for j, name in enumerate(feat_order[n_cont:]):\n", " raw_coef[name] = coef[n_cont + j]\n", "\n", " table3 = pd.DataFrame([\n", " {\"Variable\": n,\n", " \"Coef_standardised\": round(float(coef[i]), 4),\n", " \"Coef_raw_units\": round(float(raw_coef[n]), 4),\n", " \"Type\": \"continuous\" if n in cont else \"indicator\"}\n", " for i, n in enumerate(feat_order)\n", " ])\n", " table3.loc[len(table3)] = {\"Variable\": \"Intercept\",\n", " \"Coef_standardised\": round(intercept, 4),\n", " \"Coef_raw_units\": round(raw_intercept, 4),\n", " \"Type\": \"intercept\"}\n", " table3.to_csv(out / \"table3_coefficients.csv\", index=False)\n", " print(\"TABLE 3 -- coefficients\\n\" + table3.to_string(index=False), \"\\n\")\n", "\n", " eq = [f\"CS = {raw_intercept:.4f}\"] + [\n", " f\" {'+' if raw_coef[n] >= 0 else '-'} {abs(raw_coef[n]):.4f} x {n}\"\n", " for n in feat_order]\n", " eq_txt = (\n", " \"Equation in RAW units (directly applicable, reproduces Table 4):\\n\\n\"\n", " + \"\\n\".join(eq)\n", " + f\"\\n\\nReference category: {ref_used} (coefficient = 0 by construction)\"\n", " + f\"\\nSelected alpha: {alpha:.6g}\"\n", " + \"\\n\\nStandardisation applied to continuous predictors before fitting:\\n\"\n", " + \"\\n\".join(f\" {n}: mean={scaler.mean_[i]:.4f}, sd={scaler.scale_[i]:.4f}\"\n", " for i, n in enumerate(cont))\n", " + \"\\n\\nState this standardisation explicitly in Section 2.4.\\n\"\n", " )\n", " (out / \"equation.txt\").write_text(eq_txt)\n", " print(eq_txt)\n", "\n", " # ---- Table 4 ------------------------------------------------------------\n", " def scenario(scm, days, sulph, pct=None):\n", " r = {c: 0.0 for c in X.columns}\n", " r[\"exposure\"] = float(days)\n", " r[\"sulphate\"] = float(sulph)\n", " if pct is None:\n", " sub = df[df[\"scm_type\"] == scm][\"scm_pct\"]\n", " pct = float(sub.median()) if len(sub) else 0.0\n", " r[\"scm_pct\"] = pct\n", " col = f\"SCM_{scm}\"\n", " if col in X.columns:\n", " r[col] = 1.0\n", " return pd.DataFrame([r])[X.columns]\n", "\n", " t4 = []\n", " obs = df[(df[\"scm_type\"] == TABLE4_SCM) &\n", " np.isclose(df[\"sulphate\"], TABLE4_SULPHATE, atol=1e-3)]\n", " for d in TABLE4_DAYS:\n", " a = obs[obs[\"exposure\"] == d][\"strength\"]\n", " t4.append({\n", " \"Exposure_days\": d,\n", " \"Actual_MPa\": round(float(a.mean()), 2) if len(a) else np.nan,\n", " \"n_obs\": len(a),\n", " \"Predicted_MPa\": round(float(\n", " ridge.predict(scenario(TABLE4_SCM, d, TABLE4_SULPHATE))[0]), 2),\n", " \"In_training_range\": \"yes\",\n", " })\n", " for d in TABLE4_EXTRAPOLATE:\n", " t4.append({\n", " \"Exposure_days\": d, \"Actual_MPa\": np.nan, \"n_obs\": 0,\n", " \"Predicted_MPa\": round(float(\n", " ridge.predict(scenario(TABLE4_SCM, d, TABLE4_SULPHATE))[0]), 2),\n", " \"In_training_range\": \"NO -- extrapolation\",\n", " })\n", " table4 = pd.DataFrame(t4)\n", " table4.to_csv(out / \"table4_cpa_predictions.csv\", index=False)\n", " print(f\"TABLE 4 -- {TABLE4_SCM} at {TABLE4_SULPHATE} mol/L\\n\"\n", " + table4.to_string(index=False), \"\\n\")\n", "\n", " # ==========================================================================\n", " # FIGURES -- all from this same fitted model\n", " # ==========================================================================\n", " p_te = ridge.predict(X_te)\n", "\n", " # Figure 3\n", " fig, ax = plt.subplots(figsize=(5, 5))\n", " ax.scatter(y_te, p_te, s=45, alpha=.8, edgecolor=\"k\", linewidth=.5,\n", " zorder=3, label=\"Test set\")\n", " lim = [min(y.min(), p_te.min()) - 3, max(y.max(), p_te.max()) + 3]\n", " ax.plot(lim, lim, \"r--\", lw=1.2, label=\"Actual = predicted\")\n", " ax.set(xlim=lim, ylim=lim,\n", " xlabel=\"Actual compressive strength (MPa)\",\n", " ylabel=\"Predicted compressive strength (MPa)\",\n", " title=f\"Actual vs predicted (ridge, $\\\\alpha$={alpha:.3g})\")\n", " ax.legend(); fig.savefig(out / \"figure3_actual_vs_pred.png\"); plt.close(fig)\n", "\n", " # Figure 4\n", " fig, ax = plt.subplots(figsize=(6, 4))\n", " if len(obs):\n", " ax.scatter(obs[\"exposure\"], obs[\"strength\"], marker=\"x\", c=\"tab:blue\",\n", " s=60, label=\"Measured\", zorder=3)\n", " grid = np.linspace(min(TABLE4_DAYS), max(TABLE4_DAYS), 200)\n", " ax.plot(grid, [ridge.predict(scenario(TABLE4_SCM, d, TABLE4_SULPHATE))[0]\n", " for d in grid], \"r-\", lw=1.6, label=\"Ridge prediction\")\n", " if TABLE4_EXTRAPOLATE:\n", " eg = np.linspace(max(TABLE4_DAYS), max(TABLE4_EXTRAPOLATE), 100)\n", " ax.plot(eg, [ridge.predict(scenario(TABLE4_SCM, d, TABLE4_SULPHATE))[0]\n", " for d in eg], \"r:\", lw=1.6, label=\"Extrapolated\")\n", " ax.axvline(max(TABLE4_DAYS), color=\"grey\", ls=\"--\", lw=.8)\n", " ax.set(xlabel=\"Exposure (days)\", ylabel=\"Compressive strength (MPa)\",\n", " title=f\"Actual vs predicted, {TABLE4_SCM} \"\n", " f\"({TABLE4_SULPHATE} mol/L sulphate)\")\n", " ax.legend(); fig.savefig(out / \"figure4_cpa_timeseries.png\"); plt.close(fig)\n", "\n", " # Figure 5\n", " fig, ax = plt.subplots(figsize=(6.5, 4))\n", " days = np.linspace(df[\"exposure\"].min(), FIG5_MAX_DAYS, 100)\n", " for scm in sorted(df[\"scm_type\"].unique()):\n", " ax.plot(days, [ridge.predict(scenario(scm, d, FIG5_SULPHATE))[0]\n", " for d in days], lw=1.3, label=scm)\n", " ax.set(xlabel=\"Exposure (days)\",\n", " ylabel=\"Predicted compressive strength (MPa)\",\n", " title=f\"Predicted strength by SCM type ({FIG5_SULPHATE} mol/L)\")\n", " ax.legend(fontsize=7, ncol=2, loc=\"best\")\n", " fig.savefig(out / \"figure5_strength_by_scm.png\"); plt.close(fig)\n", "\n", " # Figure 6 -- standardised coefficients, so bars are comparable\n", " fig, ax = plt.subplots(figsize=(6, 4.5))\n", " cd = table3[table3[\"Type\"] != \"intercept\"].copy()\n", " cd = cd.sort_values(\"Coef_standardised\")\n", " ax.barh(cd[\"Variable\"], cd[\"Coef_standardised\"],\n", " color=[\"tab:red\" if v < 0 else \"tab:blue\"\n", " for v in cd[\"Coef_standardised\"]])\n", " ax.axvline(0, color=\"k\", lw=.8)\n", " ax.set(xlabel=\"Standardised coefficient (MPa per SD, or vs reference)\",\n", " title=\"Ridge coefficients\")\n", " fig.savefig(out / \"figure6_coefficients.png\"); plt.close(fig)\n", "\n", " # ---- manifest -----------------------------------------------------------\n", " (out / \"manifest.txt\").write_text(\"\\n\".join([\n", " \"FROZEN RUN MANIFEST\",\n", " \"=\" * 50,\n", " f\"Timestamp: {datetime.now().isoformat(timespec='seconds')}\",\n", " f\"Data file: {args.data}\",\n", " f\"Rows modelled: {len(X)}\",\n", " f\"Predictors: {X.shape[1]}\",\n", " f\"Reference cat.: {ref_used}\",\n", " f\"Random seed: {SEED}\",\n", " f\"Train/test: {1-TEST_SIZE:.0%}/{TEST_SIZE:.0%}\",\n", " f\"CV folds: {CV_FOLDS}\",\n", " f\"Alpha grid: 1e-4 to 1e4, {len(ALPHA_GRID)} points\",\n", " f\"Selected alpha: {alpha:.6g}\",\n", " \"\",\n", " f\"Python: {platform.python_version()}\",\n", " f\"numpy: {np.__version__}\",\n", " f\"pandas: {pd.__version__}\",\n", " f\"scikit-learn: {sklearn.__version__}\",\n", " f\"matplotlib: {matplotlib.__version__}\",\n", " \"\",\n", " \"Cite these versions in Section 2.4. Replace the 'Python 3.0' claim.\",\n", " ]))\n", "\n", " print(f\"Done. Everything written to {out.resolve()}\")\n", " print(\"Read audit_report.txt before you touch the manuscript.\")\n", "\n", "\n", "if __name__ == \"__main__\":\n", " main()" ] }, { "cell_type": "code", "source": [ "import shutil\n", "from pathlib import Path\n", "from google.colab import files\n", "\n", "output_dir = Path(\"run_output\")\n", "zip_filename = \"run_output.zip\"\n", "\n", "if output_dir.exists() and output_dir.is_dir():\n", " shutil.make_archive(str(output_dir), 'zip', output_dir)\n", " files.download(zip_filename)\n", " print(f\"'{zip_filename}' downloaded successfully.\")\n", "else:\n", " print(f\"Error: The directory '{output_dir}' does not exist. Please ensure the analysis script has been run successfully to generate outputs.\")" ], "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 35 }, "id": "DIOFDTmVJtPm", "outputId": "e043b950-76d1-43da-d15b-fbf813428ab4" }, "execution_count": null, "outputs": [ { "output_type": "display_data", "data": { "text/plain": [ "" ], "application/javascript": [ "\n", " async function download(id, filename, size) {\n", " if (!google.colab.kernel.accessAllowed) {\n", " return;\n", " }\n", " const div = document.createElement('div');\n", " const label = document.createElement('label');\n", " label.textContent = `Downloading \"${filename}\": `;\n", " div.appendChild(label);\n", " const progress = document.createElement('progress');\n", " progress.max = size;\n", " div.appendChild(progress);\n", " document.body.appendChild(div);\n", "\n", " const buffers = [];\n", " let downloaded = 0;\n", "\n", " const channel = await google.colab.kernel.comms.open(id);\n", " // Send a message to notify the kernel that we're ready.\n", " channel.send({})\n", "\n", " for await (const message of channel.messages) {\n", " // Send a message to notify the kernel that we're ready.\n", " channel.send({})\n", " if (message.buffers) {\n", " for (const buffer of message.buffers) {\n", " buffers.push(buffer);\n", " downloaded += buffer.byteLength;\n", " progress.value = downloaded;\n", " }\n", " }\n", " }\n", " const blob = new Blob(buffers, {type: 'application/binary'});\n", " const a = document.createElement('a');\n", " a.href = window.URL.createObjectURL(blob);\n", " a.download = filename;\n", " div.appendChild(a);\n", " a.click();\n", " div.remove();\n", " }\n", " " ] }, "metadata": {} }, { "output_type": "display_data", "data": { "text/plain": [ "" ], "application/javascript": [ "download(\"download_68f8ade2-758e-4d1b-8544-b7c415d9db27\", \"run_output.zip\", 396379)" ] }, "metadata": {} }, { "output_type": "stream", "name": "stdout", "text": [ "'run_output.zip' downloaded successfully.\n" ] } ] }, { "cell_type": "code", "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 35 }, "id": "61cfe588", "outputId": "d90afd45-055d-46e3-ed8d-c3a5bb49f2bf" }, "source": [ "from pathlib import Path\n", "from google.colab import files\n", "\n", "output_dir = Path(\"run_output\")\n", "summary_file_name = \"summary_report.txt\"\n", "\n", "# List of text files to combine\n", "text_files_to_combine = [\n", " output_dir / \"audit_report.txt\",\n", " output_dir / \"equation.txt\",\n", " output_dir / \"manifest.txt\",\n", "]\n", "\n", "combined_content = []\n", "for f_path in text_files_to_combine:\n", " if f_path.exists():\n", " combined_content.append(f\"--- Content of {f_path.name} ---\\n\")\n", " combined_content.append(f_path.read_text())\n", " combined_content.append(\"\\n\\n\")\n", " else:\n", " print(f\"Warning: {f_path.name} not found and will be skipped.\")\n", "\n", "if combined_content:\n", " with open(summary_file_name, \"w\") as f:\n", " f.write(\"\".join(combined_content))\n", " files.download(summary_file_name)\n", " print(f\"'{summary_file_name}' downloaded successfully, containing combined text reports.\")\n", "else:\n", " print(\"No text reports found to combine.\")" ], "execution_count": null, "outputs": [ { "output_type": "display_data", "data": { "text/plain": [ "" ], "application/javascript": [ "\n", " async function download(id, filename, size) {\n", " if (!google.colab.kernel.accessAllowed) {\n", " return;\n", " }\n", " const div = document.createElement('div');\n", " const label = document.createElement('label');\n", " label.textContent = `Downloading \"${filename}\": `;\n", " div.appendChild(label);\n", " const progress = document.createElement('progress');\n", " progress.max = size;\n", " div.appendChild(progress);\n", " document.body.appendChild(div);\n", "\n", " const buffers = [];\n", " let downloaded = 0;\n", "\n", " const channel = await google.colab.kernel.comms.open(id);\n", " // Send a message to notify the kernel that we're ready.\n", " channel.send({})\n", "\n", " for await (const message of channel.messages) {\n", " // Send a message to notify the kernel that we're ready.\n", " channel.send({})\n", " if (message.buffers) {\n", " for (const buffer of message.buffers) {\n", " buffers.push(buffer);\n", " downloaded += buffer.byteLength;\n", " progress.value = downloaded;\n", " }\n", " }\n", " }\n", " const blob = new Blob(buffers, {type: 'application/binary'});\n", " const a = document.createElement('a');\n", " a.href = window.URL.createObjectURL(blob);\n", " a.download = filename;\n", " div.appendChild(a);\n", " a.click();\n", " div.remove();\n", " }\n", " " ] }, "metadata": {} }, { "output_type": "display_data", "data": { "text/plain": [ "" ], "application/javascript": [ "download(\"download_97272ad6-c3be-4f29-ac98-95004c222b03\", \"summary_report.txt\", 3092)" ] }, "metadata": {} }, { "output_type": "stream", "name": "stdout", "text": [ "'summary_report.txt' downloaded successfully, containing combined text reports.\n" ] } ] }, { "cell_type": "code", "source": [ "#!/usr/bin/env python3\n", "\"\"\"\n", "clean_data.py\n", "================================================================================\n", "Cleans final_cleaned_sulfate_data.csv before the analysis rerun.\n", "\n", "Does three things:\n", " 1. Reports what is actually in the file (category counts, RHA count, SF rows)\n", " 2. Removes exact duplicate rows\n", " 3. Applies label fixes you configure below\n", "\n", "Run it TWICE:\n", "\n", " Pass 1 -- report only, changes nothing:\n", " python clean_data.py --data final_cleaned_sulfate_data.csv --report-only\n", "\n", " Read the output. Decide what the bare \"SF\" rows should become.\n", "\n", " Pass 2 -- apply fixes after you have set SF_MERGE below:\n", " python clean_data.py --data final_cleaned_sulfate_data.csv\n", "\n", " Writes final_cleaned_sulfate_data_v2.csv and cleaning_log.txt\n", "================================================================================\n", "\"\"\"\n", "\n", "import argparse\n", "import sys\n", "from pathlib import Path\n", "\n", "import numpy as np\n", "import pandas as pd\n", "\n", "# ==============================================================================\n", "# CONFIGURATION\n", "# ==============================================================================\n", "\n", "COL_SOURCE = \"Source Title\"\n", "COL_SCM_TYPE = \"SCM_Type\"\n", "COL_SCM_PCT = \"SCM_%\"\n", "COL_EXPOSURE = \"Exposure_Days\"\n", "COL_SULPHATE = \"Sulfate_Conc_mol/L\"\n", "COL_STRENGTH = \"Comp_Strength_MPa\"\n", "\n", "CONTROL_LABEL = \"Control\"\n", "\n", "# Dosage each label implies. SCM_% is rewritten from this, because the label\n", "# is the more specific record. Verify against your source papers first.\n", "IMPLIED_DOSAGE = {\n", " \"SF5\": 5, \"SF10\": 10,\n", " \"NS1\": 1, \"NS3\": 3,\n", " \"FA10\": 10, \"FA30\": 30,\n", " \"GBS10\": 10, \"GBS30\": 30,\n", " CONTROL_LABEL: 0,\n", "}\n", "\n", "# Labels with no dosage in the name. Leave SCM_% as recorded in the CSV.\n", "DOSAGE_FROM_CSV = [\"CPA\", \"RHA\"]\n", "\n", "# What the bare \"SF\" rows should become. Leave as None until pass 1 tells you\n", "# what those rows contain. Then set one of:\n", "# SF_MERGE = \"SF5\" merge all SF rows into SF5\n", "# SF_MERGE = \"SF10\" merge all SF rows into SF10\n", "# SF_MERGE = \"by_pct\" use each row's SCM_% (5 -> SF5, 10 -> SF10)\n", "# SF_MERGE = \"drop\" remove the SF rows entirely\n", "SF_MERGE = \"SF10\"\n", "\n", "# Categories with fewer than this many rows get flagged as too thin to\n", "# interpret. Their coefficients are noise.\n", "MIN_ROWS_PER_CATEGORY = 6\n", "\n", "# Rewrite sulphate concentration to 0 for unexposed control rows.\n", "# Set True only if your control rows are genuinely water-cured, not immersed.\n", "ZERO_SULPHATE_FOR_CONTROL = False\n", "\n", "\n", "# ==============================================================================\n", "\n", "def load(path):\n", " df = pd.read_csv(path) if str(path).lower().endswith(\".csv\") \\\n", " else pd.read_excel(path)\n", " df.columns = [str(c).strip() for c in df.columns]\n", " need = [COL_SOURCE, COL_SCM_TYPE, COL_SCM_PCT, COL_EXPOSURE,\n", " COL_SULPHATE, COL_STRENGTH]\n", " missing = [c for c in need if c not in df.columns]\n", " if missing:\n", " sys.exit(f\"ERROR: missing columns {missing}\\nFile has: {list(df.columns)}\")\n", " df[COL_SCM_TYPE] = (df[COL_SCM_TYPE].astype(str).str.strip()\n", " .replace({\"nan\": CONTROL_LABEL, \"\": CONTROL_LABEL,\n", " \"None\": CONTROL_LABEL}))\n", " for c in [COL_SCM_PCT, COL_EXPOSURE, COL_SULPHATE, COL_STRENGTH]:\n", " df[c] = pd.to_numeric(df[c], errors=\"coerce\")\n", " return df\n", "\n", "\n", "def report(df, log):\n", " def P(s=\"\"):\n", " print(s); log.append(s)\n", "\n", " P(\"=\" * 78)\n", " P(\"DATASET REPORT\")\n", " P(\"=\" * 78)\n", " P(f\"Rows: {len(df)}\")\n", " P(\"\")\n", "\n", " P(\"-- Rows per SCM category (this answers the RHA question) --\")\n", " P(f\"{'Category':<12}{'n':>5}{'sources':>9} {'strength range':<22}\"\n", " f\"{'SCM_% values':<16}{'exposure days'}\")\n", " P(\"-\" * 78)\n", " thin = []\n", " for scm, g in df.groupby(COL_SCM_TYPE):\n", " pcts = sorted(g[COL_SCM_PCT].dropna().unique().tolist())\n", " days = sorted(g[COL_EXPOSURE].dropna().unique().tolist())\n", " rng = f\"{g[COL_STRENGTH].min():.2f}-{g[COL_STRENGTH].max():.2f} MPa\"\n", " flag = \" <-- THIN\" if len(g) < MIN_ROWS_PER_CATEGORY else \"\"\n", " if len(g) < MIN_ROWS_PER_CATEGORY:\n", " thin.append((scm, len(g)))\n", " P(f\"{scm:<12}{len(g):>5}{g[COL_SOURCE].nunique():>9} \"\n", " f\"{rng:<22}{str(pcts):<16}{days}{flag}\")\n", " P(\"\")\n", " if thin:\n", " P(\" THIN CATEGORIES -- their coefficients are noise, not findings:\")\n", " for scm, n in thin:\n", " P(f\" {scm}: {n} rows\")\n", " P(\" Either merge them into a broader class (e.g. all agro-ashes\")\n", " P(\" together) or exclude them and say so in the Limitations.\")\n", " P(\"\")\n", "\n", " P(\"-- Bare 'SF' rows (decide what these become) --\")\n", " sf = df[df[COL_SCM_TYPE] == \"SF\"]\n", " if len(sf):\n", " P(f\" {len(sf)} rows:\")\n", " cols = [COL_SOURCE, COL_SCM_PCT, COL_EXPOSURE, COL_SULPHATE,\n", " COL_STRENGTH]\n", " for line in sf[cols].to_string(index=False).split(\"\\n\"):\n", " P(\" \" + line)\n", " P(\" Set SF_MERGE at the top of this script, then rerun.\")\n", " else:\n", " P(\" none present\")\n", " P(\"\")\n", "\n", " P(\"-- Exact duplicate rows --\")\n", " dup = df[df.duplicated(keep=False)]\n", " P(f\" {df.duplicated().sum()} removable duplicate(s), \"\n", " f\"{len(dup)} rows involved\")\n", " if len(dup):\n", " for line in dup.sort_values(list(df.columns)).to_string(\n", " index=False).split(\"\\n\")[:20]:\n", " P(\" \" + line)\n", " P(\"\")\n", "\n", " P(\"-- SCM_% disagreeing with the label --\")\n", " bad = 0\n", " for lbl, dose in IMPLIED_DOSAGE.items():\n", " g = df[df[COL_SCM_TYPE] == lbl]\n", " if len(g) and not np.allclose(g[COL_SCM_PCT].dropna(), dose):\n", " vals = sorted(g[COL_SCM_PCT].dropna().unique().tolist())\n", " P(f\" {lbl}: label says {dose}%, CSV has {vals} ({len(g)} rows)\")\n", " bad += 1\n", " P(\" all consistent\" if not bad else\n", " \" These will be rewritten from the label.\")\n", " P(\"\")\n", "\n", " P(\"-- Control rows --\")\n", " ctrl = df[df[COL_SCM_TYPE] == CONTROL_LABEL]\n", " P(f\" {len(ctrl)} control rows, \"\n", " f\"{(ctrl[COL_SULPHATE] > 0).sum()} with non-zero sulphate\")\n", " P(f\" ZERO_SULPHATE_FOR_CONTROL is currently {ZERO_SULPHATE_FOR_CONTROL}\")\n", " P(\"\")\n", "\n", " P(\"-- Sulphate concentration by source (the VIF=744 problem) --\")\n", " for src, g in df.groupby(COL_SOURCE):\n", " P(f\" {src:<24} {sorted(g[COL_SULPHATE].dropna().unique().tolist())}\")\n", " P(\" If each source used one fixed concentration, sulphate is a proxy for\")\n", " P(\" study identity and its coefficient cannot be read as a chemical\")\n", " P(\" effect. Say this in the Discussion regardless of what you drop.\")\n", " P(\"\")\n", "\n", " P(\"-- Exposure coverage --\")\n", " ct = pd.crosstab(df[COL_SOURCE], df[COL_EXPOSURE])\n", " for line in ct.to_string().split(\"\\n\"):\n", " P(\" \" + line)\n", " P(\"\")\n", " return log\n", "\n", "\n", "def clean(df, log):\n", " def P(s=\"\"):\n", " print(s); log.append(s)\n", "\n", " P(\"=\" * 78)\n", " P(\"APPLYING FIXES\")\n", " P(\"=\" * 78)\n", " n0 = len(df)\n", "\n", " # 1 -- duplicates\n", " df = df.drop_duplicates().reset_index(drop=True)\n", " P(f\"1. Duplicates removed: {n0 - len(df)} ({n0} -> {len(df)} rows)\")\n", "\n", " # 2 -- SF merge\n", " n_sf = (df[COL_SCM_TYPE] == \"SF\").sum()\n", " if n_sf:\n", " if SF_MERGE is None:\n", " sys.exit(f\"\\nSTOP: {n_sf} bare 'SF' rows and SF_MERGE is None.\\n\"\n", " f\"Run with --report-only, read the SF rows, then set\\n\"\n", " f\"SF_MERGE at the top of this script.\\n\")\n", " if SF_MERGE == \"drop\":\n", " df = df[df[COL_SCM_TYPE] != \"SF\"].reset_index(drop=True)\n", " P(f\"2. SF rows dropped: {n_sf}\")\n", " elif SF_MERGE == \"by_pct\":\n", " m = df[COL_SCM_TYPE] == \"SF\"\n", " df.loc[m, COL_SCM_TYPE] = df.loc[m, COL_SCM_PCT].map(\n", " lambda p: \"SF5\" if p == 5 else (\"SF10\" if p == 10 else \"SF\"))\n", " left = (df[COL_SCM_TYPE] == \"SF\").sum()\n", " P(f\"2. SF mapped by SCM_%; {n_sf - left} reassigned, {left} unmapped\")\n", " if left:\n", " P(\" <-- unmapped SF rows remain. Their SCM_% is neither 5 \"\n", " \"nor 10.\")\n", " else:\n", " df.loc[df[COL_SCM_TYPE] == \"SF\", COL_SCM_TYPE] = SF_MERGE\n", " P(f\"2. SF -> {SF_MERGE}: {n_sf} rows\")\n", " else:\n", " P(\"2. No bare SF rows\")\n", "\n", " # 3 -- dosage from label\n", " fixed = 0\n", " for lbl, dose in IMPLIED_DOSAGE.items():\n", " m = df[COL_SCM_TYPE] == lbl\n", " if m.any():\n", " n = int((df.loc[m, COL_SCM_PCT] != dose).sum())\n", " if n:\n", " df.loc[m, COL_SCM_PCT] = dose\n", " P(f\" {lbl}: {n} rows rewritten to {dose}%\")\n", " fixed += n\n", " P(f\"3. SCM_% rewritten from label: {fixed} rows\")\n", " for lbl in DOSAGE_FROM_CSV:\n", " if (df[COL_SCM_TYPE] == lbl).any():\n", " vals = sorted(df.loc[df[COL_SCM_TYPE] == lbl,\n", " COL_SCM_PCT].dropna().unique().tolist())\n", " P(f\" {lbl}: SCM_% left as recorded {vals}\")\n", "\n", " # 4 -- control sulphate\n", " if ZERO_SULPHATE_FOR_CONTROL:\n", " m = df[COL_SCM_TYPE] == CONTROL_LABEL\n", " n = int((df.loc[m, COL_SULPHATE] != 0).sum())\n", " df.loc[m, COL_SULPHATE] = 0.0\n", " P(f\"4. Control sulphate zeroed: {n} rows\")\n", " else:\n", " P(\"4. Control sulphate left unchanged\")\n", "\n", " # duplicates can reappear once labels are normalised\n", " n1 = len(df)\n", " df = df.drop_duplicates().reset_index(drop=True)\n", " if n1 != len(df):\n", " P(f\"5. Duplicates created by relabelling, removed: {n1 - len(df)}\")\n", "\n", " P(\"\")\n", " P(\"-- Final category counts --\")\n", " for scm, g in df.groupby(COL_SCM_TYPE):\n", " P(f\" {scm:<12}{len(g):>4} rows\"\n", " + (\" <-- THIN\" if len(g) < MIN_ROWS_PER_CATEGORY else \"\"))\n", " P(f\"\\nFinal: {len(df)} rows ({n0 - len(df)} removed overall)\")\n", " return df\n", "\n", "\n", "def main():\n", " ap = argparse.ArgumentParser()\n", " ap.add_argument(\"--data\", required=True)\n", " ap.add_argument(\"--out\", default=None)\n", " ap.add_argument(\"--report-only\", action=\"store_true\")\n", " # Simulate command-line arguments for Colab environment\n", " # Use 'final_cleaned_sulfate_data.csv' as the default data file\n", " args = ap.parse_args(['--data', 'final_cleaned_sulfate_data.csv'])\n", "\n", " df = load(args.data)\n", " log = []\n", " report(df, log)\n", "\n", " if args.report_only:\n", " Path(\"cleaning_report.txt\").write_text(\"\\n\".join(log))\n", " print(\"\\nReport only. Nothing written to the dataset.\")\n", " print(\"Set SF_MERGE, then rerun without --report-only.\")\n", " return\n", "\n", " df = clean(df, log)\n", " out = args.out or str(Path(args.data).with_name(\n", " Path(args.data).stem + \"_v2.csv\"))\n", " df.to_csv(out, index=False)\n", " Path(\"cleaning_log.txt\").write_text(\"\\n\".join(log))\n", " print(f\"\\nWritten: {out}\")\n", " print(\"Log: cleaning_log.txt\")\n", " print(\"\\nNext: edit run_analysis.py, change\")\n", " print(' CONTINUOUS = [\"scm_pct\", \"exposure\", \"sulphate\"]')\n", " print(\"to\")\n", " print(' CONTINUOUS = [\"exposure\", \"sulphate\"]')\n", " print(f\"then run: python run_analysis.py --data {out}\")\n", "\n", "\n", "if __name__ == \"__main__\":\n", " main()" ], "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "p0B-9ERWTTJx", "outputId": "0a3403d0-f86b-494b-f718-ac90e74a95a6" }, "execution_count": null, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "==============================================================================\n", "DATASET REPORT\n", "==============================================================================\n", "Rows: 96\n", "\n", "-- Rows per SCM category (this answers the RHA question) --\n", "Category n sources strength range SCM_% values exposure days\n", "------------------------------------------------------------------------------\n", "CPA 40 1 6.12-21.50 MPa [20] [7, 14, 28, 56, 90]\n", "Control 4 1 48.00-52.00 MPa [0] [7, 28, 180, 500] <-- THIN\n", "FA10 4 1 52.00-57.00 MPa [10] [7, 28, 180, 500] <-- THIN\n", "FA30 4 1 47.00-52.00 MPa [10] [7, 28, 180, 500] <-- THIN\n", "GBS10 4 1 55.00-60.00 MPa [10] [7, 28, 180, 500] <-- THIN\n", "GBS30 4 1 49.00-54.00 MPa [10] [7, 28, 180, 500] <-- THIN\n", "NS1 4 1 50.00-56.00 MPa [10] [7, 28, 180, 500] <-- THIN\n", "NS3 4 1 52.00-58.00 MPa [10] [7, 28, 180, 500] <-- THIN\n", "RHA 10 1 43.00-64.00 MPa [10] [0, 14, 42, 70, 140]\n", "SF 10 1 43.00-62.00 MPa [10] [0, 14, 42, 70, 140]\n", "SF10 4 1 58.00-66.00 MPa [10] [7, 28, 180, 500] <-- THIN\n", "SF5 4 1 55.00-63.00 MPa [10] [7, 28, 180, 500] <-- THIN\n", "\n", " THIN CATEGORIES -- their coefficients are noise, not findings:\n", " Control: 4 rows\n", " FA10: 4 rows\n", " FA30: 4 rows\n", " GBS10: 4 rows\n", " GBS30: 4 rows\n", " NS1: 4 rows\n", " NS3: 4 rows\n", " SF10: 4 rows\n", " SF5: 4 rows\n", " Either merge them into a broader class (e.g. all agro-ashes\n", " together) or exclude them and say so in the Limitations.\n", "\n", "-- Bare 'SF' rows (decide what these become) --\n", " 10 rows:\n", " Source Title SCM_% Exposure_Days Sulfate_Conc_mol/L Comp_Strength_MPa\n", " Diego+et+al 10 0 0.704 43.0\n", " Diego+et+al 10 14 0.704 54.0\n", " Diego+et+al 10 42 0.704 55.0\n", " Diego+et+al 10 70 0.704 55.0\n", " Diego+et+al 10 140 0.704 56.0\n", " Diego+et+al 10 0 0.704 43.0\n", " Diego+et+al 10 14 0.704 56.0\n", " Diego+et+al 10 42 0.704 59.0\n", " Diego+et+al 10 70 0.704 61.0\n", " Diego+et+al 10 140 0.704 62.0\n", " Set SF_MERGE at the top of this script, then rerun.\n", "\n", "-- Exact duplicate rows --\n", " 2 removable duplicate(s), 4 rows involved\n", " Source Title Cement_Type SCM_Type SCM_% Exposure_Days Sulfate_Exposure Sulfate_Conc_mol/L Comp_Strength_MPa SCM_Sulfate_Interaction Exposure_Period\n", " Diego+et+al CPV-ARI RHA 10 0 Sulfate 0.704 43.0 7.04 0-1m\n", " Diego+et+al CPV-ARI RHA 10 0 Sulfate 0.704 43.0 7.04 0-1m\n", " Diego+et+al CPV-ARI SF 10 0 Sulfate 0.704 43.0 7.04 0-1m\n", " Diego+et+al CPV-ARI SF 10 0 Sulfate 0.704 43.0 7.04 0-1m\n", "\n", "-- SCM_% disagreeing with the label --\n", " SF5: label says 5%, CSV has [10] (4 rows)\n", " NS1: label says 1%, CSV has [10] (4 rows)\n", " NS3: label says 3%, CSV has [10] (4 rows)\n", " FA30: label says 30%, CSV has [10] (4 rows)\n", " GBS30: label says 30%, CSV has [10] (4 rows)\n", " These will be rewritten from the label.\n", "\n", "-- Control rows --\n", " 4 control rows, 4 with non-zero sulphate\n", " ZERO_SULPHATE_FOR_CONTROL is currently False\n", "\n", "-- Sulphate concentration by source (the VIF=744 problem) --\n", " Abdulwahab+et+al [0.352, 0.367]\n", " Diego+et+al [0.704]\n", " Wang+et+al [0.415]\n", " If each source used one fixed concentration, sulphate is a proxy for\n", " study identity and its coefficient cannot be read as a chemical\n", " effect. Say this in the Discussion regardless of what you drop.\n", "\n", "-- Exposure coverage --\n", " Exposure_Days 0 7 14 28 42 56 70 90 140 180 500\n", " Source Title \n", " Abdulwahab+et+al 0 8 8 8 0 8 0 8 0 0 0\n", " Diego+et+al 4 0 4 0 4 0 4 0 4 0 0\n", " Wang+et+al 0 9 0 9 0 0 0 0 0 9 9\n", "\n", "==============================================================================\n", "APPLYING FIXES\n", "==============================================================================\n", "1. Duplicates removed: 2 (96 -> 94 rows)\n", "2. SF -> SF10: 9 rows\n", " SF5: 4 rows rewritten to 5%\n", " NS1: 4 rows rewritten to 1%\n", " NS3: 4 rows rewritten to 3%\n", " FA30: 4 rows rewritten to 30%\n", " GBS30: 4 rows rewritten to 30%\n", "3. SCM_% rewritten from label: 20 rows\n", " CPA: SCM_% left as recorded [20]\n", " RHA: SCM_% left as recorded [10]\n", "4. Control sulphate left unchanged\n", "\n", "-- Final category counts --\n", " CPA 40 rows\n", " Control 4 rows <-- THIN\n", " FA10 4 rows <-- THIN\n", " FA30 4 rows <-- THIN\n", " GBS10 4 rows <-- THIN\n", " GBS30 4 rows <-- THIN\n", " NS1 4 rows <-- THIN\n", " NS3 4 rows <-- THIN\n", " RHA 9 rows\n", " SF10 13 rows\n", " SF5 4 rows <-- THIN\n", "\n", "Final: 94 rows (2 removed overall)\n", "\n", "Written: final_cleaned_sulfate_data_v2.csv\n", "Log: cleaning_log.txt\n", "\n", "Next: edit run_analysis.py, change\n", " CONTINUOUS = [\"scm_pct\", \"exposure\", \"sulphate\"]\n", "to\n", " CONTINUOUS = [\"exposure\", \"sulphate\"]\n", "then run: python run_analysis.py --data final_cleaned_sulfate_data_v2.csv\n" ] } ] }, { "cell_type": "code", "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 35 }, "id": "33727ffa", "outputId": "cbb31b0a-6e7d-40f3-9747-3404841aaaa6" }, "source": [ "import os\n", "import zipfile\n", "from google.colab import files\n", "\n", "def zip_directory(path, zip_file_handle):\n", " for root, dirs, files_in_dir in os.walk(path):\n", " for file in files_in_dir:\n", " zip_file_handle.write(os.path.join(root, file), os.path.relpath(os.path.join(root, file), os.path.join(path, '..')))\n", "\n", "zip_filename = 'all_files.zip'\n", "with zipfile.ZipFile(zip_filename, 'w', zipfile.ZIP_DEFLATED) as zipf:\n", " zip_directory('.', zipf)\n", "\n", "files.download(zip_filename)\n", "print(f\"All files zipped into '{zip_filename}' and ready for download.\")" ], "execution_count": null, "outputs": [ { "output_type": "display_data", "data": { "text/plain": [ "" ], "application/javascript": [ "\n", " async function download(id, filename, size) {\n", " if (!google.colab.kernel.accessAllowed) {\n", " return;\n", " }\n", " const div = document.createElement('div');\n", " const label = document.createElement('label');\n", " label.textContent = `Downloading \"${filename}\": `;\n", " div.appendChild(label);\n", " const progress = document.createElement('progress');\n", " progress.max = size;\n", " div.appendChild(progress);\n", " document.body.appendChild(div);\n", "\n", " const buffers = [];\n", " let downloaded = 0;\n", "\n", " const channel = await google.colab.kernel.comms.open(id);\n", " // Send a message to notify the kernel that we're ready.\n", " channel.send({})\n", "\n", " for await (const message of channel.messages) {\n", " // Send a message to notify the kernel that we're ready.\n", " channel.send({})\n", " if (message.buffers) {\n", " for (const buffer of message.buffers) {\n", " buffers.push(buffer);\n", " downloaded += buffer.byteLength;\n", " progress.value = downloaded;\n", " }\n", " }\n", " }\n", " const blob = new Blob(buffers, {type: 'application/binary'});\n", " const a = document.createElement('a');\n", " a.href = window.URL.createObjectURL(blob);\n", " a.download = filename;\n", " div.appendChild(a);\n", " a.click();\n", " div.remove();\n", " }\n", " " ] }, "metadata": {} }, { "output_type": "display_data", "data": { "text/plain": [ "" ], "application/javascript": [ "download(\"download_cc343f87-fe74-4b42-9d3c-d4b82a0144f5\", \"all_files.zip\", 7870969)" ] }, "metadata": {} }, { "output_type": "stream", "name": "stdout", "text": [ "All files zipped into 'all_files.zip' and ready for download.\n" ] } ] }, { "cell_type": "code", "source": [ "#!/usr/bin/env python3\n", "\"\"\"\n", "make_figures.py\n", "================================================================================\n", "Generates all six manuscript figures for\n", "\n", " \"Between-study variance inflates apparent accuracy in literature-compiled\n", " machine learning models of sulphate resistance in SCM-blended cement\"\n", "\n", " python make_figures.py --data final_cleaned_sulfate_data_v2.csv\n", "\n", "Outputs to ./figures/ :\n", " Fig1.tif dataset composition by source and SCM category\n", " Fig2.tif predicted vs measured, random test partition\n", " Fig3.tif standardised ridge coefficients\n", " Fig4.tif R2 under random partitioning vs leave-one-source-out\n", " Fig5.tif strength distribution by source\n", " Fig6.tif CPA measured vs predicted over exposure\n", "\n", "Journal requirements applied throughout:\n", " - No titles inside the figure. Titles are captions in the manuscript.\n", " - Arial lettering, 8-12 pt, minimal size variation within a figure.\n", " - TIFF, 300 DPI (halftone/combination). Set --eps for vector output instead.\n", " - Line widths at or above 0.3 pt.\n", " - Patterns and markers as well as colour, for colour-blind readers.\n", " - Figure parts labelled with lowercase letters where a figure has parts.\n", "\n", "Everything is fitted with the same seed, alpha grid and preparation as\n", "run_analysis.py, so the numbers here match Tables 3-7.\n", "================================================================================\n", "\"\"\"\n", "\n", "import argparse\n", "from pathlib import Path\n", "\n", "import numpy as np\n", "import pandas as pd\n", "import matplotlib\n", "matplotlib.use(\"Agg\")\n", "import matplotlib.pyplot as plt\n", "from matplotlib.lines import Line2D\n", "\n", "from sklearn.compose import ColumnTransformer\n", "from sklearn.linear_model import Ridge, RidgeCV\n", "from sklearn.model_selection import KFold, LeaveOneGroupOut, train_test_split\n", "from sklearn.metrics import r2_score\n", "from sklearn.pipeline import Pipeline\n", "from sklearn.preprocessing import StandardScaler\n", "\n", "# ==============================================================================\n", "# CONFIGURATION -- must match run_analysis.py\n", "# ==============================================================================\n", "\n", "SEED = 42\n", "TEST_SIZE = 0.20\n", "CV_FOLDS = 5\n", "ALPHA_GRID = np.logspace(-4, 4, 200)\n", "\n", "COLUMN_MAP = {\n", " \"source\": \"Source Title\",\n", " \"scm_type\": \"SCM_Type\",\n", " \"scm_pct\": \"SCM_%\",\n", " \"exposure\": \"Exposure_Days\",\n", " \"sulphate\": \"Sulfate_Conc_mol/L\",\n", " \"strength\": \"Comp_Strength_MPa\",\n", "}\n", "REFERENCE_SCM = \"Control\"\n", "CONTINUOUS = [\"exposure\", \"sulphate\"] # scm_pct dropped (see Section 2.3)\n", "\n", "# Display names for the three sources. Replace once citations are resolved.\n", "SOURCE_LABELS = {\n", " \"Abdulwahab+et+al\": \"Source A\",\n", " \"Diego+et+al\": \"Source B\",\n", " \"Wang+et+al\": \"Source C\",\n", "}\n", "SOURCE_ORDER = [\"Source A\", \"Source B\", \"Source C\"]\n", "\n", "# Figure 6 scenario\n", "F6_SCM = \"CPA\"\n", "F6_SULPHATE = 0.367\n", "F6_OBS_MAX = 90 # end of observed range\n", "F6_EXTRAP_MAX = 180 # dotted beyond this point\n", "\n", "OUTDIR = Path(\"figures\")\n", "\n", "# ==============================================================================\n", "# STYLE\n", "# ==============================================================================\n", "\n", "def set_style(fmt):\n", " plt.rcParams.update({\n", " \"font.family\": \"sans-serif\",\n", " \"font.sans-serif\": [\"Arial\", \"Helvetica\", \"DejaVu Sans\"],\n", " \"font.size\": 9, # body lettering\n", " \"axes.titlesize\": 9,\n", " \"axes.labelsize\": 9,\n", " \"xtick.labelsize\": 8,\n", " \"ytick.labelsize\": 8,\n", " \"legend.fontsize\": 8,\n", " \"axes.linewidth\": 0.6, # ~0.4 pt, above the 0.3 pt floor\n", " \"grid.linewidth\": 0.5,\n", " \"lines.linewidth\": 1.2,\n", " \"xtick.major.width\": 0.6,\n", " \"ytick.major.width\": 0.6,\n", " \"axes.grid\": True,\n", " \"grid.alpha\": 0.25,\n", " \"grid.linestyle\": \"-\",\n", " \"axes.spines.top\": False,\n", " \"axes.spines.right\": False,\n", " \"figure.dpi\": 300,\n", " \"savefig.dpi\": 300,\n", " \"savefig.bbox\": \"tight\",\n", " \"savefig.pad_inches\": 0.02,\n", " })\n", "\n", "# Greyscale-safe palette. Distinguishable in colour and when printed mono.\n", "C_DARK = \"#1f3b5c\"\n", "C_MID = \"#5b8db8\"\n", "C_LIGHT = \"#a8c4dc\"\n", "C_NEG = \"#8c2d2d\"\n", "C_POS = \"#2f5f8f\"\n", "C_LINE = \"#b03030\"\n", "SOURCE_FILL = {\"Source A\": C_DARK, \"Source B\": C_MID, \"Source C\": C_LIGHT}\n", "SOURCE_HATCH = {\"Source A\": \"///\", \"Source B\": \"...\", \"Source C\": \"\"}\n", "SOURCE_MARK = {\"Source A\": \"o\", \"Source B\": \"s\", \"Source C\": \"^\"}\n", "\n", "# Single-column width 78 mm, double-column 119 mm, per journal guidance.\n", "W1 = 78 / 25.4\n", "W2 = 119 / 25.4\n", "\n", "\n", "def save(fig, name, fmt, outdir):\n", " outdir.mkdir(exist_ok=True)\n", " path = outdir / f\"{name}.{fmt}\"\n", " if fmt == \"tif\":\n", " fig.savefig(path, format=\"tiff\", pil_kwargs={\"compression\": \"tiff_lzw\"})\n", " else:\n", " fig.savefig(path, format=fmt)\n", " plt.close(fig)\n", " print(f\" wrote {path}\")\n", "\n", "\n", "# ==============================================================================\n", "# DATA AND MODEL\n", "# ==============================================================================\n", "\n", "def load(path):\n", " df = pd.read_csv(path) if str(path).lower().endswith(\".csv\") else pd.read_excel(path)\n", " df.columns = [str(c).strip() for c in df.columns]\n", " out = pd.DataFrame({k: df[v] for k, v in COLUMN_MAP.items()})\n", " out[\"scm_type\"] = (out[\"scm_type\"].astype(str).str.strip()\n", " .replace({\"nan\": REFERENCE_SCM, \"\": REFERENCE_SCM,\n", " \"None\": REFERENCE_SCM}))\n", " for c in [\"scm_pct\", \"exposure\", \"sulphate\", \"strength\"]:\n", " out[c] = pd.to_numeric(out[c], errors=\"coerce\")\n", " out[\"source\"] = (out[\"source\"].astype(str).str.strip()\n", " .map(lambda s: SOURCE_LABELS.get(s, s)))\n", " return out.dropna(subset=CONTINUOUS + [\"strength\"]).reset_index(drop=True)\n", "\n", "\n", "def build_design(df):\n", " dummies = pd.get_dummies(df[\"scm_type\"], prefix=\"SCM\", dtype=float)\n", " ref = f\"SCM_{REFERENCE_SCM}\"\n", " dummies = dummies.drop(columns=[ref]) if ref in dummies else dummies.iloc[:, 1:]\n", " X = pd.concat([df[CONTINUOUS], dummies], axis=1)\n", " return X, df[\"strength\"].astype(float)\n", "\n", "\n", "def make_pipe(model, X):\n", " cont = [c for c in CONTINUOUS if c in X.columns]\n", " dums = [c for c in X.columns if c not in cont]\n", " pre = ColumnTransformer([(\"scale\", StandardScaler(), cont),\n", " (\"pass\", \"passthrough\", dums)])\n", " return Pipeline([(\"pre\", pre), (\"model\", model)]), cont + dums\n", "\n", "\n", "def fit_everything(df):\n", " X, y = build_design(df)\n", " groups = df[\"source\"]\n", " X_tr, X_te, y_tr, y_te = train_test_split(\n", " X, y, test_size=TEST_SIZE, random_state=SEED)\n", "\n", " cv = KFold(CV_FOLDS, shuffle=True, random_state=SEED)\n", " rcv, feat_order = make_pipe(RidgeCV(alphas=ALPHA_GRID, cv=cv), X)\n", " rcv.fit(X_tr, y_tr)\n", " alpha = float(rcv.named_steps[\"model\"].alpha_)\n", "\n", " ridge, _ = make_pipe(Ridge(alpha=alpha), X)\n", " ridge.fit(X_tr, y_tr)\n", "\n", " loso = []\n", " for tr_i, te_i in LeaveOneGroupOut().split(X, y, groups):\n", " held = groups.iloc[te_i].unique()[0]\n", " p, _ = make_pipe(Ridge(alpha=alpha), X)\n", " p.fit(X.iloc[tr_i], y.iloc[tr_i])\n", " loso.append({\"source\": held,\n", " \"r2\": r2_score(y.iloc[te_i], p.predict(X.iloc[te_i]))})\n", "\n", " return dict(X=X, y=y, X_te=X_te, y_te=y_te, ridge=ridge, alpha=alpha,\n", " feat_order=feat_order,\n", " r2_random=r2_score(y_te, ridge.predict(X_te)),\n", " loso=pd.DataFrame(loso))\n", "\n", "\n", "def scenario(X, df, scm, days, sulph):\n", " r = {c: 0.0 for c in X.columns}\n", " r[\"exposure\"] = float(days)\n", " r[\"sulphate\"] = float(sulph)\n", " col = f\"SCM_{scm}\"\n", " if col in X.columns:\n", " r[col] = 1.0\n", " return pd.DataFrame([r])[X.columns]\n", "\n", "\n", "# ==============================================================================\n", "# FIGURES\n", "# ==============================================================================\n", "\n", "def fig1(df, fmt, outdir):\n", " \"\"\"Dataset composition: observations per source, grouped by SCM category.\"\"\"\n", " ct = pd.crosstab(df[\"scm_type\"], df[\"source\"])\n", " ct = ct.reindex(columns=[s for s in SOURCE_ORDER if s in ct.columns],\n", " fill_value=0)\n", " ct = ct.loc[ct.sum(axis=1).sort_values(ascending=True).index]\n", "\n", " fig, ax = plt.subplots(figsize=(W2, 3.4))\n", " left = np.zeros(len(ct))\n", " for src in ct.columns:\n", " vals = ct[src].values\n", " ax.barh(ct.index, vals, left=left, height=0.68,\n", " color=SOURCE_FILL[src], hatch=SOURCE_HATCH[src],\n", " edgecolor=\"white\", linewidth=0.6, label=src)\n", " left += vals\n", "\n", " for i, total in enumerate(ct.sum(axis=1).values):\n", " ax.text(total + 0.7, i, str(int(total)), va=\"center\", fontsize=8)\n", "\n", " ax.set_xlabel(\"Number of observations\")\n", " ax.set_ylabel(\"SCM category\")\n", " ax.set_xlim(0, ct.sum(axis=1).max() * 1.13)\n", " ax.grid(axis=\"y\", visible=False)\n", " ax.legend(frameon=False, loc=\"lower right\")\n", " save(fig, \"Fig1\", fmt, outdir)\n", "\n", "\n", "def fig2(m, fmt, outdir):\n", " \"\"\"Predicted vs measured, random 20% test partition, with 1:1 line.\"\"\"\n", " pred = m[\"ridge\"].predict(m[\"X_te\"])\n", " obs = m[\"y_te\"].values\n", "\n", " fig, ax = plt.subplots(figsize=(W1 * 1.35, W1 * 1.35))\n", " lo = min(obs.min(), pred.min()) - 4\n", " hi = max(obs.max(), pred.max()) + 4\n", " ax.plot([lo, hi], [lo, hi], ls=\"--\", color=C_LINE, lw=1.0, zorder=1)\n", " ax.scatter(obs, pred, s=34, facecolor=C_MID, edgecolor=C_DARK,\n", " linewidth=0.6, zorder=3)\n", " ax.set_xlim(lo, hi); ax.set_ylim(lo, hi)\n", " ax.set_aspect(\"equal\")\n", " ax.set_xlabel(\"Measured compressive strength (MPa)\")\n", " ax.set_ylabel(\"Predicted compressive strength (MPa)\")\n", " ax.legend(handles=[\n", " Line2D([], [], marker=\"o\", ls=\"none\", markerfacecolor=C_MID,\n", " markeredgecolor=C_DARK, markersize=5, label=\"Test observation\"),\n", " Line2D([], [], ls=\"--\", color=C_LINE, label=\"Perfect agreement\")],\n", " frameon=False, loc=\"upper left\")\n", " save(fig, \"Fig2\", fmt, outdir)\n", "\n", "\n", "def fig3(m, fmt, outdir):\n", " \"\"\"Standardised ridge coefficients, sorted, zero line marked.\"\"\"\n", " coef = m[\"ridge\"].named_steps[\"model\"].coef_\n", " names = [n.replace(\"SCM_\", \"\").replace(\"exposure\", \"Exposure duration\")\n", " .replace(\"sulphate\", \"Sulphate concentration\")\n", " for n in m[\"feat_order\"]]\n", " s = pd.Series(coef, index=names).sort_values()\n", "\n", " fig, ax = plt.subplots(figsize=(W2, 3.6))\n", " colors = [C_NEG if v < 0 else C_POS for v in s.values]\n", " hatches = [\"///\" if v < 0 else \"\" for v in s.values]\n", " bars = ax.barh(s.index, s.values, height=0.68, color=colors,\n", " edgecolor=\"black\", linewidth=0.5)\n", " for b, h in zip(bars, hatches):\n", " b.set_hatch(h)\n", " ax.axvline(0, color=\"black\", lw=0.8)\n", "\n", " pad = max(abs(s.min()), abs(s.max())) * 0.03\n", " for i, v in enumerate(s.values):\n", " ax.text(v + (pad if v >= 0 else -pad), i, f\"{v:.2f}\",\n", " va=\"center\", ha=\"left\" if v >= 0 else \"right\", fontsize=7.5)\n", "\n", " ax.set_xlabel(\"Standardised coefficient (MPa)\")\n", " ax.set_xlim(s.min() * 1.22, s.max() * 1.28)\n", " ax.grid(axis=\"y\", visible=False)\n", " save(fig, \"Fig3\", fmt, outdir)\n", "\n", "\n", "def fig4(m, fmt, outdir):\n", " \"\"\"R2 under random partitioning vs leave-one-source-out.\n", "\n", " Two panels sharing a category axis. A single linear axis cannot show\n", " 0.968 and -157.92 together; a symlog axis compresses the difference the\n", " figure exists to convey. Panel (a) shows the near-zero region, panel (b)\n", " the full range, with the panel (a) window marked on panel (b).\n", " \"\"\"\n", " loso = m[\"loso\"].set_index(\"source\").reindex(SOURCE_ORDER)\n", " labels = [\"Random\\nsplit\"] + [f\"LOSO\\n{s}\" for s in SOURCE_ORDER]\n", " vals = [m[\"r2_random\"]] + list(loso[\"r2\"].values)\n", " colors = [C_POS] + [C_NEG] * 3\n", " hatches = [\"\"] + [\"///\"] * 3\n", "\n", " fig, (ax1, ax2) = plt.subplots(\n", " 2, 1, figsize=(W2, 4.6), sharex=True,\n", " gridspec_kw={\"height_ratios\": [1, 1.5], \"hspace\": 0.18})\n", "\n", " for ax in (ax1, ax2):\n", " bars = ax.bar(labels, vals, width=0.6, color=colors,\n", " edgecolor=\"black\", linewidth=0.5)\n", " for b, h in zip(bars, hatches):\n", " b.set_hatch(h)\n", " ax.axhline(0, color=\"black\", lw=0.8)\n", " ax.grid(axis=\"x\", visible=False)\n", "\n", " ax1.set_ylim(-20, 2)\n", " ax2.set_ylim(-172, 12)\n", " ax1.text(0.012, 0.90, \"a\", transform=ax1.transAxes, fontweight=\"bold\", fontsize=10)\n", " ax2.text(0.012, 0.93, \"b\", transform=ax2.transAxes, fontweight=\"bold\", fontsize=10)\n", "\n", " for i, v in enumerate(vals):\n", " if v > -20:\n", " ax1.text(i, v + (0.8 if v >= 0 else -1.9), f\"{v:.3f}\" if v > 0 else f\"{v:.2f}\",\n", " ha=\"center\", fontsize=8)\n", " ax2.text(i, v - 8 if v < 0 else v + 4, f\"{v:.2f}\" if v < 0 else f\"{v:.3f}\",\n", " ha=\"center\", fontsize=8)\n", "\n", " ax2.axhspan(-20, 2, color=\"black\", alpha=0.06, zorder=0)\n", " ax2.annotate(\"range shown in (a)\", xy=(3.42, -10), xytext=(3.42, -46),\n", " fontsize=7, ha=\"center\", color=\"0.35\",\n", " arrowprops=dict(arrowstyle=\"->\", lw=0.6, color=\"0.5\"))\n", "\n", " ax2.set_xlabel(\"Validation scheme\", labelpad=6)\n", " fig.supylabel(\"Coefficient of determination, $R^2$\", fontsize=9, x=0.005)\n", " save(fig, \"Fig4\", fmt, outdir)\n", "\n", "\n", "def fig5(df, fmt, outdir):\n", " \"\"\"Strength distribution by source, showing the unpopulated interval.\"\"\"\n", " fig, (ax_h, ax_b) = plt.subplots(\n", " 1, 2, figsize=(W2, 3.2), gridspec_kw={\"width_ratios\": [2, 1], \"wspace\": 0.28})\n", "\n", " bins = np.arange(0, 70, 2.5)\n", " for src in SOURCE_ORDER:\n", " s = df.loc[df[\"source\"] == src, \"strength\"]\n", " if len(s):\n", " ax_h.hist(s, bins=bins, color=SOURCE_FILL[src], alpha=0.85,\n", " hatch=SOURCE_HATCH[src], edgecolor=\"white\",\n", " linewidth=0.5, label=src)\n", " ax_h.axvspan(21.50, 43.00, color=\"0.5\", alpha=0.16, zorder=0)\n", " ax_h.text(32.25, ax_h.get_ylim()[1] * 0.55, \"no\\nobservations\",\n", " ha=\"center\", va=\"center\", fontsize=7.5, color=\"0.3\")\n", " ax_h.set_xlabel(\"Compressive strength (MPa)\")\n", " ax_h.set_ylabel(\"Number of observations\")\n", " ax_h.legend(frameon=False, loc=\"upper right\", handlelength=1.4)\n", " ax_h.text(0.015, 0.965, \"a\", transform=ax_h.transAxes,\n", " fontweight=\"bold\", fontsize=10)\n", "\n", " data = [df.loc[df[\"source\"] == s, \"strength\"].values for s in SOURCE_ORDER]\n", " # Updated 'labels' to 'tick_labels' to fix Matplotlib deprecation warning\n", " bp = ax_b.boxplot(data, tick_labels=SOURCE_ORDER, widths=0.55, patch_artist=True,\n", " medianprops=dict(color=\"black\", linewidth=1.1),\n", " flierprops=dict(marker=\"o\", markersize=3.2,\n", " markerfacecolor=\"white\",\n", " markeredgecolor=\"black\",\n", " markeredgewidth=0.6))\n", " for patch, src in zip(bp[\"boxes\"], SOURCE_ORDER):\n", " patch.set_facecolor(SOURCE_FILL[src])\n", " patch.set_hatch(SOURCE_HATCH[src])\n", " patch.set_edgecolor(\"black\")\n", " patch.set_linewidth(0.6)\n", " ax_b.axhspan(21.50, 43.00, color=\"0.5\", alpha=0.16, zorder=0)\n", " ax_b.set_ylabel(\"Compressive strength (MPa)\")\n", " ax_b.grid(axis=\"x\", visible=False)\n", " ax_b.tick_params(axis=\"x\", rotation=20)\n", " ax_b.text(0.03, 0.965, \"b\", transform=ax_b.transAxes,\n", " fontweight=\"bold\", fontsize=10)\n", "\n", " save(fig, \"Fig5\", fmt, outdir)\n", "\n", "\n", "def fig6(df, m, fmt, outdir):\n", " \"\"\"CPA measured vs predicted; extrapolated segment dotted.\"\"\"\n", " obs = df[(df[\"scm_type\"] == F6_SCM) &\n", " np.isclose(df[\"sulphate\"], F6_SULPHATE, atol=1e-3)]\n", "\n", " fig, ax = plt.subplots(figsize=(W2 * 0.86, 3.2))\n", "\n", " if len(obs):\n", " ax.scatter(obs[\"exposure\"], obs[\"strength\"], s=32, marker=\"o\",\n", " facecolor=\"white\", edgecolor=C_DARK, linewidth=0.9,\n", " zorder=3, label=\"Measured specimen\")\n", " g = obs.groupby(\"exposure\")[\"strength\"].mean()\n", " ax.scatter(g.index, g.values, s=46, marker=\"D\", color=C_DARK,\n", " zorder=4, label=\"Measured mean\")\n", "\n", " grid_in = np.linspace(obs[\"exposure\"].min() if len(obs) else 7, F6_OBS_MAX, 200)\n", " ax.plot(grid_in, [m[\"ridge\"].predict(\n", " scenario(m[\"X\"], df, F6_SCM, d, F6_SULPHATE))[0] for d in grid_in],\n", " color=C_LINE, lw=1.4, zorder=2, label=\"Ridge prediction\")\n", "\n", " grid_out = np.linspace(F6_OBS_MAX, F6_EXTRAP_MAX, 120)\n", " ax.plot(grid_out, [m[\"ridge\"].predict(\n", " scenario(m[\"X\"], df, F6_SCM, d, F6_SULPHATE))[0] for d in grid_out],\n", " color=C_LINE, lw=1.4, ls=\":\", zorder=2, label=\"Extrapolated\")\n", " ax.axvline(F6_OBS_MAX, color=\"0.55\", ls=\"--\", lw=0.7)\n", " ax.text(F6_OBS_MAX + 3, ax.get_ylim()[0] + 0.6,\n", " \"observed range ends\", fontsize=7, color=\"0.35\")\n", "\n", " ax.set_xlabel(\"Exposure duration (days)\")\n", " ax.set_ylabel(\"Compressive strength (MPa)\")\n", " ax.legend(frameon=False, loc=\"upper left\", ncol=2, columnspacing=1.1)\n", " save(fig, \"Fig6\", fmt, outdir)\n", "\n", "\n", "# ==============================================================================\n", "\n", "def main():\n", " ap = argparse.ArgumentParser()\n", " ap.add_argument(\"--data\", required=True)\n", " ap.add_argument(\"--outdir\", default=str(OUTDIR))\n", " ap.add_argument(\"--format\", default=\"tif\", choices=[\"tif\", \"eps\", \"pdf\", \"png\"],\n", " help=\"tif (default) for submission; png for quick checks\")\n", " # Simulate command-line arguments for Colab environment\n", " a = ap.parse_args(['--data', 'final_cleaned_sulfate_data_v2.csv'])\n", "\n", " set_style(a.format)\n", " outdir = Path(a.outdir)\n", "\n", " df = load(a.data)\n", " print(f\"Loaded {len(df)} observations, {df['source'].nunique()} sources\")\n", "\n", " m = fit_everything(df)\n", " print(f\"alpha = {m['alpha']:.7g} random-split R2 = {m['r2_random']:.4f}\")\n", " for _, r in m[\"loso\"].iterrows():\n", " print(f\" LOSO {r['source']}: R2 = {r['r2']:.2f}\")\n", " print()\n", "\n", " fig1(df, a.format, outdir)\n", " fig2(m, a.format, outdir)\n", " fig3(m, a.format, outdir)\n", " fig4(m, a.format, outdir)\n", " fig5(df, a.format, outdir)\n", " fig6(df, m, a.format, outdir)\n", "\n", " print(f\"\\nDone. Six figures in {outdir.resolve()}\")\n", " print(\"Check that the printed alpha and R2 above match Tables 3 and 5.\")\n", " print(\"If Arial is not installed, matplotlib falls back to DejaVu Sans and\")\n", " print(\"prints a warning. Install Arial or msttcorefonts before final export.\")\n", "\n", "\n", "if __name__ == \"__main__\":\n", " main()" ], "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "bf05k5nWQOD8", "outputId": "df573e51-dda2-4381-c2e2-524af373baa0" }, "execution_count": null, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "Loaded 94 observations, 3 sources\n", "alpha = 0.07149429 random-split R2 = 0.9676\n", " LOSO Source A: R2 = -157.92\n", " LOSO Source B: R2 = -1.83\n", " LOSO Source C: R2 = -14.55\n", "\n", " wrote figures/Fig1.tif\n", " wrote figures/Fig2.tif\n", " wrote figures/Fig3.tif\n", " wrote figures/Fig4.tif\n", " wrote figures/Fig5.tif\n", " wrote figures/Fig6.tif\n", "\n", "Done. Six figures in /content/figures\n", "Check that the printed alpha and R2 above match Tables 3 and 5.\n", "If Arial is not installed, matplotlib falls back to DejaVu Sans and\n", "prints a warning. Install Arial or msttcorefonts before final export.\n" ] } ] }, { "cell_type": "code", "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 35 }, "id": "9cc36f47", "outputId": "84501e70-a365-47ac-957b-16f50de247de" }, "source": [ "import shutil\n", "from pathlib import Path\n", "from google.colab import files\n", "\n", "# Define the directory and zip filename\n", "fig_dir = Path('figures')\n", "zip_name = 'manuscript_figures.zip'\n", "\n", "if fig_dir.exists() and fig_dir.is_dir():\n", " # Create the zip archive\n", " shutil.make_archive('manuscript_figures', 'zip', fig_dir)\n", " # Trigger the download\n", " files.download(zip_name)\n", " print(f'Successfully zipped {fig_dir} into {zip_name} and initiated download.')\n", "else:\n", " print(f'Error: The directory {fig_dir} was not found. Please run the figure generation cell first.')" ], "execution_count": null, "outputs": [ { "output_type": "display_data", "data": { "text/plain": [ "" ], "application/javascript": [ "\n", " async function download(id, filename, size) {\n", " if (!google.colab.kernel.accessAllowed) {\n", " return;\n", " }\n", " const div = document.createElement('div');\n", " const label = document.createElement('label');\n", " label.textContent = `Downloading \"${filename}\": `;\n", " div.appendChild(label);\n", " const progress = document.createElement('progress');\n", " progress.max = size;\n", " div.appendChild(progress);\n", " document.body.appendChild(div);\n", "\n", " const buffers = [];\n", " let downloaded = 0;\n", "\n", " const channel = await google.colab.kernel.comms.open(id);\n", " // Send a message to notify the kernel that we're ready.\n", " channel.send({})\n", "\n", " for await (const message of channel.messages) {\n", " // Send a message to notify the kernel that we're ready.\n", " channel.send({})\n", " if (message.buffers) {\n", " for (const buffer of message.buffers) {\n", " buffers.push(buffer);\n", " downloaded += buffer.byteLength;\n", " progress.value = downloaded;\n", " }\n", " }\n", " }\n", " const blob = new Blob(buffers, {type: 'application/binary'});\n", " const a = document.createElement('a');\n", " a.href = window.URL.createObjectURL(blob);\n", " a.download = filename;\n", " div.appendChild(a);\n", " a.click();\n", " div.remove();\n", " }\n", " " ] }, "metadata": {} }, { "output_type": "display_data", "data": { "text/plain": [ "" ], "application/javascript": [ "download(\"download_d72865a2-fc62-4dd9-9c56-6e27a49e9fc7\", \"manuscript_figures.zip\", 926166)" ] }, "metadata": {} }, { "output_type": "stream", "name": "stdout", "text": [ "Successfully zipped figures into manuscript_figures.zip and initiated download.\n" ] } ] } ] }