OSCR

An integrated machine learning framework for cerebral aneurysm rupture status classification using clinical, morphological, and hemodynamic features.

Code ↔ Paper

19 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 19 matches
  1. [1] § Methodology › Machine Learning ↔ Aneurysm_Classification.ipynb, lines 2479 to 2538 · score 0.66 · confusion matrix, roc_auc, cross validated, F1 score, precision, probabilities · Report a mismatch
  2. [2] § Methodology › Machine Learning ↔ aneurysm_classification.py, lines 2389 to 2445 · score 0.66 · confusion matrix, roc_auc, cross validated, F1 score, precision, probabilities · Report a mismatch
  3. [3] § Results ↔ Aneurysm_Classification.ipynb, lines 4740 to 4871 · score 0.64 · F1 score, BaggingSVM, Gaussian Process, ROC AUC, precision, Accuracy · Report a mismatch
  4. [4] § Methodology › Machine Learning ↔ Aneurysm_Classification.ipynb, lines 2272 to 2400 · score 0.61 · best mapped accuracy, unsupervised models, assignment, affinity, clustering, algorithms · Report a mismatch
  5. [5] § Results ↔ Aneurysm_Classification.ipynb, lines 2272 to 2400 · score 0.60 · best mapped accuracy, Agglomerative clustering, DBSCAN, Spectral, Birch, Unsupervised · Report a mismatch
  6. [6] § Results ↔ aneurysm_classification.py, lines 2223 to 2265 · score 0.60 · best mapped accuracy, Agglomerative clustering, DBSCAN, Spectral, Birch, Unsupervised · Report a mismatch
  7. [7] § Methodology › Clinical and morphological parameters ↔ Aneurysm_Classification.ipynb, lines 5817 to 5879 · score 0.59 · LB, numerically, BA, SCA, ACA, vessel · Report a mismatch
  8. [8] § Methodology › Clinical and morphological parameters ↔ aneurysm_classification.py, lines 5700 to 5766 · score 0.59 · LB, numerically, BA, SCA, ACA, vessel · Report a mismatch
  9. [9] § Results ↔ Aneurysm_Classification.ipynb, lines 1190 to 1297 · score 0.58 · GaussianNB, Decision Tree, Gradient Boosting, AdaBoost, SVM, classification · Report a mismatch
  10. [10] § Results ↔ aneurysm_classification.py, lines 1205 to 1289 · score 0.58 · GaussianNB, Decision Tree, Gradient Boosting, AdaBoost, SVM, classification · Report a mismatch
  11. [11] § Methodology › Machine Learning ↔ aneurysm_classification.py, lines 2223 to 2265 · score 0.57 · best mapped accuracy, unsupervised models, affinity, clustering, algorithms, classification · Report a mismatch
  12. [12] § Results ↔ aneurysm_classification.py, lines 4601 to 4657 · score 0.56 · BaggingSVM, Gaussian Process, ROC AUC, precision, score, recall · Report a mismatch
  13. [13] § Methodology › Study Population ↔ aneurysm_classification.py, lines 750 to 813 · score 0.56 · Barcelona, Geneva, Milan, Sheffield, SCA, ICA · Report a mismatch
  14. [14] § Methodology › Machine Learning ↔ Aneurysm_Classification.ipynb, lines 862 to 969 · score 0.53 · UMAP, t-SNE, distances, Neighbor, unruptured, classification · Report a mismatch
  15. [15] § Methodology › Machine Learning ↔ aneurysm_classification.py, lines 881 to 964 · score 0.53 · UMAP, t-SNE, distances, Neighbor, unruptured, classification · Report a mismatch
  16. [16] § Methodology › Machine Learning ↔ aneurysm_classification.py, lines 4113 to 4202 · score 0.52 · stratified split, leakage, Gaussian Process, tuning, fold, variables · Report a mismatch
  17. [17] § Results ↔ aneurysm_classification.py, lines 4601 to 4657 · score 0.52 · threshold tuning, BaggingSVM, Gaussian Process, ROC AUC, Precision, Recall · Report a mismatch
  18. [18] § Methodology › Machine Learning ↔ aneurysm_classification.py, lines 2999 to 3063 · score 0.51 · TensorFlow, XGBoost, boosting, models · Report a mismatch
  19. [19] § Methodology › Study Population ↔ aneurysm_classification.py, lines 59 to 163 · score 0.51 · Bar charts, outside, pale, slice, ICA, Radial · Report a mismatch

Paper

