From 0bc7003cdb8a162ff5e08665942ac9ac0f9fa3db Mon Sep 17 00:00:00 2001 From: asistradition Date: Mon, 26 Apr 2021 15:05:48 -0400 Subject: [PATCH 01/14] Add a quick & dirty memory profiler --- inferelator/utils/profiler.py | 35 +++++++++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) create mode 100644 inferelator/utils/profiler.py diff --git a/inferelator/utils/profiler.py b/inferelator/utils/profiler.py new file mode 100644 index 00000000..db892a56 --- /dev/null +++ b/inferelator/utils/profiler.py @@ -0,0 +1,35 @@ +import os, psutil, argparse, time, csv + +def main(): + """ + Module that produces a profile of process memory usage and saves it to a file + """ + + ap = argparse.ArgumentParser(description="Create a prior from a genome, TF motifs, and an optional BED file") + + # REQUIRED ARGUMENTS ############################################################################################### + + ap.add_argument("-p", dest="pid", help="Process ID", type=int, required=True) + ap.add_argument("-o", dest="out", help="Output TSV", metavar="FILE", required=True) + ap.add_argument("-t", dest="timer", help="Time interval", type=int, default=1) + args = ap.parse_args() + + pid, t_step = args.pid, args.timer + + time_start = time.time() + with open(args.out, mode="w") as out_fh: + csv_handler = csv.writer(out_fh, delimiter="\t", lineterminator="\n", quoting=csv.QUOTE_NONE) + csv_handler.writerow(["Time", "Memory"]) + + while psutil.pid_exists(pid): + csv_handler.writerow(['%.3f' % (time.time() - time_start), int(get_current_memory(pid) / 1e6)]) + time.sleep(t_step) + + +def get_current_memory(proc_id): + proc = psutil.Process(proc_id) + return proc.memory_info().rss + + +if __name__ == '__main__': + main() From 415fe1cab889cbc21835456eeb3363f205603555 Mon Sep 17 00:00:00 2001 From: asistradition Date: Mon, 26 Apr 2021 18:08:44 -0400 Subject: [PATCH 02/14] Add flag to skip grid search without sampling --- inferelator/crossvalidation_workflow.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/inferelator/crossvalidation_workflow.py b/inferelator/crossvalidation_workflow.py index 154b0480..c5d90fc8 100644 --- a/inferelator/crossvalidation_workflow.py +++ b/inferelator/crossvalidation_workflow.py @@ -58,6 +58,7 @@ class CrossValidationManager(object): size_sample_stratified_column = None size_sample_with_replacement = False size_sample_seed = None + size_sample_only = False # Workflow storage _baseline_workflow = None @@ -180,7 +181,8 @@ def add_grouping_dropin(self, metadata_column_name, group_size=None, seed=42): self.dropin_max_size = group_size self.dropin_seed = seed - def add_size_subsampling(self, size_vector, stratified_column_name=None, with_replacement=False, seed=42): + def add_size_subsampling(self, size_vector, stratified_column_name=None, with_replacement=False, seed=42, + size_sample_only=None): """ Resample expression data to a ratio of the original data. @@ -206,6 +208,7 @@ def add_size_subsampling(self, size_vector, stratified_column_name=None, with_re self.size_sample_stratified_column = stratified_column_name self.size_sample_with_replacement = with_replacement self.size_sample_seed = seed + self.size_sample_only = size_sample_only if size_sample_only is not None else self.size_sample_only def run(self): """ @@ -222,7 +225,10 @@ def run(self): self._check_grid_search_params_exist() # Run base grid search - results = self._grid_search() + if self.size_sample_only: + results = [] + else: + results = self._grid_search() # Run size sampling if self.size_sample_vector is not None: From d1f82e2027daf7592a1751694fae568bc06d83ec Mon Sep 17 00:00:00 2001 From: asistradition Date: Mon, 26 Apr 2021 18:08:52 -0400 Subject: [PATCH 03/14] Add flush() --- inferelator/utils/profiler.py | 1 + 1 file changed, 1 insertion(+) diff --git a/inferelator/utils/profiler.py b/inferelator/utils/profiler.py index db892a56..b2d287a5 100644 --- a/inferelator/utils/profiler.py +++ b/inferelator/utils/profiler.py @@ -23,6 +23,7 @@ def main(): while psutil.pid_exists(pid): csv_handler.writerow(['%.3f' % (time.time() - time_start), int(get_current_memory(pid) / 1e6)]) + out_fh.flush() time.sleep(t_step) From 2cf642d9533f55f2eeb12e07632c02701f005c7b Mon Sep 17 00:00:00 2001 From: asistradition Date: Mon, 26 Apr 2021 18:29:48 -0400 Subject: [PATCH 04/14] Fix CV for new data object --- inferelator/crossvalidation_workflow.py | 8 ++--- .../tests/test_crossvalidation_wrapper.py | 36 +++++++++++-------- 2 files changed, 26 insertions(+), 18 deletions(-) diff --git a/inferelator/crossvalidation_workflow.py b/inferelator/crossvalidation_workflow.py index c5d90fc8..560a9237 100644 --- a/inferelator/crossvalidation_workflow.py +++ b/inferelator/crossvalidation_workflow.py @@ -454,7 +454,7 @@ def _check_metadata_column_exists(self, col_name): :param col_name: str """ - if col_name in self.workflow.meta_data.columns: + if col_name in self.workflow.data.meta_data.columns: return True else: raise ValueError("Column {col} is not present in the loaded metadata".format(col=col_name)) @@ -464,7 +464,7 @@ def _dropout_cv(self): Run grid search on all data minus one group at a time """ - meta_data = self.workflow.meta_data.copy() + meta_data = self.workflow.data.meta_data.copy() col = self.dropout_column max_size = self.dropout_max_size @@ -517,7 +517,7 @@ def _dropin_cv(self): Run grid search on one group from the data at a time """ - meta_data = self.workflow.meta_data.copy() + meta_data = self.workflow.data.meta_data.copy() col = self.dropin_column max_size = self.dropin_max_size @@ -557,7 +557,7 @@ def _size_cv(self): for i, size_ratio in enumerate(self.size_sample_vector): rgen = np.random.RandomState(self.size_sample_seed + i) - meta_data = self.workflow.meta_data.copy() + meta_data = self.workflow.data.meta_data.copy() if self.size_sample_stratified_column is not None: strat_col = self.size_sample_stratified_column diff --git a/inferelator/tests/test_crossvalidation_wrapper.py b/inferelator/tests/test_crossvalidation_wrapper.py index b64faef4..27d2a5b5 100644 --- a/inferelator/tests/test_crossvalidation_wrapper.py +++ b/inferelator/tests/test_crossvalidation_wrapper.py @@ -7,9 +7,17 @@ from inferelator import crossvalidation_workflow from inferelator.workflow import WorkflowBase +from inferelator.utils.data import InferelatorData + +from numpy.random import default_rng + +fake_obsnames = list(map(str, range(1000))) fake_metadata = pd.DataFrame({"CONST": ["A"] * 1000, - "VAR": ["A"] * 100 + ["B"] * 200 + ["C"] * 1 + ["D"] * 99 + ["E"] * 500 + ["F"] * 100}) + "VAR": ["A"] * 100 + ["B"] * 200 + ["C"] * 1 + ["D"] * 99 + ["E"] * 500 + ["F"] * 100}, + index = fake_obsnames) + +fake_data_object = InferelatorData(default_rng(12345).random(size=1000).reshape((1000,1)), meta_data=fake_metadata, sample_names=fake_obsnames) TEMP_DIR = tempfile.gettempdir() TEMP_DIR_1 = os.path.join(TEMP_DIR, "test1") @@ -29,7 +37,7 @@ class FakeWorkflow(WorkflowBase): seed = 10 def __init__(self): - self.meta_data = fake_metadata.copy() + self.data = fake_data_object.copy() def run(self): return FakeResult() @@ -219,8 +227,8 @@ def test_group_dropout_no_limit(self): def test_grid_search(slf, test=None, value=None, mask_function=None): self.assertEqual(test, "dropout") - self.assertTrue(value in slf.workflow.meta_data[slf.dropout_column].unique()) - self.assertListEqual((slf.workflow.meta_data[slf.dropout_column] != value).tolist(), + self.assertTrue(value in slf.workflow.data.meta_data[slf.dropout_column].unique()) + self.assertListEqual((slf.workflow.data.meta_data[slf.dropout_column] != value).tolist(), mask_function().tolist()) self.cv._grid_search = types.MethodType(test_grid_search, self.cv) @@ -235,10 +243,10 @@ def test_group_dropout_limit(self): def test_grid_search(slf, test=None, value=None, mask_function=None): self.assertEqual(test, "dropout") - uniques = slf.workflow.meta_data[slf.dropout_column].unique() + uniques = slf.workflow.data.meta_data[slf.dropout_column].unique() mask = mask_function() - unique_counts = slf.workflow.meta_data[slf.dropout_column].value_counts() + unique_counts = slf.workflow.data.meta_data[slf.dropout_column].value_counts() unique_counts[unique_counts > slf.dropout_max_size] = slf.dropout_max_size if value == "all": @@ -247,7 +255,7 @@ def test_grid_search(slf, test=None, value=None, mask_function=None): self.assertTrue(value in uniques) unique_counts[value] = 0 self.assertEqual(unique_counts.sum(), mask.sum()) - self.assertEqual(sum((self.cv.workflow.meta_data[self.cv.dropout_column] == value)[mask]), 0) + self.assertEqual(sum((self.cv.workflow.data.meta_data[self.cv.dropout_column] == value)[mask]), 0) self.cv._grid_search = types.MethodType(test_grid_search, self.cv) @@ -261,8 +269,8 @@ def test_group_dropin_no_limit(self): def test_grid_search(slf, test=None, value=None, mask_function=None): self.assertEqual(test, "dropin") - self.assertTrue(value in slf.workflow.meta_data[slf.dropin_column].unique()) - self.assertListEqual((slf.workflow.meta_data[slf.dropin_column] == value).tolist(), + self.assertTrue(value in slf.workflow.data.meta_data[slf.dropin_column].unique()) + self.assertListEqual((slf.workflow.data.meta_data[slf.dropin_column] == value).tolist(), mask_function().tolist()) self.cv._grid_search = types.MethodType(test_grid_search, self.cv) @@ -283,9 +291,9 @@ def test_grid_search(slf, test=None, value=None, mask_function=None): if value == "all": self.assertEqual(mask.sum(), slf.dropin_max_size) else: - self.assertTrue(value in slf.workflow.meta_data[slf.dropin_column].unique()) + self.assertTrue(value in slf.workflow.data.meta_data[slf.dropin_column].unique()) - self.assertEqual(min((slf.workflow.meta_data[slf.dropin_column] == value).sum(), + self.assertEqual(min((slf.workflow.data.meta_data[slf.dropin_column] == value).sum(), slf.dropin_max_size), mask.sum()) @@ -303,7 +311,7 @@ def test_grid_search(slf, test=None, value=None, mask_function=None): self.assertEqual(test, "size") self.assertTrue(value == "0.5") - self.assertEqual(max(int(slf.workflow.meta_data.shape[0] * float(value)), 1), + self.assertEqual(max(int(slf.workflow.data.meta_data.shape[0] * float(value)), 1), mask_function().sum()) self.cv._grid_search = types.MethodType(test_grid_search, self.cv) @@ -320,8 +328,8 @@ def test_grid_search(slf, test=None, value=None, mask_function=None): self.assertTrue(value == "0.5") mask = mask_function() - for g in slf.workflow.meta_data[slf.size_sample_stratified_column].unique(): - is_group = slf.workflow.meta_data[slf.size_sample_stratified_column] == g + for g in slf.workflow.data.meta_data[slf.size_sample_stratified_column].unique(): + is_group = slf.workflow.data.meta_data[slf.size_sample_stratified_column] == g self.assertEqual(max(int(is_group.sum() * float(value)), 1), mask[is_group].sum()) From 545d7341b1218c4dcf54d93a5e716d5b95cf0181 Mon Sep 17 00:00:00 2001 From: asistradition Date: Mon, 26 Apr 2021 19:24:37 -0400 Subject: [PATCH 05/14] Fix sampling data object --- inferelator/crossvalidation_workflow.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/inferelator/crossvalidation_workflow.py b/inferelator/crossvalidation_workflow.py index 560a9237..67db53ae 100644 --- a/inferelator/crossvalidation_workflow.py +++ b/inferelator/crossvalidation_workflow.py @@ -388,8 +388,7 @@ def _grid_search(self, test=None, value=None, mask_function=None): # Drop any observations which are False in the mask (if set) if mask_function is not None: mask = mask_function() - cv_workflow.expression_matrix.drop(cv_workflow.expression_matrix.columns[~mask], axis=1, inplace=True) - cv_workflow.meta_data.drop(cv_workflow.meta_data.index[~mask], axis=0, inplace=True) + cv_workflow.data = cv_workflow.data.get_sample_data(mask.index[mask]) n_obs = mask.sum() else: n_obs = cv_workflow._num_obs From 925f6661d560233e6ecbaa37b6ed721972b31fa7 Mon Sep 17 00:00:00 2001 From: asistradition Date: Tue, 27 Apr 2021 12:41:16 -0400 Subject: [PATCH 06/14] Add sampling to the data wrapper --- inferelator/tests/artifacts/test_data.py | 5 ++- inferelator/tests/test_data_wrapper.py | 55 +++++++++++++++++++++++ inferelator/utils/data.py | 57 +++++++++++++++++++++++- 3 files changed, 113 insertions(+), 4 deletions(-) diff --git a/inferelator/tests/artifacts/test_data.py b/inferelator/tests/artifacts/test_data.py index 9dac697b..b63a5e92 100644 --- a/inferelator/tests/artifacts/test_data.py +++ b/inferelator/tests/artifacts/test_data.py @@ -22,11 +22,12 @@ class TestDataSingleCellLike(object): transpose_expression=True, meta_data=TestDataSingleCellLike.meta_data, gene_data=TestDataSingleCellLike.gene_metadata, - gene_data_idx_column=TestDataSingleCellLike.gene_list_index) + gene_data_idx_column=TestDataSingleCellLike.gene_list_index, + sample_names=list(map(str, range(10)))) TEST_DATA_SPARSE = InferelatorData(sps.csr_matrix(TestDataSingleCellLike.expression_matrix.T.values), gene_names=TestDataSingleCellLike.expression_matrix.index, - sample_names=TestDataSingleCellLike.expression_matrix.columns, + sample_names=list(map(str, range(10))), meta_data=TestDataSingleCellLike.meta_data, gene_data=TestDataSingleCellLike.gene_metadata, gene_data_idx_column=TestDataSingleCellLike.gene_list_index) diff --git a/inferelator/tests/test_data_wrapper.py b/inferelator/tests/test_data_wrapper.py index 6da7cc40..e076f932 100644 --- a/inferelator/tests/test_data_wrapper.py +++ b/inferelator/tests/test_data_wrapper.py @@ -361,5 +361,60 @@ def test_change_sparse(self): self.assertTrue(sparse.isspmatrix_csr(self.adata_sparse.expression_data)) +class TestSampling(TestWrapperSetup): + + def setUp(self): + super(TestSampling, self).setUp() + + self.adata.trim_genes() + self.adata_sparse.trim_genes() + + def test_without_replacement(self): + + new_adata = self.adata.get_random_samples(10, with_replacement=False) + + new_sample_names = new_adata.sample_names.tolist() + new_sample_names.sort() + + old_sample_names = self.adata.sample_names.tolist() + old_sample_names.sort() + + self.assertListEqual(new_sample_names, old_sample_names) + with self.assertRaises(AssertionError): + self.assertListEqual(new_adata.sample_names.tolist(), self.adata.sample_names.tolist()) + + with self.assertRaises(ValueError): + self.adata.get_random_samples(100, with_replacement=False) + + with self.assertRaises(ValueError): + self.adata.get_random_samples(0, with_replacement=False) + + self.assertEqual(self.adata.get_random_samples(2, with_replacement=False).num_obs, 2) + + def test_with_replacement(self): + + new_adata = self.adata.get_random_samples(11, with_replacement=True, fix_names=True) + self.assertEqual(new_adata.num_obs, 11) + + new_sample_names = new_adata.sample_names.tolist() + new_sample_names.sort() + + old_sample_names = self.adata.sample_names.tolist() + old_sample_names.sort() + + with self.assertRaises(AssertionError): + self.assertListEqual(new_sample_names, old_sample_names) + + with self.assertRaises(ValueError): + self.adata.get_random_samples(0, with_replacement=True) + + self.assertEqual(self.adata.get_random_samples(200, with_replacement=True).num_obs, 200) + + def test_inplace(self): + + new_adata = self.adata.get_random_samples(11, with_replacement=True, fix_names=False, inplace=True) + self.assertEqual(id(new_adata), id(self.adata)) + + if __name__ == '__main__': unittest.main() diff --git a/inferelator/utils/data.py b/inferelator/utils/data.py index 756c7310..65c2f135 100644 --- a/inferelator/utils/data.py +++ b/inferelator/utils/data.py @@ -11,7 +11,8 @@ from sklearn.preprocessing import StandardScaler import scipy.io from anndata import AnnData -from inferelator.utils import Debug, Validator +from inferelator.utils import Debug +from inferelator.utils import Validator as check def dot_product(a, b, dense=True, cast=True): @@ -598,9 +599,61 @@ def get_sample_data(self, sample_index, copy=False, force_dense=False, to_df=Fal return pd.DataFrame(x, columns=self.gene_names, index=labels) if to_df else x def get_bootstrap(self, sample_bootstrap_index): - return InferelatorData(expression_data=self._adata[sample_bootstrap_index, :].X.copy(), + return InferelatorData(expression_data=self._adata.X[sample_bootstrap_index, :].copy(), gene_names=self.gene_names) + def get_random_samples(self, num_obs, with_replacement=False, random_seed=None, random_gen=None, inplace=False, + fix_names=True): + """ + Randomly sample to a specific number of observatons from the entire data set + + :param num_obs: Number of observations to return + :type num_obs: int + :param with_replacement: Sample with replacement, defaults to False + :type with_replacement: bool, optional + :param random_seed: Seed for numpy random generator, defaults to None. Will be ignored if a generator itself is + passed to random_gen. + :type random_seed: int, optional + :param random_gen: Numpy random generator to use, defaults to None. + :type random_gen: np.random.Generator, optional + :param inplace: Change this instance of the data structure inplace and return a reference to itself + :type inplace: bool, optional + """ + + check.argument_integer(num_obs, low=1) + check.argument_integer(random_seed, allow_none=True) + + if (num_obs > self.num_obs) and not with_replacement: + _msg = "Unable to sample {x} from {y} observations without replacement".format(x=num_obs, y=self.num_obs) + raise ValueError(_msg) + + # Make a new random generator if not provided + if random_gen is None: + random_gen = np.random.default_rng() if random_seed is None else np.random.default_rng(random_seed) + + # Sample with replacement using randint + if with_replacement: + keeper_ilocs = random_gen.integers(self.num_obs, size=(num_obs,)) + + # Sample without replacement using choice + else: + keeper_ilocs = random_gen.choice(np.arange(self.num_obs), size=(num_obs,), replace=False) + + # Change this instance's _adata (explicit copy allows the old data to be dereferenced instead of held as view) + if inplace: + self._adata = self._adata[keeper_ilocs, :].copy() + return_obj = self + + # Create a new InferelatorData instance with the _adata slice + else: + return_obj = InferelatorData(self._adata[keeper_ilocs, :].copy()) + + # Fix names + return_obj._adata.obs_names_make_unique() if with_replacement and fix_names else None + + return return_obj + + def subset_copy(self, row_index=None, column_index=None): if row_index is not None and column_index is not None: From 69d3ab90cc01906161818b0f8aafff72f255bc77 Mon Sep 17 00:00:00 2001 From: asistradition Date: Tue, 27 Apr 2021 17:09:07 -0400 Subject: [PATCH 07/14] Profiling Fixes For Runtime --- inferelator/regression/amusr_regression.py | 37 +++++++++++++++------- 1 file changed, 26 insertions(+), 11 deletions(-) diff --git a/inferelator/regression/amusr_regression.py b/inferelator/regression/amusr_regression.py index 36c62f75..cf676e0a 100644 --- a/inferelator/regression/amusr_regression.py +++ b/inferelator/regression/amusr_regression.py @@ -382,28 +382,43 @@ def updateB(C, D, B, S, lamB, prior, n_tasks, n_features): reference: Liu et al, ICML 2009. Blockwise coordinate descent procedures for the multi-task lasso, with applications to neural semantic basis discovery. """ - #p = prior.min(axis=1) + + S = np.asarray(S, order="F") + # cycles through predictors for j in range(n_features): + # initialize next coefficients alphas = np.zeros(n_tasks) + # update tasks for each predictor together + d = D[:, :, j] + for k in range(n_tasks): - # get task covariance update terms - c = C[k]; d = D[k] - # get previous block-sparse and sparse coefficients - b = B[:,k]; s = S[:,k] - # set block-sparse coefficient for feature j to zero - b_tmp = deepcopy(b) - b_tmp[j] = 0. - # calculate next coefficient based on fit only - if d[j,j] == 0: + + d_kjj = d[k, j] + + if d_kjj == 0: + alphas[k] = 0 + else: - alphas[k] = (c[j]-np.sum((b_tmp+s)*d[:,j]))/d[j,j] + + # get previous block-sparse + # copies because B is C-ordered + b = B[:, k] + + # set block-sparse coefficient for feature j to zero + b[j] = 0. + + # calculate next coefficient based on fit only + alphas[k] = (C[k, j] - np.sum((b + S[:, k]) * d[k, :])) / d_kjj + + # set all tasks to zero if l1-norm less than lamB if np.linalg.norm(alphas, 1) <= lamB: B[j,:] = np.zeros(n_tasks) + # regularized update for predictors with larger l1-norm else: # find number of coefficients that would make l1-norm greater than penalty From 12c391d251f02b9b835b24e60560133909401104 Mon Sep 17 00:00:00 2001 From: asistradition Date: Tue, 27 Apr 2021 17:32:45 -0400 Subject: [PATCH 08/14] Add relative tolerance escape hatch --- inferelator/regression/amusr_regression.py | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/inferelator/regression/amusr_regression.py b/inferelator/regression/amusr_regression.py index cf676e0a..b2b395da 100644 --- a/inferelator/regression/amusr_regression.py +++ b/inferelator/regression/amusr_regression.py @@ -22,6 +22,7 @@ MAX_ITER = 1000 TOL = 1e-2 +REL_TOL = 1e-4 MIN_WEIGHT_VAL = 0.1 MIN_RSS = 1e-10 @@ -260,7 +261,7 @@ def pileup_data(self, run_data): return weights, rescaled_weights def amusr_fit(cov_C, cov_D, lambda_B=0., lambda_S=0., sparse_matrix=None, block_matrix=None, prior=None, - max_iter=MAX_ITER, tol=TOL, min_weight=MIN_WEIGHT_VAL): + max_iter=MAX_ITER, tol=TOL, rel_tol=REL_TOL, rel_tol_min_iter=10, min_weight=MIN_WEIGHT_VAL): """ Fits regression model in which the weights matrix W (predictors x tasks) is decomposed in two components: B that captures block structure across tasks @@ -321,6 +322,8 @@ def amusr_fit(cov_C, cov_D, lambda_B=0., lambda_S=0., sparse_matrix=None, block_ # Initialize weights combined_weights = sparse_matrix + block_matrix + + iter_tols = np.zeros(max_iter) for i in range(max_iter): # Keep a reference to the old combined_weights @@ -334,9 +337,18 @@ def amusr_fit(cov_C, cov_D, lambda_B=0., lambda_S=0., sparse_matrix=None, block_ combined_weights = sparse_matrix + block_matrix # If convergence tolerance reached, break loop and move on - if np.max(np.abs(combined_weights - _combined_weights_old)) < tol: + iter_tols[i] = np.max(np.abs(combined_weights - _combined_weights_old)) + + if iter_tols[i] < tol: break + # If the maximum over the last few iterations is less than the relative tolerance, break loop and move on + if i > rel_tol_min_iter: + lb_start, lb_stop = i - rel_tol_min_iter, i + iter_rel_max = iter_tols[lb_start: lb_stop] - iter_tols[lb_start - 1: lb_stop - 1] + if np.max(iter_rel_max) < rel_tol: + break + # Set small values of W to zero # Since we don't run the algorithm until update equals zero combined_weights[np.abs(combined_weights) < min_weight] = 0 @@ -439,7 +451,7 @@ def updateB(C, D, B, S, lamB, prior, n_tasks, n_features): # update current predictor B[j,:] = new_weights - return(B) + return B def _covariance_by_task(X, Y): """ From 33fb02547cb9178879db8ab4d13429110aba7c78 Mon Sep 17 00:00:00 2001 From: asistradition Date: Wed, 28 Apr 2021 16:33:46 -0400 Subject: [PATCH 09/14] Parameterize tolerances for amusr regression --- inferelator/distributed/dask_functions.py | 5 +-- inferelator/regression/amusr_regression.py | 41 +++++++++++++++++----- inferelator/tests/test_amusr.py | 6 ++-- 3 files changed, 38 insertions(+), 14 deletions(-) diff --git a/inferelator/distributed/dask_functions.py b/inferelator/distributed/dask_functions.py index 55e4695d..2fda80ad 100644 --- a/inferelator/distributed/dask_functions.py +++ b/inferelator/distributed/dask_functions.py @@ -16,7 +16,8 @@ def amusr_regress_dask(X, Y, priors, prior_weight, n_tasks, genes, tfs, G, remove_autoregulation=True, - lambda_Bs=None, lambda_Ss=None, Cs=None, Ss=None, regression_function=None): + lambda_Bs=None, lambda_Ss=None, Cs=None, Ss=None, regression_function=None, + tol=None, rel_tol=None): """ Execute multitask (AMUSR) @@ -54,7 +55,7 @@ def regression_maker(j, x_df, y_list, prior, tf): prior = format_prior(prior, gene, tasks, prior_weight, tfs=tf) return j, regression_function(x, y, tf, tasks, gene, prior, - lambda_Bs=lambda_Bs, lambda_Ss=lambda_Ss, Cs=Cs, Ss=Ss) + lambda_Bs=lambda_Bs, lambda_Ss=lambda_Ss, Cs=Cs, Ss=Ss, tol=tol, rel_tol=rel_tol) def response_maker(y_df, i): y = [] diff --git a/inferelator/regression/amusr_regression.py b/inferelator/regression/amusr_regression.py index b2b395da..4f020e53 100644 --- a/inferelator/regression/amusr_regression.py +++ b/inferelator/regression/amusr_regression.py @@ -22,13 +22,13 @@ MAX_ITER = 1000 TOL = 1e-2 -REL_TOL = 1e-4 +REL_TOL = None MIN_WEIGHT_VAL = 0.1 MIN_RSS = 1e-10 def run_regression_EBIC(X, Y, TFs, tasks, gene, prior, Cs=None, Ss=None, lambda_Bs=None, - lambda_Ss=None, scale_data=False, return_lambdas=False): + lambda_Ss=None, scale_data=False, return_lambdas=False, tol=TOL, rel_tol=REL_TOL): """ Run multitask regression. Search the regularization coefficient space and select the model with the lowest eBIC. @@ -97,7 +97,8 @@ def run_regression_EBIC(X, Y, TFs, tasks, gene, prior, Cs=None, Ss=None, lambda_ # Fit model combined_weights, sparse_matrix, block_matrix = amusr_fit(cov_C, cov_D, b, s, - sparse_matrix, block_matrix, prior) + sparse_matrix, block_matrix, prior, + tol=tol, rel_tol=rel_tol) # Score model ebic_score = ebic(X, Y, combined_weights, n_tasks, n_samples, n_preds) @@ -139,10 +140,13 @@ class AMuSR_regression(base_regression.BaseRegression): Cs = None Ss = None + tol = None + rel_tol = None + regression_function = staticmethod(run_regression_EBIC) def __init__(self, X, Y, tfs=None, genes=None, priors=None, prior_weight=1, remove_autoregulation=True, - lambda_Bs=None, lambda_Ss=None, Cs=None, Ss=None): + lambda_Bs=None, lambda_Ss=None, Cs=None, Ss=None, tol=TOL, rel_tol=REL_TOL): """ Set up a regression object for multitask regression :param X: list(InferelatorData) @@ -162,6 +166,8 @@ def __init__(self, X, Y, tfs=None, genes=None, priors=None, prior_weight=1, remo assert check.argument_type(tfs, (list, pd.Series, pd.Index), allow_none=True) assert check.argument_type(genes, (list, pd.Series, pd.Index), allow_none=True) assert check.argument_numeric(prior_weight) + assert check.argument_numeric(tol) + assert check.argument_numeric(rel_tol, allow_none=True) assert len(X) == len(Y) # Set the data into the regression object @@ -197,6 +203,10 @@ def __init__(self, X, Y, tfs=None, genes=None, priors=None, prior_weight=1, remo self.Cs = Cs self.Ss = Ss + # Set the tolerances into the regression object + self.tol = tol + self.rel_tol = rel_tol + def regress(self, regression_function=None): """ Execute multitask (AMUSR) @@ -210,7 +220,8 @@ def regress(self, regression_function=None): from inferelator.distributed.dask_functions import amusr_regress_dask return amusr_regress_dask(self.X, self.Y, self.priors, self.prior_weight, self.n_tasks, self.genes, self.tfs, self.G, remove_autoregulation=self.remove_autoregulation, - regression_function=regression_function) + regression_function=regression_function, + tol=self.tol, rel_tol=self.rel_tol) def regression_maker(j): level = 0 if j % 100 == 0 else 2 @@ -233,7 +244,8 @@ def regression_maker(j): prior = format_prior(self.priors, gene, tasks, self.prior_weight, tfs=tfs) return regression_function(x, y, tfs, tasks, gene, prior, Cs=self.Cs, Ss=self.Ss, - lambda_Bs=self.lambda_Bs, lambda_Ss=self.lambda_Ss) + lambda_Bs=self.lambda_Bs, lambda_Ss=self.lambda_Ss, + tol=self.tol, rel_tol=self.rel_tol) return MPControl.map(regression_maker, range(self.G)) @@ -301,6 +313,7 @@ def amusr_fit(cov_C, cov_D, lambda_B=0., lambda_S=0., sparse_matrix=None, block_ assert check.argument_type(max_iter, int) assert check.argument_type(tol, float) assert check.argument_type(min_weight, float) + assert check.argument_numeric(rel_tol, allow_none=True) assert cov_C.shape[0] == cov_D.shape[0] assert cov_C.shape[1] == cov_D.shape[1] assert cov_D.shape[1] == cov_D.shape[2] @@ -343,7 +356,7 @@ def amusr_fit(cov_C, cov_D, lambda_B=0., lambda_S=0., sparse_matrix=None, block_ break # If the maximum over the last few iterations is less than the relative tolerance, break loop and move on - if i > rel_tol_min_iter: + if rel_tol is not None and (i > rel_tol_min_iter): lb_start, lb_stop = i - rel_tol_min_iter, i iter_rel_max = iter_tols[lb_start: lb_stop] - iter_tols[lb_start - 1: lb_stop - 1] if np.max(iter_rel_max) < rel_tol: @@ -750,7 +763,11 @@ class AMUSRRegressionWorkflowMixin(base_regression._MultitaskRegressionWorkflowM lambda_Ss = None heuristic_Cs = None - def set_regression_parameters(self, prior_weight=None, lambda_Bs=None, lambda_Ss=None, heuristic_Cs=None): + tol = TOL + relative_tol = REL_TOL + + def set_regression_parameters(self, prior_weight=None, lambda_Bs=None, lambda_Ss=None, heuristic_Cs=None, + tol=None, relative_tol=None): """ Set regression parameters for AmUSR. @@ -771,12 +788,18 @@ def set_regression_parameters(self, prior_weight=None, lambda_Bs=None, lambda_Ss Defaults to np.logspace(np.log10(0.01), np.log10(10), 20)[::-1]. Does not have an effect if lambda_B is provided. :type heuristic_Cs: list(floats) or np.ndarray(floats) + :param tol: Convergence tolerance for amusr regression + :type tol: float + :param relative_tol: Relative convergence tolerance for amusr regression + :type relative_tol: float """ self._set_with_warning("prior_weight", prior_weight) self._set_with_warning("lambda_Bs", lambda_Bs) self._set_with_warning("lambda_Ss", lambda_Ss) self._set_with_warning("heuristic_Cs", heuristic_Cs) + self._set_without_warning("tol", tol) + self._set_without_warning("relative_tol", relative_tol) def run_bootstrap(self, bootstrap_idx): @@ -790,7 +813,7 @@ def run_bootstrap(self, bootstrap_idx): MPControl.sync_processes(pref="amusr_pre") regress = AMuSR_regression(x, y, tfs=self._regulators, genes=self._targets, priors=self._task_priors, prior_weight=self.prior_weight, lambda_Bs=self.lambda_Bs, lambda_Ss=self.lambda_Ss, - Cs=self.heuristic_Cs) + Cs=self.heuristic_Cs, tol=self.tol, rel_tol=self.relative_tol) return regress.run() diff --git a/inferelator/tests/test_amusr.py b/inferelator/tests/test_amusr.py index 7a02217e..1b2fdb85 100644 --- a/inferelator/tests/test_amusr.py +++ b/inferelator/tests/test_amusr.py @@ -347,7 +347,7 @@ def test_lamb_b(self): priors=self.workflow.priors_data, lambda_Bs=lamb_b) def is_passed(X, Y, TFs, tasks, gene, prior, Cs=None, Ss=None, lambda_Bs=None, - lambda_Ss=None, scale_data=False): + lambda_Ss=None, scale_data=False, tol=None, rel_tol=None): npt.assert_array_equal(lambda_Bs, lamb_b) @@ -369,7 +369,7 @@ def test_lamb_s(self): priors=self.workflow.priors_data, lambda_Ss=lamb_s) def is_passed(X, Y, TFs, tasks, gene, prior, Cs=None, Ss=None, lambda_Bs=None, - lambda_Ss=None, scale_data=False): + lambda_Ss=None, scale_data=False, tol=None, rel_tol=None): npt.assert_array_equal(lambda_Ss, lamb_s) @@ -391,7 +391,7 @@ def test_heuristic_c(self): priors=self.workflow.priors_data, Cs=set_Cs) def is_passed(X, Y, TFs, tasks, gene, prior, Cs=None, Ss=None, lambda_Bs=None, - lambda_Ss=None, scale_data=False): + lambda_Ss=None, scale_data=False, tol=None, rel_tol=None): npt.assert_array_equal(set_Cs, Cs) From c775a1796601c74fdf6cd8245f82ec96d4e4eb00 Mon Sep 17 00:00:00 2001 From: asistradition Date: Thu, 29 Apr 2021 09:35:26 -0400 Subject: [PATCH 10/14] Optimize UpdateS from profiling --- inferelator/regression/amusr_regression.py | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/inferelator/regression/amusr_regression.py b/inferelator/regression/amusr_regression.py index 4f020e53..cd4566d6 100644 --- a/inferelator/regression/amusr_regression.py +++ b/inferelator/regression/amusr_regression.py @@ -376,28 +376,34 @@ def updateS(C, D, B, S, lamS, prior, n_tasks, n_features): """ # update each task independently (shared penalty only) for k in range(n_tasks): + c = C[k]; d = D[k] b = B[:,k]; s = S[:,k] - p = prior[:,k] + p = prior[:, k] + p *= lamS # cycle through predictors + for j in range(n_features): # set sparse coefficient for predictor j to zero - s_tmp = deepcopy(s) - s_tmp[j] = 0. + s[j] = 0. + # calculate next coefficient based on fit only if d[j,j] == 0: alpha = 0 + continue else: - alpha = (c[j]-np.sum((b+s_tmp)*d[j]))/d[j,j] + alpha = (c[j]-np.sum((b+s)*d[j]))/d[j,j] + # lasso regularization - if abs(alpha) <= p[j]*lamS: + if abs(alpha) <= p[j]: s[j] = 0. else: - s[j] = alpha-(np.sign(alpha)*p[j]*lamS) + s[j] = alpha - (np.sign(alpha) * p[j]) + # update current task - S[:,k] = s + S[:, k] = s - return(S) + return S def updateB(C, D, B, S, lamB, prior, n_tasks, n_features): """ From 72029818bbbcaf4ddc6919f6343b20a35cb7b100 Mon Sep 17 00:00:00 2001 From: asistradition Date: Thu, 29 Apr 2021 10:38:02 -0400 Subject: [PATCH 11/14] Fix probable problem with lambda S --- inferelator/regression/amusr_regression.py | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/inferelator/regression/amusr_regression.py b/inferelator/regression/amusr_regression.py index cd4566d6..2e443ae8 100644 --- a/inferelator/regression/amusr_regression.py +++ b/inferelator/regression/amusr_regression.py @@ -377,10 +377,13 @@ def updateS(C, D, B, S, lamS, prior, n_tasks, n_features): # update each task independently (shared penalty only) for k in range(n_tasks): - c = C[k]; d = D[k] - b = B[:,k]; s = S[:,k] - p = prior[:, k] - p *= lamS + c = C[k] + d = D[k] + + b = B[:,k] + s = S[:,k] + p = prior[:, k] * lamS + # cycle through predictors for j in range(n_features): @@ -389,10 +392,9 @@ def updateS(C, D, B, S, lamS, prior, n_tasks, n_features): # calculate next coefficient based on fit only if d[j,j] == 0: - alpha = 0 - continue + alpha = 0. else: - alpha = (c[j]-np.sum((b+s)*d[j]))/d[j,j] + alpha = (c[j]- np.sum((b + s) * d[j])) / d[j,j] # lasso regularization if abs(alpha) <= p[j]: From 3fe2eef18a202211358ff4868563511ffda1aacd Mon Sep 17 00:00:00 2001 From: asistradition Date: Thu, 29 Apr 2021 10:46:15 -0400 Subject: [PATCH 12/14] Minor formatting --- inferelator/regression/amusr_regression.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/inferelator/regression/amusr_regression.py b/inferelator/regression/amusr_regression.py index 2e443ae8..0ce64d0f 100644 --- a/inferelator/regression/amusr_regression.py +++ b/inferelator/regression/amusr_regression.py @@ -380,8 +380,8 @@ def updateS(C, D, B, S, lamS, prior, n_tasks, n_features): c = C[k] d = D[k] - b = B[:,k] - s = S[:,k] + b = B[:, k] + s = S[:, k] p = prior[:, k] * lamS # cycle through predictors From f6e0923192c6424b7bf25289dfbf7c850e0725bc Mon Sep 17 00:00:00 2001 From: asistradition Date: Thu, 29 Apr 2021 14:53:53 -0400 Subject: [PATCH 13/14] Fix branch misalignment --- inferelator/amusr_workflow.py | 2 ++ inferelator/distributed/dask_functions.py | 1 + 2 files changed, 3 insertions(+) diff --git a/inferelator/amusr_workflow.py b/inferelator/amusr_workflow.py index 5d914973..0dcbe14e 100644 --- a/inferelator/amusr_workflow.py +++ b/inferelator/amusr_workflow.py @@ -45,6 +45,8 @@ class MultitaskLearningWorkflow(single_cell_workflow.SingleCellWorkflow): def _num_obs(self): if self._task_objects is not None: return sum([t if t is not None else 0 for t in map(lambda x: x._num_obs, self._task_objects)]) + elif self._task_design is not None: + return sum([t.num_obs for t in self._task_design]) else: return None diff --git a/inferelator/distributed/dask_functions.py b/inferelator/distributed/dask_functions.py index 2fda80ad..eff77de3 100644 --- a/inferelator/distributed/dask_functions.py +++ b/inferelator/distributed/dask_functions.py @@ -82,6 +82,7 @@ def response_maker(y_df, i): DaskController.client.cancel(scatter_x) DaskController.client.cancel(scatter_priors) + DaskController.client.restart() return result_list From 4be3aab9098a66aa9fc93870036feec2e24e3f78 Mon Sep 17 00:00:00 2001 From: asistradition Date: Thu, 29 Apr 2021 14:58:42 -0400 Subject: [PATCH 14/14] v0.5.5 --- docs/changelog.rst | 11 +++++++++++ docs/conf.py | 2 +- setup.py | 2 +- 3 files changed, 13 insertions(+), 2 deletions(-) diff --git a/docs/changelog.rst b/docs/changelog.rst index ccf8ee12..498ae87b 100644 --- a/docs/changelog.rst +++ b/docs/changelog.rst @@ -1,6 +1,17 @@ Change Log ========== +Inferelator v0.5.5 `April 29, 2021` +----------------------------------- + +New Functionality: + +- Added ``.set_regression_parameters(tol=None)`` to parameterize tolerances in AMuSR regression + +Code Refactoring: + +- Profiled and optimized AMuSR code + Inferelator v0.5.4 `April 23, 2021` ----------------------------------- diff --git a/docs/conf.py b/docs/conf.py index e230a7ea..09556691 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -23,7 +23,7 @@ author = 'Chris Jackson' # The full version, including alpha/beta/rc tags -release = 'v0.5.4' +release = 'v0.5.5' # -- General configuration --------------------------------------------------- diff --git a/setup.py b/setup.py index 49ad4503..0dd12e14 100644 --- a/setup.py +++ b/setup.py @@ -2,7 +2,7 @@ from setuptools import setup, find_packages # Current Inferelator Version Number -version = "0.5.4" +version = "0.5.5" # Description from README.md