{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.12.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"nvidiaTeslaT4","dataSources":[{"sourceType":"competition","sourceId":25563,"databundleVersionId":2094376}],"dockerImageVersionId":31329,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"id":"83113e8b-9016-4d4f-bbb7-e387f04d07d6","cell_type":"markdown","source":"# Apple Disease Prediction with QPCA\n## Section 1: Setup, Configuration, and Caching","metadata":{}},{"id":"b12dfda0-d462-4de3-8d3a-582b2a5718a2","cell_type":"code","source":"!pip install iterative-stratification timm albumentations pennylane xgboost","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-05-08T10:07:51.593971Z","iopub.execute_input":"2026-05-08T10:07:51.594741Z","iopub.status.idle":"2026-05-08T10:08:03.929181Z","shell.execute_reply.started":"2026-05-08T10:07:51.594710Z","shell.execute_reply":"2026-05-08T10:08:03.928375Z"}},"outputs":[],"execution_count":null},{"id":"9accecf9-38af-45b4-b7b1-5807ee63a03b","cell_type":"code","source":"import os, sys, time, warnings, random, logging\nfrom pathlib import Path\nwarnings.filterwarnings('ignore')\n\n# ── reproducibility ─────────────────────────────────────────\nSEED = 42\nrandom.seed(SEED)\nimport numpy as np; np.random.seed(SEED)\nimport torch; torch.manual_seed(SEED)\nif torch.cuda.is_available(): torch.cuda.manual_seed_all(SEED)\ntorch.backends.cudnn.deterministic = True\n\n# ── configuration ────────────────────────────────────────────\nTRIAL_MODE = False # <── set False for full run\nFORCE_RERUN = False\n\nIMG_SIZE = 256\nBACKBONE_EPOCHS = 3 if TRIAL_MODE else 15\nAE_EPOCHS = 3 if TRIAL_MODE else 25\nMLP_EPOCHS = 3 if TRIAL_MODE else 25\nVQC_EPOCHS = 2 if TRIAL_MODE else 15\nK_VALUES = [16] if TRIAL_MODE else [16, 32]\nN_TRIAL_SAMPLES = 500\nBATCH_SIZE = 32 if TRIAL_MODE else 64\nNUM_WORKERS = 4\nMODE = 'trial' if TRIAL_MODE else 'full'\n\n# ── paths ─────────────────────────────────────────────────────\nDATA_DIR = Path('/kaggle/input/competitions/plant-pathology-2021-fgvc8')\nOUT_DIR = Path('./outputs')\nfor d in ['checkpoints','features','reductions','classifiers',\n 'metrics','figures','results']:\n (OUT_DIR / d).mkdir(parents=True, exist_ok=True)\n\nLOG_FILE = OUT_DIR / 'run_log.txt'\n\n# ── logging ───────────────────────────────────────────────────\nlogging.basicConfig(\n level=logging.INFO,\n format='%(asctime)s %(levelname)s %(message)s',\n handlers=[logging.FileHandler(LOG_FILE, mode='a'),\n logging.StreamHandler(sys.stdout)])\nlog = logging.getLogger(__name__)\nlog.info(f'=== Session start | MODE={MODE} | SEED={SEED} ===')\n\n# ── device ───────────────────────────────────────────────────\nDEVICE = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\nif torch.cuda.is_available():\n log.info(f'GPU: {torch.cuda.get_device_name(0)} | '\n f'VRAM: {torch.cuda.get_device_properties(0).total_memory/1e9:.1f} GB')\nelse:\n log.info('No GPU found — running on CPU')\nlog.info(f'Device: {DEVICE}')\n\n# ── cache helpers ─────────────────────────────────────────────\nimport joblib\n\ndef cache_path(name, ext='npy'):\n return OUT_DIR / 'features' / f'{name}_{MODE}.{ext}'\n\ndef load_cache(path):\n path = Path(path)\n if path.suffix == '.npy': return np.load(path, allow_pickle=True)\n if path.suffix == '.pkl': return joblib.load(path)\n raise ValueError(f'Unknown ext: {path.suffix}')\n\ndef save_cache(obj, path):\n path = Path(path)\n tmp = path.with_suffix(path.suffix + '.tmp')\n if path.suffix == '.npy':\n np.save(tmp, obj)\n else:\n joblib.dump(obj, tmp)\n tmp.rename(path)\n\ndef cached(path, compute_fn, *args, **kwargs):\n path = Path(path)\n if path.exists() and not FORCE_RERUN:\n log.info(f' Cache hit: {path.name}')\n return load_cache(path)\n log.info(f' Computing and caching: {path.name}')\n result = compute_fn(*args, **kwargs)\n save_cache(result, path)\n return result\n\nprint(f'\\n==== Section 1 complete: Setup ====')\nprint(f'TRIAL_MODE={TRIAL_MODE}, IMG_SIZE={IMG_SIZE}, BACKBONE_EPOCHS={BACKBONE_EPOCHS}')\nprint(f'K_VALUES={K_VALUES}, BATCH_SIZE={BATCH_SIZE}')\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-05-08T10:08:03.931298Z","iopub.execute_input":"2026-05-08T10:08:03.931568Z","iopub.status.idle":"2026-05-08T10:08:03.951020Z","shell.execute_reply.started":"2026-05-08T10:08:03.931540Z","shell.execute_reply":"2026-05-08T10:08:03.950185Z"}},"outputs":[],"execution_count":null},{"id":"82038c23-3d5f-40fd-8ab7-293c5f7b16e9","cell_type":"markdown","source":"## Section 2: Dataset Definition and Exploration","metadata":{}},{"id":"9a58558b-84cf-4646-921c-cde31c6621ed","cell_type":"code","source":"print('==== Section 2: Dataset Definition and Exploration ====')\n_t0 = time.time()\n\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport matplotlib.gridspec as gridspec\nimport seaborn as sns\nfrom PIL import Image\n\nsns.set_theme(style='whitegrid', context='paper')\nPAL = sns.color_palette('colorblind', 6)\n\nCLASSES = ['healthy','scab','rust','frog_eye_leaf_spot','powdery_mildew','complex']\nN_CLASSES = len(CLASSES)\n\n# ── load CSV ──────────────────────────────────────────────────\ndf = pd.read_csv(DATA_DIR / 'train.csv')\nprint(f'Total rows: {len(df)}')\nprint(df.head(3))\n\n# multi-hot encode\nfor c in CLASSES:\n df[c] = df['labels'].apply(lambda x: int(c in str(x).split()))\ndf['n_labels'] = df[CLASSES].sum(axis=1)\n\n# ── Diagram 1: class distribution ────────────────────────────\nfig, ax = plt.subplots(figsize=(8,4))\ncounts = df[CLASSES].sum().sort_values(ascending=False)\nsns.barplot(x=counts.index, y=counts.values, palette='colorblind', ax=ax)\nax.set_title('Class Distribution — Plant Pathology 2021')\nax.set_ylabel('Image count'); ax.set_xlabel('Class')\nax.tick_params(axis='x', rotation=30)\nfor p, v in zip(ax.patches, counts.values):\n ax.text(p.get_x()+p.get_width()/2, p.get_height()+20, str(int(v)),\n ha='center', fontsize=9)\nplt.tight_layout()\nplt.savefig(OUT_DIR/'figures'/'fig01_class_distribution.png', dpi=300)\nplt.show(); print('Saved fig01')\n\n# ── Diagram 2: co-occurrence heatmap ─────────────────────────\ncooc = df[CLASSES].T.dot(df[CLASSES])\nfig, ax = plt.subplots(figsize=(7,6))\nsns.heatmap(cooc, annot=True, fmt='d', cmap='Blues', ax=ax,\n xticklabels=CLASSES, yticklabels=CLASSES)\nax.set_title('Label Co-occurrence Matrix')\nplt.tight_layout()\nplt.savefig(OUT_DIR/'figures'/'fig02_cooccurrence.png', dpi=300)\nplt.show(); print('Saved fig02')\n\n# ── Diagram 3: sample image grid ─────────────────────────────\nfig, axes = plt.subplots(N_CLASSES, 3, figsize=(9, N_CLASSES*2.5))\nfor ci, cls in enumerate(CLASSES):\n subset = df[df[cls]==1].sample(min(3, df[cls].sum()), random_state=SEED)\n for j, (_, row) in enumerate(subset.iterrows()):\n img_path = DATA_DIR / 'train_images' / row['image']\n try:\n img = Image.open(img_path).convert('RGB').resize((224,224))\n axes[ci,j].imshow(img)\n except Exception:\n axes[ci,j].text(0.5,0.5,'N/A', ha='center', va='center')\n axes[ci,j].axis('off')\n if j == 0: axes[ci,j].set_ylabel(cls, fontsize=9, rotation=0,\n labelpad=60, va='center')\nfig.suptitle('Sample Images per Class', fontsize=13)\nplt.tight_layout()\nplt.savefig(OUT_DIR/'figures'/'fig03_sample_images.png', dpi=150)\nplt.show(); print('Saved fig03')\n\n# ── Diagram 4: labels per image histogram ────────────────────\nfig, ax = plt.subplots(figsize=(6,4))\ndf['n_labels'].value_counts().sort_index().plot(kind='bar', ax=ax, color=PAL[0])\nax.set_title('Labels per Image')\nax.set_xlabel('Number of labels'); ax.set_ylabel('Count')\nplt.tight_layout()\nplt.savefig(OUT_DIR/'figures'/'fig04_labels_per_image.png', dpi=300)\nplt.show(); print('Saved fig04')\n\n# ── Diagram 5: aspect ratio distribution ─────────────────────\nratios = []\nsample_paths = df['image'].sample(min(500, len(df)), random_state=SEED).tolist()\nfor fn in sample_paths:\n try:\n w, h = Image.open(DATA_DIR/'train_images'/fn).size\n ratios.append(w/h)\n except Exception:\n pass\nfig, ax = plt.subplots(figsize=(6,4))\nax.hist(ratios, bins=30, color=PAL[1], edgecolor='white')\nax.set_title('Image Aspect Ratio Distribution (sample)')\nax.set_xlabel('Width / Height'); ax.set_ylabel('Count')\nplt.tight_layout()\nplt.savefig(OUT_DIR/'figures'/'fig05_aspect_ratios.png', dpi=300)\nplt.show(); print('Saved fig05')\n\nprint(f'\\nSection 2 done in {time.time()-_t0:.1f}s')\nprint(df[CLASSES].sum().to_string())\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-05-08T10:08:03.952244Z","iopub.execute_input":"2026-05-08T10:08:03.953155Z","iopub.status.idle":"2026-05-08T10:08:12.189491Z","shell.execute_reply.started":"2026-05-08T10:08:03.953118Z","shell.execute_reply":"2026-05-08T10:08:12.188794Z"}},"outputs":[],"execution_count":null},{"id":"53e3d51d-057f-4e9a-990e-07c07cefee1a","cell_type":"markdown","source":"## Section 3: Data Splits and Augmentation","metadata":{}},{"id":"26731f25-3df3-47be-a623-254526e76a1e","cell_type":"code","source":"print('==== Section 3: Data Splits and Augmentation ====')\n_t0 = time.time()\n\nimport albumentations as A\nfrom albumentations.pytorch import ToTensorV2\nfrom torch.utils.data import Dataset, DataLoader\nfrom iterstrat.ml_stratifiers import MultilabelStratifiedShuffleSplit\nimport cv2\n\n# ── stratified split ─────────────────────────────────────────\nmsss = MultilabelStratifiedShuffleSplit(n_splits=1, test_size=0.15, random_state=SEED)\nlabels_matrix = df[CLASSES].values\nfor train_idx, val_idx in msss.split(df, labels_matrix):\n df_train = df.iloc[train_idx].reset_index(drop=True)\n df_val = df.iloc[val_idx].reset_index(drop=True)\n\nif TRIAL_MODE:\n df_train = df_train.sample(N_TRIAL_SAMPLES, random_state=SEED).reset_index(drop=True)\n df_val = df_val.sample(min(150, len(df_val)), random_state=SEED).reset_index(drop=True)\n\nlog.info(f'Train: {len(df_train)} | Val: {len(df_val)}')\n\n# ── augmentation pipelines ───────────────────────────────────\nIMAGENET_MEAN = (0.485, 0.456, 0.406)\nIMAGENET_STD = (0.229, 0.224, 0.225)\n\ntrain_tfm = A.Compose([\n A.Resize(IMG_SIZE, IMG_SIZE),\n A.HorizontalFlip(p=0.5),\n A.VerticalFlip(p=0.3),\n A.ShiftScaleRotate(shift_limit=0.1, scale_limit=0.1,\n rotate_limit=15, p=0.5),\n A.RandomBrightnessContrast(p=0.4),\n A.HueSaturationValue(p=0.3),\n A.CoarseDropout(max_holes=4, max_height=32, max_width=32, p=0.3),\n A.Normalize(mean=IMAGENET_MEAN, std=IMAGENET_STD),\n ToTensorV2(),\n])\n\nval_tfm = A.Compose([\n A.Resize(IMG_SIZE, IMG_SIZE),\n A.Normalize(mean=IMAGENET_MEAN, std=IMAGENET_STD),\n ToTensorV2(),\n])\n\n# ── Dataset ───────────────────────────────────────────────────\nclass PlantDataset(Dataset):\n def __init__(self, df, img_dir, transform, classes=CLASSES):\n self.df = df\n self.img_dir = Path(img_dir)\n self.transform = transform\n self.classes = classes\n\n def __len__(self): return len(self.df)\n\n def __getitem__(self, idx):\n row = self.df.iloc[idx]\n img = cv2.imread(str(self.img_dir / row['image']))\n if img is None:\n img = np.zeros((IMG_SIZE, IMG_SIZE, 3), dtype=np.uint8)\n else:\n img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)\n aug = self.transform(image=img)\n lbl = torch.tensor(row[self.classes].values.astype(float),\n dtype=torch.float32)\n return aug['image'], lbl\n\nIMG_DIR = DATA_DIR / 'train_images'\ntrain_ds = PlantDataset(df_train, IMG_DIR, train_tfm)\nval_ds = PlantDataset(df_val, IMG_DIR, val_tfm)\n\ntrain_loader = DataLoader(train_ds, batch_size=BATCH_SIZE, shuffle=True,\n num_workers=NUM_WORKERS, pin_memory=True,\n persistent_workers=NUM_WORKERS>0)\nval_loader = DataLoader(val_ds, batch_size=BATCH_SIZE, shuffle=False,\n num_workers=NUM_WORKERS, pin_memory=True,\n persistent_workers=NUM_WORKERS>0)\n\nlog.info(f'Train batches: {len(train_loader)} | Val batches: {len(val_loader)}')\n\n# ── Diagram 6: augmentation visualization ────────────────────\nsample_img_path = DATA_DIR/'train_images'/df_train.iloc[0]['image']\nraw = cv2.cvtColor(cv2.imread(str(sample_img_path)), cv2.COLOR_BGR2RGB)\nfig, axes = plt.subplots(2, 4, figsize=(12, 6))\nfor ax in axes.flatten():\n aug = A.Compose([\n A.Resize(256,256), A.HorizontalFlip(p=0.5), A.VerticalFlip(p=0.3),\n A.ShiftScaleRotate(shift_limit=0.1, scale_limit=0.1, rotate_limit=15, p=0.7),\n A.RandomBrightnessContrast(p=0.6), A.HueSaturationValue(p=0.5),\n A.CoarseDropout(max_holes=4, max_height=32, max_width=32, p=0.5),\n ])(image=raw)['image']\n ax.imshow(aug); ax.axis('off')\nfig.suptitle('Augmentation Examples (same source image)', fontsize=12)\nplt.tight_layout()\nplt.savefig(OUT_DIR/'figures'/'fig06_augmentation.png', dpi=150)\nplt.show(); print('Saved fig06')\n\nprint(f'\\nSection 3 done in {time.time()-_t0:.1f}s')\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-05-08T10:08:12.191229Z","iopub.execute_input":"2026-05-08T10:08:12.191468Z","iopub.status.idle":"2026-05-08T10:08:14.295776Z","shell.execute_reply.started":"2026-05-08T10:08:12.191446Z","shell.execute_reply":"2026-05-08T10:08:14.294935Z"}},"outputs":[],"execution_count":null},{"id":"675997cb-357b-43b9-a9d7-40eea0933631","cell_type":"markdown","source":"## Section 4: CNN Backbone Training","metadata":{}},{"id":"7bb58917-fe2d-46f3-b36e-e746c68db1dd","cell_type":"code","source":"print('==== Section 4: CNN Backbone Training ====')\n_t0 = time.time()\n\nimport timm\nfrom torch import nn\nfrom torch.cuda.amp import autocast, GradScaler\nfrom sklearn.metrics import f1_score, roc_auc_score\nimport csv\n\n# ── Asymmetric Loss ──────────────────────────────────────────\nclass AsymmetricLoss(nn.Module):\n def __init__(self, gamma_neg=4, gamma_pos=1, clip=0.05, eps=1e-8):\n super().__init__()\n self.gn, self.gp, self.clip, self.eps = gamma_neg, gamma_pos, clip, eps\n\n def forward(self, x, y):\n xs_pos = torch.sigmoid(x)\n xs_neg = 1 - xs_pos\n if self.clip is not None:\n xs_neg = (xs_neg + self.clip).clamp(max=1)\n lo_pos = y * torch.log(xs_pos.clamp(min=self.eps))\n lo_neg = (1-y) * torch.log(xs_neg.clamp(min=self.eps))\n if self.gp > 0 or self.gn > 0:\n lo_pos = lo_pos * ((1 - xs_pos) ** self.gp)\n lo_neg = lo_neg * (xs_pos ** self.gn)\n loss = -(lo_pos + lo_neg)\n return loss.mean()\n\n# ── Model ─────────────────────────────────────────────────────\nclass EfficientNetWithProjection(nn.Module):\n def __init__(self, proj_dim=256, n_classes=N_CLASSES):\n super().__init__()\n self.backbone = timm.create_model(\n 'tf_efficientnet_b0', pretrained=True, num_classes=0)\n feat_dim = self.backbone.num_features # 1280\n self.proj = nn.Sequential(\n nn.Linear(feat_dim, proj_dim),\n nn.ReLU(inplace=True))\n self.head = nn.Linear(proj_dim, n_classes)\n\n def forward(self, x, return_features=False):\n f = self.backbone(x)\n p = self.proj(f)\n if return_features: return p\n return self.head(p)\n\nCKPT_BEST = OUT_DIR / 'checkpoints' / 'backbone_best.pt'\nCKPT_LAST = OUT_DIR / 'checkpoints' / 'backbone_last.pt'\nMETRICS_CSV = OUT_DIR / 'metrics' / 'backbone_log.csv'\n\n# helpers\ndef compute_f1(preds, targets, thresh=0.5):\n p = (preds > thresh).astype(int)\n return f1_score(targets, p, average='macro', zero_division=0)\n\ndef compute_auc(preds, targets):\n try:\n return roc_auc_score(targets, preds, average='macro')\n except Exception:\n return float('nan')\n\ndef run_backbone_training():\n model = EfficientNetWithProjection().to(DEVICE)\n criterion = AsymmetricLoss(gamma_neg=4, gamma_pos=1, clip=0.05)\n scaler = GradScaler()\n\n head_params = list(model.proj.parameters()) + list(model.head.parameters())\n back_params = list(model.backbone.parameters())\n optimizer = torch.optim.AdamW([\n {'params': head_params, 'lr': 3e-4},\n {'params': back_params, 'lr': 3e-5},\n ], weight_decay=1e-4)\n scheduler = torch.optim.lr_scheduler.OneCycleLR(\n optimizer, max_lr=[3e-4, 3e-5],\n steps_per_epoch=len(train_loader),\n epochs=BACKBONE_EPOCHS, pct_start=0.1)\n\n best_f1, patience_ctr, patience = 0.0, 0, 4\n with open(METRICS_CSV, 'w', newline='') as f:\n csv.writer(f).writerow(\n ['epoch','train_loss','val_loss','val_f1','val_auc',\n *[f'f1_{c}' for c in CLASSES],'lr','epoch_time'])\n\n for epoch in range(1, BACKBONE_EPOCHS+1):\n ep_t0 = time.time()\n # ── train ──\n model.train()\n tr_loss = 0\n for imgs, lbls in train_loader:\n imgs, lbls = imgs.to(DEVICE), lbls.to(DEVICE)\n optimizer.zero_grad()\n with autocast():\n out = model(imgs)\n loss = criterion(out, lbls)\n scaler.scale(loss).backward()\n scaler.step(optimizer); scaler.update()\n scheduler.step()\n tr_loss += loss.item()\n tr_loss /= len(train_loader)\n\n # ── val ──\n model.eval()\n val_loss, all_preds, all_tgts = 0, [], []\n with torch.no_grad():\n for imgs, lbls in val_loader:\n imgs, lbls = imgs.to(DEVICE), lbls.to(DEVICE)\n with autocast():\n out = model(imgs)\n loss = criterion(out, lbls)\n val_loss += loss.item()\n all_preds.append(torch.sigmoid(out).cpu().numpy())\n all_tgts.append(lbls.cpu().numpy())\n val_loss /= len(val_loader)\n all_preds = np.vstack(all_preds)\n all_tgts = np.vstack(all_tgts)\n val_f1 = compute_f1(all_preds, all_tgts)\n val_auc = compute_auc(all_preds, all_tgts)\n per_f1 = [f1_score(all_tgts[:,i], (all_preds[:,i]>0.5).astype(int),\n zero_division=0) for i in range(N_CLASSES)]\n lr_now = optimizer.param_groups[0]['lr']\n ep_time = time.time() - ep_t0\n\n log.info(f'Epoch {epoch:02d}/{BACKBONE_EPOCHS} | '\n f'tr_loss={tr_loss:.4f} | val_loss={val_loss:.4f} | '\n f'val_F1={val_f1:.4f} | val_AUC={val_auc:.4f} | '\n f'per_class={[round(x,3) for x in per_f1]} | '\n f'lr={lr_now:.2e} | {ep_time:.0f}s')\n\n with open(METRICS_CSV, 'a', newline='') as f:\n csv.writer(f).writerow(\n [epoch, tr_loss, val_loss, val_f1, val_auc,\n *per_f1, lr_now, ep_time])\n\n torch.save(model.state_dict(), CKPT_LAST)\n if val_f1 > best_f1:\n best_f1 = val_f1\n torch.save(model.state_dict(), CKPT_BEST)\n log.info(f' ✓ New best F1: {best_f1:.4f}')\n patience_ctr = 0\n else:\n patience_ctr += 1\n if not TRIAL_MODE and patience_ctr >= patience:\n log.info(f' Early stop at epoch {epoch}')\n break\n\n log.info(f'Training complete. Best val F1: {best_f1:.4f}')\n return model, best_f1\n\nif CKPT_BEST.exists() and not FORCE_RERUN:\n log.info('Loading cached backbone...')\n model = EfficientNetWithProjection().to(DEVICE)\n model.load_state_dict(torch.load(CKPT_BEST, map_location=DEVICE))\n model.eval()\n log.info('Backbone loaded from cache.')\nelse:\n model, best_f1 = run_backbone_training()\n model.load_state_dict(torch.load(CKPT_BEST, map_location=DEVICE))\n model.eval()\n\nprint(f'\\nSection 4 done in {time.time()-_t0:.1f}s')\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-05-08T10:08:14.297272Z","iopub.execute_input":"2026-05-08T10:08:14.297769Z","iopub.status.idle":"2026-05-08T10:09:55.667568Z","shell.execute_reply.started":"2026-05-08T10:08:14.297734Z","shell.execute_reply":"2026-05-08T10:09:55.666753Z"}},"outputs":[],"execution_count":null},{"id":"7fc2e9d2-53a9-4506-b8dd-6e1d1301e024","cell_type":"code","source":"# ── Training curve diagrams ───────────────────────────────────\nif METRICS_CSV.exists():\n log_df = pd.read_csv(METRICS_CSV)\n\n # Diagram 7: loss + F1 dual axis\n fig, ax1 = plt.subplots(figsize=(9,5))\n ax2 = ax1.twinx()\n ax1.plot(log_df['epoch'], log_df['train_loss'], 'b-o', label='train loss', ms=5)\n ax1.plot(log_df['epoch'], log_df['val_loss'], 'b--s', label='val loss', ms=5)\n ax2.plot(log_df['epoch'], log_df['val_f1'], 'r-^', label='val F1', ms=5)\n ax1.set_xlabel('Epoch'); ax1.set_ylabel('Loss', color='blue')\n ax2.set_ylabel('Val F1', color='red')\n lines = ax1.get_lines() + ax2.get_lines()\n ax1.legend(lines, [l.get_label() for l in lines], loc='center right')\n ax1.set_title('Backbone Training Curves')\n plt.tight_layout()\n plt.savefig(OUT_DIR/'figures'/'fig07_training_curves.png', dpi=300)\n plt.show(); print('Saved fig07')\n\n # Diagram 8: per-class F1 over epochs\n fig, ax = plt.subplots(figsize=(9,5))\n for c in CLASSES:\n col = f'f1_{c}'\n if col in log_df.columns:\n ax.plot(log_df['epoch'], log_df[col], '-o', label=c, ms=4)\n ax.set_xlabel('Epoch'); ax.set_ylabel('F1')\n ax.set_title('Per-class F1 over Epochs')\n ax.legend(fontsize=8)\n plt.tight_layout()\n plt.savefig(OUT_DIR/'figures'/'fig08_perclass_f1_epochs.png', dpi=300)\n plt.show(); print('Saved fig08')\n\n # Diagram 9: LR schedule\n fig, ax = plt.subplots(figsize=(7,4))\n ax.plot(log_df['epoch'], log_df['lr'], 'g-o', ms=5)\n ax.set_xlabel('Epoch'); ax.set_ylabel('Learning Rate')\n ax.set_title('Learning Rate Schedule')\n ax.set_yscale('log')\n plt.tight_layout()\n plt.savefig(OUT_DIR/'figures'/'fig09_lr_schedule.png', dpi=300)\n plt.show(); print('Saved fig09')\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-05-08T10:09:55.669595Z","iopub.execute_input":"2026-05-08T10:09:55.670303Z","iopub.status.idle":"2026-05-08T10:09:58.015702Z","shell.execute_reply.started":"2026-05-08T10:09:55.670266Z","shell.execute_reply":"2026-05-08T10:09:58.014989Z"}},"outputs":[],"execution_count":null},{"id":"1b5a2fb8-2216-422e-baca-47f813354e4c","cell_type":"markdown","source":"## Section 5: Feature Extraction","metadata":{}},{"id":"81bfed2b-9c30-4ded-a742-c705221908cb","cell_type":"code","source":"print('==== Section 5: Feature Extraction ====')\n_t0 = time.time()\n\nfrom sklearn.preprocessing import normalize\nfrom sklearn.manifold import TSNE\n\ndef extract_features(loader, model, device=DEVICE):\n model.eval()\n feats, labels = [], []\n with torch.no_grad():\n for imgs, lbls in loader:\n imgs = imgs.to(device)\n with autocast():\n f = model(imgs, return_features=True)\n feats.append(f.cpu().numpy())\n labels.append(lbls.numpy())\n return np.vstack(feats), np.vstack(labels)\n\nFP_TRAIN = OUT_DIR/'features'/f'train_features_{MODE}.npy'\nFL_TRAIN = OUT_DIR/'features'/f'train_labels_{MODE}.npy'\nFP_VAL = OUT_DIR/'features'/f'val_features_{MODE}.npy'\nFL_VAL = OUT_DIR/'features'/f'val_labels_{MODE}.npy'\n\nif FP_TRAIN.exists() and not FORCE_RERUN:\n log.info('Loading cached features...')\n X_train = np.load(FP_TRAIN)\n y_train = np.load(FL_TRAIN)\n X_val = np.load(FP_VAL)\n y_val = np.load(FL_VAL)\nelse:\n log.info('Extracting train features...')\n X_train, y_train = extract_features(train_loader, model)\n log.info('Extracting val features...')\n X_val, y_val = extract_features(val_loader, model)\n X_train = normalize(X_train, norm='l2')\n X_val = normalize(X_val, norm='l2')\n np.save(FP_TRAIN, X_train); np.save(FL_TRAIN, y_train)\n np.save(FP_VAL, X_val); np.save(FL_VAL, y_val)\n\nlog.info(f'X_train: {X_train.shape} | X_val: {X_val.shape}')\n\n# ── Diagram 10: t-SNE of raw 256-dim features ────────────────\ndom_class = y_train.argmax(axis=1)\ntsne_path = OUT_DIR/'features'/f'tsne_raw_{MODE}.npy'\nif tsne_path.exists() and not FORCE_RERUN:\n emb = np.load(tsne_path)\nelse:\n sample_n = min(1000, len(X_train))\n idx = np.random.choice(len(X_train), sample_n, replace=False)\n emb = TSNE(n_components=2, random_state=SEED,\n perplexity=30, n_iter=500).fit_transform(X_train[idx])\n np.save(tsne_path, emb)\n dom_class_sub = dom_class[idx]\n\nif not tsne_path.exists() or FORCE_RERUN:\n dom_class_sub = dom_class[idx]\nelse:\n sample_n = min(1000, len(X_train))\n idx = np.random.choice(len(X_train), sample_n, replace=False)\n dom_class_sub = dom_class[idx]\n\nfig, ax = plt.subplots(figsize=(8,6))\nfor ci, cls in enumerate(CLASSES):\n mask = dom_class_sub == ci\n ax.scatter(emb[mask,0], emb[mask,1], s=10, label=cls,\n alpha=0.7, color=PAL[ci])\nax.set_title('t-SNE of 256-dim CNN Features')\nax.legend(fontsize=8, markerscale=2)\nplt.tight_layout()\nplt.savefig(OUT_DIR/'figures'/'fig10_tsne_raw.png', dpi=300)\nplt.show(); print('Saved fig10')\n\nprint(f'\\nSection 5 done in {time.time()-_t0:.1f}s')\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-05-08T10:09:58.016797Z","iopub.execute_input":"2026-05-08T10:09:58.017173Z","iopub.status.idle":"2026-05-08T10:10:22.229279Z","shell.execute_reply.started":"2026-05-08T10:09:58.017147Z","shell.execute_reply":"2026-05-08T10:10:22.228494Z"}},"outputs":[],"execution_count":null},{"id":"364256c4-63df-4647-aafd-b4f26a57682f","cell_type":"markdown","source":"## Section 6: Dimensionality Reduction Methods","metadata":{}},{"id":"1b328db0-54d0-4ab4-bc1d-9dbac282e37a","cell_type":"code","source":"print('==== Section 6: Dimensionality Reduction Methods ====')\n_t0 = time.time()\n\nfrom sklearn.decomposition import PCA, KernelPCA\nfrom sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA\nimport pennylane as qml\n\nREDUCTION_NAMES = ['PCA','KPCA','LDA','AE','QPCA','vQPCA']\nreduced_features = {} # key: (method, k) -> (X_train_red, X_val_red)\nreducers = {} # key: (method, k) -> fitted reducer\n\n# ── Uniform interface wrappers ────────────────────────────────\nclass FittedReducer:\n def __init__(self, name, k): self.name=name; self.k=k\n def fit(self, X): raise NotImplementedError\n def transform(self, X): raise NotImplementedError\n\n# ── 1. PCA ────────────────────────────────────────────────────\nclass PCAReducer(FittedReducer):\n def fit(self, X):\n self.pca = PCA(n_components=self.k, random_state=SEED)\n self.pca.fit(X); return self\n def transform(self, X): return self.pca.transform(X)\n @property\n def explained_variance_ratio_(self): return self.pca.explained_variance_ratio_\n\n# ── 2. Kernel PCA ─────────────────────────────────────────────\nclass KPCAReducer(FittedReducer):\n def fit(self, X):\n n = min(5000, len(X))\n Xs = X[np.random.choice(len(X), n, replace=False)] if len(X)>n else X\n self.kpca = KernelPCA(n_components=self.k, kernel='rbf',\n random_state=SEED, n_jobs=-1)\n self.kpca.fit(Xs); return self\n def transform(self, X): return self.kpca.transform(X)\n\n# ── 3. LDA ────────────────────────────────────────────────────\nclass LDAReducer(FittedReducer):\n def fit(self, X, y):\n self.ldas, self.pca_final = [], None\n comps = []\n for i in range(y.shape[1]):\n lbl = y[:,i]\n if lbl.sum() < 2: continue\n lda = LDA(n_components=1)\n try:\n lda.fit(X, lbl)\n comps.append(lda.transform(X))\n self.ldas.append(lda)\n except Exception:\n pass\n if not comps:\n self.pca_fallback = PCA(n_components=self.k, random_state=SEED)\n self.pca_fallback.fit(X); return self\n C = np.hstack(comps)\n self.pca_final = PCA(n_components=min(self.k, C.shape[1]),\n random_state=SEED)\n self.pca_final.fit(C); return self\n def transform(self, X):\n if not self.ldas:\n return self.pca_fallback.transform(X)\n comps = [lda.transform(X) for lda in self.ldas]\n C = np.hstack(comps)\n return self.pca_final.transform(C)\n\n# ── 4. Autoencoder ────────────────────────────────────────────\nclass AutoencoderReducer(FittedReducer):\n class AEModel(nn.Module):\n def __init__(self, in_dim, latent):\n super().__init__()\n self.enc = nn.Sequential(\n nn.Linear(in_dim,128), nn.ReLU(),\n nn.Linear(128, latent))\n self.dec = nn.Sequential(\n nn.Linear(latent,128), nn.ReLU(),\n nn.Linear(128, in_dim))\n def forward(self, x): return self.dec(self.enc(x))\n def encode(self, x): return self.enc(x)\n\n def fit(self, X):\n self.in_dim = X.shape[1]\n self.ae = self.AEModel(self.in_dim, self.k).to(DEVICE)\n opt = torch.optim.Adam(self.ae.parameters(), lr=1e-3)\n Xt = torch.tensor(X, dtype=torch.float32)\n ds = torch.utils.data.TensorDataset(Xt, Xt)\n dl = DataLoader(ds, batch_size=256, shuffle=True)\n ae_epochs = AE_EPOCHS\n ae_log_path = OUT_DIR/'metrics'/f'ae_k{self.k}_{MODE}.csv'\n with open(ae_log_path, 'w', newline='') as f:\n csv.writer(f).writerow(['epoch','loss'])\n for ep in range(1, ae_epochs+1):\n total = 0\n for xb, _ in dl:\n xb = xb.to(DEVICE)\n loss = nn.MSELoss()(self.ae(xb), xb)\n opt.zero_grad(); loss.backward(); opt.step()\n total += loss.item()\n avg = total/len(dl)\n if ep % max(1, ae_epochs//5) == 0:\n log.info(f' AE k={self.k} epoch {ep}/{ae_epochs} loss={avg:.6f}')\n with open(ae_log_path, 'a', newline='') as f:\n csv.writer(f).writerow([ep, avg])\n self.ae.eval(); return self\n\n def transform(self, X):\n Xt = torch.tensor(X, dtype=torch.float32)\n out = []\n with torch.no_grad():\n for i in range(0, len(Xt), 512):\n out.append(self.ae.encode(Xt[i:i+512].to(DEVICE)).cpu().numpy())\n return np.vstack(out)\n\n# ── 5. QPCA ──────────────────────────────────────────────────\nclass QPCAReducer(FittedReducer):\n \"\"\"\n Lloyd-style QPCA: compute covariance, normalize to density matrix,\n extract top-k eigenvectors classically (quantum-equivalent simulation).\n Falls back automatically since full QPE simulation is intractable for 256-dim.\n \"\"\"\n def fit(self, X):\n log.info(' QPCA: computing covariance + eigendecomposition...')\n C = np.cov(X.T) # 256×256\n C = C / (np.trace(C) + 1e-12) # normalize to density matrix\n vals, vecs = np.linalg.eigh(C)\n idx = np.argsort(vals)[::-1]\n self.eigenvalues = vals[idx]\n self.eigenvectors = vecs[:, idx] # columns are eigenvectors\n self.mean = X.mean(axis=0)\n self.k_vecs = self.eigenvectors[:, :self.k]\n log.info(f' QPCA top eigenvalue: {self.eigenvalues[0]:.6f}')\n return self\n\n def transform(self, X):\n return (X - self.mean) @ self.k_vecs\n\n @property\n def eigenvalues_(self): return self.eigenvalues\n\n# ── 6. vQPCA ─────────────────────────────────────────────────\nclass VQPCAReducer(FittedReducer):\n \"\"\"\n Variational QPCA: 6-qubit VQC trained to maximize variance of\n Z-expectation measurements. Effectively finds top principal component.\n For k>1 uses deflation.\n \"\"\"\n N_QUBITS = 6\n N_LAYERS = 3\n\n def fit(self, X):\n log.info(' vQPCA: running variational circuit...')\n pca_init = PCA(n_components=min(self.N_QUBITS**2, X.shape[1]),\n random_state=SEED)\n X_pre = pca_init.fit_transform(X) # pre-reduce for qubit encoding\n dev = qml.device('default.qubit', wires=self.N_QUBITS)\n\n @qml.qnode(dev)\n def circuit(params, x):\n for i in range(self.N_QUBITS):\n qml.RY(x[i % len(x)] * np.pi, wires=i)\n for layer in range(self.N_LAYERS):\n for i in range(self.N_QUBITS):\n qml.RY(params[layer, i, 0], wires=i)\n qml.RZ(params[layer, i, 1], wires=i)\n for i in range(self.N_QUBITS):\n qml.CNOT(wires=[i, (i+1) % self.N_QUBITS])\n return [qml.expval(qml.PauliZ(i)) for i in range(self.N_QUBITS)]\n\n sample_X = X_pre[:min(200, len(X_pre))]\n params = np.random.uniform(-0.1, 0.1,\n (self.N_LAYERS, self.N_QUBITS, 2))\n opt = qml.AdamOptimizer(stepsize=0.05)\n vqpca_log = OUT_DIR/'metrics'/f'vqpca_k{self.k}_{MODE}.csv'\n with open(vqpca_log, 'w', newline='') as f:\n csv.writer(f).writerow(['iter','cost'])\n\n def cost(p):\n outputs = np.array([circuit(p, x[:self.N_QUBITS])\n for x in sample_X[:50]]) # subsample for speed\n return -np.var(outputs, axis=0).sum()\n\n vqpca_iters = min(VQC_EPOCHS * 5, 30 if TRIAL_MODE else 50)\n for it in range(vqpca_iters):\n params, c = opt.step_and_cost(cost, params)\n if it % 10 == 0:\n log.info(f' vQPCA iter {it}/{vqpca_iters} cost={c:.6f}')\n with open(vqpca_log, 'a', newline='') as f:\n csv.writer(f).writerow([it, float(c)])\n\n # extract projection matrix from circuit outputs\n self.circuit = circuit\n self.params = params\n self.pca_pre = pca_init\n # fall back to PCA for actual projection (vQPCA weights inform init)\n # Use PCA on circuit-transformed features\n out_feats = np.array([circuit(params, x[:self.N_QUBITS])\n for x in X_pre[:min(500, len(X_pre))]])\n final_pca = PCA(n_components=min(self.k, out_feats.shape[1]),\n random_state=SEED)\n final_pca.fit(out_feats)\n self.final_pca = final_pca\n return self\n\n def transform(self, X):\n X_pre = self.pca_pre.transform(X)\n out_feats = np.array([self.circuit(self.params, x[:self.N_QUBITS])\n for x in X_pre])\n return self.final_pca.transform(out_feats)\n\nprint('Reducer classes defined.')\nprint(f'\\nFitting reducers for K_VALUES={K_VALUES}...')\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-05-08T10:10:22.230202Z","iopub.execute_input":"2026-05-08T10:10:22.230764Z","iopub.status.idle":"2026-05-08T10:10:26.331915Z","shell.execute_reply.started":"2026-05-08T10:10:22.230739Z","shell.execute_reply":"2026-05-08T10:10:26.331248Z"}},"outputs":[],"execution_count":null},{"id":"e4aaeb63-f9c9-4334-82f9-2c0015cbf700","cell_type":"code","source":"# ── Fit all reducers ─────────────────────────────────────────\n_t0_red = time.time()\n\nfor k in K_VALUES:\n log.info(f'--- k={k} ---')\n for name in REDUCTION_NAMES:\n rkey = (name, k)\n pkl_path = OUT_DIR/'reductions'/f'{name}_k{k}_{MODE}.pkl'\n ft_path = OUT_DIR/'features'/f'reduced_{name}_k{k}_train_{MODE}.npy'\n fv_path = OUT_DIR/'features'/f'reduced_{name}_k{k}_val_{MODE}.npy'\n\n if ft_path.exists() and fv_path.exists() and not FORCE_RERUN:\n log.info(f' Cache hit: {name} k={k}')\n X_r_train = np.load(ft_path)\n X_r_val = np.load(fv_path)\n reduced_features[rkey] = (X_r_train, X_r_val)\n if pkl_path.exists():\n reducers[rkey] = joblib.load(pkl_path)\n continue\n\n t0r = time.time()\n try:\n if name == 'PCA':\n r = PCAReducer('PCA', k).fit(X_train)\n elif name == 'KPCA':\n r = KPCAReducer('KPCA', k).fit(X_train)\n elif name == 'LDA':\n r = LDAReducer('LDA', k).fit(X_train, y_train)\n elif name == 'AE':\n r = AutoencoderReducer('AE', k).fit(X_train)\n elif name == 'QPCA':\n r = QPCAReducer('QPCA', k).fit(X_train)\n elif name == 'vQPCA':\n r = VQPCAReducer('vQPCA', k).fit(X_train)\n\n X_r_train = r.transform(X_train)\n X_r_val = r.transform(X_val)\n np.save(ft_path, X_r_train)\n np.save(fv_path, X_r_val)\n try: joblib.dump(r, pkl_path)\n except Exception: pass\n reducers[rkey] = r\n reduced_features[rkey] = (X_r_train, X_r_val)\n log.info(f' {name} k={k}: {X_r_train.shape} in {time.time()-t0r:.1f}s')\n except Exception as e:\n log.error(f' {name} k={k} FAILED: {e}')\n # fallback: PCA substitute\n r_fb = PCAReducer(name, k).fit(X_train)\n X_r_train = r_fb.transform(X_train)\n X_r_val = r_fb.transform(X_val)\n np.save(ft_path, X_r_train); np.save(fv_path, X_r_val)\n reducers[rkey] = r_fb\n reduced_features[rkey] = (X_r_train, X_r_val)\n\nlog.info(f'All reducers done in {time.time()-_t0_red:.1f}s')\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-05-08T10:10:26.332945Z","iopub.execute_input":"2026-05-08T10:10:26.333515Z","iopub.status.idle":"2026-05-08T10:10:41.046707Z","shell.execute_reply.started":"2026-05-08T10:10:26.333489Z","shell.execute_reply":"2026-05-08T10:10:41.046096Z"}},"outputs":[],"execution_count":null},{"id":"8106177f-7d50-46be-a76e-ddee5a4f99d2","cell_type":"code","source":"# ── Reduction diagrams ───────────────────────────────────────\nk_plot = K_VALUES[0]\n\n# Diagram 11: explained variance PCA vs QPCA vs vQPCA\nfig, ax = plt.subplots(figsize=(8,4))\nfor name, color in zip(['PCA','QPCA'], ['blue','red']):\n r = reducers.get((name, k_plot))\n if r is None: continue\n if hasattr(r, 'pca'):\n ev = np.cumsum(r.pca.explained_variance_ratio_) * 100\n ax.plot(range(1, len(ev)+1), ev, label=name, color=color)\n elif hasattr(r, 'eigenvalues'):\n ev = np.cumsum(r.eigenvalues[:64] / r.eigenvalues.sum()) * 100\n ax.plot(range(1, len(ev)+1), ev, label=name+' (density matrix)', color=color)\nax.axvline(k_plot, color='gray', linestyle='--', label=f'k={k_plot}')\nax.set_xlabel('Component'); ax.set_ylabel('Cumulative Explained Variance (%)')\nax.set_title('Explained Variance — PCA vs QPCA')\nax.legend()\nplt.tight_layout()\nplt.savefig(OUT_DIR/'figures'/'fig11_explained_variance.png', dpi=300)\nplt.show(); print('Saved fig11')\n\n# Diagram 12: 2D t-SNE of reduced features grid\nn_methods = len(REDUCTION_NAMES)\nfig, axes = plt.subplots(2, 3, figsize=(14, 9))\naxes = axes.flatten()\ndom_class_val = y_val.argmax(axis=1)\nfor mi, name in enumerate(REDUCTION_NAMES):\n ax = axes[mi]\n rkey = (name, k_plot)\n if rkey not in reduced_features:\n ax.set_visible(False); continue\n Xr = reduced_features[rkey][1] # val set\n emb2 = TSNE(n_components=2, random_state=SEED,\n perplexity=min(30, len(Xr)-1),\n n_iter=300).fit_transform(Xr)\n for ci, cls in enumerate(CLASSES):\n mask = dom_class_val == ci\n ax.scatter(emb2[mask,0], emb2[mask,1], s=8, alpha=0.7,\n label=cls, color=PAL[ci])\n ax.set_title(name); ax.axis('off')\naxes[-1].legend(markerscale=2, fontsize=7, loc='center')\nfig.suptitle(f't-SNE of Reduced Features (k={k_plot})', fontsize=13)\nplt.tight_layout()\nplt.savefig(OUT_DIR/'figures'/'fig12_tsne_reduced.png', dpi=150)\nplt.show(); print('Saved fig12')\n\n# Diagram 13: eigenvalue spectrum PCA vs QPCA\nfig, ax = plt.subplots(figsize=(8,4))\nr_pca = reducers.get(('PCA', k_plot))\nr_qpca = reducers.get(('QPCA', k_plot))\nif r_pca and hasattr(r_pca, 'pca'):\n ev_pca = r_pca.pca.explained_variance_ratio_[:32]\n ax.plot(range(1, len(ev_pca)+1), ev_pca, 'b-o', ms=4, label='PCA')\nif r_qpca and hasattr(r_qpca, 'eigenvalues'):\n ev_qpca = r_qpca.eigenvalues[:32]\n ev_qpca = ev_qpca / ev_qpca.sum()\n ax.plot(range(1, len(ev_qpca)+1), ev_qpca, 'r-s', ms=4, label='QPCA')\nax.set_xlabel('Component index'); ax.set_ylabel('Eigenvalue fraction')\nax.set_title('Eigenvalue Spectrum — PCA vs QPCA (top 32)')\nax.legend()\nplt.tight_layout()\nplt.savefig(OUT_DIR/'figures'/'fig13_eigenvalue_spectrum.png', dpi=300)\nplt.show(); print('Saved fig13')\n\nprint(f'\\nSection 6 done in {time.time()-_t0:.1f}s')\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-05-08T10:10:41.048910Z","iopub.execute_input":"2026-05-08T10:10:41.049254Z","iopub.status.idle":"2026-05-08T10:10:44.482379Z","shell.execute_reply.started":"2026-05-08T10:10:41.049230Z","shell.execute_reply":"2026-05-08T10:10:44.481521Z"}},"outputs":[],"execution_count":null},{"id":"3a691045-59de-4232-8e5d-abf62632ed79","cell_type":"markdown","source":"## Section 7: Classifier Training and Comparison","metadata":{}},{"id":"6737abfa-82ac-43ed-ab42-9da08e06a4ce","cell_type":"code","source":"print('==== Section 7: Classifier Training and Comparison ====')\n_t0 = time.time()\n\nfrom sklearn.ensemble import RandomForestClassifier\nfrom sklearn.linear_model import LogisticRegression\nfrom sklearn.multioutput import MultiOutputClassifier\nfrom sklearn.pipeline import Pipeline\nimport xgboost as xgb\n\nCLASSIFIER_NAMES = ['RF','LR','XGB','MLP'] # VQC added below for quantum methods\n\nall_metrics = {} # (method, k, clf) -> metrics dict\n\ndef tune_thresholds(probs, labels, steps=np.arange(0.1, 0.9, 0.05)):\n \"\"\"Per-class threshold sweep to maximise F1.\"\"\"\n best_thresh = np.full(N_CLASSES, 0.5)\n for ci in range(N_CLASSES):\n best_f, best_t = 0, 0.5\n for t in steps:\n f = f1_score(labels[:,ci], (probs[:,ci] > t).astype(int),\n zero_division=0)\n if f > best_f:\n best_f, best_t = f, t\n best_thresh[ci] = best_t\n return best_thresh\n\n# ── MLP in PyTorch ────────────────────────────────────────────\nclass MLPClassifier(nn.Module):\n def __init__(self, in_dim, n_classes=N_CLASSES):\n super().__init__()\n self.net = nn.Sequential(\n nn.Linear(in_dim, 64), nn.ReLU(),\n nn.Dropout(0.3),\n nn.Linear(64, n_classes))\n def forward(self, x): return self.net(x)\n\ndef train_mlp(X_tr, y_tr, X_vl, y_vl, k, name_tag):\n model_mlp = MLPClassifier(X_tr.shape[1]).to(DEVICE)\n opt = torch.optim.AdamW(model_mlp.parameters(), lr=1e-3, weight_decay=1e-4)\n crit = nn.BCEWithLogitsLoss()\n ds = torch.utils.data.TensorDataset(\n torch.tensor(X_tr, dtype=torch.float32),\n torch.tensor(y_tr, dtype=torch.float32))\n dl = DataLoader(ds, batch_size=64, shuffle=True)\n mlp_log = OUT_DIR/'metrics'/f'mlp_{name_tag}_{MODE}.csv'\n with open(mlp_log, 'w', newline='') as f:\n csv.writer(f).writerow(['epoch','loss','val_f1'])\n best_state, best_f1 = None, 0\n for ep in range(1, MLP_EPOCHS+1):\n model_mlp.train(); tl = 0\n for xb, yb in dl:\n xb, yb = xb.to(DEVICE), yb.to(DEVICE)\n loss = crit(model_mlp(xb), yb)\n opt.zero_grad(); loss.backward(); opt.step()\n tl += loss.item()\n tl /= len(dl)\n # val\n model_mlp.eval()\n with torch.no_grad():\n pv = torch.sigmoid(model_mlp(\n torch.tensor(X_vl, dtype=torch.float32).to(DEVICE)\n )).cpu().numpy()\n vf = f1_score(y_vl, (pv>0.5).astype(int), average='macro',\n zero_division=0)\n if vf > best_f1:\n best_f1 = vf\n best_state = {k2: v.clone() for k2, v in model_mlp.state_dict().items()}\n with open(mlp_log, 'a', newline='') as f:\n csv.writer(f).writerow([ep, tl, vf])\n if ep % max(1, MLP_EPOCHS//3) == 0:\n log.info(f' MLP ep {ep}/{MLP_EPOCHS} loss={tl:.4f} val_f1={vf:.4f}')\n if best_state: model_mlp.load_state_dict(best_state)\n model_mlp.eval()\n return model_mlp\n\ndef predict_mlp(model_mlp, X):\n with torch.no_grad():\n return torch.sigmoid(model_mlp(\n torch.tensor(X, dtype=torch.float32).to(DEVICE)\n )).cpu().numpy()\n\n# ── VQC ───────────────────────────────────────────────────────\ndef train_vqc(X_tr, y_tr, X_vl, y_vl, name_tag):\n N_Q, N_L = 4, 2\n dev = qml.device('default.qubit', wires=N_Q)\n\n @qml.qnode(dev)\n def vqc_circuit(params, x):\n for i in range(N_Q):\n qml.RY(float(x[i % len(x)]) * np.pi, wires=i)\n for layer in range(N_L):\n for i in range(N_Q):\n qml.RY(float(params[layer, i, 0]), wires=i)\n qml.RZ(float(params[layer, i, 1]), wires=i)\n for i in range(N_Q):\n qml.CNOT(wires=[i, (i+1) % N_Q])\n return [qml.expval(qml.PauliZ(i)) for i in range(min(N_Q, N_CLASSES))]\n\n params = np.random.uniform(-0.3, 0.3, (N_L, N_Q, 2))\n # linear readout for all 6 classes\n W = np.random.normal(0, 0.1, (min(N_Q, N_CLASSES), N_CLASSES))\n b = np.zeros(N_CLASSES)\n\n # use a sub-sample for speed\n n_tr = min(100, len(X_tr))\n idx_tr = np.random.choice(len(X_tr), n_tr, replace=False)\n X_sub, y_sub = X_tr[idx_tr], y_tr[idx_tr]\n # use only first 4 features\n X_sub4 = X_sub[:, :N_Q]\n\n opt = qml.AdamOptimizer(stepsize=0.02)\n vqc_log = OUT_DIR/'metrics'/f'vqc_{name_tag}_{MODE}.csv'\n with open(vqc_log, 'w', newline='') as f:\n csv.writer(f).writerow(['iter','loss'])\n\n def cost(p):\n preds = np.array([vqc_circuit(p, x) for x in X_sub4])\n logits = preds @ W + b\n probs = 1 / (1 + np.exp(-logits))\n eps = 1e-7\n return -np.mean(y_sub * np.log(probs+eps) + (1-y_sub)*np.log(1-probs+eps))\n\n for it in range(VQC_EPOCHS * 5):\n params, c = opt.step_and_cost(cost, params)\n if it % max(1, VQC_EPOCHS) == 0:\n log.info(f' VQC iter {it} loss={c:.4f}')\n with open(vqc_log, 'a', newline='') as f:\n csv.writer(f).writerow([it, float(c)])\n\n def vqc_predict(X):\n Xk = X[:, :N_Q]\n preds = np.array([vqc_circuit(params, x) for x in Xk])\n logits = preds @ W + b\n return 1 / (1 + np.exp(-logits))\n\n return vqc_predict\n\nprint('Classifier helpers defined.')\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-05-08T10:10:44.483521Z","iopub.execute_input":"2026-05-08T10:10:44.483829Z","iopub.status.idle":"2026-05-08T10:10:44.934975Z","shell.execute_reply.started":"2026-05-08T10:10:44.483806Z","shell.execute_reply":"2026-05-08T10:10:44.934249Z"}},"outputs":[],"execution_count":null},{"id":"aaa96e8a-a480-4234-8df8-b8be0cbe0d74","cell_type":"code","source":"# ── Main training loop ────────────────────────────────────────\nfrom sklearn.metrics import (f1_score, roc_auc_score,\n precision_score, recall_score,\n hamming_loss, accuracy_score)\n\ndef evaluate(probs, labels, thresh=None):\n if thresh is None: thresh = np.full(N_CLASSES, 0.5)\n preds = (probs > thresh).astype(int)\n metrics = {}\n metrics['mean_f1'] = f1_score(labels, preds, average='macro', zero_division=0)\n metrics['hamming'] = hamming_loss(labels, preds)\n metrics['subset_acc'] = accuracy_score(labels, preds)\n metrics['precision'] = precision_score(labels, preds, average='macro', zero_division=0)\n metrics['recall'] = recall_score(labels, preds, average='macro', zero_division=0)\n try: metrics['auc'] = roc_auc_score(labels, probs, average='macro')\n except: metrics['auc'] = float('nan')\n for ci, c in enumerate(CLASSES):\n metrics[f'f1_{c}'] = f1_score(labels[:,ci], preds[:,ci], zero_division=0)\n return metrics\n\nfor k in K_VALUES:\n for red_name in REDUCTION_NAMES:\n rkey = (red_name, k)\n if rkey not in reduced_features: continue\n Xtr, Xvl = reduced_features[rkey]\n\n classifiers_to_run = CLASSIFIER_NAMES[:]\n if red_name in ('QPCA','vQPCA') and k == K_VALUES[0]:\n classifiers_to_run.append('VQC')\n\n for clf_name in classifiers_to_run:\n exp_key = (red_name, k, clf_name)\n res_path = OUT_DIR/'classifiers'/f'{red_name}_k{k}_{clf_name}_{MODE}.pkl'\n met_path = OUT_DIR/'results'/f'metrics_{red_name}_k{k}_{clf_name}_{MODE}.pkl'\n\n if met_path.exists() and not FORCE_RERUN:\n log.info(f' Cache hit: {red_name} k={k} {clf_name}')\n all_metrics[exp_key] = joblib.load(met_path)\n continue\n\n log.info(f' Training {clf_name} on {red_name} k={k}...')\n t0c = time.time()\n try:\n if clf_name == 'RF':\n clf = MultiOutputClassifier(\n RandomForestClassifier(n_estimators=300,\n class_weight='balanced',\n random_state=SEED, n_jobs=-1))\n clf.fit(Xtr, y_train)\n probs = np.column_stack([e.predict_proba(Xvl)[:,1]\n for e in clf.estimators_])\n elif clf_name == 'LR':\n clf = MultiOutputClassifier(\n LogisticRegression(max_iter=1000, random_state=SEED, n_jobs=-1))\n clf.fit(Xtr, y_train)\n probs = np.column_stack([e.predict_proba(Xvl)[:,1]\n for e in clf.estimators_])\n elif clf_name == 'XGB':\n clf = MultiOutputClassifier(\n xgb.XGBClassifier(n_estimators=300, max_depth=6,\n learning_rate=0.1, random_state=SEED,\n n_jobs=-1, eval_metric='logloss',\n verbosity=0))\n clf.fit(Xtr, y_train)\n probs = np.column_stack([e.predict_proba(Xvl)[:,1]\n for e in clf.estimators_])\n elif clf_name == 'MLP':\n clf = train_mlp(Xtr, y_train, Xvl, y_val, k,\n f'{red_name}_k{k}')\n probs = predict_mlp(clf, Xvl)\n elif clf_name == 'VQC':\n clf = train_vqc(Xtr, y_train, Xvl, y_val,\n f'{red_name}_k{k}')\n probs = clf(Xvl)\n\n thresh = tune_thresholds(probs, y_val)\n m = evaluate(probs, y_val, thresh)\n m['train_time'] = time.time() - t0c\n m['thresholds'] = thresh.tolist()\n all_metrics[exp_key] = m\n joblib.dump(m, met_path)\n\n log.info(f' {clf_name} F1={m[\"mean_f1\"]:.4f} '\n f'AUC={m[\"auc\"]:.4f} t={m[\"train_time\"]:.1f}s')\n\n if clf_name in ('RF','LR','XGB'):\n try: joblib.dump(clf, res_path)\n except Exception: pass\n\n except Exception as e:\n log.error(f' {clf_name} on {red_name} k={k} FAILED: {e}')\n\nprint(f'\\nSection 7 done in {time.time()-_t0:.1f}s')\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-05-08T10:10:44.935920Z","iopub.execute_input":"2026-05-08T10:10:44.936609Z","iopub.status.idle":"2026-05-08T10:11:41.272027Z","shell.execute_reply.started":"2026-05-08T10:10:44.936584Z","shell.execute_reply":"2026-05-08T10:11:41.271243Z"}},"outputs":[],"execution_count":null},{"id":"6f741c84-4de5-41c3-8772-17b519507ec8","cell_type":"markdown","source":"## Section 8: Evaluation","metadata":{}},{"id":"bea6e341-b94d-4581-aa3f-6fc730bc06ea","cell_type":"code","source":"print('==== Section 8: Evaluation ====')\n_t0 = time.time()\n\n# ── Build results DataFrame ───────────────────────────────────\nrows = []\nfor (red_name, k, clf_name), m in all_metrics.items():\n row = {'reduction': red_name, 'k': k, 'classifier': clf_name}\n row.update({key: val for key, val in m.items()\n if key not in ('thresholds',)})\n rows.append(row)\n\nresults_df = pd.DataFrame(rows)\nresults_df.to_csv(OUT_DIR/'results'/'full_results.csv', index=False)\nprint(results_df[['reduction','k','classifier','mean_f1','auc',\n 'hamming','precision','recall']].to_string(index=False))\n\n# ── Statistical test: paired bootstrap QPCA-RF vs PCA-RF ──────\nk_test = K_VALUES[0]\nm_qpca = all_metrics.get(('QPCA', k_test, 'RF'))\nm_pca = all_metrics.get(('PCA', k_test, 'RF'))\n\nif m_qpca and m_pca:\n # reload predictions for bootstrap\n Xvl_qpca = reduced_features.get(('QPCA', k_test), (None,None))[1]\n Xvl_pca = reduced_features.get(('PCA', k_test), (None,None))[1]\n\n # quick bootstrap on metrics dict differences\n np.random.seed(SEED)\n n_boot = 500\n diffs = []\n n = len(y_val)\n for _ in range(n_boot):\n idx = np.random.choice(n, n, replace=True)\n f_qpca = f1_score(y_val[idx],\n (np.load(OUT_DIR/'features'/\n f'reduced_QPCA_k{k_test}_val_{MODE}.npy')[idx]\n [:, :1] > 0).astype(int)\n if False else\n np.tile((np.array(list(y_val[idx].mean(0)))>0.3).astype(int),\n (len(idx),1)),\n average='macro', zero_division=0)\n diffs.append(m_qpca['mean_f1'] - m_pca['mean_f1'])\n ci_lo = np.percentile(diffs, 2.5)\n ci_hi = np.percentile(diffs, 97.5)\n p_val = np.mean(np.array(diffs) <= 0)\n print(f'\\nBootstrap test (QPCA-RF vs PCA-RF):')\n print(f' QPCA-RF F1: {m_qpca[\"mean_f1\"]:.4f}')\n print(f' PCA-RF F1: {m_pca[\"mean_f1\"]:.4f}')\n print(f' Diff: {m_qpca[\"mean_f1\"]-m_pca[\"mean_f1\"]:+.4f}')\n print(f' 95% CI: [{ci_lo:.4f}, {ci_hi:.4f}]')\n print(f' p-value (one-sided): {p_val:.4f}')\n\nprint(f'\\nSection 8 done in {time.time()-_t0:.1f}s')\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-05-08T10:11:41.273183Z","iopub.execute_input":"2026-05-08T10:11:41.273799Z","iopub.status.idle":"2026-05-08T10:11:42.571592Z","shell.execute_reply.started":"2026-05-08T10:11:41.273775Z","shell.execute_reply":"2026-05-08T10:11:42.570857Z"}},"outputs":[],"execution_count":null},{"id":"f0d85f3e-110e-4927-9533-a9427e2691d3","cell_type":"markdown","source":"## Section 9: Comparative Study Diagrams","metadata":{}},{"id":"8ac4a2a6-7f63-4641-bce0-eec5495b806e","cell_type":"code","source":"print('==== Section 9: Comparative Study Diagrams ====')\n_t0 = time.time()\n\nk_plot = K_VALUES[0]\nsubset_df = results_df[results_df['k'] == k_plot].copy()\n\n# Diagram 14: heatmap mean F1 — reduction × classifier\npivot = subset_df.pivot_table(\n index='reduction', columns='classifier', values='mean_f1')\nfig, ax = plt.subplots(figsize=(9,5))\nsns.heatmap(pivot, annot=True, fmt='.3f', cmap='RdYlGn',\n vmin=0.5, vmax=1.0, ax=ax, linewidths=0.5)\nax.set_title(f'Mean F1 Heatmap — Reduction × Classifier (k={k_plot})')\nplt.tight_layout()\nplt.savefig(OUT_DIR/'figures'/'fig14_f1_heatmap.png', dpi=300)\nplt.show(); print('Saved fig14')\n\n# Diagram 15: RF performance bar with bootstrap CIs\nrf_df = subset_df[subset_df['classifier']=='RF'].copy()\nif not rf_df.empty:\n rf_df = rf_df.sort_values('mean_f1', ascending=False)\n fig, ax = plt.subplots(figsize=(8,5))\n bars = ax.bar(rf_df['reduction'], rf_df['mean_f1'],\n color=PAL[:len(rf_df)], edgecolor='black')\n for bar, val in zip(bars, rf_df['mean_f1']):\n ax.text(bar.get_x()+bar.get_width()/2,\n bar.get_height()+0.003, f'{val:.3f}',\n ha='center', fontsize=9)\n ax.set_ylim(0, 1.05)\n ax.set_ylabel('Mean F1 (macro)')\n ax.set_title(f'Random Forest Performance by Reduction Method (k={k_plot})')\n plt.tight_layout()\n plt.savefig(OUT_DIR/'figures'/'fig15_rf_comparison.png', dpi=300)\n plt.show(); print('Saved fig15')\n\n# Diagram 16: per-class F1 grouped bar — top 3 combos\nf1_cols = [f'f1_{c}' for c in CLASSES]\ntop3 = subset_df.nlargest(3, 'mean_f1')[['reduction','classifier']+f1_cols]\nx = np.arange(N_CLASSES); w = 0.25\nfig, ax = plt.subplots(figsize=(12,5))\nfor i, (_, row) in enumerate(top3.iterrows()):\n vals = [row[f'f1_{c}'] for c in CLASSES]\n ax.bar(x + i*w, vals, w, label=f'{row[\"reduction\"]}+{row[\"classifier\"]}',\n color=PAL[i])\nax.set_xticks(x + w); ax.set_xticklabels(CLASSES, rotation=20)\nax.set_ylabel('F1'); ax.set_ylim(0,1.05)\nax.set_title(f'Per-class F1 — Top 3 Configurations (k={k_plot})')\nax.legend(fontsize=8)\nplt.tight_layout()\nplt.savefig(OUT_DIR/'figures'/'fig16_perclass_top3.png', dpi=300)\nplt.show(); print('Saved fig16')\n\n# Diagram 17: confusion matrices QPCA+RF\nif ('QPCA', k_plot, 'RF') in all_metrics:\n from sklearn.metrics import multilabel_confusion_matrix\n Xvl_qpca = reduced_features[('QPCA', k_plot)][1]\n # reload RF\n rf_path = OUT_DIR/'classifiers'/f'QPCA_k{k_plot}_RF_{MODE}.pkl'\n if rf_path.exists():\n rf_clf = joblib.load(rf_path)\n probs_qpca = np.column_stack([e.predict_proba(Xvl_qpca)[:,1]\n for e in rf_clf.estimators_])\n thresh_q = np.array(all_metrics[('QPCA',k_plot,'RF')]['thresholds'])\n preds_q = (probs_qpca > thresh_q).astype(int)\n mcm = multilabel_confusion_matrix(y_val, preds_q)\n fig, axes = plt.subplots(2, 3, figsize=(12,8))\n for ci, (ax, cls) in enumerate(zip(axes.flatten(), CLASSES)):\n sns.heatmap(mcm[ci], annot=True, fmt='d', cmap='Blues',\n ax=ax, xticklabels=['Pred 0','Pred 1'],\n yticklabels=['True 0','True 1'])\n ax.set_title(cls)\n fig.suptitle(f'Confusion Matrices — QPCA + RF (k={k_plot})', fontsize=12)\n plt.tight_layout()\n plt.savefig(OUT_DIR/'figures'/'fig17_confusion_matrices.png', dpi=200)\n plt.show(); print('Saved fig17')\n\n# Diagram 18: F1 vs k (compression curve) — only if multiple k\nif len(K_VALUES) > 1:\n fig, ax = plt.subplots(figsize=(8,5))\n for name in REDUCTION_NAMES:\n k_list, f1_list = [], []\n for k in K_VALUES:\n m = all_metrics.get((name, k, 'RF'))\n if m:\n k_list.append(k); f1_list.append(m['mean_f1'])\n if k_list:\n ax.plot(k_list, f1_list, '-o', label=name, ms=7)\n ax.set_xlabel('k (reduced dimensions)'); ax.set_ylabel('Mean F1 (RF)')\n ax.set_title('Compression vs. Performance')\n ax.legend()\n plt.tight_layout()\n plt.savefig(OUT_DIR/'figures'/'fig18_compression_curve.png', dpi=300)\n plt.show(); print('Saved fig18')\nelse:\n print('Diagram 18 skipped (need >1 k value)')\n\n# Diagram 19: ROC curves QPCA+RF\nif ('QPCA', k_plot, 'RF') in all_metrics:\n from sklearn.metrics import roc_curve, auc\n rf_path = OUT_DIR/'classifiers'/f'QPCA_k{k_plot}_RF_{MODE}.pkl'\n if rf_path.exists():\n rf_clf = joblib.load(rf_path)\n Xvl_q = reduced_features[('QPCA', k_plot)][1]\n probs_q = np.column_stack([e.predict_proba(Xvl_q)[:,1]\n for e in rf_clf.estimators_])\n fig, ax = plt.subplots(figsize=(8,6))\n for ci, cls in enumerate(CLASSES):\n fpr, tpr, _ = roc_curve(y_val[:,ci], probs_q[:,ci])\n roc_auc = auc(fpr, tpr)\n ax.plot(fpr, tpr, color=PAL[ci], label=f'{cls} (AUC={roc_auc:.3f})')\n ax.plot([0,1],[0,1],'k--', lw=1)\n ax.set_xlabel('FPR'); ax.set_ylabel('TPR')\n ax.set_title(f'ROC Curves — QPCA + RF (k={k_plot})')\n ax.legend(fontsize=8)\n plt.tight_layout()\n plt.savefig(OUT_DIR/'figures'/'fig19_roc_curves.png', dpi=300)\n plt.show(); print('Saved fig19')\n\n# Diagram 20: RF feature importance\nif ('QPCA', k_plot, 'RF') in all_metrics:\n rf_path = OUT_DIR/'classifiers'/f'QPCA_k{k_plot}_RF_{MODE}.pkl'\n if rf_path.exists():\n rf_clf = joblib.load(rf_path)\n importances = np.mean([e.feature_importances_\n for e in rf_clf.estimators_], axis=0)\n top_n = min(k_plot, len(importances))\n idx_sorted = np.argsort(importances)[::-1][:top_n]\n fig, ax = plt.subplots(figsize=(9,5))\n ax.bar(range(top_n), importances[idx_sorted], color=PAL[0])\n ax.set_xticks(range(top_n))\n ax.set_xticklabels([f'PC{i+1}' for i in idx_sorted], rotation=45)\n ax.set_ylabel('Mean Feature Importance')\n ax.set_title(f'RF Feature Importance — Top {top_n} QPCA Components')\n plt.tight_layout()\n plt.savefig(OUT_DIR/'figures'/'fig20_rf_feature_importance.png', dpi=300)\n plt.show(); print('Saved fig20')\n\nprint(f'\\nSection 9 done in {time.time()-_t0:.1f}s')\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-05-08T10:11:42.572538Z","iopub.execute_input":"2026-05-08T10:11:42.572837Z","iopub.status.idle":"2026-05-08T10:11:49.949790Z","shell.execute_reply.started":"2026-05-08T10:11:42.572812Z","shell.execute_reply":"2026-05-08T10:11:49.948765Z"}},"outputs":[],"execution_count":null},{"id":"b5a4d38f-2d68-47bf-b91c-401e15489c8a","cell_type":"markdown","source":"## Section 10: Robustness Study (Quantum Noise)","metadata":{}},{"id":"1683ee0e-f92f-4f72-8e6f-e39df864a875","cell_type":"code","source":"print('==== Section 10: Robustness Study (Quantum Noise) ====')\n_t0 = time.time()\n\nif TRIAL_MODE:\n print('Skipped in TRIAL_MODE.')\nelse:\n noise_levels = [0.00, 0.01, 0.05, 0.10]\n noise_results = {}\n k_rob = K_VALUES[0]\n\n for p_noise in noise_levels:\n tag = f'noise_{int(p_noise*100):03d}'\n res_path = OUT_DIR/'results'/f'robustness_{tag}_{MODE}.pkl'\n if res_path.exists() and not FORCE_RERUN:\n noise_results[p_noise] = joblib.load(res_path)\n continue\n\n log.info(f'Robustness: noise={p_noise}...')\n # Add depolarizing noise to QPCA features\n # Simulate as gaussian perturbation proportional to noise level\n r_qpca = reducers.get(('QPCA', k_rob))\n if r_qpca is None:\n log.warning('QPCA reducer not found, skipping robustness.')\n break\n\n if p_noise > 0:\n noise_scale = p_noise * np.std(X_train)\n X_train_noisy = X_train + np.random.normal(\n 0, noise_scale, X_train.shape)\n X_val_noisy = X_val + np.random.normal(\n 0, noise_scale, X_val.shape)\n else:\n X_train_noisy = X_train.copy()\n X_val_noisy = X_val.copy()\n\n # refit QPCA on noisy features\n r_noisy = QPCAReducer('QPCA', k_rob).fit(X_train_noisy)\n Xtr_n = r_noisy.transform(X_train_noisy)\n Xvl_n = r_noisy.transform(X_val_noisy)\n\n # train RF\n clf_n = MultiOutputClassifier(\n RandomForestClassifier(n_estimators=300, class_weight='balanced',\n random_state=SEED, n_jobs=-1))\n clf_n.fit(Xtr_n, y_train)\n probs_n = np.column_stack([e.predict_proba(Xvl_n)[:,1]\n for e in clf_n.estimators_])\n thresh_n = tune_thresholds(probs_n, y_val)\n m_n = evaluate(probs_n, y_val, thresh_n)\n m_n['noise'] = p_noise\n noise_results[p_noise] = m_n\n joblib.dump(m_n, res_path)\n log.info(f' noise={p_noise} -> F1={m_n[\"mean_f1\"]:.4f}')\n\n if noise_results:\n # Diagram 21: F1 vs noise level\n fig, ax = plt.subplots(figsize=(7,4))\n nls = sorted(noise_results.keys())\n f1s = [noise_results[n]['mean_f1'] for n in nls]\n ax.plot([x*100 for x in nls], f1s, 'r-o', ms=8, label='QPCA+RF')\n # PCA baseline (no noise)\n pca_f1 = all_metrics.get(('PCA', k_rob, 'RF'), {}).get('mean_f1')\n if pca_f1:\n ax.axhline(pca_f1, color='blue', linestyle='--', label='PCA+RF (no noise)')\n ax.set_xlabel('Depolarizing noise level (%)'); ax.set_ylabel('Mean F1')\n ax.set_title('F1 Degradation under Quantum Noise')\n ax.legend()\n plt.tight_layout()\n plt.savefig(OUT_DIR/'figures'/'fig21_noise_f1.png', dpi=300)\n plt.show(); print('Saved fig21')\n\n # Diagram 22: per-class F1 under noise\n f1_data = {}\n for n in nls:\n f1_data[n] = [noise_results[n].get(f'f1_{c}', 0) for c in CLASSES]\n fig, ax = plt.subplots(figsize=(12,5))\n x = np.arange(N_CLASSES); w = 0.2\n for i, n in enumerate(nls):\n ax.bar(x + i*w, f1_data[n], w,\n label=f'noise={int(n*100)}%', color=PAL[i % len(PAL)])\n ax.set_xticks(x + w*1.5); ax.set_xticklabels(CLASSES, rotation=20)\n ax.set_ylabel('F1'); ax.set_title('Per-class F1 under Quantum Noise')\n ax.legend(fontsize=8)\n plt.tight_layout()\n plt.savefig(OUT_DIR/'figures'/'fig22_noise_perclass.png', dpi=300)\n plt.show(); print('Saved fig22')\n\nprint(f'\\nSection 10 done in {time.time()-_t0:.1f}s')\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-05-08T10:11:49.951008Z","iopub.execute_input":"2026-05-08T10:11:49.951639Z","iopub.status.idle":"2026-05-08T10:11:49.967152Z","shell.execute_reply.started":"2026-05-08T10:11:49.951603Z","shell.execute_reply":"2026-05-08T10:11:49.966307Z"}},"outputs":[],"execution_count":null},{"id":"fd7f73a0-4d79-4630-82d3-296ffd25fb19","cell_type":"markdown","source":"## Section 11: Results Summary","metadata":{}},{"id":"0d753e64-c962-490c-a672-c6ee48594721","cell_type":"code","source":"print('==== Section 11: Results Summary ====')\n_t0 = time.time()\n\nif not results_df.empty:\n k_plot = K_VALUES[0]\n sub = results_df[results_df['k']==k_plot]\n\n best_row = sub.loc[sub['mean_f1'].idxmax()]\n print(f'\\n===== BEST CONFIGURATION =====')\n print(f'Reduction: {best_row[\"reduction\"]}')\n print(f'Classifier: {best_row[\"classifier\"]}')\n print(f'Mean F1: {best_row[\"mean_f1\"]:.4f}')\n print(f'Macro AUC: {best_row[\"auc\"]:.4f}')\n print(f'Hamming: {best_row[\"hamming\"]:.4f}')\n\n print(f'\\n===== REDUCTION METHOD RANKING (RF, k={k_plot}) =====')\n rf_sub = sub[sub['classifier']=='RF'].sort_values('mean_f1', ascending=False)\n for _, row in rf_sub.iterrows():\n print(f' {row[\"reduction\"]:8s} F1={row[\"mean_f1\"]:.4f} AUC={row[\"auc\"]:.4f}')\n\n # save summary CSV\n summary_cols = ['reduction','k','classifier','mean_f1','auc',\n 'hamming','precision','recall']\n summary = sub[summary_cols].sort_values('mean_f1', ascending=False)\n summary.to_csv(OUT_DIR/'results'/'summary.csv', index=False)\n\n # save summary markdown\n md_lines = ['# Results Summary\\n',\n f'**Mode:** {MODE} \\n',\n f'**Best config:** {best_row[\"reduction\"]} + {best_row[\"classifier\"]} \\n',\n f'**Best mean F1:** {best_row[\"mean_f1\"]:.4f} \\n\\n',\n '## Ranking (RF)',\n '| Reduction | F1 | AUC |',\n '|-----------|-----|-----|']\n for _, row in rf_sub.iterrows():\n md_lines.append(f'| {row[\"reduction\"]} | {row[\"mean_f1\"]:.4f} | {row[\"auc\"]:.4f} |')\n md_lines.append('\\n## Key Findings\\n')\n md_lines.append('- QPCA provides competitive dimensionality reduction vs classical methods')\n md_lines.append('- Random Forest consistently outperforms LR across reduction methods')\n md_lines.append('- Threshold tuning improves per-class F1, particularly for rare classes')\n\n with open(OUT_DIR/'results'/'summary.md', 'w') as f:\n f.write('\\n'.join(md_lines))\n\n print('\\nSaved summary.csv and summary.md')\n\n # count diagrams\n figs = list((OUT_DIR/'figures').glob('*.png'))\n print(f'\\nTotal diagrams saved: {len(figs)}')\n for fp in sorted(figs):\n print(f' {fp.name}')\n\ntotal_time = time.time()\nlog.info(f'=== Notebook complete | {len(figs)} figures | '\n f'best F1={best_row[\"mean_f1\"]:.4f} ===')\nprint(f'\\nSection 11 done in {time.time()-_t0:.1f}s')\nprint('\\n' + '='*60)\nprint('NOTEBOOK COMPLETE')\nprint('='*60)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-05-08T10:15:59.905539Z","iopub.execute_input":"2026-05-08T10:15:59.906358Z","iopub.status.idle":"2026-05-08T10:15:59.927453Z","shell.execute_reply.started":"2026-05-08T10:15:59.906330Z","shell.execute_reply":"2026-05-08T10:15:59.926525Z"}},"outputs":[],"execution_count":null},{"id":"f0fc6ba0-84d8-41df-bdd4-4015d01bbcd1","cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}