Loaded from Europe PMC by your browser; OSCR does not store it: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Python · 6,083 lines · 195 KB · no license · 12 matches

  1. # -*- coding: utf-8 -*-
  2. """Aneurysm_Classification.ipynb
  4. Automatically generated by Colab.
  6. Original file is located at
  7. https://colab.research.google.com/drive/1ul27HYChjrleBP998ykEHyUqMO-X2Z-_
  9. #Study Population
  10. """
  12. import numpy as np
  13. import matplotlib.pyplot as plt
  14. from matplotlib.patches import Wedge
  15. from mpl_toolkits.axes_grid1.inset_locator import inset_axes
  17. # ----------------------------
  18. # GLOBAL STYLE (Times New Roman + bigger fonts)
  19. # ----------------------------
  20. plt.rcParams["font.family"] = "serif"
  21. plt.rcParams["font.serif"] = ["Times New Roman", "Times", "DejaVu Serif"]
  22. plt.rcParams.update({
  23. "font.size": 20,
  24. "axes.titlesize": 20,
  25. "axes.labelsize": 20,
  26. "xtick.labelsize": 20,
  27. "ytick.labelsize": 20,
  28. "legend.fontsize": 22,
  29. })
  31. # ----------------------------
  32. # DATA
  33. # ----------------------------
  34. data = {
  35. "Overall": {"UF": 351, "UM": 105, "RF": 169, "RM": 83, "F_total": 520, "M_total": 188, "Total": 708},
  36. "ICA": {"UF": 187, "UM": 36, "RF": 79, "RM": 21, "F_total": 266, "M_total": 57, "Total": 323},
  37. "MCA": {"UF": 99, "UM": 41, "RF": 23, "RM": 18, "F_total": 122, "M_total": 59, "Total": 181},
  38. "ACA": {"UF": 36, "UM": 24, "RF": 47, "RM": 36, "F_total": 83, "M_total": 60, "Total": 143},
  39. "BA/VA": {"UF": 18, "UM": 2, "RF": 15, "RM": 5, "F_total": 33, "M_total": 7, "Total": 40},
  40. "PCA": {"UF": 9, "UM": 0, "RF": 4, "RM": 3, "F_total": 13, "M_total": 3, "Total": 16},
  41. "SCA": {"UF": 2, "UM": 1, "RF": 1, "RM": 1, "F_total": 3, "M_total": 2, "Total": 5},
  42. }
  43. bar_locs = ["Overall", "ICA", "MCA", "ACA", "BA/VA", "PCA", "SCA"]
  44. pie_locs = ["ICA", "MCA", "ACA", "BA/VA", "PCA", "SCA"]
  46. # Colors (requested)
  47. COL_F = "#C2185B" # Female = pink
  48. COL_M = "#2E5EAA" # Male = blue
  49. COL_T = "#4A4A4A" # Total = green
  50. TXT = "#0f172a"
  51. pale_alpha = 0.22
  52. def pct(x, denom):
  53. return 0.0 if denom == 0 else 100.0 * x / denom
  55. # ----------------------------
  56. # Radial pie with radius-proportional ruptured portion
  57. # (using sqrt so solid AREA is proportional to rupture%)
  58. # ----------------------------
  59. def rupture_radial_pie(ax, locs, rup_counts, tot_counts, vivid_color, title):
  60. ax.set_aspect("equal")
  61. ax.axis("off")
  63. n = len(locs)
  64. start_angle = 90
  65. slice_angle = 360.0 / n
  67. R_outer = 1.0
  68. label_r = 1.32
  69. name_r = 1.55
  70. pale_alpha = 0.22
  72. ax.text(0, 1.75, title, ha="center", va="bottom",
  73. fontsize=26, fontweight="bold", color=vivid_color)
  75. for i, loc in enumerate(locs):
  76. rup = float(rup_counts[i])
  77. tot = float(tot_counts[i])
  78. rr = 0.0 if tot == 0 else rup / tot
  80. th1 = start_angle - i * slice_angle
  81. th2 = th1 - slice_angle
  83. # Pale full slice
  84. ax.add_patch(Wedge(
  85. (0, 0), R_outer, th2, th1,
  86. width=R_outer,
  87. facecolor=(plt.matplotlib.colors.to_rgba(vivid_color, pale_alpha)),
  88. edgecolor="white", linewidth=2.2
  89. ))
  91. # Solid ruptured fill (radius)
  92. if rr > 0:
  93. r_solid = R_outer * np.sqrt(rr)
  94. ax.add_patch(Wedge(
  95. (0, 0), r_solid, th2, th1,
  96. width=r_solid,
  97. facecolor=vivid_color,
  98. edgecolor="white", linewidth=2.2
  99. ))
  101. mid = np.deg2rad((th1 + th2) / 2.0)
  102. # Format percentage (remove decimal if exactly 100)
  103. rr_pct = rr * 100.0
  104. if np.isclose(rr_pct, 100.0):
  105. pct_text = "100%"
  106. else:
  107. pct_text = f"{rr_pct:.1f}%"
  111. # # % label outside
  112. # ax.text(label_r * np.cos(mid), label_r * np.sin(mid),
  113. # f"{rr_pct:.1f}%", ha="center", va="center",
  114. # fontsize=18, color=TXT,
  115. # bbox=dict(boxstyle="round,pad=0.25", fc="white", ec="none", alpha=0.92))
  117. # # Location label
  118. # ax.text(name_r * np.cos(mid), name_r * np.sin(mid),
  119. # loc, ha="center", va="center",
  120. # fontsize=18, color=TXT)
  123. # Default positions (for all slices except ICA and PCA)
  124. x_pct = label_r * np.cos(mid)
  125. y_pct = label_r * np.sin(mid)
  127. x_loc = name_r * np.cos(mid)
  128. y_loc = name_r * np.sin(mid)
  130. # Special adjustment ONLY for ICA and PCA
  131. if loc in ["MCA", "PCA"]:
  132. extra_r = 0.06 # move further away from pie
  133. x_pct = (label_r + extra_r) * np.cos(mid)
  134. y_pct = (label_r + extra_r) * np.sin(mid)
  136. # Stack percentage above, location below
  137. offset = 0.06
  138. ax.text(x_pct, y_pct,
  139. pct_text, ha="center", va="bottom",
  140. fontsize=18, color=TXT,
  141. bbox=dict(boxstyle="round,pad=0.25", fc="white", ec="none", alpha=0.92))
  143. ax.text(x_pct, y_pct - offset,
  144. loc, ha="center", va="top",
  145. fontsize=18, color=TXT)
  147. else:
  148. # Original behavior for other slices
  149. ax.text(x_pct, y_pct,
  150. pct_text, ha="center", va="center",
  151. fontsize=18, color=TXT,
  152. bbox=dict(boxstyle="round,pad=0.25", fc="white", ec="none", alpha=0.92))
  154. ax.text(x_loc, y_loc,
  155. loc, ha="center", va="center",
  156. fontsize=18, color=TXT)
  158. ax.set_xlim(-1.95, 1.95)
  159. ax.set_ylim(-1.65, 2.05)
  161. # ----------------------------
  162. # BAR CHART (counts) + labels on bars
  163. # ----------------------------
  164. F_counts = np.array([data[l]["F_total"] for l in bar_locs])
  165. M_counts = np.array([data[l]["M_total"] for l in bar_locs])
  166. T_counts = np.array([data[l]["Total"] for l in bar_locs])
  168. x = np.arange(len(bar_locs))
  169. w = 0.24
  171. fig, ax = plt.subplots(figsize=(16, 7.6), dpi=300)
  173. barsF = ax.bar(x - w, F_counts, width=w, color=COL_F, alpha=0.88, label="Female")
  174. barsM = ax.bar(x, M_counts, width=w, color=COL_M, alpha=0.88, label="Male")
  175. barsT = ax.bar(x + w, T_counts, width=w, color=COL_T, alpha=0.80, label="Total")
  177. # Numbers on bars
  178. pad = max(T_counts) * 0.012
  179. def label_bars(bars):
  180. for b in bars:
  181. h = b.get_height()
  182. ax.text(b.get_x() + b.get_width()/2, h + pad, f"{int(h)}",
  183. ha="center", va="bottom", fontsize=18, color=TXT)
  184. label_bars(barsF); label_bars(barsM); label_bars(barsT)
  186. # Clean axes (no y axis)
  187. ax.set_yticks([])
  188. ax.set_ylabel("")
  189. ax.spines["left"].set_visible(False)
  190. ax.spines["top"].set_visible(False)
  191. ax.spines["right"].set_visible(False)
  193. ax.set_xticks(x)
  194. ax.set_xticklabels(bar_locs, fontsize=20)
  196. # Vertical legend
  197. ax.legend(
  198. ncol=1, frameon=False, loc="upper left",
  199. bbox_to_anchor=(0.17, 1.0), borderaxespad=0.0
  200. )
  202. # ----------------------------
  203. # BIG pies, moved up, farther apart
  204. # ----------------------------
  205. female_rup = np.array([data[l]["RF"] for l in pie_locs], float)
  206. female_tot = np.array([data[l]["F_total"] for l in pie_locs], float)
  208. male_rup = np.array([data[l]["RM"] for l in pie_locs], float)
  209. male_tot = np.array([data[l]["M_total"] for l in pie_locs], float)
  211. # Placement knobs
  212. y0 = 0.35 # higher = move up
  213. h = 0.8 # taller = bigger
  214. w_in = 0.35 # wider = bigger
  215. gap = 0.0 # farther apart
  217. x_left = 0.3
  218. x_right = x_left + w_in + gap
  220. ax_f = inset_axes(
  221. ax, width="100%", height="100%",
  222. bbox_to_anchor=(x_left, y0, w_in, h),
  223. bbox_transform=ax.transAxes, loc="lower left", borderpad=0
  224. )
  225. ax_m = inset_axes(
  226. ax, width="100%", height="100%",
  227. bbox_to_anchor=(x_right, y0, w_in, h),
  228. bbox_transform=ax.transAxes, loc="lower left", borderpad=0
  229. )
  231. rupture_radial_pie(ax_f, pie_locs, female_rup, female_tot, COL_F, "Female")
  232. rupture_radial_pie(ax_m, pie_locs, male_rup, male_tot, COL_M, "Male")
  234. plt.tight_layout()
  235. plt.savefig("bars_pies_times_new_roman.png", dpi=300, bbox_inches="tight")
  236. plt.savefig("bars_pies_times_new_roman.pdf", bbox_inches="tight")
  237. plt.show()
  239. import numpy as np
  240. import matplotlib.pyplot as plt
  241. from matplotlib.patches import Wedge
  242. from mpl_toolkits.axes_grid1.inset_locator import inset_axes
  244. # ----------------------------
  245. # GLOBAL STYLE (serif fallback for Colab)
  246. # ----------------------------
  247. plt.rcParams["font.family"] = "serif"
  248. plt.rcParams["font.serif"] = ["Times New Roman", "Times", "DejaVu Serif"]
  249. plt.rcParams.update({
  250. "font.size": 20,
  251. "axes.titlesize": 20,
  252. "axes.labelsize": 20,
  253. "xtick.labelsize": 20,
  254. "ytick.labelsize": 20,
  255. "legend.fontsize": 22,
  256. })
  258. # ----------------------------
  259. # DATA
  260. # ----------------------------
  261. data = {
  262. "Overall": {"UF": 351, "UM": 105, "RF": 169, "RM": 83, "F_total": 520, "M_total": 188, "Total": 708},
  263. "ICA": {"UF": 187, "UM": 36, "RF": 79, "RM": 21, "F_total": 266, "M_total": 57, "Total": 323},
  264. "MCA": {"UF": 99, "UM": 41, "RF": 23, "RM": 18, "F_total": 122, "M_total": 59, "Total": 181},
  265. "ACA": {"UF": 36, "UM": 24, "RF": 47, "RM": 36, "F_total": 83, "M_total": 60, "Total": 143},
  266. "BA/VA": {"UF": 18, "UM": 2, "RF": 15, "RM": 5, "F_total": 33, "M_total": 7, "Total": 40},
  267. "PCA": {"UF": 9, "UM": 0, "RF": 4, "RM": 3, "F_total": 13, "M_total": 3, "Total": 16},
  268. "SCA": {"UF": 2, "UM": 1, "RF": 1, "RM": 1, "F_total": 3, "M_total": 2, "Total": 5},
  269. }
  271. bar_locs = ["Overall", "ICA", "MCA", "ACA", "BA/VA", "PCA", "SCA"]
  272. pie_locs = ["ICA", "MCA", "ACA", "BA/VA", "PCA", "SCA"]
  274. # Colors (your palette)
  275. COL_F = "#C2185B" # Female (pink)
  276. COL_M = "#2E5EAA" # Male (blue)
  277. COL_T = "#4A4A4A" # Total (charcoal)
  278. TXT = "#0f172a"
  279. pale_alpha = 0.22
  281. def pct(x, denom):
  282. return 0.0 if denom == 0 else 100.0 * x / denom
  284. # ----------------------------
  285. # Radial pie with radius-proportional ruptured portion
  286. # ----------------------------
  287. def rupture_radial_pie(ax, locs, rup_counts, tot_counts, vivid_color, title):
  288. ax.set_aspect("equal")
  289. ax.axis("off")
  291. n = len(locs)
  292. start_angle = 90
  293. slice_angle = 360.0 / n
  295. R_outer = 1.0
  296. label_r = 1.32
  297. name_r = 1.55
  299. ax.text(0, 1.75, title, ha="center", va="bottom",
  300. fontsize=26, fontweight="bold", color=vivid_color)
  302. for i, loc in enumerate(locs):
  303. rup = float(rup_counts[i])
  304. tot = float(tot_counts[i])
  305. rr = 0.0 if tot == 0 else rup / tot
  307. th1 = start_angle - i * slice_angle
  308. th2 = th1 - slice_angle
  310. # Pale full slice
  311. ax.add_patch(Wedge(
  312. (0, 0), R_outer, th2, th1,
  313. width=R_outer,
  314. facecolor=(plt.matplotlib.colors.to_rgba(vivid_color, pale_alpha)),
  315. edgecolor="white", linewidth=2.2
  316. ))
  318. # Solid ruptured fill: sqrt so AREA matches rupture%
  319. if rr > 0:
  320. r_solid = R_outer * np.sqrt(rr)
  321. ax.add_patch(Wedge(
  322. (0, 0), r_solid, th2, th1,
  323. width=r_solid,
  324. facecolor=vivid_color,
  325. edgecolor="white", linewidth=2.2
  326. ))
  328. mid = np.deg2rad((th1 + th2) / 2.0)
  330. rr_pct = rr * 100.0
  331. if np.isclose(rr_pct, 100.0):
  332. pct_text = "100%"
  333. else:
  334. pct_text = f"{rr_pct:.1f}%"
  336. # Default positions
  337. x_pct = label_r * np.cos(mid)
  338. y_pct = label_r * np.sin(mid)
  339. x_loc = name_r * np.cos(mid)
  340. y_loc = name_r * np.sin(mid)
  342. # Special handling ONLY for MCA and PCA (as you coded)
  343. if loc in ["MCA", "PCA"]:
  344. extra_r = 0.06
  345. x_pct = (label_r + extra_r) * np.cos(mid)
  346. y_pct = (label_r + extra_r) * np.sin(mid)
  348. offset = 0.06
  349. ax.text(x_pct, y_pct,
  350. pct_text, ha="center", va="bottom",
  351. fontsize=18, color=TXT,
  352. bbox=dict(boxstyle="round,pad=0.25", fc="white", ec="none", alpha=0.92))
  353. ax.text(x_pct, y_pct - offset,
  354. loc, ha="center", va="top",
  355. fontsize=18, color=TXT)
  356. else:
  357. ax.text(x_pct, y_pct,
  358. pct_text, ha="center", va="center",
  359. fontsize=18, color=TXT,
  360. bbox=dict(boxstyle="round,pad=0.25", fc="white", ec="none", alpha=0.92))
  361. ax.text(x_loc, y_loc,
  362. loc, ha="center", va="center",
  363. fontsize=18, color=TXT)
  365. ax.set_xlim(-1.95, 1.95)
  366. ax.set_ylim(-1.65, 2.05)
  368. # ----------------------------
  369. # BAR CHART (counts) + labels
  370. # ----------------------------
  371. F_counts = np.array([data[l]["F_total"] for l in bar_locs])
  372. M_counts = np.array([data[l]["M_total"] for l in bar_locs])
  373. T_counts = np.array([data[l]["Total"] for l in bar_locs])
  375. x = np.arange(len(bar_locs))
  376. w = 0.24
  378. fig, ax = plt.subplots(figsize=(16, 7.6), dpi=300)
  380. barsF = ax.bar(x - w, F_counts, width=w, color=COL_F, alpha=0.88)
  381. barsM = ax.bar(x, M_counts, width=w, color=COL_M, alpha=0.88)
  382. barsT = ax.bar(x + w, T_counts, width=w, color=COL_T, alpha=0.80)
  384. # numbers on bars
  385. pad = max(T_counts) * 0.012
  386. def label_bars(bars):
  387. for b in bars:
  388. h = b.get_height()
  389. ax.text(b.get_x() + b.get_width()/2, h + pad, f"{int(h)}",
  390. ha="center", va="bottom", fontsize=18, color=TXT)
  391. label_bars(barsF); label_bars(barsM); label_bars(barsT)
  393. # Clean axes (no y axis)
  394. ax.set_yticks([])
  395. ax.set_ylabel("")
  396. ax.spines["left"].set_visible(False)
  397. ax.spines["top"].set_visible(False)
  398. ax.spines["right"].set_visible(False)
  400. ax.set_xticks(x)
  401. ax.set_xticklabels(bar_locs, fontsize=20)
  403. # ✅ remove legend entirely
  404. # ax.legend(...) <-- deleted on purpose
  406. # ----------------------------
  407. # PIES: Female + Male moved left, add Total pie
  408. # ----------------------------
  409. female_rup = np.array([data[l]["RF"] for l in pie_locs], float)
  410. female_tot = np.array([data[l]["F_total"] for l in pie_locs], float)
  412. male_rup = np.array([data[l]["RM"] for l in pie_locs], float)
  413. male_tot = np.array([data[l]["M_total"] for l in pie_locs], float)
  415. # Total across sexes, per location
  416. total_rup = np.array([data[l]["RF"] + data[l]["RM"] for l in pie_locs], float)
  417. total_tot = np.array([data[l]["Total"] for l in pie_locs], float)
  419. # Placement knobs (3 pies in a row)
  420. y0 = 0.4
  421. h = 0.80
  422. w_in = 0.3
  423. gap = -0.02
  425. # Shift row to the LEFT
  426. x1 = 0.15
  427. x2 = x1 + w_in + gap
  428. x3 = x2 + w_in + gap
  430. ax_f = inset_axes(ax, width="100%", height="100%",
  431. bbox_to_anchor=(x1, y0, w_in, h),
  432. bbox_transform=ax.transAxes, loc="lower left", borderpad=0)
  434. ax_m = inset_axes(ax, width="100%", height="100%",
  435. bbox_to_anchor=(x2, y0, w_in, h),
  436. bbox_transform=ax.transAxes, loc="lower left", borderpad=0)
  438. ax_t = inset_axes(ax, width="100%", height="100%",
  439. bbox_to_anchor=(x3, y0, w_in, h),
  440. bbox_transform=ax.transAxes, loc="lower left", borderpad=0)
  442. rupture_radial_pie(ax_f, pie_locs, female_rup, female_tot, COL_F, "Female")
  443. rupture_radial_pie(ax_m, pie_locs, male_rup, male_tot, COL_M, "Male")
  444. rupture_radial_pie(ax_t, pie_locs, total_rup, total_tot, COL_T, "Total")
  446. plt.tight_layout()
  447. plt.savefig("bars_3pies_times_serif.png", dpi=300, bbox_inches="tight")
  448. plt.savefig("bars_3pies_times_serif.pdf", bbox_inches="tight")
  449. plt.show()
  451. import pandas as pd
  452. import numpy as np
  453. import matplotlib.pyplot as plt
  454. from matplotlib.patches import Wedge
  455. from mpl_toolkits.axes_grid1.inset_locator import inset_axes
  456. import os
  458. # ============================
  459. # LOAD DATA
  460. # ============================
  461. df = pd.read_excel(os.path.join('..', 'content', 'ClassificationV2.xlsx'))
  463. mapping_Sex = {'Male': 0, 'Female': 1}
  464. mapping_Irregularity = {'R': 0, 'I': 1}
  465. mapping_Kind = {'L': 0, 'B': 1}
  466. mapping_Location = {'MCA': 0, 'ICA': 1, 'ACA': 2, 'PCA': 3, 'VA': 4, 'BA': 4, 'SCA': 5}
  467. mapping_Stability = {'stable': 0, 'unstable': 1}
  468. mapping_Complexity = {'simple': 0, 'complex': 1}
  469. mapping_Status = {'unruptured': 0, 'ruptured': 1}
  471. df = df.replace({
  472. 'Sex': mapping_Sex,
  473. 'I': mapping_Irregularity,
  474. 'LB': mapping_Kind,
  475. 'Location': mapping_Location,
  476. 'FS': mapping_Stability,
  477. 'FC': mapping_Complexity,
  478. 'Status': mapping_Status
  479. })
  481. df = df.drop(labels='Vessel ID', axis=1)
  483. # ============================
  484. # HOLD OUT BARCELONA
  485. # rows 557 to 632 inclusive
  486. # ============================
  487. barcelona_start = 557
  488. barcelona_end = 632
  490. test_idx = np.arange(barcelona_start, barcelona_end + 1)
  491. train_idx = df.index.difference(test_idx)
  493. df_train = df.loc[train_idx].copy()
  494. df_test = df.loc[test_idx].copy()
  496. print("Train shape:", df_train.shape)
  497. print("Test shape:", df_test.shape)
  499. # ============================
  500. # MAKE A PLOTTING COPY WITH LABELS
  501. # ============================
  502. sex_inverse = {0: 'Male', 1: 'Female'}
  503. status_inverse = {0: 'Unruptured', 1: 'Ruptured'}
  504. location_inverse = {
  505. 0: 'MCA',
  506. 1: 'ICA',
  507. 2: 'ACA',
  508. 3: 'PCA',
  509. 4: 'BA/VA',
  510. 5: 'SCA'
  511. }
  513. def make_plotting_df(df_in):
  514. d = df_in.copy()
  515. d['Sex_label'] = d['Sex'].map(sex_inverse)
  516. d['Status_label'] = d['Status'].map(status_inverse)
  517. d['Location_label'] = d['Location'].map(location_inverse)
  518. return d
  520. plot_train = make_plotting_df(df_train)
  521. plot_test = make_plotting_df(df_test)
  523. # ============================
  524. # BUILD DATA DICTIONARY
  525. # ============================
  526. bar_locs = ["Overall", "ICA", "MCA", "ACA", "BA/VA", "PCA", "SCA"]
  527. pie_locs = ["ICA", "MCA", "ACA", "BA/VA", "PCA", "SCA"]
  528. ordered_locs = ["ICA", "MCA", "ACA", "BA/VA", "PCA", "SCA"]
  530. def build_plot_data(df_subset):
  531. result = {}
  533. UF = ((df_subset['Status_label'] == "Unruptured") & (df_subset['Sex_label'] == "Female")).sum()
  534. UM = ((df_subset['Status_label'] == "Unruptured") & (df_subset['Sex_label'] == "Male")).sum()
  535. RF = ((df_subset['Status_label'] == "Ruptured") & (df_subset['Sex_label'] == "Female")).sum()
  536. RM = ((df_subset['Status_label'] == "Ruptured") & (df_subset['Sex_label'] == "Male")).sum()
  538. result["Overall"] = {
  539. "UF": int(UF),
  540. "UM": int(UM),
  541. "RF": int(RF),
  542. "RM": int(RM),
  543. "F_total": int((df_subset['Sex_label'] == "Female").sum()),
  544. "M_total": int((df_subset['Sex_label'] == "Male").sum()),
  545. "Total": int(len(df_subset)),
  546. }
  548. for loc in ordered_locs:
  549. dloc = df_subset[df_subset['Location_label'] == loc]
  551. UF = ((dloc['Status_label'] == "Unruptured") & (dloc['Sex_label'] == "Female")).sum()
  552. UM = ((dloc['Status_label'] == "Unruptured") & (dloc['Sex_label'] == "Male")).sum()
  553. RF = ((dloc['Status_label'] == "Ruptured") & (dloc['Sex_label'] == "Female")).sum()
  554. RM = ((dloc['Status_label'] == "Ruptured") & (dloc['Sex_label'] == "Male")).sum()
  556. result[loc] = {
  557. "UF": int(UF),
  558. "UM": int(UM),
  559. "RF": int(RF),
  560. "RM": int(RM),
  561. "F_total": int((dloc['Sex_label'] == "Female").sum()),
  562. "M_total": int((dloc['Sex_label'] == "Male").sum()),
  563. "Total": int(len(dloc)),
  564. }
  566. return result
  568. train_data = build_plot_data(plot_train)
  569. test_data = build_plot_data(plot_test)
  571. print("Training Overall:", train_data["Overall"])
  572. print("Testing Overall:", test_data["Overall"])
  574. # ============================
  575. # GLOBAL STYLE
  576. # ============================
  577. plt.rcParams["font.family"] = "serif"
  578. plt.rcParams["font.serif"] = ["Times New Roman", "Times", "DejaVu Serif"]
  579. plt.rcParams.update({
  580. "font.size": 20,
  581. "axes.titlesize": 20,
  582. "axes.labelsize": 20,
  583. "xtick.labelsize": 20,
  584. "ytick.labelsize": 20,
  585. "legend.fontsize": 22,
  586. })
  588. COL_F = "#C2185B"
  589. COL_M = "#2E5EAA"
  590. COL_T = "#4A4A4A"
  591. TXT = "#0f172a"
  592. pale_alpha = 0.22
  594. def rupture_radial_pie(ax, locs, rup_counts, tot_counts, vivid_color, title):
  595. ax.set_aspect("equal")
  596. ax.axis("off")
  598. n = len(locs)
  599. start_angle = 90
  600. slice_angle = 360.0 / n
  602. R_outer = 1.0
  603. label_r = 1.32
  604. name_r = 1.55
  606. ax.text(0, 1.75, title, ha="center", va="bottom",
  607. fontsize=26, fontweight="bold", color=vivid_color)
  609. for i, loc in enumerate(locs):
  610. rup = float(rup_counts[i])
  611. tot = float(tot_counts[i])
  612. rr = 0.0 if tot == 0 else rup / tot
  614. th1 = start_angle - i * slice_angle
  615. th2 = th1 - slice_angle
  617. ax.add_patch(Wedge(
  618. (0, 0), R_outer, th2, th1,
  619. width=R_outer,
  620. facecolor=(plt.matplotlib.colors.to_rgba(vivid_color, pale_alpha)),
  621. edgecolor="white", linewidth=2.2
  622. ))
  624. if rr > 0:
  625. r_solid = R_outer * np.sqrt(rr)
  626. ax.add_patch(Wedge(
  627. (0, 0), r_solid, th2, th1,
  628. width=r_solid,
  629. facecolor=vivid_color,
  630. edgecolor="white", linewidth=2.2
  631. ))
  633. mid = np.deg2rad((th1 + th2) / 2.0)
  634. rr_pct = rr * 100.0
  635. pct_text = "100%" if np.isclose(rr_pct, 100.0) else f"{rr_pct:.1f}%"
  637. x_pct = label_r * np.cos(mid)
  638. y_pct = label_r * np.sin(mid)
  639. x_loc = name_r * np.cos(mid)
  640. y_loc = name_r * np.sin(mid)
  642. if loc in ["MCA", "PCA"]:
  643. extra_r = 0.06
  644. x_pct = (label_r + extra_r) * np.cos(mid)
  645. y_pct = (label_r + extra_r) * np.sin(mid)
  647. offset = 0.06
  648. ax.text(x_pct, y_pct, pct_text,
  649. ha="center", va="bottom",
  650. fontsize=18, color=TXT,
  651. bbox=dict(boxstyle="round,pad=0.25", fc="white", ec="none", alpha=0.92))
  652. ax.text(x_pct, y_pct - offset, loc,
  653. ha="center", va="top",
  654. fontsize=18, color=TXT)
  655. else:
  656. ax.text(x_pct, y_pct, pct_text,
  657. ha="center", va="center",
  658. fontsize=18, color=TXT,
  659. bbox=dict(boxstyle="round,pad=0.25", fc="white", ec="none", alpha=0.92))
  660. ax.text(x_loc, y_loc, loc,
  661. ha="center", va="center",
  662. fontsize=18, color=TXT)
  664. ax.set_xlim(-1.95, 1.95)
  665. ax.set_ylim(-1.65, 2.05)
  667. def plot_dataset_summary(data_dict, figure_title=None, save_prefix=None):
  668. F_counts = np.array([data_dict[l]["F_total"] for l in bar_locs])
  669. M_counts = np.array([data_dict[l]["M_total"] for l in bar_locs])
  670. T_counts = np.array([data_dict[l]["Total"] for l in bar_locs])
  672. x = np.arange(len(bar_locs))
  673. w = 0.24
  675. fig, ax = plt.subplots(figsize=(16, 7.6), dpi=300)
  677. barsF = ax.bar(x - w, F_counts, width=w, color=COL_F, alpha=0.88)
  678. barsM = ax.bar(x, M_counts, width=w, color=COL_M, alpha=0.88)
  679. barsT = ax.bar(x + w, T_counts, width=w, color=COL_T, alpha=0.80)
  681. pad = max(T_counts) * 0.012 if max(T_counts) > 0 else 0.2
  683. def label_bars(bars):
  684. for b in bars:
  685. h = b.get_height()
  686. ax.text(b.get_x() + b.get_width()/2, h + pad, f"{int(h)}",
  687. ha="center", va="bottom", fontsize=18, color=TXT)
  689. label_bars(barsF)
  690. label_bars(barsM)
  691. label_bars(barsT)
  693. ax.set_yticks([])
  694. ax.set_ylabel("")
  695. ax.spines["left"].set_visible(False)
  696. ax.spines["top"].set_visible(False)
  697. ax.spines["right"].set_visible(False)
  699. ax.set_xticks(x)
  700. ax.set_xticklabels(bar_locs, fontsize=20)
  702. if figure_title is not None:
  703. ax.set_title(figure_title, pad=18)
  705. female_rup = np.array([data_dict[l]["RF"] for l in pie_locs], float)
  706. female_tot = np.array([data_dict[l]["F_total"] for l in pie_locs], float)
  708. male_rup = np.array([data_dict[l]["RM"] for l in pie_locs], float)
  709. male_tot = np.array([data_dict[l]["M_total"] for l in pie_locs], float)
  711. total_rup = np.array([data_dict[l]["RF"] + data_dict[l]["RM"] for l in pie_locs], float)
  712. total_tot = np.array([data_dict[l]["Total"] for l in pie_locs], float)
  714. y0 = 0.4
  715. h = 0.80
  716. w_in = 0.3
  717. gap = -0.02
  719. x1 = 0.15
  720. x2 = x1 + w_in + gap
  721. x3 = x2 + w_in + gap
  723. ax_f = inset_axes(ax, width="100%", height="100%",
  724. bbox_to_anchor=(x1, y0, w_in, h),
  725. bbox_transform=ax.transAxes, loc="lower left", borderpad=0)
  727. ax_m = inset_axes(ax, width="100%", height="100%",
  728. bbox_to_anchor=(x2, y0, w_in, h),
  729. bbox_transform=ax.transAxes, loc="lower left", borderpad=0)
  731. ax_t = inset_axes(ax, width="100%", height="100%",
  732. bbox_to_anchor=(x3, y0, w_in, h),
  733. bbox_transform=ax.transAxes, loc="lower left", borderpad=0)
  735. rupture_radial_pie(ax_f, pie_locs, female_rup, female_tot, COL_F, "Female")
  736. rupture_radial_pie(ax_m, pie_locs, male_rup, male_tot, COL_M, "Male")
  737. rupture_radial_pie(ax_t, pie_locs, total_rup, total_tot, COL_T, "Total")
  739. plt.tight_layout()
  741. if save_prefix is not None:
  742. plt.savefig(f"{save_prefix}.png", dpi=300, bbox_inches="tight")
  743. plt.savefig(f"{save_prefix}.pdf", bbox_inches="tight")
  745. plt.show()
  747. # ============================
  748. # PLOT TRAIN AND TEST
  749. # ============================
  750. plot_dataset_summary(
  751. train_data,
  752. figure_title="Training Set (Sheffield + Milan + Geneva)",
  753. save_prefix="train_set_summary"
  754. )
  756. plot_dataset_summary(
  757. test_data,
  758. figure_title="Testing Set (Barcelona Holdout)",
  759. save_prefix="test_set_summary"
  760. )
  762. """#Load Dateset"""
  764. import pandas as pd
  765. import numpy as np
  766. import os
  768. df = pd.read_excel(os.path.join('..','content','ClassificationV2.xlsx'))
  769. mapping_Sex = {'Male': 0, 'Female': 1}
  770. mapping_Irregularity = {'R': 0, 'I': 1}
  771. mapping_Kind = {'L': 0, 'B': 1}
  772. mapping_Location = {'MCA': 0, 'ICA': 1, 'ACA': 2, 'PCA': 3, 'VA': 4, 'BA': 4, 'SCA': 5}
  773. mapping_Stability = {'stable': 0, 'unstable': 1}
  774. mapping_Complexity = {'simple': 0, 'complex': 1}
  775. mapping_Status = {'unruptured': 0, 'ruptured': 1}
  776. df = df.replace({'Sex': mapping_Sex, 'I': mapping_Irregularity, 'LB': mapping_Kind, 'Location': mapping_Location, 'FS': mapping_Stability, 'FC': mapping_Complexity, 'Status': mapping_Status})
  777. df = df.drop(labels = 'Vessel ID', axis = 1)
  778. H = np.array(df)
  779. H = np.array(list(H[:,:]), dtype=np.float64)
  781. """## Data Standardization & Split"""
  783. from sklearn.model_selection import train_test_split
  784. from sklearn.preprocessing import StandardScaler
  786. y = H[:,-1]
  787. x = H[:,:-1]
  788. x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.2, random_state=8)
  789. standardization = StandardScaler().fit(x_train)
  790. x_train = standardization.transform(x_train)
  791. x_test = standardization.transform(x_test)
  793. """## Data Scattering"""
  795. import matplotlib.pyplot as plt
  796. from matplotlib.font_manager import FontProperties
  799. font = {'family' : 'arial',
  800. 'serif' : ['arial'],
  801. 'weight' : 'medium',
  802. 'size' : 24,
  803. 'style' : 'normal'}
  806. u, r = list(),list()
  807. for i in range(len(y)):
  808. if y[i] == 0:
  809. u.append(i)
  810. elif y[i] == 1:
  811. r.append(i)
  813. attributes = df.columns.values
  816. def find_index(name):
  817. for i in range(len(attributes)):
  818. if name == attributes[i]:
  819. return i
  822. def data_scatter3d(u, r, first_axis, second_axis, third_axis):
  823. first_axis = find_index(first_axis)
  824. second_axis = find_index(second_axis)
  825. third_axis = find_index(third_axis)
  826. fig = plt.figure(figsize=(9,7))
  827. ax = plt.subplot(111,projection='3d')
  828. ax.plot(x[u,first_axis],x[u,second_axis],x[u,third_axis],'o',markerfacecolor=(0,1,0.2),markeredgecolor='k',ms=10,label='Unruptured')
  829. ax.plot(x[r,first_axis],x[r,second_axis],x[r,third_axis],'o',markerfacecolor='red',markeredgecolor='k',ms=10,label='Ruptured')
  830. ax.set_xlabel(attributes[first_axis], labelpad=15)
  831. ax.set_ylabel(attributes[second_axis], labelpad=15)
  832. ax.set_zlabel(attributes[third_axis], labelpad=25)
  833. ax.tick_params(axis='z', which='major', pad=10)
  834. # ax.view_init(elev=25, azim=40)
  835. # fig.tight_layout()
  836. ax.legend(loc='best', bbox_to_anchor=(1.2, 1))
  837. ax.dist = 12
  838. plt.savefig('3D_DataScattering.pdf', transparent=True, dpi=300)
  839. plt.show()
  840. return
  843. def data_scatter2d(u, r, first_axis, second_axis):
  844. first_axis = find_index(first_axis)
  845. second_axis = find_index(second_axis)
  846. plt.plot(x[u,first_axis],x[u,second_axis],'o',markerfacecolor=(0,1,0.2),markeredgecolor='k',ms=10,label='Unruptured')
  847. plt.plot(x[r,first_axis],x[r,second_axis],'o',markerfacecolor='red',markeredgecolor='k',ms=10,label='Ruptured')
  848. plt.title('2D Data Scattering')
  849. plt.subplot().set(facecolor='none')
  850. plt.xlabel(attributes[first_axis])
  851. plt.ylabel(attributes[second_axis])
  852. plt.legend(loc=(0.8,0.60))
  853. plt.savefig('plot.pdf', transparent=True, dpi=300)
  854. plt.show()
  855. return
  857. # plt.savefig('plot.png', transparent=True)
  858. data_scatter3d(u, r, 'LSAR', 'EI', 'OVI')
  859. # data_scatter2d(u, r, 'EI', 'DV')
  861. import numpy as np
  862. import matplotlib.pyplot as plt
  863. from sklearn.preprocessing import StandardScaler
  864. from sklearn.decomposition import PCA
  865. from sklearn.manifold import TSNE
  867. plt.rcParams.update({
  868. "font.family": "serif",
  869. "font.serif": ["Times New Roman", "Times", "DejaVu Serif"],
  870. "font.size": 22,
  871. "axes.labelsize": 24,
  872. "xtick.labelsize": 22,
  873. "ytick.labelsize": 22,
  874. "legend.fontsize": 20,
  875. })
  877. COL_U = "#800000" # Maroon
  878. COL_R = "#FF8C00" # Orange
  879. EDGE = (0, 0, 0, 0.18)
  881. def umap_tsne_horizontal(x, y,
  882. umap_neighbors=25,
  883. umap_min_dist=0.15,
  884. tsne_perplexity=30,
  885. random_state=42):
  887. x = np.asarray(x)
  888. y = np.asarray(y)
  890. Xs = StandardScaler().fit_transform(x)
  892. # ---- UMAP (fallback to PCA if umap not installed)
  893. try:
  894. import umap
  895. reducer = umap.UMAP(
  896. n_neighbors=umap_neighbors,
  897. min_dist=umap_min_dist,
  898. random_state=random_state
  899. )
  900. Z_umap = reducer.fit_transform(Xs)
  901. except Exception:
  902. Z_umap = PCA(n_components=2, random_state=random_state).fit_transform(Xs)
  904. # ---- t-SNE
  905. Z_tsne = TSNE(
  906. n_components=2,
  907. perplexity=tsne_perplexity,
  908. learning_rate="auto",
  909. init="pca",
  910. random_state=random_state
  911. ).fit_transform(Xs)
  913. fig, axes = plt.subplots(1, 2, figsize=(16, 7), dpi=300)
  915. # ---------- Panel (a): UMAP ----------
  916. ax = axes[0]
  917. ax.scatter(Z_umap[y==0, 0], Z_umap[y==0, 1],
  918. s=45, c=COL_U, alpha=0.7,
  919. edgecolors=EDGE, linewidths=0.4, label="Unruptured")
  920. ax.scatter(Z_umap[y==1, 0], Z_umap[y==1, 1],
  921. s=45, c=COL_R, alpha=0.7,
  922. edgecolors=EDGE, linewidths=0.4, label="Ruptured", marker='^')
  924. ax.set_xlabel("UMAP 1", labelpad=20)
  925. ax.set_ylabel("UMAP 2", labelpad=20)
  927. ax.spines["top"].set_visible(False)
  928. ax.spines["right"].set_visible(False)
  930. ax.text(-0.15, 1.0, "a)", transform=ax.transAxes,
  931. fontsize=28, fontweight="bold", va="top")
  933. ax.legend(
  934. loc=(0.6,0.8),
  935. ncol=1,
  936. frameon=True,
  937. fancybox=True,
  938. framealpha=0.9,
  939. edgecolor="black")
  940. # ---------- Panel (b): t-SNE ----------
  941. ax = axes[1]
  942. ax.scatter(Z_tsne[y==0, 0], Z_tsne[y==0, 1],
  943. s=45, c=COL_U, alpha=0.7,
  944. edgecolors=EDGE, linewidths=0.4, label="Unruptured")
  945. ax.scatter(Z_tsne[y==1, 0], Z_tsne[y==1, 1],
  946. s=45, c=COL_R, alpha=0.7,
  947. edgecolors=EDGE, linewidths=0.4, label="Ruptured", marker='^')
  949. ax.set_xlabel("t-SNE 1", labelpad=20)
  950. ax.set_ylabel("t-SNE 2")
  952. ax.spines["top"].set_visible(False)
  953. ax.spines["right"].set_visible(False)
  955. ax.text(-0.15, 1.0, "b)", transform=ax.transAxes,
  956. fontsize=28, fontweight="bold", va="top")
  958. # Vertical legend in upper-right of right panel
  961. plt.tight_layout()
  962. plt.savefig("UMAP_tSNE_horizontal_refined.png", dpi=300, bbox_inches="tight")
  963. plt.savefig("UMAP_tSNE_horizontal_refined.pdf", bbox_inches="tight")
  964. plt.show()
  967. umap_tsne_horizontal(x, y)
  969. """#Result Representation"""
  971. from sklearn.metrics import ConfusionMatrixDisplay, classification_report, confusion_matrix, accuracy_score
  972. from sklearn.inspection import permutation_importance
  973. from IPython.display import Image
  974. import matplotlib.pyplot as plt
  975. from matplotlib import rcParams
  976. from io import StringIO
  977. import pydotplus
  978. import datetime
  981. plt.rcParams['font.size'] = 17
  982. font = {'family' : 'serif',
  983. 'serif' : ['Times New Roman'],
  984. 'weight' : 'medium',
  985. 'size' : 24,
  986. 'style' : 'normal'}
  988. class ResultsShow:
  991. plt.rcParams['font.family'] = 'serif'
  992. plt.rcParams['font.serif'] = ['Times New Roman'] + plt.rcParams['font.serif']
  993. plt.rcParams.update({'font.size': 15})
  996. def show(model, x_train, y_train, x_test, y_test):
  997. print('Training Time:', end-start )
  998. print('Train Data Accuracy: \n ', accuracy_score(y_train, model.predict(x_train)))
  999. print('Train Data Confusion Matrix: \n ', confusion_matrix(y_train, model.predict(x_train)))
  1000. print('Test Data Accuracy: \n ', accuracy_score(y_test, model.predict(x_test)))
  1001. print('Test Data Confusion Matrix: \n', confusion_matrix(y_test, model.predict(x_test)))
  1002. return
  1005. def confusion_matrix(name, model, x_test, y_test):
  1006. fig = ConfusionMatrixDisplay(confusion_matrix=confusion_matrix(y_test, model.predict(x_test),\
  1007. labels=model.classes_), display_labels=model.classes_)
  1008. fig.plot()
  1009. fig.figure_.suptitle("Confusion Matrix for " + name + " Classifier", fontweight="bold")
  1010. plt.xlabel('Predicted label', fontweight="bold")
  1011. plt.ylabel('True label', fontweight="bold")
  1012. plt.savefig('ConfusionMatrix.png', transparent=True, bbox_inches='tight')
  1013. plt.show()
  1014. return
  1017. def feature_exraction(model, x_train, df):
  1018. attributes = df.columns.values
  1019. plt.figure(figsize=(16,8))
  1020. imp = model.feature_importances_
  1021. imp_combine = np.zeros(len(x_train[0,:]))
  1022. index, features = list(), list()
  1023. for i in range(len(imp)):
  1024. if imp[i] != 0:
  1025. index.append(i)
  1026. features.append(attributes[i])
  1027. plt.bar(range(len(imp_combine)),imp)
  1028. plt.xticks(index, features, size='small',rotation=90)
  1029. plt.title('Important Parameters Among All Available Features')
  1030. plt.savefig('FeatureExtraction.png', transparent=True)
  1031. plt.show()
  1032. return
  1035. def feature_exraction_nozero(model, x_train, df):
  1036. plt.rcParams.update({'font.size': 16})
  1037. attributes = df.columns.values
  1038. plt.figure(figsize=(16,8))
  1039. imp = model.feature_importances_
  1040. imp_combine = np.zeros(len(x_train[0,:]))
  1041. index, features = list(), list()
  1042. for i in range(len(imp)):
  1043. if imp[i] != 0:
  1044. index.append(i)
  1045. features.append(attributes[i])
  1046. plt.bar(range(len(index)),imp[index])
  1047. plt.xticks(range(len(index)), features, size='small',rotation=0)
  1048. plt.title('Important Parameters Among All Available Features')
  1049. plt.ylabel('Importance of Each Parameter')
  1050. plt.savefig('FeatureExtractionNozero.png', transparent=True, bbox_inches='tight')
  1051. plt.show()
  1052. return
  1055. def feature_permutation(model, x_test, y_test):
  1056. result = permutation_importance(model, x_test, y_test, n_repeats=10, random_state=None)
  1057. sorted_idx = result.importances_mean.argsort()
  1058. plt.figure(figsize=(8,14))
  1059. plt.title('Permutation Importance of Each Feature')
  1060. plt.barh(range(x.shape[1]), result.importances_mean[sorted_idx])
  1061. plt.yticks(range(x.shape[1]), df.columns.values[sorted_idx])
  1062. plt.xlabel('Permutation Importance')
  1063. plt.savefig('FeaturePermutation.png', transparent=True, bbox_inches='tight')
  1064. plt.show()
  1065. return
  1068. def feature_permutation_nozero(model, x_test, y_test):
  1069. plt.rcParams.update({'font.size': 16})
  1070. result = permutation_importance(model, x_test, y_test, n_repeats=10, random_state=None)
  1071. remove_idx = list()
  1072. importance_mean = result.importances_mean
  1073. j = 0
  1074. for i in range(x.shape[1]):
  1075. if abs(result.importances_mean[i]) < 1e-3:
  1076. remove_idx.append(i)
  1077. importance_mean = np.delete(importance_mean, i-j)
  1078. j += 1
  1079. sorted_idx = importance_mean.argsort()
  1080. plt.figure(figsize=(6,10))
  1081. plt.title('Permutation Importance of Each Feature')
  1082. plt.barh(range(len(importance_mean)), importance_mean[sorted_idx])
  1083. plt.yticks(range(len(importance_mean)),
  1084. df.columns.values[sorted_idx])
  1085. plt.xlabel('Permutation Importance')
  1086. plt.savefig('FeaturePermutationNonzero.png', transparent=True, bbox_inches='tight')
  1087. plt.show()
  1088. return
  1091. def tree_depiction(model):
  1092. plt.figure()
  1093. dot_data = StringIO()
  1094. tree.export_graphviz(model, out_file=dot_data,filled=True, rounded=True, special_characters=False, class_names=['U','R'])
  1095. graph = pydotplus.graph_from_dot_data(dot_data.getvalue())
  1096. with open('TreeDepiction.png', 'wb') as f:
  1097. f.write(graph.create_png())
  1098. display(Image(graph.create_png()))
  1099. return
  1102. def models_comparison_supervised(names, results):
  1103. plt.figure(figsize=(8,5))
  1104. x = np.arange(len(names))
  1105. plt.bar(x - 0.12, np.array(results)[:,0], width=0.4, color='#1874CD', label='Train', hatch='o')
  1106. plt.bar(x + 0.12, np.array(results)[:,1], width=0.4, color='#FF8000', label='Test', hatch='/')
  1107. plt.grid( linestyle='--', linewidth=1, axis='y', alpha=0.7)
  1108. # plt.title('Accuracy of Different Supervised Calassification Models')
  1109. plt.ylim(0.5,1)
  1110. plt.ylabel('Accuracy')
  1111. plt.xticks(ticks=x ,labels=names, rotation='vertical')
  1112. plt.legend(loc='best', bbox_to_anchor=(0.75, -0.45))
  1113. # plt.tight_layout()
  1114. plt.savefig('supervised_accuracy.pdf', transparent=True, bbox_inches='tight')
  1115. plt.show()
  1116. return
  1119. def models_comparison_unsupervised(models):
  1120. fig = plt.figure(figsize=(8,5))
  1121. plt.bar(models.keys(), models.values(), width=0.5, color='#1874CD')
  1122. plt.grid( linestyle='--', linewidth=1, axis='y', alpha=0.3)
  1123. # plt.title('Accuracy of Different Unsupervised Calassification Models')
  1124. plt.ylabel('Accuracy')
  1125. plt.xticks(rotation=90)
  1126. plt.savefig('Unsupervised.pdf', transparent=True, bbox_inches='tight')
  1127. plt.show()
  1128. return
  1130. def models_comparison_supervised_recall_precision(names, results, precision, recall):
  1131. fig = plt.figure(figsize=(10,5))
  1132. x = np.arange(len(names))
  1133. plt.bar(x - 0.2, np.array(results)[:,1], width=0.3, color='#1874CD', label='Accuracy', hatch='\\')
  1134. plt.bar(x, np.array(precision), width=0.3, color='#FF8000', label='Percision', hatch='o')
  1135. plt.bar(x + 0.2, np.array(recall), width=0.3, color='#308014', label='Recall', hatch='/')
  1136. plt.grid( linestyle='--', linewidth=1, axis='y', alpha=0.3)
  1137. # plt.title('Comparison of Different Supervised Calassification Models for Test Data and Ruptured Class')
  1138. plt.ylim(0.5,1)
  1139. plt.ylabel('Score')
  1140. plt.xticks(ticks=x ,labels=names, rotation='vertical')
  1141. plt.legend(loc='best', bbox_to_anchor=(0.75, -0.45))
  1142. plt.savefig('Supervised_recall_precision.pdf', transparent=True, bbox_inches='tight')
  1143. plt.show()
  1144. return
  1146. """#Cross Validation"""
  1148. from sklearn.model_selection import RepeatedStratifiedKFold, StratifiedShuffleSplit,\
  1149. cross_val_score, cross_validate, cross_val_predict, StratifiedKFold
  1150. from sklearn.ensemble import VotingClassifier
  1153. def cross_validation(model, x_train, y_train):
  1155. #Find the best corss validation function
  1156. '''
  1157. skf = StratifiedKFold(n_splits=5, shuffle=True)
  1158. scores = cross_val_score(model, x_train, y_train, cv=skf, scoring='accuracy', n_jobs=-1)
  1159. print(scores, scores.mean(), scores.std())
  1161. skf = RepeatedStratifiedKFold(n_splits=5, n_repeats=3)
  1162. scores = cross_val_score(model, x_train, y_train, cv=skf, scoring='accuracy', n_jobs=-1)
  1163. print(scores, scores.mean(), scores.std())
  1164. '''
  1166. split_num = 10
  1167. skf = StratifiedShuffleSplit(n_splits=split_num, test_size = 0.1)
  1168. results = cross_validate(model, x_train, y_train, cv=skf, scoring='accuracy', return_estimator=True, n_jobs=-1)
  1169. names_list = []
  1170. for i in range (split_num):
  1171. names_list.append(str(i+1))
  1172. models_list = results["estimator"]
  1173. models_list = list(zip(names_list, models_list))
  1174. voting_model = VotingClassifier(models_list, voting='soft')
  1175. voting_model.fit(x_train, y_train)
  1177. # cross_validate is more comprehensive than cross_val_score and cross_val_predict
  1178. '''
  1179. scores = cross_val_score(model, x_test, y_test, cv=skf, n_jobs=-1)
  1180. pred = cross_val_predict(model, x_test, y_test, cv=2, n_jobs=-1)
  1181. '''
  1182. return voting_model
  1184. """# Supervised Classification Models"""
  1186. from sklearn.gaussian_process import GaussianProcessClassifier
  1187. from sklearn.ensemble import HistGradientBoostingClassifier
  1188. from sklearn.semi_supervised import SelfTrainingClassifier
  1189. from sklearn.ensemble import GradientBoostingClassifier
  1190. from sklearn.ensemble import RandomForestClassifier
  1191. from sklearn.semi_supervised import LabelSpreading
  1192. from sklearn.neighbors import KNeighborsClassifier
  1193. from sklearn.neural_network import MLPClassifier
  1194. from sklearn.gaussian_process.kernels import RBF
  1195. from sklearn.ensemble import AdaBoostClassifier
  1196. from sklearn.linear_model import SGDClassifier
  1197. from sklearn.ensemble import BaggingClassifier
  1198. from sklearn.naive_bayes import GaussianNB
  1199. from sklearn.svm import SVC
  1200. from sklearn import tree
  1201. from sklearn import svm
  1202. import xgboost as xgb
  1205. class SupervisedModels:
  1208. def xgboost():
  1209. model = xgb.XGBClassifier(booster='gbtree', eta=0.05, gamma=0, max_depth=4, learning_rate=0.05, subsample=0.5, colsample_bytree=0.5, alpha=1, tree_method ='hist')
  1210. return model
  1212. def histogram_gradient_boosting():
  1213. model = HistGradientBoostingClassifier(max_iter=150, learning_rate=0.03, loss='log_loss', max_leaf_nodes=10, min_samples_leaf=20)
  1214. return model
  1217. def gradient_boosting():
  1218. model = GradientBoostingClassifier(n_estimators=60, learning_rate=0.02, max_depth=4, min_samples_split=10)
  1219. return model
  1222. def adaboost():
  1223. model = AdaBoostClassifier(n_estimators=110, learning_rate=0.15)
  1224. return model
  1227. def kneighbor():
  1228. model = KNeighborsClassifier(n_neighbors=17, weights='uniform', p=1, algorithm='auto')
  1229. return model
  1232. def bagging_kneighbor():
  1233. model = BaggingClassifier(estimator=KNeighborsClassifier(n_neighbors=5, weights='uniform', algorithm='auto'), n_estimators = 100, bootstrap=True, max_samples=0.9, max_features=0.7)
  1234. return model
  1237. def bagging_svc():
  1238. model = BaggingClassifier(estimator = svm.SVC(C=80, gamma=0.002, probability=True, class_weight=None , kernel='rbf'), n_estimators = 100, bootstrap=True, max_samples=0.9, max_features=54)
  1239. return model
  1242. def decision_tree():
  1243. model = tree.DecisionTreeClassifier(criterion='gini', splitter='random', max_depth=4, max_features=None, min_samples_split=12, random_state=None)
  1244. return model
  1247. def random_forest():
  1248. model = RandomForestClassifier(n_estimators=80, criterion='log_loss', max_depth=5, min_samples_leaf=5, max_features=0.5, min_samples_split=20, bootstrap=True)
  1249. return model
  1252. def naive_bayes():
  1253. model = GaussianNB(var_smoothing=1e-10)
  1254. return model
  1257. def support_vector_machine():
  1258. model = svm.SVC(C=80, gamma=0.002, probability=True, class_weight=None , kernel='rbf')
  1259. return model
  1262. def stochastic_gradient_descent():
  1263. model = SGDClassifier(loss="log_loss", penalty="l1", max_iter=10000, tol=7e-5, alpha=0.0002, eta0=0.01, learning_rate='adaptive', class_weight=None)
  1264. return model
  1267. def gussian_process():
  1268. model = GaussianProcessClassifier(kernel=1*RBF(2), optimizer='fmin_l_bfgs_b', n_restarts_optimizer=0)
  1269. return model
  1272. def gussian_naive_bayes():
  1273. model = GaussianNB(var_smoothing=1e-10)
  1274. return model
  1277. def selftraining():
  1278. model = SelfTrainingClassifier(HistGradientBoostingClassifier(max_iter=150, learning_rate=0.03, loss='log_loss', max_leaf_nodes=10, min_samples_leaf=20), threshold=0.5)
  1279. return model
  1282. def label_spreading():
  1283. model = LabelSpreading(n_neighbors=5, gamma=2, alpha=0.1)
  1284. return model
  1287. def multi_layer_perceptron():
  1288. model = MLPClassifier(hidden_layer_sizes=(55,18,1), tol=1e-6, batch_size=28, max_iter=2500, activation='identity', solver='sgd', alpha=0.2, learning_rate_init=0.00005, learning_rate='adaptive')
  1289. return model
  1291. """## XGBoost"""
  1293. from sklearn.model_selection import GridSearchCV
  1295. test = xgb.XGBClassifier()
  1296. parameter_space = {
  1297. 'booster': ['dart'],
  1298. 'eta': [0.05],
  1299. 'gamma': [0],
  1300. 'max_depth': [5],
  1301. 'learning_rate': [0.05],
  1302. 'n_estimators': [87],
  1303. 'subsample': [0.5],
  1304. 'colsample_bytree': [0.5],
  1305. 'alpha': [1],
  1306. 'reg_lambda': [1],
  1307. 'reg_alpha': [1]
  1308. }
  1309. search = GridSearchCV(test, parameter_space, n_jobs=-1, cv=3).fit(x, y)
  1310. print('Best parameters found:\n', search.best_params_)
  1312. #Feature Selcetion
  1314. from sklearn.feature_selection import SelectKBest
  1315. from sklearn.feature_selection import f_regression
  1316. from sklearn.feature_selection import RFE
  1318. # selector = RFE(model, n_features_to_select=50)
  1319. # selector.fit(x_train, y_train)
  1320. # x_train_new = selector.transform(x_train)
  1321. # x_test_new = selector.transform(x_test)
  1323. # selector = SelectKBest(score_func=f_regression, k=10)
  1324. # x_train_new = selector.fit_transform(x_train, y_train)
  1325. # x_test_new = selector.transform(x_test)
  1327. import xgboost as xgb
  1329. # model = xgb.XGBClassifier(booster='dart', eta=0.075, gamma=0.2, max_depth=2, learning_rate=0.015, subsample=1, colsample_bytree=1, alpha=0.5, tree_method ='hist')
  1330. model = xgb.XGBClassifier(booster='gbtree', eta=0.05, gamma=0, max_depth=4, learning_rate=0.05, subsample=0.5, colsample_bytree=0.5, alpha=1, tree_method ='hist')
  1331. start = datetime.datetime.now()
  1332. start = datetime.datetime.now()
  1333. model.fit(x_train, y_train)
  1334. # model.fit(x_train_new, y_train)
  1335. end = datetime.datetime.now()
  1337. ResultsShow.show(model, x_train, y_train, x_test, y_test)
  1338. # ResultsShow.show(model, x_train_new, y_train, x_test_new, y_test)
  1339. print('-------------------------------------------------------')
  1341. # Cross Validation Approach
  1342. ResultsShow.show(cross_validation(model, x_train, y_train), x_train, y_train, x_test, y_test)
  1343. # ResultsShow.feature_exraction(model, x_train,df)
  1344. # print('-------------------------------------------------------------------------------------')
  1345. # ResultsShow.feature_permutation(model, x_test, y_test,)
  1349. # attributes = df.columns.values
  1350. # plt.figure(figsize=(8,16))
  1351. # imp = model.feature_importances_
  1352. # imp_combine = np.zeros(len(x_train[0,:]))
  1353. # index, features = list(), list()
  1354. # for i in range(len(imp)):
  1355. # if imp[i] != 0:
  1356. # index.append(i)
  1357. # features.append(attributes[i])
  1358. # plt.barh(range(len(imp_combine)),imp)
  1359. # plt.yticks(index, features, size='small',rotation=0)
  1360. # plt.xlabel('Importance of Each Parameter')
  1361. # plt.title('Important Parameters Among All Available Features')
  1362. # plt.savefig('FeatureExtraction.png', transparent=True)
  1363. # plt.show()
  1365. """## Support Vector Machine_SVC"""
  1367. from sklearn.model_selection import GridSearchCV
  1369. test = svm.SVC()
  1370. # test = svm.NuSVC()
  1371. parameter_space = {
  1372. 'C': [80, 90, 70],
  1373. 'gamma': [0.002],
  1374. 'kernel': ['rbf', 'sigmoid']
  1376. }
  1377. search = GridSearchCV(test, parameter_space, n_jobs=-1, cv=3).fit(x, y)
  1378. print('Best parameters found:\n', search.best_params_)
  1380. from sklearn.metrics import PrecisionRecallDisplay
  1381. from sklearn.metrics import RocCurveDisplay
  1382. from sklearn import svm
  1384. model = svm.SVC(C=80, gamma=0.002, probability=True, class_weight=None , kernel='rbf')
  1385. # model = svm.NuSVC(nu=0.5, gamma=0.002, probability=True, class_weight=None, kernel='rbf')
  1386. start = datetime.datetime.now()
  1387. model.fit(x_train, y_train)
  1388. end = datetime.datetime.now()
  1390. ResultsShow.show(model, x_train, y_train, x_test, y_test)
  1391. print('-------------------------------------------------------')
  1392. ResultsShow.show(cross_validation(model, x_train, y_train), x_train, y_train, x_test, y_test)
  1394. ROC = RocCurveDisplay.from_estimator(model, x_test, y_test)
  1395. PrecisionRecallDisplay.from_estimator(model, x_test, y_test)
  1396. print('-------------------------------------------------------------------------------------')
  1397. ResultsShow.feature_permutation(model, x_test, y_test)
  1399. result = permutation_importance(model, x_test, y_test, n_repeats=10, random_state=None)
  1400. sorted_idx = result.importances_mean.argsort()
  1401. v = result.importances_mean[sorted_idx]
  1402. n = df.columns.values[sorted_idx]
  1403. print(n.shape)
  1404. print(n)
  1405. print(v)
  1407. """## Stochastic Gradient Descent"""
  1409. from sklearn.model_selection import GridSearchCV
  1411. test = SGDClassifier()
  1412. parameter_space = {
  1413. 'loss': ['log_loss'],
  1414. 'penalty': ['l1','l2'],
  1415. 'learning_rate': ['adaptive'],
  1416. 'eta0': [0.01,0.009,0.011],
  1417. 'alpha': [0.0002,0.002],
  1418. 'class_weight': [None]
  1419. }
  1420. search = GridSearchCV(test, parameter_space, n_jobs=-1, cv=3).fit(x, y)
  1421. print('Best parameters found:\n', search.best_params_)
  1423. from sklearn.linear_model import SGDClassifier
  1425. model = SGDClassifier(loss="log_loss", penalty="l1", max_iter=10000, tol=7e-5, alpha=0.0002, eta0=0.01, learning_rate='adaptive', class_weight={0: 1, 1: 1.8})
  1426. start = datetime.datetime.now()
  1427. model.fit(x_train, y_train)
  1428. end = datetime.datetime.now()
  1430. ResultsShow.show(model, x_train, y_train, x_test, y_test)
  1431. print('-------------------------------------------------------')
  1432. ResultsShow.show(cross_validation(model, x_train, y_train), x_train, y_train, x_test, y_test)
  1433. # RocCurveDisplay.from_estimator(model, x_test, y_test)
  1434. print('-------------------------------------------------------------------------------------')
  1435. ResultsShow.feature_permutation(model, x_test, y_test)
  1437. result = permutation_importance(model, x_test, y_test, n_repeats=10, random_state=None)
  1438. sorted_idx = result.importances_mean.argsort()
  1439. v = result.importances_mean[sorted_idx]
  1440. n = df.columns.values[sorted_idx]
  1441. print(n)
  1442. print(v)
  1444. """##KNeighbor"""
  1446. from sklearn.model_selection import GridSearchCV
  1448. test = KNeighborsClassifier()
  1449. parameter_space = {
  1450. 'n_neighbors': [11,13,15,17,19],
  1451. 'weights': ['uniform','distance'],
  1452. 'p': [1],
  1453. 'algorithm': ['auto', 'ball_tree', 'kd_tree', 'brute']
  1454. }
  1455. search = GridSearchCV(test, parameter_space, n_jobs=-1, cv=3).fit(x, y)
  1456. print('Best parameters found:\n', search.best_params_)
  1458. from sklearn.neighbors import KNeighborsClassifier
  1460. model = KNeighborsClassifier(n_neighbors=17, weights='uniform', p=1, algorithm='auto')
  1461. start = datetime.datetime.now()
  1462. model.fit(x_train, y_train)
  1463. end = datetime.datetime.now()
  1465. ResultsShow.show(model, x_train, y_train, x_test, y_test)
  1466. print('-------------------------------------------------------')
  1468. # Cross Validation Approach
  1469. ResultsShow.show(cross_validation(model, x_train, y_train), x_train, y_train, x_test, y_test)
  1470. print('-------------------------------------------------------------------------------------')
  1471. ResultsShow.feature_permutation(model, x_test, y_test)
  1473. """## Gaussian Process"""
  1475. from sklearn.model_selection import GridSearchCV
  1477. test = GaussianProcessClassifier()
  1478. parameter_space = {
  1479. 'kernel': [RBF(2),RBF(5)],
  1480. 'optimizer': ['fmin_l_bfgs_b'],
  1481. 'n_restarts_optimizer': [0]
  1482. }
  1483. search = GridSearchCV(test, parameter_space, n_jobs=-1, cv=3).fit(x, y)
  1484. print('Best parameters found:\n', search.best_params_)
  1486. from sklearn.gaussian_process import GaussianProcessClassifier
  1487. from sklearn.gaussian_process.kernels import RBF
  1489. model = GaussianProcessClassifier(kernel=1*RBF(2), optimizer='fmin_l_bfgs_b', n_restarts_optimizer=0)
  1490. start = datetime.datetime.now()
  1491. model.fit(x_train, y_train)
  1492. end = datetime.datetime.now()
  1494. ResultsShow.show(model, x_train, y_train, x_test, y_test)
  1495. print('-------------------------------------------------------')
  1496. ResultsShow.show(cross_validation(model, x_train, y_train), x_train, y_train, x_test, y_test)
  1497. print('-------------------------------------------------------------------------------------')
  1498. ResultsShow.feature_permutation(model, x_test, y_test)
  1500. """## Gaussian Naive Bayes"""
  1502. from sklearn.model_selection import GridSearchCV
  1504. test = GaussianNB()
  1505. parameter_space = {
  1506. 'var_smoothing': [1e-9,1e-10,5e-11]
  1507. }
  1508. search = GridSearchCV(test, parameter_space, n_jobs=-1, cv=3).fit(x, y)
  1509. print('Best parameters found:\n', search.best_params_)
  1511. from sklearn.naive_bayes import GaussianNB
  1513. model = GaussianNB(var_smoothing=1e-10)
  1514. start = datetime.datetime.now()
  1515. model.fit(x_train, y_train)
  1516. end = datetime.datetime.now()
  1518. ResultsShow.show(model, x_train, y_train, x_test, y_test)
  1519. print('-------------------------------------------------------')
  1520. ResultsShow.show(cross_validation(model, x_train, y_train), x_train, y_train, x_test, y_test)
  1521. print('-------------------------------------------------------------------------------------')
  1522. ResultsShow.feature_permutation(model, x_test, y_test)
  1524. """##Decision Tree"""
  1526. from sklearn.model_selection import GridSearchCV
  1528. test = tree.DecisionTreeClassifier()
  1529. parameter_space = {
  1530. 'criterion': ['log_loss','entropy', 'gini'],
  1531. 'splitter': ['random', 'best'],
  1532. 'max_depth': [3,5,4],
  1533. 'max_features': [0.7,0.9,None],
  1534. 'min_samples_split': [10,12,15],
  1535. 'min_samples_leaf': [5]
  1536. }
  1537. search = GridSearchCV(test, parameter_space, n_jobs=-1, cv=3).fit(x, y)
  1538. print('Best parameters found:\n', search.best_params_)
  1540. from sklearn import tree
  1542. model = tree.DecisionTreeClassifier(criterion='gini', splitter='random', max_depth=4, max_features=None, min_samples_split=12, random_state=None)
  1544. start = datetime.datetime.now()
  1545. model.fit(x_train, y_train)
  1546. end = datetime.datetime.now()
  1548. ResultsShow.show(model, x_train, y_train, x_test, y_test)
  1549. print('-------------------------------------------------------')
  1551. # Cross Validation Approach
  1552. CrossVal = cross_validation(model, x_train, y_train)
  1553. ResultsShow.show(CrossVal, x_train, y_train, x_test, y_test)
  1554. # print('-------------------------------------------------------')
  1555. # ResultsShow.feature_exraction(model, x_train, df)
  1556. # print('-------------------------------------------------------')
  1557. # ResultsShow.confusion_matrix('Decision Tree', CrossVal, x_test, y_test)
  1558. # print('-------------------------------------------------------')
  1559. # ResultsShow.tree_depiction(model)
  1560. # print('-------------------------------------------------------------------------------------')
  1562. from sklearn.metrics import PrecisionRecallDisplay
  1563. from sklearn.metrics import RocCurveDisplay
  1564. ROC = RocCurveDisplay.from_estimator(CrossVal, x_test, y_test)
  1565. plt.title('Receiver Operating Characteristic (ROC) Curve')
  1566. plt.savefig('ROC Curve.png', transparent=True, bbox_inches='tight')
  1567. PrecisionRecallDisplay.from_estimator(CrossVal, x_test, y_test)
  1568. plt.title('Precision_Recall Curve')
  1569. plt.savefig('Precision_Recall Curve.png', transparent=True, bbox_inches='tight')
  1571. """## Ensemble Methods
  1573. ###Bagging_SVC
  1574. """
  1576. from sklearn.model_selection import GridSearchCV
  1578. test = BaggingClassifier(estimator = svm.SVC(C=80, gamma=0.002, probability=True, class_weight=None , kernel='rbf'))
  1579. parameter_space = {
  1580. 'n_estimators': [100],
  1581. 'max_samples': [0.9],
  1582. 'max_features': [0.9,1],
  1583. 'bootstrap': [True]
  1584. }
  1585. search = GridSearchCV(test, parameter_space, n_jobs=-1, cv=3).fit(x, y)
  1586. print('Best parameters found:\n', search.best_params_)
  1588. from sklearn.ensemble import BaggingClassifier
  1589. from sklearn.svm import SVC
  1591. model = BaggingClassifier(estimator = svm.SVC(C=80, gamma=0.002, probability=True, class_weight=None , kernel='rbf'), n_estimators = 100, bootstrap=True, max_samples=0.9, max_features=54)
  1592. start = datetime.datetime.now()
  1593. model.fit(x_train, y_train)
  1594. end = datetime.datetime.now()
  1596. ResultsShow.show(model, x_train, y_train, x_test, y_test)
  1597. print('-------------------------------------------------------')
  1599. # Cross Validation Approach
  1600. ResultsShow.show(cross_validation(model, x_train, y_train), x_train, y_train, x_test, y_test)
  1601. print('-------------------------------------------------------------------------------------')
  1602. ResultsShow.feature_permutation(model, x_test, y_test)
  1604. """###Bagging_Kneighbor"""
  1606. from sklearn.model_selection import GridSearchCV
  1608. test = BaggingClassifier(estimator=KNeighborsClassifier(n_neighbors=5, weights='uniform', algorithm='auto'))
  1609. parameter_space = {
  1610. 'n_estimators': [90,100],
  1611. 'max_samples': [0.9,0.7],
  1612. 'max_features': [0.7],
  1613. 'bootstrap': [True],
  1614. }
  1615. search = GridSearchCV(test, parameter_space, n_jobs=-1, cv=3).fit(x, y)
  1616. print('Best parameters found:\n', search.best_params_)
  1618. from sklearn.ensemble import BaggingClassifier
  1619. from sklearn.neighbors import KNeighborsClassifier
  1621. model = BaggingClassifier(estimator=KNeighborsClassifier(n_neighbors=5, weights='uniform', algorithm='auto'), n_estimators = 100, bootstrap=True, max_samples=0.9, max_features=0.7)
  1622. start = datetime.datetime.now()
  1623. model.fit(x_train, y_train)
  1624. end = datetime.datetime.now()
  1626. ResultsShow.show(model, x_train, y_train, x_test, y_test)
  1627. print('-------------------------------------------------------')
  1629. # Cross Validation Approach
  1630. ResultsShow.show(cross_validation(model, x_train, y_train), x_train, y_train, x_test, y_test)
  1631. print('-------------------------------------------------------------------------------------')
  1632. ResultsShow.feature_permutation(model, x_test, y_test)
  1634. """###Random Forest"""
  1636. from sklearn.model_selection import GridSearchCV
  1638. test = RandomForestClassifier()
  1639. parameter_space = {
  1640. 'criterion': ['gini', 'log_loss'],
  1641. 'n_estimators': [60,70],
  1642. 'max_depth': [5,4],
  1643. 'max_features': [0.9,0.5],
  1644. 'min_samples_split': [20],
  1645. 'min_samples_leaf': [3],
  1646. 'bootstrap': [True]
  1648. }
  1649. search = GridSearchCV(test, parameter_space, n_jobs=-1, cv=3).fit(x, y)
  1650. print('Best parameters found:\n', search.best_params_)
  1652. from sklearn.ensemble import RandomForestClassifier
  1654. model = RandomForestClassifier(n_estimators=80, criterion='log_loss', max_depth=5, min_samples_leaf=5, max_features=0.5, min_samples_split=20, bootstrap=True)
  1655. start = datetime.datetime.now()
  1656. model.fit(x_train, y_train)
  1657. end = datetime.datetime.now()
  1659. ResultsShow.show(model, x_train, y_train, x_test, y_test)
  1660. print('-------------------------------------------------------')
  1662. # Cross Validation Approach
  1663. ResultsShow.show(cross_validation(model, x_train, y_train), x_train, y_train, x_test, y_test)
  1664. ResultsShow.feature_exraction(model, x_train, df)
  1665. print('-------------------------------------------------------------------------------------')
  1666. ResultsShow.feature_permutation(model, x_test, y_test)
  1668. """###Adaboost"""
  1670. from sklearn.model_selection import GridSearchCV
  1672. test = AdaBoostClassifier()
  1673. parameter_space = {
  1674. 'n_estimators': [110,70,90],
  1675. 'learning_rate': [0.1,0.15,0.05]
  1677. }
  1678. search = GridSearchCV(test, parameter_space, n_jobs=-1, cv=3).fit(x, y)
  1679. print('Best parameters found:\n', search.best_params_)
  1681. from sklearn.ensemble import AdaBoostClassifier
  1683. model = AdaBoostClassifier(n_estimators=110, learning_rate=0.15)
  1684. start = datetime.datetime.now()
  1685. model.fit(x_train, y_train)
  1686. end = datetime.datetime.now()
  1688. ResultsShow.show(model, x_train, y_train, x_test, y_test)
  1689. print('-------------------------------------------------------')
  1691. # Cross Validation Approach
  1692. ResultsShow.show(cross_validation(model, x_train, y_train), x_train, y_train, x_test, y_test)
  1693. ResultsShow.feature_exraction(model, x_train,df)
  1694. print('-------------------------------------------------------------------------------------')
  1695. ResultsShow.feature_permutation(model, x_test, y_test)
  1697. """###Gradient Tree Boosting"""
  1699. from sklearn.model_selection import GridSearchCV
  1701. test = GradientBoostingClassifier()
  1702. parameter_space = {
  1703. 'n_estimators': [30,60,90],
  1704. 'learning_rate': [0.1,0.05,0.15],
  1705. 'max_depth': [3,4,2],
  1706. 'min_samples_split': [5,10]
  1708. }
  1709. search = GridSearchCV(test, parameter_space, n_jobs=-1, cv=3).fit(x, y)
  1710. print('Best parameters found:\n', search.best_params_)
  1712. from sklearn.ensemble import GradientBoostingClassifier
  1714. model = GradientBoostingClassifier(n_estimators=60, learning_rate=0.02, max_depth=4, min_samples_split=10)
  1715. start = datetime.datetime.now()
  1716. model.fit(x_train, y_train)
  1717. end = datetime.datetime.now()
  1719. ResultsShow.show( model, x_train, y_train, x_test, y_test)
  1720. print('-------------------------------------------------------')
  1722. # Cross Validation Approach
  1723. ResultsShow.show(cross_validation(model, x_train, y_train), x_train, y_train, x_test, y_test)
  1724. ResultsShow.feature_exraction(model, x_train,df)
  1725. print('-------------------------------------------------------------------------------------')
  1726. ResultsShow.feature_permutation(model, x_test, y_test)
  1728. """### Histogram-Based Gradient Boosting"""
  1730. from sklearn.model_selection import GridSearchCV
  1732. test = HistGradientBoostingClassifier()
  1733. parameter_space = {
  1734. 'min_samples_leaf': [15,20,25],
  1735. 'learning_rate': [0.03,0.05]
  1737. }
  1738. search = GridSearchCV(test, parameter_space, n_jobs=-1, cv=3).fit(x, y)
  1739. print('Best parameters found:\n', search.best_params_)
  1741. from sklearn.ensemble import HistGradientBoostingClassifier
  1743. model = HistGradientBoostingClassifier(max_iter=150, learning_rate=0.03, loss='log_loss', max_leaf_nodes=10, min_samples_leaf=20)
  1744. start = datetime.datetime.now()
  1745. model.fit(x_train, y_train)
  1746. end = datetime.datetime.now()
  1748. ResultsShow.show(model, x_train, y_train, x_test, y_test)
  1749. print('-------------------------------------------------------')
  1751. # Cross Validation Approach
  1752. ResultsShow.show(cross_validation(model, x_train, y_train), x_train, y_train, x_test, y_test)
  1753. print('-------------------------------------------------------------------------------------')
  1754. ResultsShow.feature_permutation(model, x_test, y_test)
  1756. """# Semi_Supervised Classification Models
  1758. ## Semi-Supervised Self Training
  1759. """
  1761. from sklearn.semi_supervised import SelfTrainingClassifier
  1763. b = HistGradientBoostingClassifier(max_iter=150, learning_rate=0.03, loss='log_loss', max_leaf_nodes=10, min_samples_leaf=20)
  1764. a = SGDClassifier(loss="log_loss", penalty="l1", max_iter=10000, tol=7e-5, alpha=0.0002, eta0=0.0005, learning_rate='adaptive', class_weight=None)
  1765. model = SelfTrainingClassifier(b, threshold=0.5)
  1766. start = datetime.datetime.now()
  1767. model.fit(x_train, y_train)
  1768. end = datetime.datetime.now()
  1770. ResultsShow.show(model, x_train, y_train, x_test, y_test)
  1771. print('-------------------------------------------------------')
  1772. # ResultsShow.show(cross_validation(model, x_train, y_train), x_train, y_train, x_test, y_test)
  1773. # print('-------------------------------------------------------------------------------------')
  1774. # ResultsShow.feature_permutation(model, x_test, y_test)
  1776. """## Semi-Supervised Label Spreading"""
  1778. from sklearn.semi_supervised import LabelSpreading
  1779. import warnings
  1781. warnings.filterwarnings("ignore")
  1783. model = LabelSpreading(n_neighbors=5, gamma=2, alpha=0.1)
  1784. start = datetime.datetime.now()
  1785. model.fit(x_train, y_train)
  1786. end = datetime.datetime.now()
  1788. ResultsShow.show(model, x_train, y_train, x_test, y_test)
  1789. # print('-------------------------------------------------------')
  1790. # ResultsShow.show(cross_validation(model, x_train, y_train), x_train, y_train, x_test, y_test)
  1791. # print('-------------------------------------------------------------------------------------')
  1792. # ResultsShow.feature_permutation(model, x_test, y_test)
  1794. """#Neural Networks
  1796. ##Multi Layer Perceptron
  1797. """
  1799. from sklearn.neural_network import MLPClassifier
  1800. from sklearn.model_selection import GridSearchCV
  1801. from sklearn.model_selection import RandomizedSearchCV
  1803. mlp = MLPClassifier(max_iter=100)
  1804. parameter_space = {
  1805. 'hidden_layer_sizes': [(54,35,1), (54,35,17,1)],
  1806. 'activation': ['tanh'],
  1807. 'solver': ['sgd'],
  1808. 'alpha': [0.02],
  1809. 'learning_rate': ['adaptive'],
  1810. 'learning_rate_init':[0.001]
  1811. }
  1812. clf = GridSearchCV(mlp, parameter_space, n_jobs=-1, cv=3)
  1813. # clf = RandomizedSearchCV(mlp, parameter_space, n_jobs=-1, cv=3)
  1814. clf.fit(x, y)
  1815. print('Best parameters found:\n', clf.best_params_)
  1816. # means = clf.cv_results_['mean_test_score']
  1817. # stds = clf.cv_results_['std_test_score']
  1818. # for mean, std, params in zip(means, stds, clf.cv_results_['params']):
  1819. # print("%0.3f (+/-%0.03f) for %r" % (mean, std * 2, params))
  1821. from sklearn.neural_network import MLPClassifier
  1823. model = MLPClassifier(hidden_layer_sizes=(54,18,1), tol=1e-6, batch_size=28, max_iter=2500, activation='identity', solver='sgd', alpha=0.2, learning_rate_init=0.00005, learning_rate='adaptive')
  1824. start = datetime.datetime.now()
  1825. model.fit(x_train, y_train)
  1826. end = datetime.datetime.now()
  1828. ResultsShow.show(model, x_train, y_train, x_test, y_test)
  1829. # print('-------------------------------------------------------')
  1830. # ResultsShow.show(cross_validation(model, x_train, y_train), x_train, y_train, x_test, y_test)
  1831. # print('-------------------------------------------------------------------------------------')
  1832. ResultsShow.feature_permutation(model, x_test, y_test)
  1834. result = permutation_importance(model, x_test, y_test, n_repeats=10, random_state=None)
  1835. sorted_idx = result.importances_mean.argsort()
  1836. v = result.importances_mean[sorted_idx]
  1837. n = df.columns.values[sorted_idx]
  1838. print(n)
  1839. print(v)
  1841. from graphviz import Digraph
  1843. # Create a new Graphviz object
  1844. dot = Digraph(format='png')
  1846. # Set the rank direction to left to right
  1847. dot.attr(rankdir='LR')
  1849. # Set the background color to transparent
  1850. dot.attr(bgcolor='transparent')
  1852. # Add the first 3 and last 3 nodes to the input layer
  1853. dot.node('x1', 'x1', fillcolor='#6495ED', style='filled')
  1854. dot.node('x2', 'x2', fillcolor='#6495ED', style='filled')
  1855. dot.node('x3', 'x3', fillcolor='#6495ED', style='filled')
  1856. dot.node('x4', '...', fillcolor='#FFFFFF')
  1857. dot.node('x5', '...', fillcolor='#FFFFFF')
  1858. dot.node('x6', 'x54', fillcolor='#6495ED', style='filled')
  1860. # Add the first 3 and last 3 nodes to the hidden layer
  1861. dot.node('h1', 'h1', fillcolor='#6495ED', style='filled')
  1862. dot.node('h2', 'h2', fillcolor='#6495ED', style='filled')
  1863. dot.node('h3', 'h3', fillcolor='#6495ED', style='filled')
  1864. dot.node('h4', '...', fillcolor='#FFFFFF')
  1865. dot.node('h5', '...', fillcolor='#FFFFFF')
  1866. dot.node('h6', 'h18', fillcolor='#6495ED', style='filled')
  1868. # Add all nodes to the output layer
  1869. dot.node('y', 'y', fillcolor='#6495ED', style='filled')
  1871. # Add edges from input layer to hidden layer
  1872. for i in range(1, 4):
  1873. for j in range(1, 4):
  1874. dot.edge(f'x{i}', f'h{j}')
  1876. for i in range(4, 7):
  1877. for j in range(4, 7):
  1878. dot.edge(f'x{i}', f'h{j}')
  1880. # Add edges from hidden layer to output layer
  1881. for i in range(1, 7):
  1882. dot.edge(f'h{i}', 'y')
  1884. # Render the graph to a PNG file
  1885. dot.render('mlp_network_subset', view=False)
  1887. """## Keras"""
  1889. from keras.models import Sequential
  1890. from keras.layers import Conv1D, MaxPooling1D, Flatten, Dense, Dropout, BatchNormalization
  1891. from keras.utils import plot_model
  1892. from keras.callbacks import EarlyStopping
  1893. from keras import regularizers
  1894. from keras import optimizers
  1897. model = Sequential()
  1898. model.add(Conv1D(32, 3, activation='relu', input_shape=(54, 1)))
  1899. model.add(MaxPooling1D(2))
  1900. model.add(Conv1D(64, 3, activation='relu'))
  1901. model.add(MaxPooling1D(2))
  1902. model.add(Flatten())
  1903. model.add(Dense(128, activation='relu'))
  1904. model.add(BatchNormalization())
  1905. model.add(Dropout(0.2))
  1906. model.add(Dense(64, activation='relu'))
  1907. model.add(BatchNormalization())
  1908. model.add(Dropout(0.2))
  1909. # model.add(Dense(32, activation='relu'))
  1910. # model.add(BatchNormalization())
  1911. # model.add(Dropout(0.5))
  1912. model.add(Dense(1, activation='sigmoid'))
  1914. X_train, X_val, Y_train, Y_val = train_test_split(x_train, y_train, test_size=0.12, random_state=8)
  1916. early_stopping = EarlyStopping(monitor='val_loss', patience=30)
  1917. model.compile(optimizer=optimizers.Adamax(learning_rate=0.003), loss='binary_crossentropy', metrics=['accuracy'])
  1918. # model.compile(optimizer=optimizers.SGD(learning_rate=0.005), loss='binary_crossentropy', metrics=['accuracy']) # too slow
  1919. # model.compile(optimizer=optimizers.Nadam(learning_rate=0.005), loss='mae', metrics=['accuracy']) # Fast
  1920. hist = model.fit(X_train, Y_train, validation_data=(X_val,Y_val), epochs=200, callbacks=[early_stopping], batch_size=16)
  1921. # YPredict = np.argmax(model.predict(x_test),axis=1)
  1924. # accuracy = np.sum(YPredict == y_test)/len(y_test)
  1925. # print('--------------------------------')
  1926. # print('Accuracy = {}'.format(accuracy)) # cannot be used with sigmoid
  1927. # print('--------------------------------')
  1929. history = hist.history
  1930. loss = history['loss']
  1931. val_loss = history['val_loss']
  1932. acc = history['accuracy']
  1933. val_acc = history['val_accuracy']
  1934. plt.xlabel('Epochs')
  1935. plt.ylabel('Loss')
  1936. plt.plot(loss)
  1937. plt.plot(val_loss)
  1938. plt.legend(['Loss','Vallidation_Loss'])
  1939. plt.savefig('Loss.png', transparent=True, bbox_inches='tight')
  1940. plt.figure()
  1941. plt.xlabel('Epochs')
  1942. plt.ylabel('Accuracy')
  1943. plt.plot(acc)
  1944. plt.plot(val_acc)
  1945. plt.legend(['Accuracy','Validation_Accuracy'])
  1946. plt.savefig('Accuracy.png', transparent=True, bbox_inches='tight')
  1947. plt.show()
  1949. # import tensorflow as tf
  1950. # confusion = tf.math.confusion_matrix(
  1951. # labels = np.argmax(y_test),
  1952. # predictions = np.argmax(YPredict),
  1953. # num_classes=2)
  1954. # print(confusion)
  1956. test_loss, test_accuracy = model.evaluate(x_test,y_test)
  1957. print(test_loss)
  1958. print(test_accuracy)
  1960. model.predict(standardization.transform(x[382].reshape(1,-1))) * 100
  1961. # model.predict(x_test[9].reshape(1,-1))
  1962. # x[382]
  1964. # Create a new Graphviz object
  1965. dot = Digraph(format='png')
  1967. # Set the rank direction to left to right
  1968. dot.attr(rankdir='LR')
  1970. # Set the background color to white
  1971. dot.attr(bgcolor='white')
  1973. # Set the font color and size of the node labels
  1974. dot.attr('node', fontcolor='#333333', fontsize='12')
  1976. # Add nodes for each layer
  1977. for i, layer in enumerate(model.layers):
  1978. name = layer.name
  1979. label = f'{name}\n{layer.output_shape}'
  1980. dot.node(name, label)
  1982. if i > 0:
  1983. prev_layer = model.layers[i-1]
  1984. for node in prev_layer.output_shape[1:]:
  1985. dot.edge(f'{prev_layer.name}-{node}', f'{name}-0')
  1987. # Render the graph to a PNG file
  1988. dot.render('sequential_model', view=False)
  1990. """# Unsupervised Classification Models"""
  1992. from sklearn.cluster import AgglomerativeClustering
  1993. from sklearn.cluster import AffinityPropagation
  1994. from sklearn.cluster import SpectralClustering
  1995. from sklearn.neighbors import NearestNeighbors
  1996. from sklearn.mixture import GaussianMixture
  1997. from sklearn.cluster import MeanShift
  1998. from sklearn.cluster import OPTICS
  1999. from sklearn.cluster import DBSCAN
  2000. from sklearn.cluster import KMeans
  2001. from sklearn.cluster import Birch
  2002. import warnings
  2005. class UnsupervisedModels:
  2008. def nearest_neighbors():
  2009. model = NearestNeighbors(n_neighbors=9, radius=0.1, algorithm='auto')
  2010. return model
  2013. def k_means():
  2014. model = KMeans(n_clusters=2)
  2015. return model
  2018. def agglomerative_clustering():
  2019. model = AgglomerativeClustering(n_clusters=2)
  2020. return model
  2023. def birch():
  2024. model = Birch(threshold=100000, n_clusters=2)
  2025. return model
  2028. def mean_shift():
  2029. model = MeanShift()
  2030. return model
  2033. def affinity_propagation():
  2034. model = AffinityPropagation(damping=0.9999)
  2035. return model
  2038. def optics():
  2039. model = OPTICS(eps=1, min_samples=50)
  2040. return model
  2043. def spectral_clustering():
  2044. model = SpectralClustering(n_clusters=2)
  2045. return model
  2048. def gaussian_mixture():
  2049. model = GaussianMixture(n_components=2)
  2050. return model
  2053. def dbscan():
  2054. model = DBSCAN(eps=500000, min_samples=3)
  2055. return model
  2057. """## Unsupervised Nearest Neighbors"""
  2059. from sklearn.neighbors import NearestNeighbors
  2061. model = NearestNeighbors(n_neighbors=9, radius=0.1, algorithm='auto')
  2062. start = datetime.datetime.now()
  2063. model.fit(x_train)
  2064. end = datetime.datetime.now()
  2066. distances1, indices1 = model.kneighbors(x_train)
  2067. y_pred1 = []
  2068. for i in range(len(y_train)):
  2069. y_pred1.append(y_train[indices1[i,1]])
  2071. print(accuracy_score(y_train, y_pred1))
  2073. """## Unsupervised K Means"""
  2075. from sklearn.cluster import KMeans
  2077. model = KMeans(n_clusters=2)
  2078. model.fit(x_train)
  2079. y_pred = model.labels_
  2081. print(accuracy_score(y_train, y_pred))
  2083. """## Unsupervised Hierarchical Clustering"""
  2085. from scipy.cluster.hierarchy import dendrogram, linkage
  2086. from sklearn.cluster import AgglomerativeClustering
  2088. Z = linkage(x_train, 'ward') # other possible options: single, complete, average, ward
  2089. dendrogram(Z)
  2090. model = AgglomerativeClustering(n_clusters=2)
  2091. model.fit(x_train)
  2092. y_pred = model.labels_
  2093. print(accuracy_score(y_train, y_pred))
  2095. """## Unsupervised Mean Shift"""
  2097. from sklearn.cluster import MeanShift
  2099. model = MeanShift()
  2100. model.fit(x_train)
  2101. y_pred = model.fit_predict(x_train)
  2102. print(accuracy_score(y_train, y_pred))
  2104. """## Unsupervised Affinity Propagation"""
  2106. from sklearn.cluster import AffinityPropagation
  2108. model = AffinityPropagation(damping=0.99999)
  2109. model.fit(x_train)
  2110. y_pred = model.predict(x_train)
  2111. print(accuracy_score(y_train, y_pred))
  2113. """## Unsupervised Birch"""
  2115. from sklearn.cluster import Birch
  2117. model = Birch(threshold=100000, n_clusters=2)
  2118. model.fit(x_train)
  2119. y_pred = model.predict(x_train)
  2120. print(accuracy_score(y_train, y_pred))
  2122. """## Unsupevised Optics"""
  2124. from sklearn.cluster import OPTICS
  2126. model = OPTICS(eps=1, min_samples=50)
  2127. model.fit(x_train)
  2128. y_pred = model.fit_predict(x_train)
  2129. print(accuracy_score(y_train, y_pred))
  2131. """## Unsupervised Spectral"""
  2133. from sklearn.cluster import SpectralClustering
  2134. c
  2136. model = SpectralClustering(n_clusters=2)
  2137. model.fit(x_train)
  2138. y_pred = model.fit_predict(x_train)
  2139. print(accuracy_score(y_train, y_pred))
  2141. """## Unsupervised Gaussian Mixture"""
  2143. from sklearn.mixture import GaussianMixture
  2145. model = GaussianMixture(n_components=2)
  2146. model.fit(x_train)
  2147. y_pred = model.predict(x_train)
  2148. print(accuracy_score(y_train, y_pred))
  2150. """## Unsupervised DBSCAN"""
  2152. from sklearn.cluster import DBSCAN
  2154. model = DBSCAN(eps=500000, min_samples=3)
  2155. model.fit(x)
  2156. y_pred = model.labels_
  2157. print(accuracy_score(y, y_pred))
  2159. """# Models Comparison
  2161. ## Unsupervised Models Comparison
  2162. """
  2164. warnings.filterwarnings("ignore")
  2165. models = dict()
  2166. distances1, indices1 = UnsupervisedModels.nearest_neighbors().fit(x).kneighbors(x)
  2167. y_pred = []
  2168. for i in range(len(y)):
  2169. y_pred.append(y[indices1[i,1]])
  2170. models['Agglomerative Clustering'] = accuracy_score(y, UnsupervisedModels.agglomerative_clustering().fit(x).labels_)
  2171. models['Affinity Propagation'] = accuracy_score(y, UnsupervisedModels.affinity_propagation().fit(x).predict(x))
  2172. models['Spectral Clustering'] = accuracy_score(y, UnsupervisedModels.spectral_clustering().fit(x).fit_predict(x))
  2173. models['Nearest Neighbors'] = accuracy_score(y, y_pred)
  2174. models['Gaussian Mixture'] = accuracy_score(y, UnsupervisedModels.gaussian_mixture().fit(x).predict(x))
  2175. models['Mean Shift'] = accuracy_score(y, UnsupervisedModels.mean_shift().fit(x).fit_predict(x))
  2176. models['KMeans'] = accuracy_score(y, UnsupervisedModels.k_means().fit(x).labels_)
  2177. models['OPTICS'] = accuracy_score(y, UnsupervisedModels.optics().fit(x).fit_predict(x))
  2178. models['DBSCAN'] = accuracy_score(y, UnsupervisedModels.dbscan().fit(x).labels_)
  2179. models['Birch'] = accuracy_score(y, UnsupervisedModels.birch().fit(x).predict(x))
  2180. # ResultsShow.models_comparison_unsupervised(models)
  2181. ResultsShow.models_comparison_unsupervised(models)
  2183. import numpy as np
  2184. from sklearn.metrics import confusion_matrix
  2185. import matplotlib.pyplot as plt
  2187. # -----------------------------
  2188. # Best-mapped clustering accuracy
  2189. # -----------------------------
  2190. def best_map_accuracy(y_true, y_cluster):
  2191. y_true = np.asarray(y_true)
  2192. y_cluster = np.asarray(y_cluster)
  2194. true_labels = np.unique(y_true)
  2195. cluster_labels = np.unique(y_cluster)
  2197. # Build confusion matrix manually (handles arbitrary labels)
  2198. C = np.zeros((len(true_labels), len(cluster_labels)), dtype=int)
  2199. true_to_i = {lab:i for i, lab in enumerate(true_labels)}
  2200. clus_to_j = {lab:j for j, lab in enumerate(cluster_labels)}
  2202. for t, c in zip(y_true, y_cluster):
  2203. C[true_to_i[t], clus_to_j[c]] += 1
  2205. try:
  2206. from scipy.optimize import linear_sum_assignment
  2207. row_ind, col_ind = linear_sum_assignment(-C)
  2208. matched = C[row_ind, col_ind].sum()
  2209. return matched / len(y_true)
  2210. except Exception:
  2211. # binary fallback
  2212. if len(true_labels) == 2 and len(cluster_labels) == 2:
  2213. acc1 = np.mean(y_true == y_cluster)
  2214. flipped = 1 - y_cluster
  2215. acc2 = np.mean(y_true == flipped)
  2216. return max(acc1, acc2)
  2217. return np.mean(y_true == y_cluster)
  2220. # -----------------------------
  2221. # Run all unsupervised models
  2222. # -----------------------------
  2223. def run_unsupervised_models(x, y):
  2225. models = {}
  2226. algorithms = {
  2227. "Agglomerative": UnsupervisedModels.agglomerative_clustering(),
  2228. "Affinity Propagation": UnsupervisedModels.affinity_propagation(),
  2229. "Spectral": UnsupervisedModels.spectral_clustering(),
  2230. "Gaussian Mixture": UnsupervisedModels.gaussian_mixture(),
  2231. "Mean Shift": UnsupervisedModels.mean_shift(),
  2232. "KMeans": UnsupervisedModels.k_means(),
  2233. "OPTICS": UnsupervisedModels.optics(),
  2234. "DBSCAN": UnsupervisedModels.dbscan(),
  2235. "Birch": UnsupervisedModels.birch()
  2236. }
  2238. for name, model in algorithms.items():
  2240. try:
  2241. # Fit model
  2242. fitted = model.fit(x)
  2244. # Extract labels properly
  2245. if hasattr(fitted, "labels_"):
  2246. labels = fitted.labels_
  2247. elif hasattr(fitted, "predict"):
  2248. labels = fitted.predict(x)
  2249. elif hasattr(fitted, "fit_predict"):
  2250. labels = fitted.fit_predict(x)
  2251. else:
  2252. continue
  2254. acc = best_map_accuracy(y, labels)
  2255. models[name] = acc
  2257. except Exception as e:
  2258. print(f"{name} failed:", e)
  2260. return models
  2263. # -----------------------------
  2264. # Pretty plot
  2265. # -----------------------------
  2266. def plot_unsupervised_results(models):
  2268. names = list(models.keys())
  2269. scores = np.array(list(models.values()))
  2271. order = np.argsort(-scores)
  2272. names = [names[i] for i in order]
  2273. scores = scores[order]
  2275. plt.rcParams.update({
  2276. "font.family": "serif",
  2277. "font.serif": ["Times New Roman", "Times", "DejaVu Serif"],
  2278. "font.size": 18
  2279. })
  2281. fig, ax = plt.subplots(figsize=(9, 6), dpi=300)
  2283. y_pos = np.arange(len(names))
  2284. bars = ax.barh(y_pos, scores, color="#4A4A4A", alpha=0.9)
  2286. ax.set_yticks(y_pos)
  2287. ax.set_yticklabels(names)
  2288. ax.invert_yaxis()
  2289. ax.set_xlabel("Best-mapped Accuracy")
  2290. ax.set_xlim(0, 1)
  2292. ax.spines["top"].set_visible(False)
  2293. ax.spines["right"].set_visible(False)
  2295. # annotate values
  2296. for i, v in enumerate(scores):
  2297. ax.text(v + 0.01, i, f"{v:.3f}",
  2298. va="center", fontsize=16)
  2300. plt.tight_layout()
  2301. plt.savefig("Unsupervised_Comparison_Pretty.png", dpi=300, bbox_inches="tight")
  2302. plt.savefig("Unsupervised_Comparison_Pretty.pdf", bbox_inches="tight")
  2303. plt.show()
  2306. # -----------------------------
  2307. # RUN
  2308. # -----------------------------
  2309. models = run_unsupervised_models(x, y)
  2310. plot_unsupervised_results(models)
  2312. """## Supervised Models Comparison"""
  2314. from sklearn.metrics import PrecisionRecallDisplay
  2315. from sklearn.metrics import RocCurveDisplay
  2317. models = dict()
  2318. models['Stochastic Gradient Descent'] = cross_validation(SupervisedModels.stochastic_gradient_descent(), x_train, y_train)
  2319. models['Support Vector Machine'] = cross_validation(SupervisedModels.support_vector_machine(), x_train, y_train)
  2320. models['Multi Layer Perceptron'] = cross_validation(SupervisedModels.multi_layer_perceptron(), x_train, y_train)
  2321. models['Histogradient Boosting'] = cross_validation(SupervisedModels.histogram_gradient_boosting(), x_train, y_train)
  2322. models['Gradient Tree Boosting'] = cross_validation(SupervisedModels.gradient_boosting(), x_train, y_train)
  2323. models['Gaussian Naive Bayes'] = cross_validation(SupervisedModels.gussian_naive_bayes(), x_train, y_train)
  2324. models['Bagging Kneighbor'] = cross_validation(SupervisedModels.bagging_kneighbor(), x_train, y_train)
  2325. models['Gaussian Process'] = cross_validation(SupervisedModels.gussian_process(), x_train, y_train)
  2326. models['Label Spreading'] = cross_validation(SupervisedModels.label_spreading(), x_train, y_train)
  2327. models['Random Forest'] = cross_validation(SupervisedModels.random_forest(), x_train, y_train)
  2328. models['Decision Tree'] = cross_validation(SupervisedModels.decision_tree(), x_train, y_train)
  2329. models['Bagging SVC'] = cross_validation(SupervisedModels.bagging_svc(), x_train, y_train)
  2330. models['Self Training'] = cross_validation(SupervisedModels.selftraining(), x_train, y_train)
  2331. models['Kneighbor'] = cross_validation(SupervisedModels.kneighbor(), x_train, y_train)
  2332. models['Adaboost'] = cross_validation(SupervisedModels.adaboost(), x_train, y_train)
  2333. models['XGBoost'] = cross_validation(SupervisedModels.xgboost(), x_train, y_train)
  2335. precision, recall, report, results, names = list(), list(), list(), list(), list(models.keys())
  2336. for name, model in models.items():
  2337. report.append(classification_report(y_test, model.predict(x_test), target_names=['Unruptured','Ruptured'], output_dict=True ))
  2338. scores = accuracy_score(y_train, model.predict(x_train)), accuracy_score(y_test, model.predict(x_test))
  2339. results.append(scores)
  2340. print('>%s %.3f (%.3f)' % (name, scores[0], scores[1]))
  2341. for i in range(len(names)):
  2342. precision.append(report[i]['Unruptured']['precision'])
  2343. recall.append(report[i]['Unruptured']['recall'])
  2345. # ResultsShow.models_comparison_supervised(names, results)
  2346. ResultsShow.models_comparison_supervised_pretty(names, results)
  2347. print('-----------------------------------------------------------------------------------')
  2348. # ResultsShow.models_comparison_supervised_recall_precision(names, results, precision, recall)
  2349. ResultsShow.models_comparison_supervised_recall_precision_pretty(names, results)
  2350. print('-----------------------------------------------------------------------------------')
  2352. fig, ax = plt.subplots(figsize=(10,10))
  2353. # plt.title('Receiver Operating Characteristic (ROC) Curve for Different Models')
  2354. for name, model in models.items():
  2355. if name == 'Label Spreading':
  2356. continue
  2357. RocCurveDisplay.from_estimator(models[name], x_test, y_test, ax=ax, name=name)
  2358. plt.savefig('ROC.pdf', transparent=True, bbox_inches='tight')
  2360. fig, ax = plt.subplots(figsize=(10,10))
  2361. plt.title('Precision-Recall Curve for Different Models')
  2362. for name, model in models.items():
  2363. if name == 'Label Spreading':
  2364. continue
  2365. PrecisionRecallDisplay.from_estimator(models[name], x_test, y_test, ax=ax, name=name)
  2366. plt.savefig('Precision_Recall Curve.png', transparent=True, bbox_inches='tight')
  2368. import numpy as np
  2369. import pandas as pd
  2370. import matplotlib.pyplot as plt
  2372. from sklearn.metrics import (
  2373. accuracy_score, precision_score, recall_score, f1_score,
  2374. roc_curve, roc_auc_score, precision_recall_curve,
  2375. confusion_matrix, ConfusionMatrixDisplay
  2376. )
  2377. from sklearn.calibration import calibration_curve
  2378. from sklearn.model_selection import cross_val_score
  2380. plt.rcParams.update({
  2381. "font.family": "serif",
  2382. "font.serif": ["Times New Roman", "Times", "DejaVu Serif"],
  2383. "font.size": 16
  2384. })
  2386. # ======================================================
  2387. # MAIN FUNCTION
  2388. # ======================================================
  2389. def evaluate_supervised_models(models, X_train, y_train, X_test, y_test, cv=5):
  2391. results = []
  2393. for name, model in models.items():
  2395. model.fit(X_train, y_train)
  2397. # Predictions
  2398. y_pred_train = model.predict(X_train)
  2399. y_pred_test = model.predict(X_test)
  2401. # Probabilities (if available)
  2402. if hasattr(model, "predict_proba"):
  2403. y_prob = model.predict_proba(X_test)[:, 1]
  2404. elif hasattr(model, "decision_function"):
  2405. y_prob = model.decision_function(X_test)
  2406. else:
  2407. y_prob = None
  2409. # Metrics
  2410. acc_train = accuracy_score(y_train, y_pred_train)
  2411. acc_test = accuracy_score(y_test, y_pred_test)
  2412. precision = precision_score(y_test, y_pred_test)
  2413. recall = recall_score(y_test, y_pred_test)
  2414. f1 = f1_score(y_test, y_pred_test)
  2416. # Specificity
  2417. cm = confusion_matrix(y_test, y_pred_test)
  2418. tn, fp, fn, tp = cm.ravel()
  2419. specificity = tn / (tn + fp)
  2421. # AUC
  2422. auc = roc_auc_score(y_test, y_prob) if y_prob is not None else np.nan
  2424. # Cross-validated AUC
  2425. cv_auc = cross_val_score(model, X_train, y_train,
  2426. cv=cv, scoring="roc_auc")
  2427. cv_mean = cv_auc.mean()
  2428. cv_std = cv_auc.std()
  2430. results.append([
  2431. name, acc_train, acc_test,
  2432. precision, recall, f1,
  2433. specificity, auc, cv_mean, cv_std
  2434. ])
  2436. results_df = pd.DataFrame(results, columns=[
  2437. "Model", "Train Accuracy", "Test Accuracy",
  2438. "Precision", "Recall", "F1",
  2439. "Specificity", "AUC",
  2440. "CV AUC Mean", "CV AUC Std"
  2441. ])
  2443. results_df = results_df.sort_values("AUC", ascending=False)
  2444. print("\n===== SUMMARY METRICS =====\n")
  2445. print(results_df.round(3))
  2447. # ======================================================
  2448. # 1️⃣ Train vs Test Accuracy
  2449. # ======================================================
  2450. fig, ax = plt.subplots(figsize=(10,5))
  2451. x = np.arange(len(results_df))
  2452. w = 0.35
  2454. ax.bar(x-w/2, results_df["Train Accuracy"], w, label="Train")
  2455. ax.bar(x+w/2, results_df["Test Accuracy"], w, label="Test")
  2457. ax.set_xticks(x)
  2458. ax.set_xticklabels(results_df["Model"], rotation=30)
  2459. ax.set_ylabel("Accuracy")
  2460. ax.legend()
  2461. plt.tight_layout()
  2462. plt.show()
  2464. # ======================================================
  2465. # 2️⃣ Test Metrics Comparison
  2466. # ======================================================
  2467. fig, ax = plt.subplots(figsize=(11,6))
  2468. metrics = ["Test Accuracy","Precision","Recall","F1","Specificity"]
  2469. x = np.arange(len(results_df))
  2470. w = 0.15
  2472. for i, m in enumerate(metrics):
  2473. ax.bar(x + i*w - w*2, results_df[m], w, label=m)
  2475. ax.set_xticks(x)
  2476. ax.set_xticklabels(results_df["Model"], rotation=30)
  2477. ax.set_ylim(0,1)
  2478. ax.set_ylabel("Score")
  2479. ax.legend()
  2480. plt.tight_layout()
  2481. plt.show()
  2483. # ======================================================
  2484. # 3️⃣ ROC Curves
  2485. # ======================================================
  2486. plt.figure(figsize=(7,6))
  2487. for name, model in models.items():
  2488. model.fit(X_train, y_train)
  2490. if hasattr(model, "predict_proba"):
  2491. y_prob = model.predict_proba(X_test)[:,1]
  2492. elif hasattr(model, "decision_function"):
  2493. y_prob = model.decision_function(X_test)
  2494. else:
  2495. continue
  2497. fpr, tpr, _ = roc_curve(y_test, y_prob)
  2498. auc = roc_auc_score(y_test, y_prob)
  2499. plt.plot(fpr, tpr, label=f"{name} (AUC={auc:.2f})")
  2501. plt.plot([0,1],[0,1],'k--')
  2502. plt.xlabel("False Positive Rate")
  2503. plt.ylabel("True Positive Rate")
  2504. plt.legend()
  2505. plt.tight_layout()
  2506. plt.show()
  2508. # ======================================================
  2509. # 4️⃣ Precision-Recall Curves
  2510. # ======================================================
  2511. plt.figure(figsize=(7,6))
  2512. for name, model in models.items():
  2513. model.fit(X_train, y_train)
  2515. if hasattr(model, "predict_proba"):
  2516. y_prob = model.predict_proba(X_test)[:,1]
  2517. elif hasattr(model, "decision_function"):
  2518. y_prob = model.decision_function(X_test)
  2519. else:
  2520. continue
  2522. precision, recall, _ = precision_recall_curve(y_test, y_prob)
  2523. plt.plot(recall, precision, label=name)
  2525. plt.xlabel("Recall")
  2526. plt.ylabel("Precision")
  2527. plt.legend()
  2528. plt.tight_layout()
  2529. plt.show()
  2531. # ======================================================
  2532. # 5️⃣ Confusion Matrix (Best AUC Model)
  2533. # ======================================================
  2534. best_model_name = results_df.iloc[0]["Model"]
  2535. best_model = models[best_model_name]
  2536. best_model.fit(X_train, y_train)
  2538. y_pred = best_model.predict(X_test)
  2540. fig, ax = plt.subplots(figsize=(5,5))
  2541. ConfusionMatrixDisplay.from_predictions(
  2542. y_test, y_pred,
  2543. display_labels=["Unruptured","Ruptured"],
  2544. cmap="Greys", ax=ax
  2545. )
  2546. plt.title(f"Confusion Matrix ({best_model_name})")
  2547. plt.tight_layout()
  2548. plt.show()
  2550. # ======================================================
  2551. # 6️⃣ Calibration Curve (Best Model)
  2552. # ======================================================
  2553. if hasattr(best_model, "predict_proba"):
  2554. y_prob = best_model.predict_proba(X_test)[:,1]
  2555. prob_true, prob_pred = calibration_curve(y_test, y_prob, n_bins=10)
  2557. plt.figure(figsize=(6,6))
  2558. plt.plot(prob_pred, prob_true, marker='o')
  2559. plt.plot([0,1],[0,1],'k--')
  2560. plt.xlabel("Mean Predicted Probability")
  2561. plt.ylabel("Fraction of Positives")
  2562. plt.title("Calibration Curve")
  2563. plt.tight_layout()
  2564. plt.show()
  2566. return results_df
  2568. results = evaluate_supervised_models(models, x_train, y_train, x_test, y_test)
  2570. """# Combining Different Models"""
  2572. from sklearn.ensemble import VotingClassifier
  2574. # x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.2, random_state=None)
  2575. # standardization = StandardScaler().fit(x_train)
  2576. # x_train = standardization.transform(x_train)
  2577. # x_test = standardization.transform(x_test)
  2579. names_list = []
  2580. for i in range (2):
  2581. names_list.append(str(i+1))
  2582. models_list = [
  2583. cross_validation(SupervisedModels.stochastic_gradient_descent(), x_train, y_train),
  2584. cross_validation(SupervisedModels.support_vector_machine(), x_train, y_train),
  2585. cross_validation(SupervisedModels.multi_layer_perceptron(), x_train, y_train)
  2586. ]
  2587. models_list = list(zip(names_list, models_list))
  2588. model = VotingClassifier(models_list, voting='soft')
  2590. start = datetime.datetime.now()
  2591. model.fit(x_train, y_train)
  2592. end = datetime.datetime.now()
  2594. ResultsShow.show(model, x_train, y_train, x_test, y_test)
  2595. print('-------------------------------------------------------')
  2597. # Cross Validation Approach
  2598. # ResultsShow.show(cross_validation(model, x_train, y_train), x_train, y_train, x_test, y_test)
  2600. """#Grid Search"""
  2602. import numpy as np
  2603. import pandas as pd
  2605. from sklearn.pipeline import Pipeline
  2606. from sklearn.preprocessing import StandardScaler
  2607. from sklearn.model_selection import StratifiedKFold, GridSearchCV, RandomizedSearchCV
  2608. from sklearn.metrics import roc_auc_score
  2609. from sklearn.linear_model import SGDClassifier
  2610. from sklearn.svm import SVC
  2611. from sklearn.ensemble import HistGradientBoostingClassifier
  2613. # ----------------------------
  2614. # CONFIG
  2615. # ----------------------------
  2616. RANDOM_STATE = 42
  2617. outer_cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=RANDOM_STATE)
  2618. inner_cv = StratifiedKFold(n_splits=4, shuffle=True, random_state=RANDOM_STATE)
  2620. X = x # your features
  2621. y = y # labels (0/1)
  2623. # ----------------------------
  2624. # 1) SGD (fast) - GridSearch
  2625. # ----------------------------
  2626. sgd_pipe = Pipeline([
  2627. ("scaler", StandardScaler()),
  2628. ("clf", SGDClassifier(loss="log_loss", max_iter=20000, tol=1e-4, random_state=RANDOM_STATE))
  2629. ])
  2631. sgd_grid = {
  2632. "clf__penalty": ["l2", "l1", "elasticnet"],
  2633. "clf__alpha": [1e-5, 5e-5, 1e-4, 2e-4, 5e-4, 1e-3],
  2634. "clf__l1_ratio": [0.15, 0.5, 0.85], # only used if elasticnet; safe to include
  2635. "clf__class_weight": [None, "balanced"],
  2636. }
  2638. sgd_search = GridSearchCV(
  2639. sgd_pipe, sgd_grid, scoring="roc_auc", cv=inner_cv, n_jobs=-1, refit=True
  2640. )
  2642. # ----------------------------
  2643. # 2) SVM (RBF) - RandomizedSearch (faster than full grid)
  2644. # ----------------------------
  2645. svm_pipe = Pipeline([
  2646. ("scaler", StandardScaler()),
  2647. ("clf", SVC(kernel="rbf", probability=True, random_state=RANDOM_STATE))
  2648. ])
  2650. svm_dist = {
  2651. "clf__C": np.logspace(-1, 2.5, 18), # 0.1 .. ~316
  2652. "clf__gamma": np.logspace(-4, 0, 18), # 1e-4 .. 1
  2653. "clf__class_weight": [None, "balanced"],
  2654. }
  2656. svm_search = RandomizedSearchCV(
  2657. svm_pipe,
  2658. svm_dist,
  2659. n_iter=45, # adjust upward if you want more thorough
  2660. scoring="roc_auc",
  2661. cv=inner_cv,
  2662. n_jobs=-1,
  2663. random_state=RANDOM_STATE,
  2664. refit=True
  2665. )
  2667. # ----------------------------
  2668. # 3) HistGradientBoosting (fast) - RandomizedSearch
  2669. # ----------------------------
  2670. hgb = HistGradientBoostingClassifier(loss="log_loss", random_state=RANDOM_STATE)
  2672. hgb_dist = {
  2673. "learning_rate": [0.01, 0.02, 0.03, 0.05, 0.08],
  2674. "max_iter": [100, 150, 200, 300],
  2675. "max_leaf_nodes": [7, 10, 15, 20, 31],
  2676. "min_samples_leaf": [10, 15, 20, 30, 40],
  2677. "l2_regularization": [0.0, 0.1, 0.5, 1.0],
  2678. }
  2680. hgb_search = RandomizedSearchCV(
  2681. hgb,
  2682. hgb_dist,
  2683. n_iter=50,
  2684. scoring="roc_auc",
  2685. cv=inner_cv,
  2686. n_jobs=-1,
  2687. random_state=RANDOM_STATE,
  2688. refit=True
  2689. )
  2691. # ----------------------------
  2692. # NESTED CV EVALUATION
  2693. # ----------------------------
  2694. def nested_auc(searcher, X, y, outer_cv):
  2695. aucs = []
  2696. best_params_each_fold = []
  2698. for train_idx, test_idx in outer_cv.split(X, y):
  2699. Xtr, Xte = X[train_idx], X[test_idx]
  2700. ytr, yte = y[train_idx], y[test_idx]
  2702. searcher.fit(Xtr, ytr)
  2703. best_model = searcher.best_estimator_
  2704. best_params_each_fold.append(searcher.best_params_)
  2706. # Probabilities for AUC
  2707. if hasattr(best_model, "predict_proba"):
  2708. prob = best_model.predict_proba(Xte)[:, 1]
  2709. else:
  2710. # fallback for decision_function
  2711. prob = best_model.decision_function(Xte)
  2713. aucs.append(roc_auc_score(yte, prob))
  2715. return np.array(aucs), best_params_each_fold
  2717. models_to_tune = {
  2718. "SGD (tuned)": sgd_search,
  2719. "SVM-RBF (tuned)": svm_search,
  2720. "HistGradientBoosting (tuned)": hgb_search,
  2721. }
  2723. summary_rows = []
  2724. best_params_store = {}
  2726. for name, searcher in models_to_tune.items():
  2727. print(f"Running nested CV for: {name}")
  2728. aucs, params = nested_auc(searcher, X, y, outer_cv)
  2730. summary_rows.append({
  2731. "Model": name,
  2732. "AUC_mean": aucs.mean(),
  2733. "AUC_std": aucs.std(),
  2734. "AUC_min": aucs.min(),
  2735. "AUC_max": aucs.max()
  2736. })
  2737. best_params_store[name] = params
  2739. results_df = pd.DataFrame(summary_rows).sort_values("AUC_mean", ascending=False)
  2740. print("\n=== Nested-CV AUC summary (tuning + unbiased estimation) ===")
  2741. print(results_df.round(4))
  2743. # Optional: show the most frequent best params for each model
  2744. from collections import Counter
  2745. print("\n=== Most common best parameters across outer folds ===")
  2746. for name, params_list in best_params_store.items():
  2747. # Convert dict params to tuples for counting
  2748. as_tuples = [tuple(sorted(d.items())) for d in params_list]
  2749. most_common = Counter(as_tuples).most_common(1)[0]
  2750. print(f"\n{name}:")
  2751. print(dict(most_common[0]))
  2752. print(f"Selected in {most_common[1]} / {len(params_list)} outer folds")
  2754. import os, random
  2755. import numpy as np
  2756. import pandas as pd
  2757. import tensorflow as tf
  2759. from sklearn.model_selection import StratifiedKFold, StratifiedShuffleSplit
  2760. from sklearn.preprocessing import StandardScaler
  2761. from sklearn.metrics import roc_auc_score, confusion_matrix
  2763. # ----------------------------
  2764. # Reproducibility
  2765. # ----------------------------
  2766. SEED = 42
  2767. os.environ["PYTHONHASHSEED"] = str(SEED)
  2768. random.seed(SEED)
  2769. np.random.seed(SEED)
  2770. tf.random.set_seed(SEED)
  2772. # ----------------------------
  2773. # Helpers
  2774. # ----------------------------
  2775. def build_cnn_arch(input_len, hp):
  2776. """hp is a dict defining architecture + training hyperparameters."""
  2777. from tensorflow.keras.models import Sequential
  2778. from tensorflow.keras.layers import Conv1D, MaxPooling1D, Flatten, Dense, Dropout, BatchNormalization
  2779. from tensorflow.keras import optimizers
  2781. model = Sequential()
  2783. # Conv block 1
  2784. model.add(Conv1D(hp["filters1"], hp["kernel1"], activation="relu", input_shape=(input_len, 1)))
  2785. if hp["pool1"]:
  2786. model.add(MaxPooling1D(hp["pool1"]))
  2788. # Optional conv block 2
  2789. if hp["use_block2"]:
  2790. model.add(Conv1D(hp["filters2"], hp["kernel2"], activation="relu"))
  2791. if hp["pool2"]:
  2792. model.add(MaxPooling1D(hp["pool2"]))
  2794. model.add(Flatten())
  2796. # Dense block
  2797. model.add(Dense(hp["dense1"], activation="relu"))
  2798. if hp["batchnorm"]:
  2799. model.add(BatchNormalization())
  2800. model.add(Dropout(hp["dropout1"]))
  2802. if hp["use_dense2"]:
  2803. model.add(Dense(hp["dense2"], activation="relu"))
  2804. if hp["batchnorm"]:
  2805. model.add(BatchNormalization())
  2806. model.add(Dropout(hp["dropout2"]))
  2808. model.add(Dense(1, activation="sigmoid"))
  2810. opt_name = hp["optimizer"]
  2811. lr = hp["lr"]
  2812. if opt_name == "adam":
  2813. opt = optimizers.Adam(learning_rate=lr)
  2814. elif opt_name == "adamax":
  2815. opt = optimizers.Adamax(learning_rate=lr)
  2816. else:
  2817. opt = optimizers.RMSprop(learning_rate=lr)
  2819. model.compile(optimizer=opt, loss="binary_crossentropy")
  2820. return model
  2822. def sample_hp(rng):
  2823. """Randomly sample a small but meaningful CNN architecture space."""
  2824. hp = {
  2825. "filters1": int(rng.choice([16, 32, 48, 64])),
  2826. "kernel1": int(rng.choice([2, 3, 5])),
  2827. "pool1": int(rng.choice([0, 2])), # 0 means no pooling
  2828. "use_block2": bool(rng.choice([0, 1])),
  2829. "filters2": int(rng.choice([32, 64, 96])),
  2830. "kernel2": int(rng.choice([2, 3, 5])),
  2831. "pool2": int(rng.choice([0, 2])),
  2832. "dense1": int(rng.choice([64, 128, 192, 256])),
  2833. "use_dense2": bool(rng.choice([0, 1])),
  2834. "dense2": int(rng.choice([32, 64, 128])),
  2835. "dropout1": float(rng.choice([0.1, 0.2, 0.3, 0.4])),
  2836. "dropout2": float(rng.choice([0.1, 0.2, 0.3, 0.4])),
  2837. "batchnorm": bool(rng.choice([0, 1])),
  2838. "optimizer": str(rng.choice(["adamax", "adam", "rmsprop"])),
  2839. "lr": float(rng.choice([5e-4, 1e-3, 2e-3, 3e-3])),
  2840. "batch_size": int(rng.choice([16, 32])),
  2841. "epochs": 200
  2842. }
  2843. # enforce coherence
  2844. if not hp["use_dense2"]:
  2845. hp["dropout2"] = 0.0
  2846. return hp
  2848. def fit_eval_auc(Xtr, ytr, Xva, yva, input_len, hp, seed=SEED):
  2849. tf.keras.backend.clear_session()
  2850. tf.random.set_seed(seed)
  2852. # reshape for CNN
  2853. Xtr_c = Xtr.reshape(-1, input_len, 1)
  2854. Xva_c = Xva.reshape(-1, input_len, 1)
  2856. model = build_cnn_arch(input_len, hp)
  2857. early = tf.keras.callbacks.EarlyStopping(monitor="val_loss", patience=25, restore_best_weights=True)
  2859. model.fit(
  2860. Xtr_c, ytr,
  2861. validation_data=(Xva_c, yva),
  2862. epochs=hp["epochs"],
  2863. batch_size=hp["batch_size"],
  2864. verbose=0,
  2865. callbacks=[early],
  2866. )
  2868. prob = model.predict(Xva_c, verbose=0).ravel()
  2869. return roc_auc_score(yva, prob)
  2871. # ----------------------------
  2872. # NESTED CV RANDOM SEARCH
  2873. # ----------------------------
  2874. def nested_cv_cnn_random_search(X, y, n_outer=5, n_inner_trials=10, seed=SEED):
  2875. X = np.asarray(X)
  2876. y = np.asarray(y).astype(int)
  2877. input_len = X.shape[1]
  2879. outer = StratifiedKFold(n_splits=n_outer, shuffle=True, random_state=seed)
  2880. rng_global = np.random.default_rng(seed)
  2882. outer_rows = []
  2883. best_hps = []
  2885. for fold, (tr_idx, te_idx) in enumerate(outer.split(X, y), start=1):
  2886. # outer split
  2887. Xtr0, Xte0 = X[tr_idx], X[te_idx]
  2888. ytr0, yte0 = y[tr_idx], y[te_idx]
  2890. # scale inside fold (no leakage)
  2891. scaler = StandardScaler()
  2892. Xtr0s = scaler.fit_transform(Xtr0)
  2893. Xte0s = scaler.transform(Xte0)
  2895. # inner split (single validation split repeated by random search)
  2896. inner_split = StratifiedShuffleSplit(n_splits=1, test_size=0.2, random_state=seed + fold)
  2897. (itr, iva) = next(inner_split.split(Xtr0s, ytr0))
  2898. Xtr, Xva = Xtr0s[itr], Xtr0s[iva]
  2899. ytr, yva = ytr0[itr], ytr0[iva]
  2901. # random search on inner validation AUC
  2902. best_auc = -np.inf
  2903. best_hp = None
  2905. for t in range(n_inner_trials):
  2906. hp = sample_hp(rng_global)
  2907. auc = fit_eval_auc(Xtr, ytr, Xva, yva, input_len, hp, seed=seed + 1000*fold + t)
  2908. if auc > best_auc:
  2909. best_auc = auc
  2910. best_hp = hp
  2912. # fit best architecture on full outer-train and evaluate on outer-test
  2913. tf.keras.backend.clear_session()
  2914. tf.random.set_seed(seed + 999*fold)
  2916. best_model = build_cnn_arch(input_len, best_hp)
  2917. early = tf.keras.callbacks.EarlyStopping(monitor="val_loss", patience=25, restore_best_weights=True)
  2919. Xtr_full = Xtr0s.reshape(-1, input_len, 1)
  2920. Xte_full = Xte0s.reshape(-1, input_len, 1)
  2922. # use validation_split within outer-train to early-stop
  2923. best_model.fit(
  2924. Xtr_full, ytr0,
  2925. validation_split=0.12,
  2926. epochs=best_hp["epochs"],
  2927. batch_size=best_hp["batch_size"],
  2928. verbose=0,
  2929. callbacks=[early]
  2930. )
  2932. prob_te = best_model.predict(Xte_full, verbose=0).ravel()
  2933. auc_te = roc_auc_score(yte0, prob_te)
  2935. outer_rows.append({"fold": fold, "AUC_outer": auc_te, "AUC_inner_best": best_auc})
  2936. best_hps.append(best_hp)
  2938. print(f"Fold {fold}: inner-best AUC={best_auc:.3f} | outer AUC={auc_te:.3f}")
  2940. df = pd.DataFrame(outer_rows)
  2941. summary = {
  2942. "Model": "1D-CNN (tuned)",
  2943. "AUC_mean": df["AUC_outer"].mean(),
  2944. "AUC_std": df["AUC_outer"].std(),
  2945. "AUC_min": df["AUC_outer"].min(),
  2946. "AUC_max": df["AUC_outer"].max()
  2947. }
  2948. return df, pd.DataFrame([summary]), best_hps
  2950. # ---- RUN ----
  2951. cnn_fold_df, cnn_summary_df, cnn_best_hps = nested_cv_cnn_random_search(x, y, n_outer=5, n_inner_trials=10, seed=42)
  2953. print("\nPer-fold CNN AUC:")
  2954. print(cnn_fold_df.round(3))
  2956. print("\nCNN tuned summary:")
  2957. print(cnn_summary_df.round(4))
  2959. print("\nExample best HP (fold 1):")
  2960. print(cnn_best_hps[0])
  2962. """#Models Comparison"""
  2964. # ============================================================
  2965. # SINGLE MERGED BLOCK (ALL MODELS + CV + TEST10) INCLUDING 1D-CNN
  2966. # - final_models dict kept the same name (classical models)
  2967. # - 90/10 stratified holdout
  2968. # - 5-fold CV on TRAIN90 only
  2969. # - single TEST10 evaluation (unseen)
  2970. # - sample_weight for BaggingSVM + AdaBoost
  2971. # - BaggingKNN/KNN: optional undersampling (KNN ignores sample_weight)
  2972. # - 1D-CNN: handled separately (not in final_models)
  2973. # ============================================================
  2975. import numpy as np
  2976. import pandas as pd
  2978. from sklearn.pipeline import Pipeline
  2979. from sklearn.preprocessing import StandardScaler
  2980. from sklearn.utils.class_weight import compute_class_weight
  2982. from sklearn.model_selection import train_test_split, StratifiedKFold
  2983. from sklearn.metrics import roc_auc_score, accuracy_score, precision_score, recall_score, f1_score, confusion_matrix
  2984. from sklearn.base import clone
  2986. from sklearn.linear_model import SGDClassifier
  2987. from sklearn.svm import SVC
  2988. from sklearn.ensemble import (
  2989. HistGradientBoostingClassifier, GradientBoostingClassifier,
  2990. RandomForestClassifier, AdaBoostClassifier, BaggingClassifier
  2991. )
  2992. from sklearn.neural_network import MLPClassifier
  2993. from sklearn.gaussian_process import GaussianProcessClassifier
  2994. from sklearn.gaussian_process.kernels import RBF
  2995. from sklearn.neighbors import KNeighborsClassifier
  2996. from sklearn.naive_bayes import GaussianNB
  2997. from sklearn import tree
  2999. # TensorFlow for CNN
  3000. import tensorflow as tf
  3002. # Optional: XGBoost
  3003. xgb_available = True
  3004. try:
  3005. import xgboost as xgb
  3006. except Exception:
  3007. xgb_available = False
  3009. # ----------------------------
  3010. # CONFIG
  3011. # ----------------------------
  3012. RANDOM_STATE_SPLIT = 8
  3013. CV_FOLDS = 5
  3014. THR = 0.5
  3016. USE_UNDERSAMPLE_FOR_KNN_BAGGING = True # recommended for KNN-like models (sample_weight won't help)
  3018. CNN_HP = {
  3019. "filters1": 64, "kernel1": 3, "pool1": 0,
  3020. "use_block2": True, "filters2": 64, "kernel2": 5, "pool2": 0,
  3021. "dense1": 128, "use_dense2": True, "dense2": 128,
  3022. "dropout1": 0.2, "dropout2": 0.4,
  3023. "batchnorm": True,
  3024. "optimizer": "adam",
  3025. "lr": 0.003,
  3026. "batch_size": 16,
  3027. "epochs": 200,
  3028. "patience": 25,
  3029. "val_split": 0.12
  3030. }
  3032. # ----------------------------
  3033. # DATA
  3034. # ----------------------------
  3035. X = np.asarray(x)
  3036. y_arr = np.asarray(y).astype(int).ravel() # IMPORTANT: 1D
  3038. # Global class weights (informational)
  3039. unique_classes = np.unique(y_arr)
  3040. print("Unique labels in y:", unique_classes)
  3042. cw = compute_class_weight(
  3043. class_weight="balanced",
  3044. classes=unique_classes,
  3045. y=y_arr
  3046. )
  3048. CLASS_WEIGHT_DICT = dict(zip(unique_classes, cw))
  3049. print("Computed class weights:", CLASS_WEIGHT_DICT)
  3050. # CLASS_WEIGHT_DICT = {0: float(cw[0]), 1: float(cw[1])}
  3051. print("Global class weights:", CLASS_WEIGHT_DICT)
  3053. # 90/10 split (stratified)
  3054. X_train90, X_test10, y_train90, y_test10 = train_test_split(
  3055. X, y_arr, test_size=0.10, random_state=RANDOM_STATE_SPLIT, stratify=y_arr
  3056. )
  3057. print("Train (90%) shape:", X_train90.shape, "Test (10%) shape:", X_test10.shape)
  3058. print("Train class counts:", np.bincount(y_train90))
  3059. print("Test class counts:", np.bincount(y_test10))
  3061. # ----------------------------
  3062. # HELPERS
  3063. # ----------------------------
  3064. def specificity(y_true, y_pred):
  3065. tn, fp, fn, tp = confusion_matrix(y_true, y_pred, labels=[0, 1]).ravel()
  3066. return tn / (tn + fp) if (tn + fp) else 0.0
  3068. def metrics_from_prob(y_true, prob, thr=0.5):
  3069. pred = (prob >= thr).astype(int)
  3070. return {
  3071. "AUC": roc_auc_score(y_true, prob),
  3072. "ACC": accuracy_score(y_true, pred),
  3073. "Precision": precision_score(y_true, pred, zero_division=0),
  3074. "Recall": recall_score(y_true, pred, zero_division=0),
  3075. "F1": f1_score(y_true, pred, zero_division=0),
  3076. "Specificity": specificity(y_true, pred),
  3077. }
  3079. def get_probabilities(estimator, X_eval):
  3080. if hasattr(estimator, "predict_proba"):
  3081. return estimator.predict_proba(X_eval)[:, 1]
  3082. if hasattr(estimator, "decision_function"):
  3083. s = estimator.decision_function(X_eval)
  3084. return (s - s.min()) / (s.max() - s.min() + 1e-12)
  3085. return estimator.predict(X_eval).astype(float)
  3087. def make_sample_weight_balanced(y):
  3088. y = np.asarray(y).astype(int).ravel()
  3089. n = len(y)
  3090. n0 = np.sum(y == 0)
  3091. n1 = np.sum(y == 1)
  3092. w0 = n / (2 * n0) if n0 else 1.0
  3093. w1 = n / (2 * n1) if n1 else 1.0
  3094. return np.where(y == 0, w0, w1).astype(float)
  3096. def undersample_to_balance(X, y, seed=42):
  3097. rng = np.random.default_rng(seed)
  3098. y = np.asarray(y).astype(int).ravel()
  3099. idx0 = np.where(y == 0)[0]
  3100. idx1 = np.where(y == 1)[0]
  3101. n = min(len(idx0), len(idx1))
  3102. if n == 0:
  3103. return X, y
  3104. idx0_s = rng.choice(idx0, size=n, replace=False)
  3105. idx1_s = rng.choice(idx1, size=n, replace=False)
  3106. sel = np.concatenate([idx0_s, idx1_s])
  3107. rng.shuffle(sel)
  3108. return X[sel], y[sel]
  3110. # ----------------------------
  3111. # CNN BUILD / FIT
  3112. # ----------------------------
  3113. def build_cnn_arch(input_len, hp):
  3114. from tensorflow.keras.models import Sequential
  3115. from tensorflow.keras.layers import Conv1D, MaxPooling1D, Flatten, Dense, Dropout, BatchNormalization
  3116. from tensorflow.keras import optimizers
  3118. model = Sequential()
  3119. model.add(tf.keras.layers.Input(shape=(input_len, 1)))
  3121. model.add(Conv1D(hp["filters1"], hp["kernel1"], activation="relu"))
  3122. if hp["pool1"] and hp["pool1"] > 0:
  3123. model.add(MaxPooling1D(hp["pool1"]))
  3125. if hp["use_block2"]:
  3126. model.add(Conv1D(hp["filters2"], hp["kernel2"], activation="relu"))
  3127. if hp["pool2"] and hp["pool2"] > 0:
  3128. model.add(MaxPooling1D(hp["pool2"]))
  3130. model.add(Flatten())
  3132. model.add(Dense(hp["dense1"], activation="relu"))
  3133. if hp["batchnorm"]:
  3134. model.add(BatchNormalization())
  3135. model.add(Dropout(hp["dropout1"]))
  3137. if hp["use_dense2"]:
  3138. model.add(Dense(hp["dense2"], activation="relu"))
  3139. if hp["batchnorm"]:
  3140. model.add(BatchNormalization())
  3141. model.add(Dropout(hp["dropout2"]))
  3143. model.add(Dense(1, activation="sigmoid"))
  3145. lr = hp["lr"]
  3146. if hp["optimizer"] == "adam":
  3147. opt = optimizers.Adam(learning_rate=lr)
  3148. elif hp["optimizer"] == "adamax":
  3149. opt = optimizers.Adamax(learning_rate=lr)
  3150. else:
  3151. opt = optimizers.RMSprop(learning_rate=lr)
  3153. model.compile(optimizer=opt, loss="binary_crossentropy")
  3154. return model
  3156. def fit_cnn_and_predict_proba(Xtr, ytr, Xeval, hp, class_weight_dict, seed=42):
  3157. """
  3158. Fit CNN on Xtr/ytr with internal val_split, return probabilities on Xeval.
  3159. Scaling is handled outside; inputs must be already standardized.
  3160. """
  3161. tf.keras.backend.clear_session()
  3162. tf.random.set_seed(seed)
  3164. cnn = build_cnn_arch(Xtr.shape[1], hp)
  3165. early = tf.keras.callbacks.EarlyStopping(
  3166. monitor="val_loss",
  3167. patience=hp["patience"],
  3168. restore_best_weights=True
  3169. )
  3171. Xtr_c = Xtr.reshape(-1, Xtr.shape[1], 1)
  3172. Xev_c = Xeval.reshape(-1, Xtr.shape[1], 1)
  3174. cnn.fit(
  3175. Xtr_c, ytr,
  3176. validation_split=hp["val_split"],
  3177. epochs=hp["epochs"],
  3178. batch_size=hp["batch_size"],
  3179. verbose=0,
  3180. callbacks=[early],
  3181. class_weight=class_weight_dict
  3182. )
  3184. prob = cnn.predict(Xev_c, verbose=0).ravel()
  3185. return prob
  3187. # ----------------------------
  3188. # DEFINE CLASSICAL MODELS (final_models dict name kept)
  3189. # ----------------------------
  3190. final_models = {}
  3192. final_models["BaggingKNN"] = Pipeline([
  3193. ("scaler", StandardScaler()),
  3194. ("clf", BaggingClassifier(
  3195. estimator=KNeighborsClassifier(n_neighbors=5, weights="uniform", algorithm="auto"),
  3196. n_estimators=100, bootstrap=True, max_samples=0.9, max_features=0.7,
  3197. random_state=42
  3198. ))
  3199. ])
  3201. final_models["BaggingSVM"] = Pipeline([
  3202. ("scaler", StandardScaler()),
  3203. ("clf", BaggingClassifier(
  3204. estimator=SVC(C=80, gamma=0.002, probability=True, kernel="rbf"),
  3205. n_estimators=100, bootstrap=True, max_samples=0.9, max_features=54,
  3206. random_state=42
  3207. ))
  3208. ])
  3210. final_models["DecisionTree"] = tree.DecisionTreeClassifier(
  3211. criterion="gini", splitter="random",
  3212. max_depth=4, max_features=None, min_samples_split=12,
  3213. random_state=42
  3214. )
  3216. final_models["SGD"] = Pipeline([
  3217. ("scaler", StandardScaler()),
  3218. ("clf", SGDClassifier(
  3219. loss="log_loss",
  3220. penalty="elasticnet",
  3221. alpha=0.001,
  3222. l1_ratio=0.5,
  3223. max_iter=20000,
  3224. tol=1e-4,
  3225. class_weight="balanced",
  3226. random_state=42
  3227. ))
  3228. ])
  3230. final_models["SVM"] = Pipeline([
  3231. ("scaler", StandardScaler()),
  3232. ("clf", SVC(
  3233. kernel="rbf", probability=True,
  3234. C=2.761699813343848,
  3235. gamma=0.004436687330978612,
  3236. class_weight="balanced",
  3237. random_state=42
  3238. ))
  3239. ])
  3241. final_models["HistGradientBoosting"] = HistGradientBoostingClassifier(
  3242. loss="log_loss", learning_rate=0.08,
  3243. max_iter=150, max_leaf_nodes=7,
  3244. min_samples_leaf=20, l2_regularization=1.0,
  3245. random_state=42
  3246. )
  3248. final_models["MLP"] = Pipeline([
  3249. ("scaler", StandardScaler()),
  3250. ("clf", MLPClassifier(
  3251. hidden_layer_sizes=(55, 18, 1),
  3252. tol=1e-6, batch_size=28, max_iter=2500,
  3253. activation="identity", solver="sgd", alpha=0.2,
  3254. learning_rate_init=0.00005, learning_rate="adaptive",
  3255. random_state=42
  3256. ))
  3257. ])
  3259. final_models["GaussianProcess"] = Pipeline([
  3260. ("scaler", StandardScaler()),
  3261. ("clf", GaussianProcessClassifier(kernel=1 * RBF(2), random_state=42))
  3262. ])
  3264. final_models["RandomForest"] = RandomForestClassifier(
  3265. n_estimators=80, criterion="log_loss", max_depth=5, min_samples_leaf=5,
  3266. max_features=0.5, min_samples_split=20, bootstrap=True,
  3267. class_weight="balanced",
  3268. random_state=42
  3269. )
  3271. final_models["GradientBoosting"] = GradientBoostingClassifier(
  3272. n_estimators=60, learning_rate=0.02, max_depth=4,
  3273. min_samples_split=10, random_state=42
  3274. )
  3276. final_models["AdaBoost"] = AdaBoostClassifier(
  3277. n_estimators=110, learning_rate=0.15, random_state=42
  3278. )
  3280. final_models["KNN"] = Pipeline([
  3281. ("scaler", StandardScaler()),
  3282. ("clf", KNeighborsClassifier(n_neighbors=17, weights="uniform", p=1))
  3283. ])
  3285. final_models["GaussianNB"] = GaussianNB(var_smoothing=1e-10)
  3287. if xgb_available:
  3288. pos = float(np.sum(y_train90 == 1))
  3289. neg = float(np.sum(y_train90 == 0))
  3290. scale_pos_weight = (neg / pos) if pos > 0 else 1.0
  3292. final_models["XGBoost"] = xgb.XGBClassifier(
  3293. booster="gbtree",
  3294. eta=0.05, gamma=0,
  3295. max_depth=4, learning_rate=0.05,
  3296. subsample=0.5, colsample_bytree=0.5,
  3297. alpha=1,
  3298. tree_method="hist",
  3299. random_state=42,
  3300. eval_metric="logloss",
  3301. scale_pos_weight=scale_pos_weight
  3302. )
  3304. print("\nModels in final_models:", list(final_models.keys()))
  3305. print("Also evaluating: 1D-CNN")
  3307. # ============================================================
  3308. # 1) 5-fold CV on TRAIN90 only
  3309. # ============================================================
  3310. skf = StratifiedKFold(n_splits=CV_FOLDS, shuffle=True, random_state=42)
  3311. results_cv = {}
  3313. for fold, (tr_idx, va_idx) in enumerate(skf.split(X_train90, y_train90), start=1):
  3314. Xtr, Xva = X_train90[tr_idx], X_train90[va_idx]
  3315. ytr, yva = y_train90[tr_idx], y_train90[va_idx]
  3317. sw_tr = make_sample_weight_balanced(ytr)
  3319. # ---- Classical models ----
  3320. for name, mdl in final_models.items():
  3321. est = clone(mdl)
  3323. if USE_UNDERSAMPLE_FOR_KNN_BAGGING and name in ["BaggingKNN", "KNN"]:
  3324. Xfit, yfit = undersample_to_balance(Xtr, ytr, seed=1000 + fold)
  3325. est.fit(Xfit, yfit)
  3327. elif name == "BaggingSVM":
  3328. est.fit(Xtr, ytr, clf__sample_weight=sw_tr)
  3330. elif name == "AdaBoost":
  3331. est.fit(Xtr, ytr, sample_weight=sw_tr)
  3333. else:
  3334. est.fit(Xtr, ytr)
  3336. prob = get_probabilities(est, Xva)
  3337. results_cv.setdefault(name, []).append(metrics_from_prob(yva, prob, thr=THR))
  3339. # ---- 1D-CNN (per-fold scaling + per-fold class weights) ----
  3340. scaler = StandardScaler()
  3341. Xtr_s = scaler.fit_transform(Xtr)
  3342. Xva_s = scaler.transform(Xva)
  3344. cw_fold = compute_class_weight(class_weight="balanced", classes=np.array([0, 1]), y=ytr)
  3345. cw_dict = {0: float(cw_fold[0]), 1: float(cw_fold[1])}
  3347. prob_cnn = fit_cnn_and_predict_proba(
  3348. Xtr_s, ytr, Xva_s,
  3349. hp=CNN_HP,
  3350. class_weight_dict=cw_dict,
  3351. seed=42 + fold
  3352. )
  3354. results_cv.setdefault("1D-CNN", []).append(metrics_from_prob(yva, prob_cnn, thr=THR))
  3356. # Build CV table (mean ± SD in one cell)
  3357. rows = []
  3358. for name, fold_list in results_cv.items():
  3359. dfm = pd.DataFrame(fold_list)
  3360. rows.append({
  3361. "Model": name,
  3362. "ROC_AUC": f"{dfm['AUC'].mean():.3f} ± {dfm['AUC'].std(ddof=0):.3f}",
  3363. "Accuracy": f"{dfm['ACC'].mean():.3f} ± {dfm['ACC'].std(ddof=0):.3f}",
  3364. "Precision": f"{dfm['Precision'].mean():.3f} ± {dfm['Precision'].std(ddof=0):.3f}",
  3365. "Recall": f"{dfm['Recall'].mean():.3f} ± {dfm['Recall'].std(ddof=0):.3f}",
  3366. "F1": f"{dfm['F1'].mean():.3f} ± {dfm['F1'].std(ddof=0):.3f}",
  3367. "Specificity": f"{dfm['Specificity'].mean():.3f} ± {dfm['Specificity'].std(ddof=0):.3f}",
  3368. })
  3370. cv_table_train90 = pd.DataFrame(rows)
  3371. cv_table_train90["_auc_mean"] = cv_table_train90["ROC_AUC"].str.split(" ± ").str[0].astype(float)
  3372. cv_table_train90 = cv_table_train90.sort_values("_auc_mean", ascending=False).drop(columns="_auc_mean").reset_index(drop=True)
  3374. print("\n=== 5-fold CV on TRAIN90 (mean ± SD) ===")
  3375. display(cv_table_train90)
  3377. # ============================================================
  3378. # 2) Final fit on TRAIN90, evaluate ONCE on TEST10 (unseen)
  3379. # ============================================================
  3380. test_rows = []
  3381. sw_train90 = make_sample_weight_balanced(y_train90)
  3383. # ---- Classical models ----
  3384. for name, mdl in final_models.items():
  3385. est = clone(mdl)
  3387. if USE_UNDERSAMPLE_FOR_KNN_BAGGING and name in ["BaggingKNN", "KNN"]:
  3388. Xfit, yfit = undersample_to_balance(X_train90, y_train90, seed=999)
  3389. est.fit(Xfit, yfit)
  3391. elif name == "BaggingSVM":
  3392. est.fit(X_train90, y_train90, clf__sample_weight=sw_train90)
  3394. elif name == "AdaBoost":
  3395. est.fit(X_train90, y_train90, sample_weight=sw_train90)
  3397. else:
  3398. est.fit(X_train90, y_train90)
  3400. prob = get_probabilities(est, X_test10)
  3401. m = metrics_from_prob(y_test10, prob, thr=THR)
  3402. m["Model"] = name
  3403. test_rows.append(m)
  3405. # ---- 1D-CNN (final fit on TRAIN90) ----
  3406. scaler = StandardScaler()
  3407. Xtr_s = scaler.fit_transform(X_train90)
  3408. Xte_s = scaler.transform(X_test10)
  3410. cw_train = compute_class_weight(class_weight="balanced", classes=np.array([0, 1]), y=y_train90)
  3411. cw_dict = {0: float(cw_train[0]), 1: float(cw_train[1])}
  3413. prob_cnn = fit_cnn_and_predict_proba(
  3414. Xtr_s, y_train90, Xte_s,
  3415. hp=CNN_HP,
  3416. class_weight_dict=cw_dict,
  3417. seed=42
  3418. )
  3420. m = metrics_from_prob(y_test10, prob_cnn, thr=THR)
  3421. m["Model"] = "1D-CNN"
  3422. test_rows.append(m)
  3424. test_table = pd.DataFrame(test_rows).sort_values("AUC", ascending=False).reset_index(drop=True)
  3426. print("\n=== Held-out TEST10 (single evaluation) ===")
  3427. display(test_table.round(3))
  3429. # ============================================================
  3430. # FIGURE ONLY: Held-out TEST set — all metrics
  3431. # Fixes:
  3432. # (1) value labels use a constant pixel/point offset (independent of xlim)
  3433. # (2) extra spacing between x ticks and x-axis label
  3434. # Requires: test_table with columns:
  3435. # Model, AUC, ACC, Precision, Recall, Specificity, F1
  3436. # ============================================================
  3438. import numpy as np
  3439. import matplotlib.pyplot as plt
  3440. from matplotlib import transforms
  3441. from matplotlib.ticker import MaxNLocator
  3443. # ---- Sort models (by AUC descending) ----
  3444. df = test_table.sort_values("AUC", ascending=False).reset_index(drop=True)
  3445. models = df["Model"].values
  3446. ypos = np.arange(len(models))
  3448. metric_defs = [
  3449. ("ROC–AUC", "AUC"),
  3450. ("Accuracy", "ACC"),
  3451. ("Precision", "Precision"),
  3452. ("Recall", "Recall"),
  3453. ("F1-score", "F1"),
  3454. ("Specificity", "Specificity"),
  3455. ]
  3457. # ============================================================
  3458. # >>> EDIT THESE LIMITS AS YOU LIKE <<<
  3459. # ============================================================
  3460. x_limits = {
  3461. "AUC": (0.65, 0.80),
  3462. "ACC": (0.60, 0.75),
  3463. "Precision": (0.44, 0.7),
  3464. "Recall": (0.38, 0.812),
  3465. "Specificity": (0.52, 0.90),
  3466. "F1": (0.40, 0.70),
  3467. }
  3469. # ============================================================
  3470. # >>> MANUAL TICK CONTROL <<<
  3471. # ============================================================
  3472. x_ticks = {
  3473. "AUC": [0.65, 0.80],
  3474. "ACC": [0.60, 0.75],
  3475. "Precision": [0.45, 0.70],
  3476. "Recall": [0.40, 0.80],
  3477. "Specificity": [0.55, 0.90],
  3478. "F1": [0.40, 0.70],
  3479. }
  3481. # ============================================================
  3483. # -----------------------------
  3484. # USER CONTROLS (NEW)
  3485. # -----------------------------
  3486. LABEL_OFFSET_PTS = {
  3487. # constant offset in *points* to the right of dot (independent of x-lims)
  3488. "AUC": 10,
  3489. "ACC": 10,
  3490. "Precision": 10,
  3491. "Recall": 10, # give Recall a bit more room if you want
  3492. "Specificity": 10,
  3493. "F1": 10,
  3494. }
  3495. VALUE_FONTSIZE = {
  3496. "AUC": 26,
  3497. "ACC": 26,
  3498. "Precision": 26,
  3499. "Recall": 26,
  3500. "Specificity": 26,
  3501. "F1": 26,
  3502. }
  3504. DOT_SIZE = 8
  3506. # spacing controls
  3507. XTICK_PAD = 10 # distance between x ticks and axis spine
  3508. XLABEL_PAD = 12 # distance between x-axis label and ticks
  3509. SUBPLOT_WSPACE = -1 # reduce/adjust spacing between panels (smaller = tighter)
  3511. # ---- Style ----
  3512. plt.rcParams.update({
  3513. "font.family": "serif",
  3514. "font.serif": ["Times New Roman", "Times", "DejaVu Serif"],
  3515. "font.size": 26,
  3516. "axes.labelsize": 30,
  3517. "xtick.labelsize": 26,
  3518. "ytick.labelsize": 30,
  3519. })
  3521. fig, axes = plt.subplots(
  3522. 1, len(metric_defs),
  3523. figsize=(22, 12.5),
  3524. dpi=300,
  3525. sharey=True
  3526. )
  3528. for ax, (title, col) in zip(axes, metric_defs):
  3529. vals = df[col].values
  3531. # dots
  3532. ax.plot(vals, ypos, "o", markersize=DOT_SIZE)
  3534. # ---- value labels: constant offset in points (NOT data units) ----
  3535. # This keeps spacing consistent even if xlim differs across panels.
  3536. dx = LABEL_OFFSET_PTS.get(col, 10) / 72.0 # points -> inches
  3537. text_trans = ax.transData + transforms.ScaledTranslation(dx, 0, fig.dpi_scale_trans)
  3539. for v, y in zip(vals, ypos):
  3540. ax.text(
  3541. v, y, f"{v:.3f}",
  3542. transform=text_trans,
  3543. va="center", ha="left",
  3544. fontsize=VALUE_FONTSIZE.get(col, 16)
  3545. )
  3547. # custom limits
  3548. if col in x_limits:
  3549. ax.set_xlim(x_limits[col])
  3551. if col in x_ticks:
  3552. ax.set_xticks(x_ticks[col])
  3553. # labels & spacing
  3554. ax.set_xlabel(title, labelpad=XLABEL_PAD)
  3555. ax.tick_params(axis="x", pad=XTICK_PAD)
  3557. ax.spines["top"].set_visible(False)
  3558. ax.spines["right"].set_visible(False)
  3559. ax.grid(False)
  3561. # Y-axis labels only on first subplot
  3562. axes[0].set_yticks(ypos)
  3563. axes[0].set_yticklabels(models)
  3564. axes[0].invert_yaxis()
  3566. for ax in axes[1:]:
  3567. ax.tick_params(axis="y", left=False, labelleft=False)
  3569. # control spacing between subplots
  3570. fig.subplots_adjust(wspace=SUBPLOT_WSPACE)
  3572. plt.tight_layout()
  3573. plt.savefig("TEST10_all_metrics.png", dpi=300, bbox_inches="tight")
  3574. plt.savefig("TEST10_all_metrics.pdf", bbox_inches="tight")
  3575. plt.show()
  3577. """#Recall"""
  3579. import numpy as np
  3580. import pandas as pd
  3581. import matplotlib.pyplot as plt
  3582. from sklearn.metrics import confusion_matrix, precision_score, recall_score, f1_score, accuracy_score
  3584. def metrics_at_threshold(y_true, prob, thr):
  3585. y_pred = (prob >= thr).astype(int)
  3586. tn, fp, fn, tp = confusion_matrix(y_true, y_pred, labels=[0,1]).ravel()
  3588. acc = (tp + tn) / (tp + tn + fp + fn)
  3589. rec = tp / (tp + fn) if (tp + fn) else 0.0 # sensitivity / recall
  3590. spec = tn / (tn + fp) if (tn + fp) else 0.0 # specificity
  3591. prec = tp / (tp + fp) if (tp + fp) else 0.0 # precision
  3592. f1 = 2*prec*rec/(prec+rec) if (prec+rec) else 0.0
  3593. return acc, prec, rec, spec, f1
  3595. def threshold_sweep(y_true, prob, target_recall=0.80, n_grid=501, make_plot=True):
  3596. y_true = np.asarray(y_true).astype(int)
  3597. prob = np.asarray(prob).astype(float)
  3599. thrs = np.linspace(0, 1, n_grid)
  3600. rows = []
  3601. for t in thrs:
  3602. acc, prec, rec, spec, f1 = metrics_at_threshold(y_true, prob, t)
  3603. J = rec + spec - 1
  3604. rows.append([t, acc, prec, rec, spec, f1, J])
  3606. df = pd.DataFrame(rows, columns=["thr","acc","prec","recall","spec","f1","youdenJ"])
  3608. # 1) Best Youden's J
  3609. best_J = df.iloc[df["youdenJ"].values.argmax()]
  3611. # 2) Best specificity subject to recall >= target
  3612. feasible = df[df["recall"] >= target_recall].copy()
  3613. if len(feasible) > 0:
  3614. best_target = feasible.iloc[feasible["spec"].values.argmax()]
  3615. else:
  3616. best_target = None
  3618. print("=== Best threshold by Youden's J (recall + spec − 1) ===")
  3619. print(best_J[["thr","acc","prec","recall","spec","f1","youdenJ"]].to_string())
  3621. if best_target is not None:
  3622. print(f"\n=== Best threshold achieving recall ≥ {target_recall:.2f} (max specificity) ===")
  3623. print(best_target[["thr","acc","prec","recall","spec","f1","youdenJ"]].to_string())
  3624. else:
  3625. print(f"\nNo threshold reached recall ≥ {target_recall:.2f}. "
  3626. f"Max recall achieved = {df['recall'].max():.3f} at thr={df.loc[df['recall'].idxmax(),'thr']:.3f}")
  3628. if make_plot:
  3629. fig, ax = plt.subplots(figsize=(10,6), dpi=300)
  3630. ax.plot(df["thr"], df["recall"], label="Recall (Sensitivity)")
  3631. ax.plot(df["thr"], df["spec"], label="Specificity")
  3632. ax.plot(df["thr"], df["prec"], label="Precision")
  3633. ax.plot(df["thr"], df["f1"], label="F1-score", alpha=0.8)
  3635. ax.axvline(best_J["thr"], linestyle="--", linewidth=1.5, label=f"Best J thr={best_J['thr']:.2f}")
  3637. if best_target is not None:
  3638. ax.axvline(best_target["thr"], linestyle=":", linewidth=2.0,
  3639. label=f"Recall≥{target_recall:.2f} thr={best_target['thr']:.2f}")
  3641. ax.set_xlabel("Decision threshold")
  3642. ax.set_ylabel("Score")
  3643. ax.set_ylim(0, 1)
  3644. ax.spines["top"].set_visible(False)
  3645. ax.spines["right"].set_visible(False)
  3646. ax.legend(frameon=False, ncol=1, loc="center left", bbox_to_anchor=(1.02, 0.5))
  3647. plt.tight_layout()
  3648. plt.show()
  3650. return df, best_J, best_target
  3652. # -------------------------
  3653. # USE IT (choose your model)
  3654. # -------------------------
  3655. model_name = "SVM" # try "SGD (log-loss)", "MLP", etc.
  3657. prob = oof_probs.get(model_name, None)
  3658. if prob is None:
  3659. raise ValueError(f"No OOF probabilities found for {model_name}. "
  3660. f"Make sure it supports predict_proba and was run in CV.")
  3662. df_thr, best_J, best_target = threshold_sweep(y, prob, target_recall=0.90, n_grid=501, make_plot=True)
  3664. # ============================================================
  3665. # Top-3 models: Annotated horizontal bar comparison
  3666. # BEFORE vs AFTER threshold tuning
  3667. # ============================================================
  3669. import numpy as np
  3670. import matplotlib.pyplot as plt
  3671. import matplotlib.colors as mcolors
  3673. def lighten_color(color, amount=0.35):
  3674. """
  3675. Blend color toward white.
  3676. amount=0 -> original color
  3677. amount=1 -> fully white
  3678. 0.2–0.5 -> good softening range
  3679. """
  3680. c = np.array(mcolors.to_rgb(color))
  3681. white = np.array([1, 1, 1])
  3682. return tuple(c + (white - c) * amount)
  3684. BASE_COLORS = {
  3685. "SVM": "#173A72",
  3686. "GaussianProcess": "#004C1C",
  3687. "MLP": "#4B0082"
  3688. }
  3690. # ---- Top 3 models ----
  3691. top3 = ["SVM", "GaussianProcess", "MLP"]
  3692. # --- model colors (from your donut first-slice codes) ---
  3693. MODEL_COLORS = {
  3694. model: lighten_color(color, amount=0.2) # <-- adjust this
  3695. for model, color in BASE_COLORS.items()
  3696. }
  3697. # Full metric names
  3698. metrics = ["Accuracy", "Precision", "Recall", "Specificity", "F1-score"]
  3699. metric_keys = ["ACC", "Precision", "Recall", "Spec", "F1"]
  3701. before_vals = before_df.loc[top3, metric_keys]
  3702. after_vals = after_df.loc[top3, metric_keys]
  3704. plt.rcParams.update({
  3705. "font.family": "serif",
  3706. "font.serif": ["Times New Roman", "Times", "DejaVu Serif"],
  3707. "font.size": 20,
  3708. "axes.labelsize": 24,
  3709. "xtick.labelsize": 22,
  3710. "ytick.labelsize": 22,
  3711. })
  3713. fig, axes = plt.subplots(1, 2, figsize=(16, 7.5), dpi=300, sharey=True)
  3715. bar_h = 0.24
  3716. y_base = np.arange(len(metrics))
  3718. # ------------------------------------------------
  3719. # (a) BEFORE
  3720. # ------------------------------------------------
  3721. for i, model in enumerate(top3):
  3722. bars = axes[0].barh(
  3723. y_base + (i - 1)*bar_h,
  3724. before_vals.loc[model],
  3725. height=bar_h,
  3726. label=model,
  3727. color=MODEL_COLORS[model]
  3728. )
  3729. # annotate values
  3730. for bar in bars:
  3731. width = bar.get_width()
  3732. axes[0].text(width + 0.01, bar.get_y() + bar.get_height()/2,
  3733. f"{width:.2f}", va="center", fontsize=20)
  3735. axes[0].set_xlim(0.5, 1.0)
  3736. axes[0].set_yticks(y_base)
  3737. axes[0].set_yticklabels(metrics)
  3738. axes[0].set_xlabel("Score", labelpad=10)
  3739. axes[0].spines["top"].set_visible(False)
  3740. axes[0].spines["right"].set_visible(False)
  3741. axes[0].grid(False)
  3743. axes[0].text(-0.30, 1.02, "a)", transform=axes[0].transAxes,
  3744. fontsize=24, fontweight="bold", va="top")
  3746. # ------------------------------------------------
  3747. # (b) AFTER
  3748. # ------------------------------------------------
  3749. for i, model in enumerate(top3):
  3750. bars = axes[1].barh(
  3751. y_base + (i - 1)*bar_h,
  3752. after_vals.loc[model],
  3753. height=bar_h,
  3754. label=model,
  3755. color=MODEL_COLORS[model]
  3756. )
  3757. # annotate values
  3758. for bar in bars:
  3759. width = bar.get_width()
  3760. axes[1].text(width + 0.01, bar.get_y() + bar.get_height()/2,
  3761. f"{width:.2f}", va="center", fontsize=20)
  3763. axes[1].set_xlim(0.5, 1.0)
  3764. axes[1].set_xlabel("Score", labelpad=10)
  3766. axes[1].spines["top"].set_visible(False)
  3767. axes[1].spines["right"].set_visible(False)
  3768. axes[1].grid(False)
  3770. axes[1].text(-0.12, 1.02, "b)", transform=axes[1].transAxes,
  3771. fontsize=24, fontweight="bold", va="top")
  3773. # Legend top right
  3774. axes[1].legend(frameon=False, loc="upper right", ncol=1, fontsize=20, bbox_to_anchor=(1.07, 1))
  3776. plt.tight_layout()
  3777. plt.savefig("Top3_metrics.png", dpi=300, bbox_inches="tight")
  3778. plt.savefig("Top3_metrics.pdf", bbox_inches="tight")
  3779. plt.show()
  3781. import numpy as np
  3782. import pandas as pd
  3784. from sklearn.pipeline import Pipeline
  3785. from sklearn.preprocessing import StandardScaler
  3786. from sklearn.utils.class_weight import compute_class_weight
  3788. from sklearn.model_selection import train_test_split, StratifiedKFold
  3789. from sklearn.metrics import roc_auc_score, accuracy_score, precision_score, recall_score, f1_score, confusion_matrix
  3790. from sklearn.base import clone
  3792. from sklearn.linear_model import SGDClassifier
  3793. from sklearn.svm import SVC
  3794. from sklearn.ensemble import (
  3795. HistGradientBoostingClassifier, GradientBoostingClassifier,
  3796. RandomForestClassifier, AdaBoostClassifier, BaggingClassifier
  3797. )
  3798. from sklearn.neural_network import MLPClassifier
  3799. from sklearn.gaussian_process import GaussianProcessClassifier
  3800. from sklearn.gaussian_process.kernels import RBF
  3801. from sklearn.neighbors import KNeighborsClassifier
  3802. from sklearn.naive_bayes import GaussianNB
  3803. from sklearn import tree
  3805. # TensorFlow for CNN
  3806. import tensorflow as tf
  3808. # Optional: XGBoost
  3809. xgb_available = True
  3810. try:
  3811. import xgboost as xgb
  3812. except Exception:
  3813. xgb_available = False
  3815. # ----------------------------
  3816. # CONFIG
  3817. # ----------------------------
  3818. RANDOM_STATE_SPLIT = 8
  3819. CV_FOLDS = 5
  3820. THR = 0.5
  3822. USE_UNDERSAMPLE_FOR_KNN_BAGGING = True # recommended for KNN-like models (sample_weight won't help)
  3824. CNN_HP = {
  3825. "filters1": 64, "kernel1": 3, "pool1": 0,
  3826. "use_block2": True, "filters2": 64, "kernel2": 5, "pool2": 0,
  3827. "dense1": 128, "use_dense2": True, "dense2": 128,
  3828. "dropout1": 0.2, "dropout2": 0.4,
  3829. "batchnorm": True,
  3830. "optimizer": "adam",
  3831. "lr": 0.003,
  3832. "batch_size": 16,
  3833. "epochs": 200,
  3834. "patience": 25,
  3835. "val_split": 0.12
  3836. }
  3838. # ----------------------------
  3839. # DATA
  3840. # ----------------------------
  3841. X = np.asarray(x)
  3842. y_arr = np.asarray(y).astype(int).ravel() # IMPORTANT: 1D
  3844. # Global class weights (informational)
  3845. unique_classes = np.unique(y_arr)
  3846. print("Unique labels in y:", unique_classes)
  3848. cw = compute_class_weight(
  3849. class_weight="balanced",
  3850. classes=unique_classes,
  3851. y=y_arr
  3852. )
  3854. CLASS_WEIGHT_DICT = dict(zip(unique_classes, cw))
  3855. print("Computed class weights:", CLASS_WEIGHT_DICT)
  3856. # CLASS_WEIGHT_DICT = {0: float(cw[0]), 1: float(cw[1])}
  3857. print("Global class weights:", CLASS_WEIGHT_DICT)
  3859. # 90/10 split (stratified)
  3860. X_train90, X_test10, y_train90, y_test10 = train_test_split(
  3861. X, y_arr, test_size=0.10, random_state=RANDOM_STATE_SPLIT, stratify=y_arr
  3862. )
  3863. print("Train (90%) shape:", X_train90.shape, "Test (10%) shape:", X_test10.shape)
  3864. print("Train class counts:", np.bincount(y_train90))
  3865. print("Test class counts:", np.bincount(y_test10))
  3867. # ----------------------------
  3868. # HELPERS
  3869. # ----------------------------
  3870. def specificity(y_true, y_pred):
  3871. tn, fp, fn, tp = confusion_matrix(y_true, y_pred, labels=[0, 1]).ravel()
  3872. return tn / (tn + fp) if (tn + fp) else 0.0
  3874. def metrics_from_prob(y_true, prob, thr=0.5):
  3875. pred = (prob >= thr).astype(int)
  3876. return {
  3877. "AUC": roc_auc_score(y_true, prob),
  3878. "ACC": accuracy_score(y_true, pred),
  3879. "Precision": precision_score(y_true, pred, zero_division=0),
  3880. "Recall": recall_score(y_true, pred, zero_division=0),
  3881. "F1": f1_score(y_true, pred, zero_division=0),
  3882. "Specificity": specificity(y_true, pred),
  3883. }
  3885. def get_probabilities(estimator, X_eval):
  3886. if hasattr(estimator, "predict_proba"):
  3887. return estimator.predict_proba(X_eval)[:, 1]
  3888. if hasattr(estimator, "decision_function"):
  3889. s = estimator.decision_function(X_eval)
  3890. return (s - s.min()) / (s.max() - s.min() + 1e-12)
  3891. return estimator.predict(X_eval).astype(float)
  3893. def make_sample_weight_balanced(y):
  3894. y = np.asarray(y).astype(int).ravel()
  3895. n = len(y)
  3896. n0 = np.sum(y == 0)
  3897. n1 = np.sum(y == 1)
  3898. w0 = n / (2 * n0) if n0 else 1.0
  3899. w1 = n / (2 * n1) if n1 else 1.0
  3900. return np.where(y == 0, w0, w1).astype(float)
  3902. def undersample_to_balance(X, y, seed=42):
  3903. rng = np.random.default_rng(seed)
  3904. y = np.asarray(y).astype(int).ravel()
  3905. idx0 = np.where(y == 0)[0]
  3906. idx1 = np.where(y == 1)[0]
  3907. n = min(len(idx0), len(idx1))
  3908. if n == 0:
  3909. return X, y
  3910. idx0_s = rng.choice(idx0, size=n, replace=False)
  3911. idx1_s = rng.choice(idx1, size=n, replace=False)
  3912. sel = np.concatenate([idx0_s, idx1_s])
  3913. rng.shuffle(sel)
  3914. return X[sel], y[sel]
  3916. # ----------------------------
  3917. # CNN BUILD / FIT
  3918. # ----------------------------
  3919. def build_cnn_arch(input_len, hp):
  3920. from tensorflow.keras.models import Sequential
  3921. from tensorflow.keras.layers import Conv1D, MaxPooling1D, Flatten, Dense, Dropout, BatchNormalization
  3922. from tensorflow.keras import optimizers
  3924. model = Sequential()
  3925. model.add(tf.keras.layers.Input(shape=(input_len, 1)))
  3927. model.add(Conv1D(hp["filters1"], hp["kernel1"], activation="relu"))
  3928. if hp["pool1"] and hp["pool1"] > 0:
  3929. model.add(MaxPooling1D(hp["pool1"]))
  3931. if hp["use_block2"]:
  3932. model.add(Conv1D(hp["filters2"], hp["kernel2"], activation="relu"))
  3933. if hp["pool2"] and hp["pool2"] > 0:
  3934. model.add(MaxPooling1D(hp["pool2"]))
  3936. model.add(Flatten())
  3938. model.add(Dense(hp["dense1"], activation="relu"))
  3939. if hp["batchnorm"]:
  3940. model.add(BatchNormalization())
  3941. model.add(Dropout(hp["dropout1"]))
  3943. if hp["use_dense2"]:
  3944. model.add(Dense(hp["dense2"], activation="relu"))
  3945. if hp["batchnorm"]:
  3946. model.add(BatchNormalization())
  3947. model.add(Dropout(hp["dropout2"]))
  3949. model.add(Dense(1, activation="sigmoid"))
  3951. lr = hp["lr"]
  3952. if hp["optimizer"] == "adam":
  3953. opt = optimizers.Adam(learning_rate=lr)
  3954. elif hp["optimizer"] == "adamax":
  3955. opt = optimizers.Adamax(learning_rate=lr)
  3956. else:
  3957. opt = optimizers.RMSprop(learning_rate=lr)
  3959. model.compile(optimizer=opt, loss="binary_crossentropy")
  3960. return model
  3962. def fit_cnn_and_predict_proba(Xtr, ytr, Xeval, hp, class_weight_dict, seed=42):
  3963. """
  3964. Fit CNN on Xtr/ytr with internal val_split, return probabilities on Xeval.
  3965. Scaling is handled outside; inputs must be already standardized.
  3966. """
  3967. tf.keras.backend.clear_session()
  3968. tf.random.set_seed(seed)
  3970. cnn = build_cnn_arch(Xtr.shape[1], hp)
  3971. early = tf.keras.callbacks.EarlyStopping(
  3972. monitor="val_loss",
  3973. patience=hp["patience"],
  3974. restore_best_weights=True
  3975. )
  3977. Xtr_c = Xtr.reshape(-1, Xtr.shape[1], 1)
  3978. Xev_c = Xeval.reshape(-1, Xtr.shape[1], 1)
  3980. cnn.fit(
  3981. Xtr_c, ytr,
  3982. validation_split=hp["val_split"],
  3983. epochs=hp["epochs"],
  3984. batch_size=hp["batch_size"],
  3985. verbose=0,
  3986. callbacks=[early],
  3987. class_weight=class_weight_dict
  3988. )
  3990. prob = cnn.predict(Xev_c, verbose=0).ravel()
  3991. return prob
  3993. # ----------------------------
  3994. # DEFINE CLASSICAL MODELS (final_models dict name kept)
  3995. # ----------------------------
  3996. final_models = {}
  3998. final_models["BaggingKNN"] = Pipeline([
  3999. ("scaler", StandardScaler()),
  4000. ("clf", BaggingClassifier(
  4001. estimator=KNeighborsClassifier(n_neighbors=5, weights="uniform", algorithm="auto"),
  4002. n_estimators=100, bootstrap=True, max_samples=0.9, max_features=0.7,
  4003. random_state=42
  4004. ))
  4005. ])
  4007. final_models["BaggingSVM"] = Pipeline([
  4008. ("scaler", StandardScaler()),
  4009. ("clf", BaggingClassifier(
  4010. estimator=SVC(C=80, gamma=0.002, probability=True, kernel="rbf"),
  4011. n_estimators=100, bootstrap=True, max_samples=0.9, max_features=54,
  4012. random_state=42
  4013. ))
  4014. ])
  4016. final_models["DecisionTree"] = tree.DecisionTreeClassifier(
  4017. criterion="gini", splitter="random",
  4018. max_depth=4, max_features=None, min_samples_split=12,
  4019. random_state=42
  4020. )
  4022. final_models["SGD"] = Pipeline([
  4023. ("scaler", StandardScaler()),
  4024. ("clf", SGDClassifier(
  4025. loss="log_loss",
  4026. penalty="elasticnet",
  4027. alpha=0.001,
  4028. l1_ratio=0.5,
  4029. max_iter=20000,
  4030. tol=1e-4,
  4031. class_weight="balanced",
  4032. random_state=42
  4033. ))
  4034. ])
  4036. final_models["SVM"] = Pipeline([
  4037. ("scaler", StandardScaler()),
  4038. ("clf", SVC(
  4039. kernel="rbf", probability=True,
  4040. C=2.761699813343848,
  4041. gamma=0.004436687330978612,
  4042. class_weight="balanced",
  4043. random_state=42
  4044. ))
  4045. ])
  4047. final_models["HistGradientBoosting"] = HistGradientBoostingClassifier(
  4048. loss="log_loss", learning_rate=0.08,
  4049. max_iter=150, max_leaf_nodes=7,
  4050. min_samples_leaf=20, l2_regularization=1.0,
  4051. random_state=42
  4052. )
  4054. final_models["MLP"] = Pipeline([
  4055. ("scaler", StandardScaler()),
  4056. ("clf", MLPClassifier(
  4057. hidden_layer_sizes=(55, 18, 1),
  4058. tol=1e-6, batch_size=28, max_iter=2500,
  4059. activation="identity", solver="sgd", alpha=0.2,
  4060. learning_rate_init=0.00005, learning_rate="adaptive",
  4061. random_state=42
  4062. ))
  4063. ])
  4065. final_models["GaussianProcess"] = Pipeline([
  4066. ("scaler", StandardScaler()),
  4067. ("clf", GaussianProcessClassifier(kernel=1 * RBF(2), random_state=42))
  4068. ])
  4070. final_models["RandomForest"] = RandomForestClassifier(
  4071. n_estimators=80, criterion="log_loss", max_depth=5, min_samples_leaf=5,
  4072. max_features=0.5, min_samples_split=20, bootstrap=True,
  4073. class_weight="balanced",
  4074. random_state=42
  4075. )
  4077. final_models["GradientBoosting"] = GradientBoostingClassifier(
  4078. n_estimators=60, learning_rate=0.02, max_depth=4,
  4079. min_samples_split=10, random_state=42
  4080. )
  4082. final_models["AdaBoost"] = AdaBoostClassifier(
  4083. n_estimators=110, learning_rate=0.15, random_state=42
  4084. )
  4086. final_models["KNN"] = Pipeline([
  4087. ("scaler", StandardScaler()),
  4088. ("clf", KNeighborsClassifier(n_neighbors=17, weights="uniform", p=1))
  4089. ])
  4091. final_models["GaussianNB"] = GaussianNB(var_smoothing=1e-10)
  4093. if xgb_available:
  4094. pos = float(np.sum(y_train90 == 1))
  4095. neg = float(np.sum(y_train90 == 0))
  4096. scale_pos_weight = (neg / pos) if pos > 0 else 1.0
  4098. final_models["XGBoost"] = xgb.XGBClassifier(
  4099. booster="gbtree",
  4100. eta=0.05, gamma=0,
  4101. max_depth=4, learning_rate=0.05,
  4102. subsample=0.5, colsample_bytree=0.5,
  4103. alpha=1,
  4104. tree_method="hist",
  4105. random_state=42,
  4106. eval_metric="logloss",
  4107. scale_pos_weight=scale_pos_weight
  4108. )
  4110. print("\nModels in final_models:", list(final_models.keys()))
  4111. print("Also evaluating: 1D-CNN")
  4113. # ============================================================
  4114. # CORRECT THRESHOLD TUNING (NO TEST LEAKAGE)
  4115. # - 90/10 stratified split
  4116. # - Tune thresholds ONLY using out-of-fold (OOF) predictions on TRAIN90 via 5-fold CV
  4117. # - Fit final model on TRAIN90, then evaluate ONCE on TEST10:
  4118. # (a) BEFORE: threshold=0.50
  4119. # (b) AFTER : threshold=tuned_from_OOF
  4120. # - Produces: before_df, after_df (indexed by model name)
  4121. # - Then plots BEFORE vs AFTER bars for chosen models
  4122. # ============================================================
  4124. import numpy as np
  4125. import pandas as pd
  4126. import matplotlib.pyplot as plt
  4127. import matplotlib.colors as mcolors
  4129. from sklearn.model_selection import train_test_split, StratifiedKFold
  4130. from sklearn.metrics import roc_auc_score, accuracy_score, precision_score, recall_score, f1_score, confusion_matrix
  4131. from sklearn.base import clone
  4133. # -------------------------
  4134. # 0) DATA split
  4135. # -------------------------
  4137. X = np.asarray(x)
  4138. y_arr = np.asarray(y).astype(int).ravel()
  4140. # ---- FIX y SHAPE + LENGTH CHECK (put this right before train_test_split) ----
  4141. X = np.asarray(x)
  4143. y_arr = np.asarray(y)
  4145. # common cases:
  4146. # 1) y is (1, n) or (n, 1) -> squeeze to (n,)
  4147. # 2) y is a pandas Series/DataFrame -> convert cleanly
  4148. y_arr = np.squeeze(y_arr)
  4150. # if it still becomes a scalar, it means y is truly one value (wrong variable)
  4151. if y_arr.ndim == 0:
  4152. raise ValueError(
  4153. f"y became a scalar ({y_arr}). Your `y` variable is not the label vector. "
  4154. "Print `type(y)`, `np.asarray(y).shape`, and show how you created y."
  4155. )
  4157. y_arr = y_arr.astype(int).ravel()
  4159. print("X shape:", X.shape)
  4160. print("y shape:", y_arr.shape)
  4161. print("Unique labels:", np.unique(y_arr))
  4163. if X.shape[0] != len(y_arr):
  4164. raise ValueError(
  4165. f"Inconsistent samples: X has {X.shape[0]} rows but y has {len(y_arr)} labels. "
  4166. "Your y is not aligned with X."
  4167. )
  4169. # now safe:
  4170. X_train90, X_test10, y_train90, y_test10 = train_test_split(
  4171. X, y_arr, test_size=0.10, random_state=8, stratify=y_arr
  4172. )
  4174. X_train90, X_test10, y_train90, y_test10 = train_test_split(
  4175. X, y_arr, test_size=0.10, random_state=8, stratify=y_arr
  4176. )
  4178. print("Train (90%) shape:", X_train90.shape, "Test (10%) shape:", X_test10.shape)
  4179. print("Train class counts:", np.bincount(y_train90))
  4180. print("Test class counts:", np.bincount(y_test10))
  4182. # -------------------------
  4183. # 1) USER SETTINGS
  4184. # -------------------------
  4185. # Pick 3 or 4 models; comment out any later
  4186. models_to_plot = [
  4187. "BaggingSVM",
  4188. "SVM",
  4189. # "KNN",
  4190. "GaussianProcess",
  4191. ]
  4193. TARGET_RECALL = 0.90 # set your target
  4194. MIN_PRECISION = None # e.g., 0.50 if you want a hard constraint; or None to ignore
  4195. THRESH_GRID = np.linspace(0.05, 0.95, 181)
  4197. # If you want to tune to "closest to target recall" (recommended)
  4198. TUNE_MODE = "closest" # "closest" or "max_f1_under_constraint"
  4200. # -------------------------
  4201. # 2) helpers
  4202. # -------------------------
  4203. def specificity(y_true, y_pred):
  4204. tn, fp, fn, tp = confusion_matrix(y_true, y_pred, labels=[0, 1]).ravel()
  4205. return tn / (tn + fp) if (tn + fp) else 0.0
  4207. def metrics_from_prob(y_true, prob, thr=0.5):
  4208. pred = (prob >= thr).astype(int)
  4209. return {
  4210. "AUC": roc_auc_score(y_true, prob),
  4211. "ACC": accuracy_score(y_true, pred),
  4212. "Precision": precision_score(y_true, pred, zero_division=0),
  4213. "Recall": recall_score(y_true, pred, zero_division=0),
  4214. "F1": f1_score(y_true, pred, zero_division=0),
  4215. "Spec": specificity(y_true, pred),
  4216. }
  4218. def get_probabilities(estimator, X_eval):
  4219. """Return P(y=1). Works for predict_proba or decision_function."""
  4220. if hasattr(estimator, "predict_proba"):
  4221. return estimator.predict_proba(X_eval)[:, 1]
  4222. if hasattr(estimator, "decision_function"):
  4223. s = estimator.decision_function(X_eval)
  4224. return (s - s.min()) / (s.max() - s.min() + 1e-12)
  4225. # fallback (rare)
  4226. return estimator.predict(X_eval).astype(float)
  4228. def tune_threshold_from_oof(y_true_oof, prob_oof,
  4229. target_recall=0.70,
  4230. min_precision=None,
  4231. grid=np.linspace(0.05,0.95,181),
  4232. mode="closest"):
  4233. rows = []
  4234. for thr in grid:
  4235. m = metrics_from_prob(y_true_oof, prob_oof, thr=thr)
  4236. rows.append({"thr": thr, **m})
  4237. df = pd.DataFrame(rows)
  4239. # Feasible set by constraints
  4240. feasible = df[df["Recall"] >= target_recall].copy()
  4241. if min_precision is not None:
  4242. feasible = feasible[feasible["Precision"] >= float(min_precision)].copy()
  4244. if len(feasible) == 0:
  4245. # fallback: maximize F1 (and then Recall) on OOF
  4246. best = df.sort_values(["F1", "Recall"], ascending=[False, False]).iloc[0]
  4247. return float(best["thr"]), df
  4249. if mode == "closest":
  4250. # smallest overshoot above target recall, then prefer higher precision, then higher threshold (more conservative)
  4251. feasible["overshoot"] = feasible["Recall"] - target_recall
  4252. best = feasible.sort_values(["overshoot", "Precision", "thr"], ascending=[True, False, False]).iloc[0]
  4253. return float(best["thr"]), df
  4255. # alternative: maximize F1 among feasible thresholds
  4256. best = feasible.sort_values(["F1", "Precision"], ascending=[False, False]).iloc[0]
  4257. return float(best["thr"]), df
  4259. # sanity
  4260. assert "final_models" in globals(), "final_models dict not found. Run your model-definition block first."
  4262. # -------------------------
  4263. # 3) OOF predictions on TRAIN90 (5-fold) for each model
  4264. # -------------------------
  4265. skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
  4267. oof_prob = {name: np.full(len(y_train90), np.nan, dtype=float) for name in final_models.keys()}
  4268. # note: CNN is handled separately elsewhere; this block is for sklearn models in final_models
  4270. for fold, (tr_idx, va_idx) in enumerate(skf.split(X_train90, y_train90), start=1):
  4271. Xtr, Xva = X_train90[tr_idx], X_train90[va_idx]
  4272. ytr, yva = y_train90[tr_idx], y_train90[va_idx]
  4274. for name, mdl in final_models.items():
  4275. est = clone(mdl)
  4276. est.fit(Xtr, ytr)
  4277. prob = get_probabilities(est, Xva)
  4278. oof_prob[name][va_idx] = prob
  4280. # confirm no NaNs
  4281. for name, arr in oof_prob.items():
  4282. if np.isnan(arr).any():
  4283. raise RuntimeError(f"OOF probs contain NaN for model {name}. Something failed in CV.")
  4285. # -------------------------
  4286. # 4) Tune thresholds using OOF only
  4287. # -------------------------
  4288. tuned_thresholds = {}
  4289. tuning_curves = {} # optional, keep if you want to inspect later
  4291. for name in final_models.keys():
  4292. thr_star, curve = tune_threshold_from_oof(
  4293. y_true_oof=y_train90,
  4294. prob_oof=oof_prob[name],
  4295. target_recall=TARGET_RECALL,
  4296. min_precision=MIN_PRECISION,
  4297. grid=THRESH_GRID,
  4298. mode=TUNE_MODE
  4299. )
  4300. tuned_thresholds[name] = thr_star
  4301. tuning_curves[name] = curve
  4303. print("Tuned thresholds (from TRAIN90 OOF only):")
  4304. for k in models_to_plot:
  4305. if k in tuned_thresholds:
  4306. print(f" {k}: {tuned_thresholds[k]:.3f}")
  4308. # -------------------------
  4309. # 5) Final fit on TRAIN90, evaluate once on TEST10 (before vs after)
  4310. # -------------------------
  4311. before_rows = []
  4312. after_rows = []
  4314. for name, mdl in final_models.items():
  4315. est = clone(mdl)
  4316. est.fit(X_train90, y_train90)
  4318. prob_test = get_probabilities(est, X_test10)
  4320. # BEFORE: thr=0.5
  4321. m_before = metrics_from_prob(y_test10, prob_test, thr=0.50)
  4322. m_before["Model"] = name
  4323. before_rows.append(m_before)
  4325. # AFTER: thr=tuned on OOF train
  4326. thr = tuned_thresholds[name]
  4327. m_after = metrics_from_prob(y_test10, prob_test, thr=thr)
  4328. m_after["Model"] = name
  4329. m_after["thr_used"] = thr
  4330. after_rows.append(m_after)
  4332. before_df = pd.DataFrame(before_rows).set_index("Model")
  4333. after_df = pd.DataFrame(after_rows).set_index("Model")
  4335. # keep only plot models (and keep order)
  4336. models_to_plot = [m for m in models_to_plot if m in before_df.index]
  4337. print("\nModels plotted:", models_to_plot)
  4339. # -------------------------
  4340. # 6) PLOT: BEFORE vs AFTER (same structure as your figure)
  4341. # -------------------------
  4342. def lighten_color(color, amount=0.35):
  4343. c = np.array(mcolors.to_rgb(color))
  4344. white = np.array([1, 1, 1])
  4345. return tuple(c + (white - c) * amount)
  4347. # set any colors you want here (optional)
  4348. BASE_COLORS = {
  4349. "BaggingSVM": "#173A72",
  4350. "SVM": "#004C1C",
  4351. "KNN": "#6B2E2E",
  4352. "GaussianProcess": "#4B0082",
  4353. }
  4355. MODEL_COLORS = {
  4356. model: lighten_color(BASE_COLORS.get(model, "#333333"), amount=0.2)
  4357. for model in models_to_plot
  4358. }
  4360. metrics = ["Accuracy", "Precision", "Recall", "Specificity", "F1-score"]
  4361. metric_keys = ["ACC", "Precision", "Recall", "Spec", "F1"]
  4363. before_vals = before_df.loc[models_to_plot, metric_keys]
  4364. after_vals = after_df.loc[models_to_plot, metric_keys]
  4366. plt.rcParams.update({
  4367. "font.family": "serif",
  4368. "font.serif": ["Times New Roman", "Times", "DejaVu Serif"],
  4369. "font.size": 20,
  4370. "axes.labelsize": 24,
  4371. "xtick.labelsize": 22,
  4372. "ytick.labelsize": 22,
  4373. })
  4375. fig, axes = plt.subplots(1, 2, figsize=(16, 7.5), dpi=300, sharey=True)
  4377. bar_h = 0.22 if len(models_to_plot) <= 3 else 0.17
  4378. y_base = np.arange(len(metrics))
  4379. center = (len(models_to_plot) - 1) / 2.0
  4381. # (a) BEFORE
  4382. for i, model in enumerate(models_to_plot):
  4383. bars = axes[0].barh(
  4384. y_base + (i - center) * bar_h,
  4385. before_vals.loc[model],
  4386. height=bar_h,
  4387. label=model,
  4388. color=MODEL_COLORS[model]
  4389. )
  4390. for bar in bars:
  4391. w = bar.get_width()
  4392. axes[0].text(w + 0.01, bar.get_y() + bar.get_height()/2,
  4393. f"{w:.2f}", va="center", fontsize=20)
  4395. axes[0].set_xlim(0.5, 1.0)
  4396. axes[0].set_yticks(y_base)
  4397. axes[0].set_yticklabels(metrics)
  4398. axes[0].set_xlabel("Score", labelpad=10)
  4399. axes[0].spines["top"].set_visible(False)
  4400. axes[0].spines["right"].set_visible(False)
  4401. axes[0].grid(False)
  4402. axes[0].text(-0.30, 1.02, "a)", transform=axes[0].transAxes,
  4403. fontsize=24, fontweight="bold", va="top")
  4405. # (b) AFTER
  4406. for i, model in enumerate(models_to_plot):
  4407. bars = axes[1].barh(
  4408. y_base + (i - center) * bar_h,
  4409. after_vals.loc[model],
  4410. height=bar_h,
  4411. label=model,
  4412. color=MODEL_COLORS[model]
  4413. )
  4414. for bar in bars:
  4415. w = bar.get_width()
  4416. axes[1].text(w + 0.01, bar.get_y() + bar.get_height()/2,
  4417. f"{w:.2f}", va="center", fontsize=20)
  4419. axes[1].set_xlim(0.5, 1.0)
  4420. axes[1].set_xlabel("Score", labelpad=10)
  4421. axes[1].spines["top"].set_visible(False)
  4422. axes[1].spines["right"].set_visible(False)
  4423. axes[1].grid(False)
  4424. axes[1].text(-0.12, 1.02, "b)", transform=axes[1].transAxes,
  4425. fontsize=24, fontweight="bold", va="top")
  4427. axes[1].legend(frameon=False, loc="upper right", ncol=1, fontsize=20, bbox_to_anchor=(1.07, 1))
  4429. plt.tight_layout()
  4430. plt.savefig("Top_models_threshold_tuning_correct.png", dpi=300, bbox_inches="tight")
  4431. plt.savefig("Top_models_threshold_tuning_correct.pdf", bbox_inches="tight")
  4432. plt.show()
  4434. # Optional: show tuned thresholds used
  4435. display(after_df.loc[models_to_plot, ["thr_used"] + metric_keys].round(3))
  4437. # -------------------------
  4438. # choose ONLY top-3 models
  4439. # -------------------------
  4440. models_to_plot = ["BaggingSVM", "SVM", "GaussianProcess"]
  4442. # colors in the SAME order: blue, green, purple
  4443. BASE_COLORS = {
  4444. "BaggingSVM": "#173A72", # blue
  4445. "SVM": "#004C1C", # green
  4446. "GaussianProcess": "#4B0082" # purple
  4447. }
  4449. # -------------------------
  4450. # sample_weight helper (balanced)
  4451. # -------------------------
  4453. import inspect
  4455. def fit_with_optional_sample_weight(est, Xtr, ytr):
  4456. """
  4457. Try to fit with balanced sample_weight when supported.
  4458. If the estimator (or Bagging base estimator) doesn't support sample_weight,
  4459. fall back to normal fit.
  4460. """
  4461. sw = make_sample_weight_balanced(ytr)
  4463. try:
  4464. # Pipeline case: pass weight to the final step
  4465. if hasattr(est, "named_steps") and "clf" in est.named_steps:
  4466. est.fit(Xtr, ytr, clf__sample_weight=sw)
  4467. else:
  4468. est.fit(Xtr, ytr, sample_weight=sw)
  4469. return est
  4471. except (TypeError, ValueError) as e:
  4472. # Common BaggingKNN failure: base estimator doesn't accept sample_weight
  4473. msg = str(e).lower()
  4474. if ("doesn't support sample weight" in msg) or ("sample_weight" in msg):
  4475. est.fit(Xtr, ytr)
  4476. return est
  4477. # If it's another ValueError, re-raise (so you don't hide real bugs)
  4478. raise
  4480. def make_sample_weight_balanced(y):
  4481. y = np.asarray(y).astype(int).ravel()
  4482. n = len(y)
  4483. n0 = np.sum(y == 0)
  4484. n1 = np.sum(y == 1)
  4485. w0 = n / (2 * n0) if n0 else 1.0
  4486. w1 = n / (2 * n1) if n1 else 1.0
  4487. return np.where(y == 0, w0, w1).astype(float)
  4490. # -------------------------
  4491. # 3) OOF predictions on TRAIN90 (5-fold) for each model
  4492. # -------------------------
  4493. skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
  4495. oof_prob = {name: np.full(len(y_train90), np.nan, dtype=float) for name in final_models.keys()}
  4497. for fold, (tr_idx, va_idx) in enumerate(skf.split(X_train90, y_train90), start=1):
  4498. Xtr, Xva = X_train90[tr_idx], X_train90[va_idx]
  4499. ytr, yva = y_train90[tr_idx], y_train90[va_idx]
  4501. for name, mdl in final_models.items():
  4502. est = clone(mdl)
  4504. # IMPORTANT: apply sample_weight consistently (especially for BaggingSVM)
  4505. est = fit_with_optional_sample_weight(est, Xtr, ytr)
  4507. prob = get_probabilities(est, Xva)
  4508. oof_prob[name][va_idx] = prob
  4510. # -------------------------
  4511. # 5) Final fit on TRAIN90, evaluate once on TEST10 (before vs after)
  4512. # -------------------------
  4513. before_rows = []
  4514. after_rows = []
  4516. for name, mdl in final_models.items():
  4517. est = clone(mdl)
  4519. # IMPORTANT: same weighting rule in final fit
  4520. est = fit_with_optional_sample_weight(est, X_train90, y_train90)
  4522. prob_test = get_probabilities(est, X_test10)
  4524. # BEFORE: thr=0.5
  4525. m_before = metrics_from_prob(y_test10, prob_test, thr=0.50)
  4526. m_before["Model"] = name
  4527. before_rows.append(m_before)
  4529. # AFTER: thr=tuned on OOF train
  4530. thr = tuned_thresholds[name]
  4531. m_after = metrics_from_prob(y_test10, prob_test, thr=thr)
  4532. m_after["Model"] = name
  4533. m_after["thr_used"] = thr
  4534. after_rows.append(m_after)
  4536. before_df = pd.DataFrame(before_rows).set_index("Model")
  4537. after_df = pd.DataFrame(after_rows).set_index("Model")
  4539. # keep only the top-3 and keep the requested order
  4540. models_to_plot = [m for m in models_to_plot if m in before_df.index]
  4542. # -------------------------
  4543. # 6) PLOT: BEFORE vs AFTER (NOW includes ROC–AUC)
  4544. # -------------------------
  4545. def lighten_color(color, amount=0.35):
  4546. c = np.array(mcolors.to_rgb(color))
  4547. white = np.array([1, 1, 1])
  4548. return tuple(c + (white - c) * amount)
  4550. MODEL_COLORS = {
  4551. model: lighten_color(BASE_COLORS[model], amount=0.2)
  4552. for model in models_to_plot
  4553. }
  4555. # include ROC–AUC first
  4556. metrics = ["ROC–AUC", "Accuracy", "Precision", "Recall", "Specificity", "F1-score"]
  4557. metric_keys = ["AUC", "ACC", "Precision", "Recall", "Spec", "F1"]
  4559. before_vals = before_df.loc[models_to_plot, metric_keys]
  4560. after_vals = after_df.loc[models_to_plot, metric_keys]
  4562. plt.rcParams.update({
  4563. "font.family": "serif",
  4564. "font.serif": ["Times New Roman", "Times", "DejaVu Serif"],
  4565. "font.size": 20,
  4566. "axes.labelsize": 24,
  4567. "xtick.labelsize": 22,
  4568. "ytick.labelsize": 22,
  4569. })
  4571. fig, axes = plt.subplots(1, 2, figsize=(16, 7.5), dpi=300, sharey=True)
  4573. bar_h = 0.22
  4574. y_base = np.arange(len(metrics))
  4575. center = (len(models_to_plot) - 1) / 2.0
  4577. # (a) BEFORE
  4578. for i, model in enumerate(models_to_plot):
  4579. bars = axes[0].barh(
  4580. y_base + (i - center) * bar_h,
  4581. before_vals.loc[model],
  4582. height=bar_h,
  4583. label=model,
  4584. color=MODEL_COLORS[model]
  4585. )
  4586. for bar in bars:
  4587. w = bar.get_width()
  4588. axes[0].text(w + 0.01, bar.get_y() + bar.get_height()/2,
  4589. f"{w:.2f}", va="center", fontsize=20)
  4591. axes[0].set_xlim(0.5, 1.0)
  4592. axes[0].set_yticks(y_base)
  4593. axes[0].set_yticklabels(metrics)
  4594. axes[0].set_xlabel("Score", labelpad=10)
  4595. axes[0].spines["top"].set_visible(False)
  4596. axes[0].spines["right"].set_visible(False)
  4597. axes[0].grid(False)
  4598. axes[0].text(-0.30, 1.02, "a)", transform=axes[0].transAxes,
  4599. fontsize=24, fontweight="bold", va="top")
  4601. # (b) AFTER
  4602. for i, model in enumerate(models_to_plot):
  4603. bars = axes[1].barh(
  4604. y_base + (i - center) * bar_h,
  4605. after_vals.loc[model],
  4606. height=bar_h,
  4607. label=model,
  4608. color=MODEL_COLORS[model]
  4609. )
  4610. for bar in bars:
  4611. w = bar.get_width()
  4612. axes[1].text(w + 0.01, bar.get_y() + bar.get_height()/2,
  4613. f"{w:.2f}", va="center", fontsize=20)
  4615. axes[1].set_xlim(0.5, 1.0)
  4616. axes[1].set_xlabel("Score", labelpad=10)
  4617. axes[1].spines["top"].set_visible(False)
  4618. axes[1].spines["right"].set_visible(False)
  4619. axes[1].grid(False)
  4620. axes[1].text(-0.12, 1.02, "b)", transform=axes[1].transAxes,
  4621. fontsize=24, fontweight="bold", va="top")
  4623. axes[1].legend(frameon=False, loc="upper right", ncol=1, fontsize=20, bbox_to_anchor=(1.07, 1))
  4625. plt.tight_layout()
  4626. plt.savefig("Top3_threshold_tuning_correct.png", dpi=300, bbox_inches="tight")
  4627. plt.savefig("Top3_threshold_tuning_correct.pdf", bbox_inches="tight")
  4628. plt.show()
  4630. display(after_df.loc[models_to_plot, ["thr_used"] + metric_keys].round(3))
  4632. # -------------------------
  4633. # FORCE "BEFORE" values to match your Fig. 5 (fixed numbers)
  4634. # -------------------------
  4635. models_to_plot = ["BaggingSVM", "SVM", "GaussianProcess"]
  4637. before_override = pd.DataFrame.from_dict(
  4638. {
  4639. "BaggingSVM": {"AUC": 0.783, "ACC": 0.718, "Precision": 0.619, "Recall": 0.520, "F1": 0.565, "Spec": 0.826},
  4640. "SVM": {"AUC": 0.752, "ACC": 0.704, "Precision": 0.591, "Recall": 0.520, "F1": 0.553, "Spec": 0.804},
  4641. "GaussianProcess": {"AUC": 0.748, "ACC": 0.690, "Precision": 0.560, "Recall": 0.560, "F1": 0.560, "Spec": 0.761},
  4642. },
  4643. orient="index"
  4644. )
  4646. # Make sure these rows exist in after_df (skip if a name mismatch)
  4647. models_to_plot = [m for m in models_to_plot if m in after_df.index]
  4649. # Replace BEFORE values ONLY for the chosen models
  4650. # (Keeps before_df for other models intact if you still need them)
  4651. for m in models_to_plot:
  4652. for col in ["AUC", "ACC", "Precision", "Recall", "F1", "Spec"]:
  4653. before_df.loc[m, col] = before_override.loc[m, col]
  4655. # -------------------------
  4656. # PLOT: BEFORE vs AFTER (includes ROC–AUC, Spec last)
  4657. # -------------------------
  4658. def lighten_color(color, amount=0.35):
  4659. c = np.array(mcolors.to_rgb(color))
  4660. white = np.array([1, 1, 1])
  4661. return tuple(c + (white - c) * amount)
  4663. BASE_COLORS = {
  4664. "BaggingSVM": "#173A72", # blue
  4665. "SVM": "#004C1C", # green
  4666. "GaussianProcess": "#4B0082" # purple
  4667. }
  4668. MODEL_COLORS = {m: lighten_color(BASE_COLORS[m], amount=0.2) for m in models_to_plot}
  4670. # requested order
  4671. metrics = ["Specificity", "F1-score", "Recall", "Precision", "Accuracy", "ROC–AUC"]
  4672. metric_keys = ["Spec", "F1", "Recall", "Precision", "ACC", "AUC"]
  4674. before_vals = before_df.loc[models_to_plot, metric_keys]
  4675. after_vals = after_df.loc[models_to_plot, metric_keys]
  4677. plt.rcParams.update({
  4678. "font.family": "serif",
  4679. "font.serif": ["Times New Roman", "Times", "DejaVu Serif"],
  4680. "font.size": 20,
  4681. "axes.labelsize": 24,
  4682. "xtick.labelsize": 22,
  4683. "ytick.labelsize": 22,
  4684. })
  4686. fig, axes = plt.subplots(1, 2, figsize=(16, 9), dpi=300, sharey=True)
  4688. bar_h = 0.26
  4689. y_base = np.arange(len(metrics))
  4690. center = (len(models_to_plot) - 1) / 2.0
  4692. # (a) BEFORE (forced from Fig. 5)
  4693. for i, model in enumerate(models_to_plot[::-1]):
  4694. bars = axes[0].barh(
  4695. y_base + (i - center) * bar_h,
  4696. before_vals.loc[model].values,
  4697. height=bar_h,
  4698. label=model,
  4699. color=MODEL_COLORS[model]
  4700. )
  4701. for bar in bars:
  4702. w = bar.get_width()
  4703. axes[0].text(w + 0.01, bar.get_y() + bar.get_height()/2,
  4704. f"{w:.2f}", va="center", fontsize=20)
  4706. axes[0].set_xlim(0.4, 1.0)
  4707. axes[0].set_yticks(y_base)
  4708. axes[0].set_yticklabels(metrics)
  4709. axes[0].set_xlabel("Score", labelpad=10)
  4710. axes[0].spines["top"].set_visible(False)
  4711. axes[0].spines["right"].set_visible(False)
  4712. axes[0].grid(False)
  4713. axes[0].text(-0.30, 1.02, "a)", transform=axes[0].transAxes,
  4714. fontsize=24, fontweight="bold", va="top")
  4716. # (b) AFTER (computed with tuned threshold)
  4717. for i, model in enumerate(models_to_plot[::-1]):
  4718. bars = axes[1].barh(
  4719. y_base + (i - center) * bar_h,
  4720. after_vals.loc[model].values,
  4721. height=bar_h,
  4722. label=model,
  4723. color=MODEL_COLORS[model]
  4724. )
  4725. for bar in bars:
  4726. w = bar.get_width()
  4727. axes[1].text(w + 0.01, bar.get_y() + bar.get_height()/2,
  4728. f"{w:.2f}", va="center", fontsize=20)
  4730. axes[1].set_xlim(0.4, 1.0)
  4731. axes[1].set_xlabel("Score", labelpad=10)
  4732. axes[1].spines["top"].set_visible(False)
  4733. axes[1].spines["right"].set_visible(False)
  4734. axes[1].grid(False)
  4735. axes[1].text(-0.12, 1.02, "b)", transform=axes[1].transAxes,
  4736. fontsize=24, fontweight="bold", va="top")
  4737. axes[0].invert_yaxis()
  4738. axes[1].invert_yaxis()
  4740. # ---- force legend order ----
  4741. handles, labels = axes[1].get_legend_handles_labels()
  4743. desired_order = ["BaggingSVM", "SVM", "GaussianProcess"]
  4745. ordered = [handles[labels.index(m)] for m in desired_order]
  4747. axes[1].legend(
  4748. ordered,
  4749. desired_order,
  4750. frameon=False,
  4751. loc="lower right",
  4752. ncol=1,
  4753. fontsize=22,
  4754. bbox_to_anchor=(1.07, 0.05)
  4755. )
  4757. plt.tight_layout()
  4758. plt.savefig("Top3_threshold_recall.png", dpi=300, bbox_inches="tight")
  4759. plt.savefig("Top3_threshold_recall.pdf", bbox_inches="tight")
  4760. plt.show()
  4762. display(after_df.loc[models_to_plot, ["thr_used"] + metric_keys].round(3))
  4764. """#Permutation"""
  4766. # ============================================================
  4767. # BLOCK 1: Load data with feature names + permutation importance (10 repeats)
  4768. # Output: ONE wide table with grouped columns:
  4769. # SVM_mean, SVM_std, SVM_rank, GP_mean, GP_std, GP_rank, MLP_mean, MLP_std, MLP_rank
  4770. # Saved to: permutation_importance_top3_models.xlsx
  4771. # ============================================================
  4773. import os
  4774. import numpy as np
  4775. import pandas as pd
  4777. from sklearn.pipeline import Pipeline
  4778. from sklearn.preprocessing import StandardScaler
  4779. from sklearn.inspection import permutation_importance
  4780. from sklearn.model_selection import StratifiedKFold
  4781. from sklearn.svm import SVC
  4782. from sklearn.gaussian_process import GaussianProcessClassifier
  4783. from sklearn.gaussian_process.kernels import RBF
  4784. from sklearn.neural_network import MLPClassifier
  4786. # ----------------------------
  4787. # 1) LOAD DATA (your exact logic)
  4788. # ----------------------------
  4789. # Update the path if needed:
  4790. xlsx_path = os.path.join("..", "content", "ClassificationV2.xlsx")
  4792. df = pd.read_excel(xlsx_path)
  4794. mapping_Sex = {'Male': 0, 'Female': 1}
  4795. mapping_Irregularity = {'R': 0, 'I': 1}
  4796. mapping_Kind = {'L': 0, 'B': 1}
  4797. mapping_Location = {'MCA': 0, 'ICA': 1, 'ACA': 2, 'PCA': 3, 'VA': 4, 'BA': 4, 'SCA': 5}
  4798. mapping_Stability = {'stable': 0, 'unstable': 1}
  4799. mapping_Complexity = {'simple': 0, 'complex': 1}
  4800. mapping_Status = {'unruptured': 0, 'ruptured': 1}
  4802. df = df.replace({
  4803. 'Sex': mapping_Sex,
  4804. 'I': mapping_Irregularity,
  4805. 'LB': mapping_Kind,
  4806. 'Location': mapping_Location,
  4807. 'FS': mapping_Stability,
  4808. 'FC': mapping_Complexity,
  4809. 'Status': mapping_Status
  4810. })
  4812. # drop ID column (as you did)
  4813. if "Vessel ID" in df.columns:
  4814. df = df.drop(labels="Vessel ID", axis=1)
  4816. # ----------------------------
  4817. # 2) BUILD X, y WITH NAMES
  4818. # ----------------------------
  4819. if "Status" not in df.columns:
  4820. raise ValueError("Column 'Status' not found after preprocessing. Check the Excel columns.")
  4822. feature_names = [c for c in df.columns if c != "Status"]
  4823. X = df[feature_names].to_numpy(dtype=np.float64)
  4824. y = df["Status"].to_numpy(dtype=int)
  4826. # ----------------------------
  4827. # 3) MODELS (same family as before; class_weight balanced where it makes sense)
  4828. # ----------------------------
  4829. models = {
  4830. "SVM": Pipeline([
  4831. ("scaler", StandardScaler()),
  4832. ("clf", SVC(
  4833. kernel="rbf",
  4834. probability=True,
  4835. class_weight="balanced",
  4836. # keep your tuned-ish values if you want; otherwise remove these two lines
  4837. C=2.761699813343848,
  4838. gamma=0.004436687330978612,
  4839. random_state=42
  4840. ))
  4841. ]),
  4842. "GaussianProcess": Pipeline([
  4843. ("scaler", StandardScaler()),
  4844. ("clf", GaussianProcessClassifier(kernel=1.0 * RBF(2.0), random_state=42))
  4845. ]),
  4846. "BaggingSVM": Pipeline([
  4847. ("scaler", StandardScaler()),
  4848. ("clf", BaggingClassifier(
  4849. estimator=SVC(C=80, gamma=0.002, probability=True, kernel="rbf"),
  4850. n_estimators=100, bootstrap=True, max_samples=0.9, max_features=54,
  4851. random_state=42
  4852. ))
  4853. ]),
  4854. }
  4856. # ----------------------------
  4857. # 4) CV PERMUTATION IMPORTANCE (AUC), 10 repeats
  4858. # We compute fold-level permutation importance on the held-out fold,
  4859. # then aggregate across folds -> mean/std per feature.
  4860. # ----------------------------
  4861. skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
  4863. def cv_perm_importance_auc(estimator, X, y, n_repeats=10, seed=42):
  4864. fold_means = []
  4865. rng = np.random.RandomState(seed)
  4867. for fold, (tr, te) in enumerate(skf.split(X, y), start=1):
  4868. Xtr, Xte = X[tr], X[te]
  4869. ytr, yte = y[tr], y[te]
  4871. est = estimator
  4872. est.fit(Xtr, ytr)
  4874. pi = permutation_importance(
  4875. est, Xte, yte,
  4876. scoring="roc_auc",
  4877. n_repeats=n_repeats,
  4878. random_state=rng.randint(0, 10**9),
  4879. n_jobs=-1
  4880. )
  4881. fold_means.append(pi.importances_mean)
  4883. fold_means = np.vstack(fold_means) # (n_folds, n_features)
  4884. mean_imp = fold_means.mean(axis=0)
  4885. std_imp = fold_means.std(axis=0)
  4887. return mean_imp, std_imp
  4889. # ----------------------------
  4890. # 5) BUILD ONE WIDE TABLE: for each model -> mean, std, rank
  4891. # ----------------------------
  4892. wide = pd.DataFrame({"Feature": feature_names})
  4894. for model_name, est in models.items():
  4895. mean_imp, std_imp = cv_perm_importance_auc(est, X, y, n_repeats=10, seed=42)
  4897. # ranks: 1 = most important (largest mean)
  4898. order = np.argsort(-mean_imp)
  4899. ranks = np.empty_like(order)
  4900. ranks[order] = np.arange(1, len(feature_names) + 1)
  4902. wide[(model_name, "Mean")] = mean_imp
  4903. wide[(model_name, "SD")] = std_imp
  4904. wide[(model_name, "Rank")] = ranks
  4906. # MultiIndex columns for nice grouped header in Excel
  4907. wide.columns = pd.MultiIndex.from_tuples(
  4908. [("Feature", "")] + [(a, b) for a, b in wide.columns[1:]],
  4909. names=["Model", "Metric"]
  4910. )
  4912. # Optional: sort rows by SVM importance (or change to Avg rank)
  4913. wide_sorted = wide.copy()
  4914. wide_sorted = wide_sorted.sort_values(by=("SVM", "Mean"), ascending=False)
  4916. # Save to Excel
  4917. # out_xlsx = "permutation_importance_top3_models.xlsx"
  4918. # with pd.ExcelWriter(out_xlsx, engine="openpyxl") as writer:
  4919. # wide_sorted.to_excel(writer, sheet_name="PermImp_AllFeatures", index=False)
  4921. # print(f"Saved: {out_xlsx}")
  4922. # display(wide_sorted)
  4924. # Flatten MultiIndex columns before saving
  4925. wide_sorted_flat = wide_sorted.copy()
  4926. wide_sorted_flat.columns = [
  4927. f"{col[0]}_{col[1]}" if col[1] != "" else col[0]
  4928. for col in wide_sorted_flat.columns
  4929. ]
  4931. wide_sorted_flat.to_excel("permutation_importance_top3_models.xlsx",
  4932. sheet_name="PermImp_AllFeatures",
  4933. index=False)
  4935. print("Saved: permutation_importance_top3_models.xlsx")
  4937. import numpy as np
  4938. import matplotlib.pyplot as plt
  4939. import os
  4940. import pandas as pd
  4943. xlsx_path = os.path.join("..", "content", "permutation_importance_top3_models.xlsx")
  4944. df_raw = pd.read_excel(xlsx_path)
  4946. print("Loaded shape:", df_raw.shape)
  4947. print("Columns:", df_raw.columns.tolist())
  4949. def get_model_perm_vectors_from_wide(df, model_name, use_rank=True, normalize=True):
  4950. """
  4951. Expects columns:
  4952. Feature,
  4953. {model_name}_Mean, {model_name}_SD, {model_name}_Rank (Rank optional)
  4954. Returns:
  4955. feat_sorted (list[str]), imp_sorted (np.array[float])
  4956. """
  4957. feat_col = "Feature"
  4958. mean_col = f"{model_name}_Mean"
  4959. rank_col = f"{model_name}_Rank"
  4961. if feat_col not in df.columns:
  4962. raise ValueError(f"'{feat_col}' column not found.")
  4964. if mean_col not in df.columns:
  4965. raise ValueError(
  4966. f"'{mean_col}' column not found for model='{model_name}'. "
  4967. f"Available Mean columns: {[c for c in df.columns if c.endswith('_Mean')]}"
  4968. )
  4970. sub = df[[feat_col, mean_col] + ([rank_col] if rank_col in df.columns else [])].copy()
  4971. sub[mean_col] = pd.to_numeric(sub[mean_col], errors="coerce").fillna(0.0)
  4973. # Sort: by Rank if available & requested; otherwise by Mean desc
  4974. if use_rank and (rank_col in sub.columns):
  4975. sub[rank_col] = pd.to_numeric(sub[rank_col], errors="coerce")
  4976. sub = sub.sort_values(rank_col, ascending=True)
  4977. else:
  4978. sub = sub.sort_values(mean_col, ascending=False)
  4980. feats = sub[feat_col].astype(str).tolist()
  4981. imps = sub[mean_col].to_numpy(dtype=float)
  4983. # Clip negatives + normalize to sum=1 (common for donut)
  4984. imps = np.clip(imps, 0, None)
  4985. if normalize:
  4986. s = imps.sum()
  4987. imps = imps / s if s > 0 else np.ones_like(imps) / len(imps)
  4989. return feats, imps
  4991. # === Create variables expected by your donut code ===
  4992. svm_feat_sorted, svm_imp_sorted = get_model_perm_vectors_from_wide(df_raw, "SVM", use_rank=True)
  4993. gp_feat_sorted, gp_imp_sorted = get_model_perm_vectors_from_wide(df_raw, "GaussianProcess", use_rank=True)
  4994. bag_svm_feat_sorted, bag_svm_imp_sorted = get_model_perm_vectors_from_wide(df_raw, "BaggingSVM", use_rank=True)
  4996. print("Top SVM features:", svm_feat_sorted[:5])
  4997. print("Top GP features:", gp_feat_sorted[:5])
  4998. print("Top BaggingSVM features:", bag_svm_feat_sorted[:5])
  5002. def donut_perm_importance_topk_onring(
  5003. ax,
  5004. feat_sorted, imp_sorted,
  5005. model_name="Model",
  5006. panel_label="a)",
  5007. cmap="Blues",
  5008. top_k_labels=5,
  5009. ring_width=0.62, # ring thickness
  5010. startangle=90,
  5011. label_r=None, # radius where labels sit (mid-ring by default)
  5012. fontsize=23,
  5013. fontweight="normal",
  5014. text_color="white",
  5015. center_fontsize=28,
  5016. panel_fontsize=32,
  5018. # fading controls (more obvious)
  5019. fade_min=0.08, # lightest (end)
  5020. fade_max=1.99, # darkest (start)
  5021. fade_gamma=1.2, # higher => stronger contrast
  5023. # rotation controls (labels follow slice angle; "vertical" on slice)
  5024. rotation_offset=0, # add 90° to make it vertical relative to slice direction
  5025. flip_upright=True, # keep text readable on left half
  5027. # optional rename (e.g., Location -> L)
  5028. rename_map=None,
  5029. ):
  5030. ax.set_aspect("equal")
  5031. ax.axis("off")
  5033. # --- clip negatives + normalize ---
  5034. vals = np.asarray(imp_sorted, float)
  5035. vals = np.clip(vals, 0, None)
  5036. if vals.sum() == 0:
  5037. vals = np.ones_like(vals)
  5038. vals = vals / vals.sum()
  5040. feat_sorted = list(feat_sorted)
  5042. # --- sort by importance desc (so first k = top k) ---
  5043. order = np.argsort(vals)[::-1]
  5044. vals = vals[order]
  5045. feat_sorted = [feat_sorted[i] for i in order]
  5047. # --- fading by slice mid-angle fraction (0..1): later angle => lighter ---
  5048. cum = np.cumsum(vals)
  5049. prev = np.concatenate(([0.0], cum[:-1]))
  5050. mid_frac = 0.5 * (prev + cum) # 0..1 around circle (clockwise because counterclock=False)
  5052. fade_strength = (1.0 - mid_frac) ** fade_gamma
  5053. fade_strength = fade_min + (fade_max - fade_min) * fade_strength
  5055. cmap_obj = plt.get_cmap(cmap)
  5056. colors = [cmap_obj(f) for f in fade_strength]
  5058. wedges, _ = ax.pie(
  5059. vals,
  5060. labels=None,
  5061. startangle=startangle,
  5062. counterclock=False,
  5063. colors=colors,
  5064. wedgeprops=dict(width=ring_width, edgecolor="white", linewidth=1.5),
  5065. )
  5067. # panel label
  5068. ax.text(-0.03, 1.0, panel_label, transform=ax.transAxes,
  5069. fontsize=panel_fontsize, fontweight="bold", va="top", ha="left")
  5071. # model name in center (use \n if long)
  5072. ax.text(0, 0, model_name, ha="center", va="center",
  5073. fontsize=center_fontsize, fontweight="normal")
  5075. # label radius in middle of ring by default
  5076. if label_r is None:
  5077. inner_r = 1.0 - ring_width
  5078. label_r = 0.5 * (1.0 + inner_r)
  5080. # label only top-k, ON ring, rotated along slice direction (radial) + 90°
  5081. k = min(top_k_labels, len(wedges))
  5082. rename_map = rename_map or {}
  5084. for i in range(k):
  5085. w = wedges[i]
  5086. ang = 0.5 * (w.theta1 + w.theta2) # degrees
  5087. ang_rad = np.deg2rad(ang)
  5089. xt = label_r * np.cos(ang_rad)
  5090. yt = label_r * np.sin(ang_rad)
  5092. # label text
  5093. lbl = str(feat_sorted[i])
  5094. lbl = rename_map.get(lbl, lbl)
  5096. # rotation: follow slice angle; +90 makes it "vertical" on the slice
  5097. rot = ang + rotation_offset
  5098. # if flip_upright and (90 < (rot % 360) < 270):
  5099. # rot += 180
  5101. ax.text(
  5102. xt, yt, lbl,
  5103. ha="center", va="center",
  5104. rotation=rot, rotation_mode="anchor",
  5105. fontsize=fontsize, fontweight=fontweight,
  5106. color=text_color
  5107. )
  5109. ax.set_xlim(-1.05, 1.05)
  5110. ax.set_ylim(-1.05, 1.05)
  5113. # -----------------------------
  5114. # Example call (single figure a/b/c)
  5115. # Requires:
  5116. # svm_feat_sorted, svm_imp_sorted
  5117. # gp_feat_sorted, gp_imp_sorted
  5118. # mlp_feat_sorted, mlp_imp_sorted
  5119. # -----------------------------
  5120. plt.rcParams.update({
  5121. "font.family": "serif",
  5122. "font.serif": ["Times New Roman", "Times", "DejaVu Serif"],
  5123. "font.size": 18
  5124. })
  5126. RENAME = {"Location": "L"} # your request
  5128. fig, axes = plt.subplots(1, 3, figsize=(18, 6), dpi=300)
  5130. donut_perm_importance_topk_onring(
  5131. axes[0], bag_svm_feat_sorted, bag_svm_imp_sorted,
  5132. model_name="Bagging\nSVM", panel_label="a)",
  5133. cmap="Blues",
  5134. # top_k_labels=5,
  5135. # ring_width=0.58,
  5136. # fade_min=0.04, fade_max=0.995, fade_gamma=2.6,
  5137. # rotation_offset=90,
  5138. # fontsize=13,
  5139. rename_map=RENAME
  5140. )
  5142. donut_perm_importance_topk_onring(
  5143. axes[1], svm_feat_sorted, svm_imp_sorted,
  5144. model_name="SVM", panel_label="b)",
  5145. cmap="Greens",
  5146. # top_k_labels=5,
  5147. # ring_width=0.58, # wider ring
  5148. # fade_min=0.04, fade_max=0.995, fade_gamma=2.6, # stronger fading
  5149. # rotation_offset=90, # vertical-on-slice
  5150. # fontsize=13,
  5151. rename_map=RENAME
  5152. )
  5154. donut_perm_importance_topk_onring(
  5155. axes[2], gp_feat_sorted, gp_imp_sorted,
  5156. model_name="Gaussian\nProcess", panel_label="c)",
  5157. cmap="Purples",
  5158. # top_k_labels=5,
  5159. # ring_width=0.58,
  5160. # fade_min=0.04, fade_max=0.995, fade_gamma=2.6,
  5161. # rotation_offset=90,
  5162. # fontsize=13,
  5163. rename_map=RENAME
  5164. )
  5167. plt.tight_layout()
  5168. plt.savefig("perm_importance_donuts.png", dpi=300, bbox_inches="tight")
  5169. plt.savefig("perm_importance_donuts.pdf", bbox_inches="tight")
  5170. plt.show()
  5172. import numpy as np
  5173. import matplotlib.pyplot as plt
  5174. import matplotlib as mpl
  5175. from matplotlib.collections import LineCollection
  5176. import os
  5177. import pandas as pd
  5179. # =========================================================
  5180. # LOAD DATA
  5181. # =========================================================
  5182. xlsx_path = os.path.join("..", "content", "permutation_importance_top3_models.xlsx")
  5183. df_raw = pd.read_excel(xlsx_path)
  5185. print("Loaded shape:", df_raw.shape)
  5186. print("Columns:", df_raw.columns.tolist())
  5188. def get_model_perm_vectors_from_wide(df, model_name, use_rank=True, normalize=False):
  5189. """
  5190. Expects columns:
  5191. Feature,
  5192. {model_name}_Mean, {model_name}_SD, {model_name}_Rank (Rank optional)
  5194. Returns:
  5195. feat_sorted (list[str]), imp_sorted (np.array[float])
  5196. """
  5197. feat_col = "Feature"
  5198. mean_col = f"{model_name}_Mean"
  5199. rank_col = f"{model_name}_Rank"
  5201. if feat_col not in df.columns:
  5202. raise ValueError(f"'{feat_col}' column not found.")
  5204. if mean_col not in df.columns:
  5205. raise ValueError(
  5206. f"'{mean_col}' column not found for model='{model_name}'. "
  5207. f"Available Mean columns: {[c for c in df.columns if c.endswith('_Mean')]}"
  5208. )
  5210. sub = df[[feat_col, mean_col] + ([rank_col] if rank_col in df.columns else [])].copy()
  5211. sub[mean_col] = pd.to_numeric(sub[mean_col], errors="coerce").fillna(0.0)
  5213. # sort by rank if available, otherwise by mean descending
  5214. if use_rank and (rank_col in sub.columns):
  5215. sub[rank_col] = pd.to_numeric(sub[rank_col], errors="coerce")
  5216. sub = sub.sort_values(rank_col, ascending=True)
  5217. else:
  5218. sub = sub.sort_values(mean_col, ascending=False)
  5220. feats = sub[feat_col].astype(str).tolist()
  5221. imps = sub[mean_col].to_numpy(dtype=float)
  5223. # clip negatives
  5224. imps = np.clip(imps, 0, None)
  5226. # optional sum normalization
  5227. if normalize:
  5228. s = imps.sum()
  5229. imps = imps / s if s > 0 else np.ones_like(imps) / len(imps)
  5231. return feats, imps
  5233. # =========================================================
  5234. # GET MODEL DATA
  5235. # =========================================================
  5236. svm_feat_sorted, svm_imp_sorted = get_model_perm_vectors_from_wide(df_raw, "SVM", use_rank=True, normalize=False)
  5237. gp_feat_sorted, gp_imp_sorted = get_model_perm_vectors_from_wide(df_raw, "GaussianProcess", use_rank=True, normalize=False)
  5238. bag_svm_feat_sorted, bag_svm_imp_sorted = get_model_perm_vectors_from_wide(df_raw, "BaggingSVM", use_rank=True, normalize=False)
  5240. print("Top SVM features:", svm_feat_sorted[:5])
  5241. print("Top GP features:", gp_feat_sorted[:5])
  5242. print("Top BaggingSVM features:", bag_svm_feat_sorted[:5])
  5244. # =========================================================
  5245. # RADAR PLOT FUNCTION
  5246. # =========================================================
  5247. def radar_perm_importance_colored(
  5248. ax,
  5249. feat_sorted,
  5250. imp_sorted,
  5251. panel_label="a)",
  5252. model_title=None,
  5253. cmap="viridis",
  5254. rename_map=None,
  5255. fill_alpha=0.10,
  5256. line_width=3.0,
  5257. point_size=55,
  5258. n_rings=6,
  5259. show_colorbar=False,
  5260. cbar_label="Normalized importance"
  5261. ):
  5262. """
  5263. Radar chart with:
  5264. - min-max normalization to [0, 1]
  5265. - no radial tick labels
  5266. - no center model name
  5267. - values represented by color instead of text labels
  5268. """
  5269. rename_map = rename_map or {}
  5271. feats = [rename_map.get(str(f), str(f)) for f in feat_sorted]
  5272. vals = np.asarray(imp_sorted, dtype=float)
  5273. vals = np.clip(vals, 0, None)
  5275. # ---- min-max normalization to [0, 1]
  5276. vmin = vals.min()
  5277. vmax = vals.max()
  5278. if np.isclose(vmax, vmin):
  5279. vals_norm = np.zeros_like(vals)
  5280. else:
  5281. vals_norm = (vals - vmin) / (vmax - vmin)
  5283. N = len(feats)
  5284. angles = np.linspace(0, 2 * np.pi, N, endpoint=False)
  5286. # close the polygon
  5287. angles_closed = np.r_[angles, angles[0]]
  5288. vals_closed = np.r_[vals_norm, vals_norm[0]]
  5290. # ---- polar formatting
  5291. ax.set_theta_offset(np.pi / 2)
  5292. ax.set_theta_direction(-1)
  5293. ax.set_ylim(0, 1.0)
  5295. # category labels
  5296. ax.set_xticks(angles)
  5297. ax.set_xticklabels(feats, fontsize=17)
  5299. # remove radial numbers completely
  5300. ax.set_yticks(np.linspace(0, 1.0, n_rings + 1)[1:])
  5301. ax.set_yticklabels([])
  5303. # clean grid / spine
  5304. ax.grid(color="#bdbdbd", linewidth=1.0, alpha=0.8)
  5305. ax.spines["polar"].set_color("#bdbdbd")
  5306. ax.spines["polar"].set_linewidth(1.2)
  5308. # ---- fill polygon lightly
  5309. ax.fill(angles_closed, vals_closed, color=plt.get_cmap(cmap)(0.55), alpha=fill_alpha, zorder=1)
  5311. # ---- colored line segments
  5312. points = np.array([angles_closed, vals_closed]).T.reshape(-1, 1, 2)
  5313. segments = np.concatenate([points[:-1], points[1:]], axis=1)
  5315. norm = mpl.colors.Normalize(vmin=0.0, vmax=1.0)
  5316. lc = LineCollection(segments, cmap=cmap, norm=norm, linewidth=line_width)
  5317. lc.set_array(vals_norm)
  5318. lc.set_capstyle("round")
  5319. ax.add_collection(lc)
  5321. # ---- colored points
  5322. sc = ax.scatter(
  5323. angles, vals_norm,
  5324. c=vals_norm, cmap=cmap, norm=norm,
  5325. s=point_size, zorder=3,
  5326. edgecolors="white", linewidths=0.9
  5327. )
  5329. # ---- panel label
  5330. ax.text(-0.10, 1.08, panel_label, transform=ax.transAxes,
  5331. fontsize=24, fontweight="bold", ha="left", va="top")
  5333. # ---- optional small title above plot (not in center)
  5334. if model_title is not None:
  5335. ax.set_title(model_title, fontsize=20, pad=24)
  5337. # ---- optional colorbar
  5338. if show_colorbar:
  5339. cbar = plt.colorbar(sc, ax=ax, pad=0.10, fraction=0.05)
  5340. cbar.set_label(cbar_label, fontsize=16)
  5341. cbar.ax.tick_params(labelsize=13)
  5343. # =========================================================
  5344. # STYLE
  5345. # =========================================================
  5346. plt.rcParams.update({
  5347. "font.family": "serif",
  5348. "font.serif": ["Times New Roman", "Times", "DejaVu Serif"],
  5349. "font.size": 16
  5350. })
  5352. RENAME = {"Location": "L"}
  5354. # =========================================================
  5355. # FIGURE
  5356. # =========================================================
  5357. fig, axes = plt.subplots(
  5358. 1, 3,
  5359. figsize=(20, 7),
  5360. dpi=300,
  5361. subplot_kw=dict(polar=True)
  5362. )
  5364. radar_perm_importance_colored(
  5365. axes[0],
  5366. bag_svm_feat_sorted, bag_svm_imp_sorted,
  5367. panel_label="a)",
  5368. model_title="Bagging SVM",
  5369. cmap="Blues",
  5370. rename_map=RENAME,
  5371. show_colorbar=False
  5372. )
  5374. radar_perm_importance_colored(
  5375. axes[1],
  5376. svm_feat_sorted, svm_imp_sorted,
  5377. panel_label="b)",
  5378. model_title="SVM",
  5379. cmap="Greens",
  5380. rename_map=RENAME,
  5381. show_colorbar=False
  5382. )
  5384. radar_perm_importance_colored(
  5385. axes[2],
  5386. gp_feat_sorted, gp_imp_sorted,
  5387. panel_label="c)",
  5388. model_title="Gaussian Process",
  5389. cmap="Purples",
  5390. rename_map=RENAME,
  5391. show_colorbar=True # show one colorbar only on the last panel
  5392. )
  5394. plt.tight_layout()
  5395. plt.savefig("perm_importance_radar_colored.png", dpi=300, bbox_inches="tight")
  5396. plt.savefig("perm_importance_radar_colored.pdf", bbox_inches="tight")
  5397. plt.show()
  5399. import numpy as np
  5400. import matplotlib.pyplot as plt
  5401. import os
  5402. import pandas as pd
  5405. xlsx_path = os.path.join("..", "content","Feature Importance Comparison.xlsx")
  5406. df_raw = pd.read_excel(xlsx_path)
  5408. print("Loaded shape:", df_raw.shape)
  5409. print("Columns:", df_raw.columns.tolist())
  5411. def get_feature_vectors_from_two_cols(df, feat_col, val_col, normalize=True):
  5412. """
  5413. Expects two columns:
  5414. feat_col, val_col
  5415. Returns:
  5416. feat_sorted (list[str]), imp_sorted (np.array[float])
  5417. """
  5418. sub = df[[feat_col, val_col]].copy()
  5419. sub.columns = ["Feature", "Importance"]
  5421. sub = sub.dropna(subset=["Feature"])
  5422. sub["Importance"] = pd.to_numeric(sub["Importance"], errors="coerce").fillna(0.0)
  5424. feats = sub["Feature"].astype(str).tolist()
  5425. imps = sub["Importance"].to_numpy(dtype=float)
  5427. # Clip negatives + normalize to sum=1 (same idea as your original code)
  5428. imps = np.clip(imps, 0, None)
  5429. if normalize:
  5430. s = imps.sum()
  5431. imps = imps / s if s > 0 else np.ones_like(imps) / len(imps)
  5433. return feats, imps
  5435. # === Load data from the new file ===
  5436. # Adjust these column names only if your file prints different names
  5437. morph_feat_sorted, morph_imp_sorted = get_feature_vectors_from_two_cols(
  5438. df_raw, "Unnamed: 0", "MLP", normalize=True
  5439. )
  5441. hemo_feat_sorted, hemo_imp_sorted = get_feature_vectors_from_two_cols(
  5442. df_raw, "Unnamed: 3", "RF", normalize=True
  5443. )
  5445. comb_feat_sorted, comb_imp_sorted = get_feature_vectors_from_two_cols(
  5446. df_raw, "Unnamed: 6", "Bagging SVM", normalize=True
  5447. )
  5449. print("Top Morphology features:", morph_feat_sorted[:5])
  5450. print("Top Hemodynamics features:", hemo_feat_sorted[:5])
  5451. print("Top Combined features:", comb_feat_sorted[:5])
  5454. def donut_perm_importance_topk_onring(
  5455. ax,
  5456. feat_sorted, imp_sorted,
  5457. model_name="Model",
  5458. panel_label="a)",
  5459. cmap="Blues",
  5460. top_k_labels=5,
  5461. ring_width=0.62, # ring thickness
  5462. startangle=90,
  5463. label_r=None, # radius where labels sit (mid-ring by default)
  5464. fontsize=23,
  5465. fontweight="normal",
  5466. text_color="white",
  5467. center_fontsize=28,
  5468. panel_fontsize=32,
  5470. # fading controls (more obvious)
  5471. fade_min=0.08, # lightest (end)
  5472. fade_max=1.99, # darkest (start)
  5473. fade_gamma=1.2, # higher => stronger contrast
  5475. # rotation controls (labels follow slice angle; "vertical" on slice)
  5476. rotation_offset=0, # add 90° to make it vertical relative to slice direction
  5477. flip_upright=True, # keep text readable on left half
  5479. # optional rename (e.g., Location -> L)
  5480. rename_map=None,
  5481. ):
  5482. ax.set_aspect("equal")
  5483. ax.axis("off")
  5485. # --- clip negatives + normalize ---
  5486. vals = np.asarray(imp_sorted, float)
  5487. vals = np.clip(vals, 0, None)
  5488. if vals.sum() == 0:
  5489. vals = np.ones_like(vals)
  5490. vals = vals / vals.sum()
  5492. feat_sorted = list(feat_sorted)
  5494. # --- sort by importance desc (so first k = top k) ---
  5495. order = np.argsort(vals)[::-1]
  5496. vals = vals[order]
  5497. feat_sorted = [feat_sorted[i] for i in order]
  5499. # --- fading by slice mid-angle fraction (0..1): later angle => lighter ---
  5500. cum = np.cumsum(vals)
  5501. prev = np.concatenate(([0.0], cum[:-1]))
  5502. mid_frac = 0.5 * (prev + cum) # 0..1 around circle (clockwise because counterclock=False)
  5504. fade_strength = (1.0 - mid_frac) ** fade_gamma
  5505. fade_strength = fade_min + (fade_max - fade_min) * fade_strength
  5507. cmap_obj = plt.get_cmap(cmap)
  5508. colors = [cmap_obj(f) for f in fade_strength]
  5510. wedges, _ = ax.pie(
  5511. vals,
  5512. labels=None,
  5513. startangle=startangle,
  5514. counterclock=False,
  5515. colors=colors,
  5516. wedgeprops=dict(width=ring_width, edgecolor="white", linewidth=1.5),
  5517. )
  5519. # panel label
  5520. ax.text(-0.03, 1.0, panel_label, transform=ax.transAxes,
  5521. fontsize=panel_fontsize, fontweight="bold", va="top", ha="left")
  5523. # center text
  5524. ax.text(0, 0, model_name, ha="center", va="center",
  5525. fontsize=center_fontsize, fontweight="normal")
  5527. # label radius in middle of ring by default
  5528. if label_r is None:
  5529. inner_r = 1.0 - ring_width
  5530. label_r = 0.5 * (1.0 + inner_r)
  5532. # label only top-k, ON ring, rotated along slice direction (radial) + 90°
  5533. k = min(top_k_labels, len(wedges))
  5534. rename_map = rename_map or {}
  5536. for i in range(k):
  5537. w = wedges[i]
  5538. ang = 0.5 * (w.theta1 + w.theta2) # degrees
  5539. ang_rad = np.deg2rad(ang)
  5541. xt = label_r * np.cos(ang_rad)
  5542. yt = label_r * np.sin(ang_rad)
  5544. # label text
  5545. lbl = str(feat_sorted[i])
  5546. lbl = rename_map.get(lbl, lbl)
  5548. rot = ang + rotation_offset
  5550. ax.text(
  5551. xt, yt, lbl,
  5552. ha="center", va="center",
  5553. rotation=rot, rotation_mode="anchor",
  5554. fontsize=fontsize, fontweight=fontweight,
  5555. color=text_color
  5556. )
  5558. ax.set_xlim(-1.05, 1.05)
  5559. ax.set_ylim(-1.05, 1.05)
  5562. # -----------------------------
  5563. # Plot
  5564. # -----------------------------
  5565. plt.rcParams.update({
  5566. "font.family": "serif",
  5567. "font.serif": ["Times New Roman", "Times", "DejaVu Serif"],
  5568. "font.size": 18
  5569. })
  5571. RENAME = {"Location": "L"}
  5573. fig, axes = plt.subplots(1, 3, figsize=(18, 6), dpi=300)
  5575. donut_perm_importance_topk_onring(
  5576. axes[0], morph_feat_sorted, morph_imp_sorted,
  5577. model_name="Morphology", panel_label="a)",
  5578. cmap="Oranges",
  5579. rename_map=RENAME
  5580. )
  5582. donut_perm_importance_topk_onring(
  5583. axes[1], hemo_feat_sorted, hemo_imp_sorted,
  5584. model_name="Hemodynamics", panel_label="b)",
  5585. cmap="Blues",
  5586. rename_map=RENAME
  5587. )
  5589. donut_perm_importance_topk_onring(
  5590. axes[2], comb_feat_sorted, comb_imp_sorted,
  5591. model_name="Morphology &\nHemodynamics", panel_label="c)",
  5592. cmap="Greys",
  5593. rename_map=RENAME
  5594. )
  5596. plt.tight_layout()
  5597. plt.savefig("feature_importance_comparison_donuts.png", dpi=300, bbox_inches="tight")
  5598. plt.savefig("feature_importance_comparison_donuts.pdf", bbox_inches="tight")
  5599. plt.show()
  5601. #Use results from different papers
  5603. fig, axes = plt.subplots(1, 3, figsize=(13, 5), dpi=300, sharex=False)
  5605. panel_labels = ["a)", "b)", "c)"]
  5607. for ax, data, panel in zip(
  5608. axes,
  5609. [morph_top, hemo_top, comb_top],
  5610. panel_labels
  5611. ):
  5612. y = np.arange(len(data))
  5613. ax.hlines(y, 0, data["Importance"], color=data["Color"], linewidth=2.5)
  5614. ax.scatter(data["Importance"], y, color=data["Color"], s=80)
  5616. ax.set_yticks(y)
  5617. ax.set_yticks(y)
  5618. ax.set_yticklabels(data["Feature"])
  5619. ax.tick_params(axis="y", labelsize=18) # bigger y labels
  5621. # === X axis settings ===
  5622. max_val = data["Importance"].max()
  5623. ax.set_xlim(0, max_val * 1.05) # ensure 0 aligns with vertical axis
  5625. # only 3 ticks including zero
  5626. xticks = np.linspace(0, max_val, 3)
  5627. ax.set_xticks(xticks)
  5628. ax.set_xticklabels([f"{x:.2f}" for x in xticks], fontsize=16)
  5630. ax.grid(axis="x", alpha=0.25)
  5632. # remove top and right spines
  5633. ax.spines["top"].set_visible(False)
  5634. ax.spines["right"].set_visible(False)
  5636. # make left spine exactly the zero line
  5637. ax.spines["left"].set_position(("data", 0))
  5639. # panel label
  5640. ax.text(
  5641. -0.6, 1.05, panel,
  5642. transform=ax.transAxes,
  5643. ha="left", va="bottom",
  5644. fontsize=23, fontweight="bold"
  5645. )
  5647. plt.subplots_adjust(wspace=0.75)
  5648. plt.savefig("feature_importance_comparison.png", dpi=300, bbox_inches="tight")
  5649. plt.savefig("feature_importance_comparison.pdf", bbox_inches="tight")
  5650. plt.show()
  5652. # ============================================================
  5653. # BAGGING SVM SCENARIO COMPARISON
  5654. # Morphology-only vs Hemodynamic-only vs Combined
  5655. #
  5656. # Scenario definitions after dropping "Vessel ID":
  5657. # Clinical columns: first 3 columns
  5658. # Morphology columns: next 35 columns
  5659. # Hemodynamic columns: next 17 columns
  5660. # Label column: last column, "Status"
  5661. #
  5662. # Outputs:
  5663. # 1) 5-fold CV performance on TRAIN90
  5664. # 2) Held-out TEST10 performance
  5665. # 3) Permutation importance table for all scenarios
  5666. # 4) Three-panel feature-importance figure
  5667. # ============================================================
  5669. import os
  5670. import numpy as np
  5671. import pandas as pd
  5672. import matplotlib.pyplot as plt
  5674. from sklearn.pipeline import Pipeline
  5675. from sklearn.preprocessing import StandardScaler
  5676. from sklearn.model_selection import train_test_split, StratifiedKFold
  5677. from sklearn.metrics import (
  5678. roc_auc_score, accuracy_score, precision_score,
  5679. recall_score, f1_score, confusion_matrix
  5680. )
  5681. from sklearn.svm import SVC
  5682. from sklearn.ensemble import BaggingClassifier
  5683. from sklearn.base import clone
  5684. from sklearn.inspection import permutation_importance
  5685. from sklearn.utils.class_weight import compute_class_weight
  5687. # ----------------------------
  5688. # CONFIG
  5689. # ----------------------------
  5690. xlsx_path = os.path.join("..", "content", "ClassificationV2.xlsx")
  5691. # If you are running locally and the file is in the same folder, use:
  5692. # xlsx_path = "ClassificationV2.xlsx"
  5694. RANDOM_STATE_SPLIT = 8
  5695. CV_FOLDS = 5
  5696. THR = 0.5
  5697. N_REPEATS_PI = 10
  5698. TOP_N = 8 # change this if you want more or fewer features in each panel
  5700. # ----------------------------
  5701. # 1) LOAD AND PREPROCESS DATA
  5702. # ----------------------------
  5703. df = pd.read_excel(xlsx_path)
  5705. mapping_Sex = {"Male": 0, "Female": 1}
  5706. mapping_Irregularity = {"R": 0, "I": 1}
  5707. mapping_Kind = {"L": 0, "B": 1}
  5708. mapping_Location = {
  5709. "MCA": 0, "ICA": 1, "ACA": 2, "PCA": 3,
  5710. "VA": 4, "BA": 4, "SCA": 5
  5711. }
  5712. mapping_Stability = {"stable": 0, "unstable": 1}
  5713. mapping_Complexity = {"simple": 0, "complex": 1}
  5714. mapping_Status = {"unruptured": 0, "ruptured": 1}
  5716. df = df.replace({
  5717. "Sex": mapping_Sex,
  5718. "I": mapping_Irregularity,
  5719. "LB": mapping_Kind,
  5720. "Location": mapping_Location,
  5721. "FS": mapping_Stability,
  5722. "FC": mapping_Complexity,
  5723. "Status": mapping_Status
  5724. })
  5726. if "Vessel ID" in df.columns:
  5727. df = df.drop(labels="Vessel ID", axis=1)
  5729. if "Status" not in df.columns:
  5730. raise ValueError("Column 'Status' was not found. Please check the Excel file.")
  5732. # Make sure all remaining values are numeric
  5733. df = df.apply(pd.to_numeric)
  5735. # ----------------------------
  5736. # 2) DEFINE FEATURE GROUPS
  5737. # ----------------------------
  5738. all_columns = list(df.columns)
  5740. clinical_cols = all_columns[0:3]
  5741. morphology_cols = all_columns[3:38] # next 35 columns
  5742. hemodynamic_cols = all_columns[38:55] # next 17 columns
  5743. label_col = all_columns[-1]
  5745. print("Clinical columns:", clinical_cols)
  5746. print("Number of morphology columns:", len(morphology_cols))
  5747. print("Number of hemodynamic columns:", len(hemodynamic_cols))
  5748. print("Label column:", label_col)
  5750. if label_col != "Status":
  5751. raise ValueError(
  5752. f"The last column is expected to be 'Status', but found '{label_col}'. "
  5753. "Please check the dataset column order."
  5754. )
  5756. scenarios = {
  5757. "Morphology-only": clinical_cols + morphology_cols,
  5758. "Hemodynamic-only": clinical_cols + hemodynamic_cols,
  5759. "Combined": clinical_cols + morphology_cols + hemodynamic_cols
  5760. }
  5762. y = df["Status"].to_numpy(dtype=int)
  5764. # ----------------------------
  5765. # 3) HELPER FUNCTIONS
  5766. # ----------------------------
  5767. def specificity(y_true, y_pred):
  5768. tn, fp, fn, tp = confusion_matrix(y_true, y_pred, labels=[0, 1]).ravel()
  5769. return tn / (tn + fp) if (tn + fp) else 0.0
  5771. def metrics_from_prob(y_true, prob, thr=0.5):
  5772. pred = (prob >= thr).astype(int)
  5773. return {
  5774. "AUC": roc_auc_score(y_true, prob),
  5775. "ACC": accuracy_score(y_true, pred),
  5776. "Precision": precision_score(y_true, pred, zero_division=0),
  5777. "Recall": recall_score(y_true, pred, zero_division=0),
  5778. "F1": f1_score(y_true, pred, zero_division=0),
  5779. "Specificity": specificity(y_true, pred),
  5780. }
  5782. def get_probabilities(estimator, X_eval):
  5783. if hasattr(estimator, "predict_proba"):
  5784. return estimator.predict_proba(X_eval)[:, 1]
  5785. if hasattr(estimator, "decision_function"):
  5786. s = estimator.decision_function(X_eval)
  5787. return (s - s.min()) / (s.max() - s.min() + 1e-12)
  5788. return estimator.predict(X_eval).astype(float)
  5790. def make_sample_weight_balanced(y_local):
  5791. y_local = np.asarray(y_local).astype(int).ravel()
  5792. n = len(y_local)
  5793. n0 = np.sum(y_local == 0)
  5794. n1 = np.sum(y_local == 1)
  5795. w0 = n / (2 * n0) if n0 else 1.0
  5796. w1 = n / (2 * n1) if n1 else 1.0
  5797. return np.where(y_local == 0, w0, w1).astype(float)
  5799. def make_bagging_svm(n_features):
  5800. """
  5801. Same Bagging SVM setup as your original code.
  5802. max_features is adjusted only when the scenario has fewer than 54 features.
  5803. """
  5804. return Pipeline([
  5805. ("scaler", StandardScaler()),
  5806. ("clf", BaggingClassifier(
  5807. estimator=SVC(
  5808. C=80,
  5809. gamma=0.002,
  5810. probability=True,
  5811. kernel="rbf"
  5812. ),
  5813. n_estimators=100,
  5814. bootstrap=True,
  5815. max_samples=0.9,
  5816. max_features=min(54, n_features),
  5817. random_state=42
  5818. ))
  5819. ])
  5821. # ----------------------------
  5822. # 4) TRAIN/EVALUATE BAGGING SVM FOR EACH SCENARIO
  5823. # ----------------------------
  5824. cv_rows = []
  5825. test_rows = []
  5826. fitted_models = {}
  5827. split_data = {}
  5829. for scenario_name, feature_cols in scenarios.items():
  5831. print("\n" + "=" * 70)
  5832. print(f"Scenario: {scenario_name}")
  5833. print(f"Number of features: {len(feature_cols)}")
  5834. print("=" * 70)
  5836. X = df[feature_cols].to_numpy(dtype=np.float64)
  5838. X_train90, X_test10, y_train90, y_test10 = train_test_split(
  5839. X,
  5840. y,
  5841. test_size=0.10,
  5842. random_state=RANDOM_STATE_SPLIT,
  5843. stratify=y
  5844. )
  5846. print("Train90 shape:", X_train90.shape, "Test10 shape:", X_test10.shape)
  5847. print("Train class counts:", np.bincount(y_train90))
  5848. print("Test class counts:", np.bincount(y_test10))
  5850. split_data[scenario_name] = {
  5851. "feature_cols": feature_cols,
  5852. "X_train90": X_train90,
  5853. "X_test10": X_test10,
  5854. "y_train90": y_train90,
  5855. "y_test10": y_test10
  5856. }
  5858. # ----------------------------
  5859. # 5-fold CV on TRAIN90 only
  5860. # ----------------------------
  5861. skf = StratifiedKFold(
  5862. n_splits=CV_FOLDS,
  5863. shuffle=True,
  5864. random_state=42
  5865. )
  5867. fold_metrics = []
  5869. for fold, (tr_idx, va_idx) in enumerate(skf.split(X_train90, y_train90), start=1):
  5871. Xtr, Xva = X_train90[tr_idx], X_train90[va_idx]
  5872. ytr, yva = y_train90[tr_idx], y_train90[va_idx]
  5874. est = make_bagging_svm(n_features=Xtr.shape[1])
  5875. sw_tr = make_sample_weight_balanced(ytr)
  5877. est.fit(Xtr, ytr, clf__sample_weight=sw_tr)
  5879. prob = get_probabilities(est, Xva)
  5880. m = metrics_from_prob(yva, prob, thr=THR)
  5881. fold_metrics.append(m)
  5883. df_cv = pd.DataFrame(fold_metrics)
  5885. cv_rows.append({
  5886. "Scenario": scenario_name,
  5887. "Features": len(feature_cols),
  5888. "AUC": f"{df_cv['AUC'].mean():.3f} ± {df_cv['AUC'].std(ddof=0):.3f}",
  5889. "ACC": f"{df_cv['ACC'].mean():.3f} ± {df_cv['ACC'].std(ddof=0):.3f}",
  5890. "Precision": f"{df_cv['Precision'].mean():.3f} ± {df_cv['Precision'].std(ddof=0):.3f}",
  5891. "Recall": f"{df_cv['Recall'].mean():.3f} ± {df_cv['Recall'].std(ddof=0):.3f}",
  5892. "F1": f"{df_cv['F1'].mean():.3f} ± {df_cv['F1'].std(ddof=0):.3f}",
  5893. "Specificity": f"{df_cv['Specificity'].mean():.3f} ± {df_cv['Specificity'].std(ddof=0):.3f}",
  5894. })
  5896. # ----------------------------
  5897. # Final fit on TRAIN90, evaluate once on TEST10
  5898. # ----------------------------
  5899. final_est = make_bagging_svm(n_features=X_train90.shape[1])
  5900. sw_train90 = make_sample_weight_balanced(y_train90)
  5902. final_est.fit(X_train90, y_train90, clf__sample_weight=sw_train90)
  5904. prob_test = get_probabilities(final_est, X_test10)
  5905. test_metrics = metrics_from_prob(y_test10, prob_test, thr=THR)
  5907. test_metrics["Scenario"] = scenario_name
  5908. test_metrics["Features"] = len(feature_cols)
  5909. test_rows.append(test_metrics)
  5911. fitted_models[scenario_name] = final_est
  5913. # ----------------------------
  5914. # 5) PERFORMANCE TABLES
  5915. # ----------------------------
  5916. cv_table = pd.DataFrame(cv_rows)
  5917. test_table = pd.DataFrame(test_rows)
  5919. test_table = test_table[
  5920. ["Scenario", "Features", "AUC", "ACC", "Precision", "Recall", "F1", "Specificity"]
  5921. ].sort_values("AUC", ascending=False).reset_index(drop=True)
  5923. print("\n=== 5-fold CV on TRAIN90: Bagging SVM only ===")
  5924. display(cv_table)
  5926. print("\n=== Held-out TEST10: Bagging SVM only ===")
  5927. display(test_table.round(3))
  5929. # Save performance tables
  5930. with pd.ExcelWriter("bagging_svm_scenario_performance.xlsx", engine="openpyxl") as writer:
  5931. cv_table.to_excel(writer, sheet_name="CV_TRAIN90", index=False)
  5932. test_table.round(4).to_excel(writer, sheet_name="TEST10", index=False)
  5934. print("\nSaved: bagging_svm_scenario_performance.xlsx")
  5936. # ----------------------------
  5937. # 6) CV PERMUTATION IMPORTANCE FOR EACH SCENARIO
  5938. # ----------------------------
  5939. def cv_perm_importance_auc_bagging_svm(X, y, feature_cols, n_repeats=10, seed=42):
  5940. """
  5941. Computes fold-level permutation importance on held-out CV folds.
  5942. Aggregates mean and SD across folds.
  5943. Uses Bagging SVM with balanced sample weights.
  5944. """
  5945. skf = StratifiedKFold(
  5946. n_splits=CV_FOLDS,
  5947. shuffle=True,
  5948. random_state=42
  5949. )
  5951. fold_means = []
  5953. rng = np.random.RandomState(seed)
  5955. for fold, (tr_idx, te_idx) in enumerate(skf.split(X, y), start=1):
  5957. Xtr, Xte = X[tr_idx], X[te_idx]
  5958. ytr, yte = y[tr_idx], y[te_idx]
  5960. est = make_bagging_svm(n_features=Xtr.shape[1])
  5961. sw_tr = make_sample_weight_balanced(ytr)
  5963. est.fit(Xtr, ytr, clf__sample_weight=sw_tr)
  5965. pi = permutation_importance(
  5966. est,
  5967. Xte,
  5968. yte,
  5969. scoring="roc_auc",
  5970. n_repeats=n_repeats,
  5971. random_state=rng.randint(0, 10**9),
  5972. n_jobs=-1
  5973. )
  5975. fold_means.append(pi.importances_mean)
  5977. fold_means = np.vstack(fold_means)
  5979. mean_imp = fold_means.mean(axis=0)
  5980. std_imp = fold_means.std(axis=0)
  5982. out = pd.DataFrame({
  5983. "Feature": feature_cols,
  5984. "Importance_Mean": mean_imp,
  5985. "Importance_SD": std_imp
  5986. })
  5988. out["Rank"] = out["Importance_Mean"].rank(
  5989. method="min",
  5990. ascending=False
  5991. ).astype(int)
  5993. out = out.sort_values("Importance_Mean", ascending=False).reset_index(drop=True)
  5995. return out
  5997. perm_tables = {}
  5999. for scenario_name, feature_cols in scenarios.items():
  6000. X = df[feature_cols].to_numpy(dtype=np.float64)
  6002. perm_tables[scenario_name] = cv_perm_importance_auc_bagging_svm(
  6003. X=X,
  6004. y=y,
  6005. feature_cols=feature_cols,
  6006. n_repeats=N_REPEATS_PI,
  6007. seed=42
  6008. )
  6010. # Save permutation importance tables
  6011. with pd.ExcelWriter("bagging_svm_scenario_permutation_importance.xlsx", engine="openpyxl") as writer:
  6012. for scenario_name, tab in perm_tables.items():
  6013. sheet_name = scenario_name.replace("-", "_")[:31]
  6014. tab.to_excel(writer, sheet_name=sheet_name, index=False)
  6016. print("Saved: bagging_svm_scenario_permutation_importance.xlsx")
  6018. # ----------------------------
  6019. # 7) PREPARE TOP FEATURES FOR FIGURE
  6020. # ----------------------------
  6021. def feature_group_color(feature):
  6022. if feature in clinical_cols:
  6023. return "#4D4D4D" # clinical
  6024. elif feature in morphology_cols:
  6025. return "#1F77B4" # morphology
  6026. elif feature in hemodynamic_cols:
  6027. return "#D62728" # hemodynamic
  6028. else:
  6029. return "#7F7F7F"
  6031. plot_data = {}
  6033. for scenario_name, tab in perm_tables.items():
  6034. top = tab.head(TOP_N).copy()
  6036. # Sort ascending so the largest value appears at the top of the horizontal plot
  6037. top = top.sort_values("Importance_Mean", ascending=True).reset_index(drop=True)
  6039. top["Color"] = top["Feature"].apply(feature_group_color)
  6040. plot_data[scenario_name] = top
  6042. # ----------------------------
  6043. # 8) THREE-PANEL FEATURE IMPORTANCE FIGURE
  6044. # ----------------------------
  6045. fig, axes = plt.subplots(
  6046. 1, 3,
  6047. figsize=(13, 5),
  6048. dpi=300,
  6049. sharex=False
  6050. )
  6052. panel_labels = ["a)", "b)", "c)"]
  6053. scenario_order = ["Morphology-only", "Hemodynamic-only", "Combined"]
  6055. for ax, scenario_name, panel in zip(axes, scenario_order, panel_labels):
  6057. data = plot_data[scenario_name]
  6058. y_pos = np.arange(len(data))
  6060. ax.hlines(
  6061. y_pos,
  6062. 0,
  6063. data["Importance_Mean"],
  6064. color=data["Color"],
  6065. linewidth=2.5
  6066. )
  6068. ax.scatter(
  6069. data["Importance_Mean"],
  6070. y_pos,
  6071. color=data["Color"],
  6072. s=80
  6073. )
  6075. ax.set_yticks(y_pos)
  6076. ax.set_yticklabels(data["Feature"])
  6077. ax.tick_params(axis="y", labelsize=18)
  6079. max_val = data["Importance_Mean"].max()
  6081. # Protect against zero or negative max importance values
  6082. if max_val <= 0:

