diff --git a/CHANGELOG.md b/CHANGELOG.md index 10ec97c6..2c3c64c0 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,9 @@ # task_predict_modality 0.2.0 +## NEW FUNCTIONALITY + +* Added `src/utils/exit_codes.py`, so components can mark themselves non-applicable for a dataset (PR #32). + ## MINOR CHANGES * `file_train_mod1`, `file_train_mod2`, `file_test_mod1`, `file_test_mod2`: Declare `uns["modality"]`, which `process_dataset` already writes and which six methods and the `run_benchmark` workflow already read (PR #29). @@ -34,6 +38,10 @@ ## BUG FIXES +* `guanlab_dengkw_pm`: Restore the consensus scheme of the original submission -- five reshuffles of the batches into two halves, ten kernel ridge models averaged. The port had replaced it with a single fixed two-way split for ADT pairs and leave-one-batch-out otherwise, so the result depended on the order the batches happened to come in. `--n_repeats` and `--seed` are now arguments; the unused `--distance_method` and `--n_pcs` are gone (PR #32). + +* `guanlab_dengkw_pm`: Only map the predictions back through the mod2 SVD when that SVD was actually fitted. When the target had fewer features than `n_mod2`, the method raised `NameError: embedder_mod2`. Unsupported modality pairs now exit 99 (non-applicable) instead of raising a `KeyError` (PR #32). + * `process_dataset`: Fall back to holding out a quarter of the batches when the dataset has no `obs["is_train"]`, rather than silently producing four empty h5ads. `obs["is_train"]` carries the NeurIPS 2021 competition split and stays optional; `obs["cell_type"]` is now declared and required (PR #28). * Fix the component paths, build paths and `rename_keys` separator in the helper scripts, which prevented `scripts/create_datasets/test_resources.sh` and both `run_test.sh` scripts from running at all (PR #22). diff --git a/src/methods/guanlab_dengkw_pm/config.vsh.yaml b/src/methods/guanlab_dengkw_pm/config.vsh.yaml index 925b5591..fb4b9c20 100644 --- a/src/methods/guanlab_dengkw_pm/config.vsh.yaml +++ b/src/methods/guanlab_dengkw_pm/config.vsh.yaml @@ -18,18 +18,22 @@ info: preferred_normalization: log_cp10k competition_submission_id: 170636 arguments: - - name: "--distance_method" - type: "string" - default: "minkowski" - description: The distance metric to use. Possible values include `euclidean` and `minkowski`. - choices: [euclidean, minkowski] - - name: "--n_pcs" + - name: "--n_repeats" type: "integer" - default: 50 - description: Number of components to use for dimensionality reduction. + default: 5 + description: | + Number of times the batches are reshuffled into two halves. Each half fits one + kernel ridge model, so the consensus averages `2 * n_repeats` predictions. + info: + test_default: 1 + - name: "--seed" + type: "integer" + default: 1000 + description: Seed for the batch reshuffling. resources: - type: python_script path: script.py + - path: /src/utils/exit_codes.py engines: - type: docker image: openproblems/base_python:1 diff --git a/src/methods/guanlab_dengkw_pm/script.py b/src/methods/guanlab_dengkw_pm/script.py index d11b7355..3a0fd029 100644 --- a/src/methods/guanlab_dengkw_pm/script.py +++ b/src/methods/guanlab_dengkw_pm/script.py @@ -1,3 +1,4 @@ +import sys import anndata as ad import numpy as np from scipy.sparse import csc_matrix @@ -10,15 +11,19 @@ 'input_train_mod1': 'resources_test/task_predict_modality/openproblems_neurips2021/bmmc_multiome/normal/train_mod1.h5ad', 'input_train_mod2': 'resources_test/task_predict_modality/openproblems_neurips2021/bmmc_multiome/normal/train_mod2.h5ad', 'input_test_mod1': 'resources_test/task_predict_modality/openproblems_neurips2021/bmmc_multiome/normal/test_mod1.h5ad', - 'output': 'output.h5ad', - 'distance_method': 'minkowski', - 'n_pcs': 50 + 'output': 'output.h5ad', + 'n_repeats': 5, + 'seed': 1000 } meta = { - 'name': 'guanlab_dengkw_pm' + 'name': 'guanlab_dengkw_pm', + 'resources_dir': 'src/utils' } ## VIASH END +sys.path.append(meta['resources_dir']) +from exit_codes import exit_non_applicable + ## Removed PCA and normalization steps, as they arr already performed with the input data print('Reading input files', flush=True) @@ -49,6 +54,11 @@ ("ATAC", "GEX"): (100, 70, 10, 0.1) } print(f"{mod1_type}, {mod2_type}", flush=True) +if (mod1_type, mod2_type) not in n_comp_dict: + exit_non_applicable( + f"Guanlab-dengkw was tuned on the NeurIPS 2021 modality pairs and has no " + f"hyperparameters for {mod1_type} -> {mod2_type}" + ) n_mod1, n_mod2, scale, alpha = n_comp_dict[(mod1_type, mod2_type)] print(f"{n_mod1}, {n_mod2}, {scale}, {alpha}", flush=True) @@ -68,6 +78,8 @@ embedder_mod2 = TruncatedSVD(n_components=n_mod2) train_gs = embedder_mod2.fit_transform(input_train_mod2.layers["normalized"]).astype(np.float32) else: + # no dimensionality reduction, so the KRR predicts the mod2 features directly + embedder_mod2 = None train_gs = input_train_mod2.to_df(layer="normalized").values.astype(np.float32) del input_train @@ -86,36 +98,36 @@ del test_matrix print('Running KRR model ...', flush=True) -if batch_len == 1: - # just in case there is only one batch - batch_subsets = [batches] -elif mod1_type == "ADT" or mod2_type == "ADT": - # two fold consensus predictions - batch_subsets = [ - batches[:batch_len//2], - batches[batch_len//2:] - ] -else: - # leave-one-batch-out consensus predictions - batch_subsets = [ - batches[:i] + batches[i+1:] - for i in range(batch_len) - ] - +# consensus over n_repeats reshuffles of the batches into two halves, as in the +# original submission y_pred = np.zeros((input_test_mod1.n_obs, input_train_mod2.n_vars), dtype=np.float32) -for batch in batch_subsets: - print(batch, flush=True) - kernel = RBF(length_scale = scale) - krr = KernelRidge(alpha=alpha, kernel=kernel) - print('Fitting KRR ... ', flush=True) - krr.fit( - train_norm[input_train_mod1.obs.batch.isin(batch)], - train_gs[input_train_mod2.obs.batch.isin(batch)] - ) - y_pred += (krr.predict(test_norm) @ embedder_mod2.components_) +np.random.seed(par['seed']) +n_models = 0 + +for _ in range(par['n_repeats']): + np.random.shuffle(batches) + for batch in [batches[:batch_len//2], batches[batch_len//2:]]: + # with a single batch one half comes out empty + if not batch: + batch = [batches[0]] + + print(batch, flush=True) + kernel = RBF(length_scale = scale) + krr = KernelRidge(alpha=alpha, kernel=kernel) + print('Fitting KRR ... ', flush=True) + krr.fit( + train_norm[input_train_mod1.obs.batch.isin(batch)], + train_gs[input_train_mod2.obs.batch.isin(batch)] + ) + y_pred_batch = krr.predict(test_norm) + if embedder_mod2 is not None: + # map the predicted components back to the mod2 feature space + y_pred_batch = y_pred_batch @ embedder_mod2.components_ + y_pred += y_pred_batch + n_models += 1 np.clip(y_pred, a_min=0, a_max=None, out=y_pred) -y_pred /= len(batch_subsets) +y_pred /= n_models # Store as sparse matrix to be efficient. # Note that this might require different classifiers/embedders before-hand. diff --git a/src/utils/exit_codes.py b/src/utils/exit_codes.py new file mode 100644 index 00000000..03f270c3 --- /dev/null +++ b/src/utils/exit_codes.py @@ -0,0 +1,6 @@ +import sys + + +def exit_non_applicable(msg): + print(f"NON-APPLICABLE ERROR: {msg}", flush=True) + sys.exit(99)