|
| 1 | +import os |
| 2 | +import math |
| 3 | +import logging |
| 4 | +from pathlib import Path |
| 5 | + |
| 6 | +import anndata as ad |
| 7 | +import numpy as np |
| 8 | + |
| 9 | +import torch |
| 10 | +import pytorch_lightning as pl |
| 11 | +from torch.utils.data import TensorDataset, DataLoader |
| 12 | +from pytorch_lightning.callbacks import ModelCheckpoint |
| 13 | +from pytorch_lightning.loggers import TensorBoardLogger,WandbLogger |
| 14 | + |
| 15 | +logging.basicConfig(level=logging.INFO) |
| 16 | + |
| 17 | +## VIASH START |
| 18 | +par = { |
| 19 | + 'input_train_mod1': 'resources_test/predict_modality/openproblems_neurips2021/bmmc_multiome/swap/train_mod1.h5ad', |
| 20 | + 'input_train_mod2': 'resources_test/predict_modality/openproblems_neurips2021/bmmc_multiome/swap/train_mod2.h5ad', |
| 21 | + 'input_test_mod1': 'resources_test/predict_modality/openproblems_neurips2021/bmmc_multiome/swap/test_mod1.h5ad', |
| 22 | + 'output': 'output/model' |
| 23 | +} |
| 24 | +meta = { |
| 25 | + 'resources_dir': 'src/tasks/predict_modality/methods/simple_mlp', |
| 26 | + 'cpus': 10 |
| 27 | +} |
| 28 | +## VIASH END |
| 29 | + |
| 30 | +resources_dir = f"{meta['resources_dir']}/resources" |
| 31 | + |
| 32 | +import sys |
| 33 | +sys.path.append(resources_dir) |
| 34 | +from models import MLP |
| 35 | +import utils |
| 36 | + |
| 37 | +def _train(X, y, Xt, yt, logger, config, num_workers): |
| 38 | + |
| 39 | + X = torch.from_numpy(X).float() |
| 40 | + y = torch.from_numpy(y).float() |
| 41 | + ymean = torch.mean(y, dim=0, keepdim=True) |
| 42 | + |
| 43 | + tr_ds = TensorDataset(X,y) |
| 44 | + tr_loader = DataLoader( |
| 45 | + tr_ds, |
| 46 | + batch_size=config.batch_size, |
| 47 | + num_workers=num_workers, |
| 48 | + shuffle=True, |
| 49 | + drop_last=True |
| 50 | + ) |
| 51 | + |
| 52 | + Xt = torch.from_numpy(Xt).float() |
| 53 | + yt = torch.from_numpy(yt).float() |
| 54 | + te_ds = TensorDataset(Xt,yt) |
| 55 | + te_loader = DataLoader( |
| 56 | + te_ds, |
| 57 | + batch_size=config.batch_size, |
| 58 | + num_workers=num_workers, |
| 59 | + shuffle=False, |
| 60 | + drop_last=False |
| 61 | + ) |
| 62 | + |
| 63 | + checkpoint_callback = ModelCheckpoint( |
| 64 | + monitor='valid_RMSE', |
| 65 | + dirpath=logger.save_dir, |
| 66 | + save_top_k=1, |
| 67 | + ) |
| 68 | + |
| 69 | + trainer = pl.Trainer( |
| 70 | + devices="auto", |
| 71 | + enable_checkpointing=True, |
| 72 | + logger=logger, |
| 73 | + max_epochs=config.epochs, |
| 74 | + callbacks=[checkpoint_callback], |
| 75 | + default_root_dir=logger.save_dir, |
| 76 | + # progress_bar_refresh_rate=5 |
| 77 | + ) |
| 78 | + |
| 79 | + net = MLP(X.shape[1], y.shape[1], ymean, config) |
| 80 | + trainer.fit(net, tr_loader, te_loader) |
| 81 | + |
| 82 | + yp = trainer.predict(net, te_loader, ckpt_path='best') |
| 83 | + yp = torch.cat(yp, dim=0) |
| 84 | + |
| 85 | + score = ((yp-yt)**2).mean()**0.5 |
| 86 | + print(f"VALID RMSE {score:.3f}") |
| 87 | + del trainer |
| 88 | + return score,yp.detach().numpy() |
| 89 | + |
| 90 | + |
| 91 | + |
| 92 | +input_train_mod1 = ad.read_h5ad(par['input_train_mod1']) |
| 93 | +input_train_mod2 = ad.read_h5ad(par['input_train_mod2']) |
| 94 | + |
| 95 | +mod_1 = input_train_mod1.uns["modality"] |
| 96 | +mod_2 = input_train_mod2.uns["modality"] |
| 97 | + |
| 98 | +task = f'{mod_1}2{mod_2}' |
| 99 | +yaml_path = f'{resources_dir}/yaml/mlp_{task}.yaml' |
| 100 | + |
| 101 | +obs_info = utils.to_site_donor(input_train_mod1) |
| 102 | +# TODO: if we want this method to work for other datasets, resolve dependence on site notation |
| 103 | +sites = obs_info.site.unique() |
| 104 | + |
| 105 | +os.makedirs(par['output'], exist_ok=True) |
| 106 | + |
| 107 | +print('Compute ymean', flush=True) |
| 108 | +ymean = np.asarray(input_train_mod2.layers["normalized"].mean(axis=0)) |
| 109 | +path = f"{par['output']}/{task}_ymean.npy" |
| 110 | +np.save(path, ymean) |
| 111 | + |
| 112 | + |
| 113 | +if task == "GEX2ATAC": |
| 114 | + logging.info(f"No training required for this task ({task}).") |
| 115 | + sys.exit(0) |
| 116 | + |
| 117 | +if not os.path.exists(yaml_path): |
| 118 | + logging.error(f"No configuration file found for task '{task}'") |
| 119 | + sys.exit(1) |
| 120 | + |
| 121 | +yaml_path = f'{resources_dir}/yaml/mlp_{task}.yaml' |
| 122 | +yps = [] |
| 123 | +scores = [] |
| 124 | + |
| 125 | +msgs = {} |
| 126 | +# TODO: if we want this method to work for other datasets, dont use hardcoded range |
| 127 | +for fold in range(3): |
| 128 | + |
| 129 | + run_name = f"{task}_fold_{fold}" |
| 130 | + save_path = f"{par['output']}/{run_name}" |
| 131 | + num_workers = meta["cpus"] or 0 |
| 132 | + |
| 133 | + Path(save_path).mkdir(parents=True, exist_ok=True) |
| 134 | + |
| 135 | + X,y,Xt,yt = utils.split(input_train_mod1, input_train_mod2, fold) |
| 136 | + |
| 137 | + logger = TensorBoardLogger(save_path, name='') |
| 138 | + |
| 139 | + config = utils.load_yaml(yaml_path) |
| 140 | + |
| 141 | + if config.batch_size > X.shape[0]: |
| 142 | + config = config._replace(batch_size=math.ceil(X.shape[0] / 2)) |
| 143 | + |
| 144 | + score, yp = _train(X, y, Xt, yt, logger, config, num_workers) |
| 145 | + yps.append(yp) |
| 146 | + scores.append(score) |
| 147 | + msg = f"{task} Fold {fold} RMSE {score:.3f}" |
| 148 | + msgs[f'Fold {fold}'] = f'{score:.3f}' |
| 149 | + print(msg) |
| 150 | + |
| 151 | +yp = np.concatenate(yps) |
| 152 | +score = np.mean(scores) |
| 153 | +msgs['Overall'] = f'{score:.3f}' |
| 154 | +print('Overall', f'{score:.3f}') |
0 commit comments