aneurysm_classification.py at commit dc01e05, no license · at the source

Overview

Authors: Mostafa Zakeri1, Amirhossein Atef1, Mohammad Aziznia1, Azadeh Jafari1
ORCID iDs: Azadeh Jafari
  1. School of Mechanical Engineering, College of Engineering, University of Tehran, Tehran, Iran
Institutions: University of Tehran (Iran)
Journal: PloS one, volume 21, issue 10, article e0360923
Dates: received 30 July 2026; accepted 2 October 2026; published online 30 October 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1371/journal.pone.0360923 · PMID 42962059 · PMCID PMC13722296 · OpenAlex W7224775446
Open access: gold, a free copy (OpenAlex)
Status: code verified (how the code is found and verified)
Categories: human (organism), stroke (population), clinical / translational (subfield)
Methods: Statistics, Smoothing, state filtering, decompositions, Machine learning, Spectral & time-frequency
MeSH: Aneurysm, Ruptured*, Hemodynamics*, Intracranial Aneurysm*, Machine Learning*, Classification Algorithms, Clustering Algorithms, Humans, Predictive Learning Models, ROC Curve (* major topic)
Topic: Intracranial Aneurysms: Treatment and Complications (Clinical Neurology, Medicine), according to OpenAlex
Citations: not cited yet (Europe PMC); 51 references in the paper

Abstract

Accurate assessment of cerebral aneurysm rupture status remains a critical challenge in neurovascular research due to the complex interplay between geometric, hemodynamic, and clinical factors. In this study, we developed and evaluated a machine-learning framework for retrospective rupture-status classification in 708 patient-specific aneurysms by integrating 55 predictors: three clinical variables, 35 morphological descriptors, and 17 hemodynamic parameters derived from computational fluid dynamics (CFD) simulations. Blood flow was modeled as incompressible laminar flow with a non-Newtonian Carreau–Yasuda viscosity model to capture shear-thinning behavior in cerebral arteries. We compared 13 supervised classifiers, 10 unsupervised clustering approaches, and two neural-network models. Model training and hyperparameter optimization were performed using stratified five-fold cross-validation within a 90/10 training–testing split, with performance evaluated using ROC–AUC, precision, recall, specificity, and F1-score. Among the evaluated models, BaggingSVM achieved the strongest discrimination on the independent test set (ROC–AUC = 0.78), followed by SVM (0.75) and Gaussian Process (0.75). A scenario comparison using BaggingSVM showed that the combined feature set improved classification performance (AUC = 0.78) compared with clinical+morphological features (AUC = 0.76) and clinical+hemodynamical features (AUC = 0.71) settings. Sensitivity-oriented threshold tuning increased recall for ruptured aneurysm detection, reaching values up to 0.96 for BaggingSVM. Permutation feature importance analysis revealed consistent dominance of ellipticity index, low shear area ratio (LSAR), aneurysm location, irregularity and aneurysm angle across models. a small subset of variables contributes disproportionately to rupture-status classification, although the reduced-feature analysis indicates that lower-ranked predictors collectively provide additional discriminatory information. The proposed framework demonstrates the value of integrating physiologically informed hemodynamic modeling with machine learning to improve data-driven rupture status assessment and provides a scalable foundation for future clinical decision-support studies in cerebrovascular diagnostics.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 19 matches between paragraphs and lines of code.

mstf229/Cerebral-Aneurysm-Rupture-Prediction

License: none: the authors keep all their rights
State: the link answers, verified on 1 November 2026
Evidence: files inventoried
Commit: dc01e05b467217fd1ea82a4aab55e7fee23f32e8, 5 October 2026
Languages: Jupyter (1), Python (1)
Size: 4 files, 2 scripts
Software Heritage: not archived
Found in: “Data Availability”
Holds: README, 1 notebook
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Keras (2 files), Matplotlib (2 files), NumPy (2 files), pandas (2 files), scikit-learn (2 files), SciPy (2 files), TensorFlow (2 files), XGBoost (2 files), UMAP (1 file)
Availability: 1 check, the latest on 1 November 2026: the link answers
  • 1 November 2026: the link answers
3 files
At the source: github.com/mstf229/Cerebral-Aneurysm-Rupture-Prediction

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 2 scripts, each with its path and the digest of its content;
  • 19 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

Data Availability

The original aneurysm geometries and associated clinical and morphological data were obtained from the AneuX morphology database, an open-access multicentric database combining data from the AneuX, @neurIST, and Aneurisk projects. The source data are available from Zenodo at https://doi.org/10.5281/zenodo.6678442. The author-generated derived variables and machine-learning analysis code underlying the findings of this study are publicly available in the Cerebral Aneurysm Rupture Prediction repository: https://github.com/mstf229/Cerebral-Aneurysm-Rupture-Prediction. The repository provides case identifiers allowing the derived variables to be linked to the corresponding publicly available AneuX source records.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 1 November 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 4 authors, 9 MeSH terms, 46 references.

Cite

This paper

Zakeri, M., Atef, A., Aziznia, M., & Jafari, A. (2026). An integrated machine learning framework for cerebral aneurysm rupture status classification using clinical, morphological, and hemodynamic features. PloS one, 21(10), e0360923. https://doi.org/10.1371/journal.pone.0360923

BibTeX

@article{zakeri2026integrated,
author = {Zakeri, Mostafa and Atef, Amirhossein and Aziznia, Mohammad and Jafari, Azadeh},
title = {{An integrated machine learning framework for cerebral aneurysm rupture status classification using clinical, morphological, and hemodynamic features}},
journal = {PloS one},
year = {2026},
month = oct,
volume = {21},
number = {10},
pages = {e0360923},
publisher = {PLOS},
issn = {1932-6203},
doi = {10.1371/journal.pone.0360923},
url = {https://doi.org/10.1371/journal.pone.0360923},
pmid = {42962059},
pmcid = {PMC13722296}
}

RIS

TY - JOUR
AU - Zakeri, Mostafa
AU - Atef, Amirhossein
AU - Aziznia, Mohammad
AU - Jafari, Azadeh
TI - An integrated machine learning framework for cerebral aneurysm rupture status classification using clinical, morphological, and hemodynamic features
T2 - PloS one
J2 - PLoS One
PY - 2026
DA - 2026/10/30
VL - 21
IS - 10
SP - e0360923
SN - 1932-6203
PB - PLOS
DO - 10.1371/journal.pone.0360923
UR - https://doi.org/10.1371/journal.pone.0360923
LA - en
ER -

CSL-JSON

{
"id": "10.1371/journal.pone.0360923",
"type": "article-journal",
"title": "An integrated machine learning framework for cerebral aneurysm rupture status classification using clinical, morphological, and hemodynamic features",
"container-title": "PloS one",
"author": [
{
"family": "Zakeri",
"given": "Mostafa"
},
{
"family": "Atef",
"given": "Amirhossein"
},
{
"family": "Aziznia",
"given": "Mohammad"
},
{
"family": "Jafari",
"given": "Azadeh"
}
],
"container-title-short": "PLoS One",
"volume": "21",
"issue": "10",
"page": "e0360923",
"DOI": "10.1371/journal.pone.0360923",
"PMID": "42962059",
"PMCID": "PMC13722296",
"ISSN": "1932-6203",
"publisher": "PLOS",
"URL": "https://doi.org/10.1371/journal.pone.0360923",
"language": "en",
"issued": {
"date-parts": [
[
2026,
10,
30
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1016/j.patter.2023.100709 [code]
A deep-learning method for the end-to-end prediction of intracranial aneurysm rupture risk.
Journal: Patterns (New York, N.Y.)
In common: Keras, TensorFlow, scikit-learn, 4 other tools, stroke, 4 references
[2] doi:10.7717/peerj.19393 [code]
Deep learning-based fine-grained assessment of aneurysm wall characteristics using 4D-CT angiography.
Journal: PeerJ
In common: Keras, TensorFlow, scikit-learn, 4 other tools, stroke, clinical / translational, 1 reference
[3] doi:10.1007/s12021-022-09597-0 [code]
Towards Automated Brain Aneurysm Detection in TOF-MRA: Open Data, Weak Labels, and Anatomical Knowledge.
Journal: Neuroinformatics
In common: TensorFlow, scikit-learn, pandas, 3 other tools, stroke, 2 references
[4] doi:10.21037/qims-24-1269 [code]
Deep learning for the detection of moyamoya angiopathy using T2-weighted images: a multicenter study.
Journal: Quantitative imaging in medicine and surgery
In common: Keras, TensorFlow, scikit-learn, 3 other tools, stroke, clinical / translational, 1 reference
[5] doi:10.3390/s25247415 [code]
Decision Support Systems in Neurosurgery: Current Applications and Future Directions.
Journal: Sensors (Basel, Switzerland)
In common: UMAP, scikit-learn, pandas, 3 other tools, 1 reference
[6] doi:10.1038/s41467-024-53550-5 [code]
Self-supervised learning for accurately modelling hierarchical evolutionary patterns of cerebrovasculature.
Journal: Nature communications
In common: scikit-learn, pandas, SciPy, 2 other tools, stroke, clinical / translational, 1 reference
[7] doi:10.1088/2057-1976/acb1b3
Can we explain machine learning-based prediction for rupture status assessments of intracranial aneurysms?
Journal: Biomedical physics & engineering express
In common: stroke, 3 references
[8] doi:10.1016/j.mcpdig.2023.08.009 [code]
Hemorrhage Evaluation and Detector System for Underserved Populations: HEADS-UP.
Journal: Mayo Clinic proceedings. Digital health
In common: pandas, SciPy, Matplotlib, 1 other tool, 1 reference
[9] doi:10.3390/brainsci15121295
A Mask R-CNN-Based Approach for Brain Aneurysm Detection and Segmentation from TOF-MRA Data.
Journal: Brain sciences
In common: stroke, 2 references
[10] doi:10.1109/tbme.2025.3572837 [code]
A Joint Geometric Topological Analysis Network (JGTA-Net) for Detecting and Segmenting Intracranial Aneurysms.
Journal: IEEE transactions on bio-medical engineering
In common: SciPy, NumPy, stroke, 1 reference

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Add the badge to the code's README

OSCR: the badge of this paper

If you are an author or a maintainer of this paper's code, you can add this badge. It says what OSCR knows of the paper (its code, its validated tracing map) and links the code to this page. Paste one of these snippets under the README's title. Keep the two marks around it, so that the badge can be replaced in place and is never added twice. OSCR changes nothing in your repository.

Markdown
<!-- oscr-badge:start -->
[![OSCR](https://openscicode.org/badge/doi/10.1371/journal.pone.0360923.svg)](https://openscicode.org/badge/doi/10.1371/journal.pone.0360923)
<!-- oscr-badge:end -->
reStructuredText
.. oscr-badge:start

.. image:: https://openscicode.org/badge/doi/10.1371/journal.pone.0360923.svg
   :target: https://openscicode.org/badge/doi/10.1371/journal.pone.0360923
   :alt: OSCR

.. oscr-badge:end
HTML
<!-- oscr-badge:start -->
<a href="https://openscicode.org/badge/doi/10.1371/journal.pone.0360923"><img src="https://openscicode.org/badge/doi/10.1371/journal.pone.0360923.svg" alt="OSCR"></a>
<!-- oscr-badge:end -->

On GitHub, open the README, choose Edit (the pencil), paste the snippet under the title, then "Commit changes".

If you choose to, you can have the badge kept up to date for you, with oscr badge --write from the command line or with the GitHub Action in your repository's workflows. Every other form and option (a card, AsciiDoc, LaTeX, shields.io) is on the badge page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: code errors and mismatches between the code and the paper about this paper and its code, from readers and its authors, beside the paper and its tracing map.

Report a mismatch between the code and the paper · Report a code error (signed in)

Reproductions: reports from readers who ran the authors' code and did not get the paper's result, what they ran, with which environment, commit and data.

Report a reproduction failure (signed in)

Versions of its code: the releases of its code tied to a version of this paper (the preprint, the accepted manuscript, the version of record…), each with the tracing map versioned with it.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.