diff --git a/maud/cli.py b/maud/cli.py index d22db928..87781f2f 100644 --- a/maud/cli.py +++ b/maud/cli.py @@ -220,13 +220,14 @@ def do_simulate(data_path, output_dir, n): .T ) print("\n\nSimulated conserved moiety pools:") - print( - idata.posterior["conc_moiety_pool_train"] - .mean(dim=["chain", "draw"]) - .to_series() - .unstack() - .T - ) + if "conc_moiety_pool_train" in idata.posterior.keys(): + print( + idata.posterior["conc_moiety_pool_train"] + .mean(dim=["chain", "draw"]) + .to_series() + .unstack() + .T + ) print("\n\nSimulated fluxes:") print( idata.posterior["flux_train"] diff --git a/maud/getting_idatas.py b/maud/getting_idatas.py index ae4f137c..399f5162 100644 --- a/maud/getting_idatas.py +++ b/maud/getting_idatas.py @@ -29,7 +29,10 @@ def get_idata(csvs: List[str], mi: MaudInput, mode: str) -> az.InferenceData: MeasurementType.ENZYME, ] ) - dependent_mics = [m.solve for m in mi.kinetic_model.conserved_moiety] + dependent_mics = ( + [m.solve for m in mi.kinetic_model.conserved_moiety] + if mi.kinetic_model.conserved_moiety is not None else [] + ) coords = { "enzymes": [e.id for e in mi.kinetic_model.enzymes], "experiments": [e.id for e in experiments], @@ -39,7 +42,8 @@ def get_idata(csvs: List[str], mi: MaudInput, mode: str) -> az.InferenceData: "mics": [m.id for m in mi.kinetic_model.mics], "conserved_moieties": [ cm.id for cm in mi.kinetic_model.conserved_moiety - ], + ] if mi.kinetic_model.conserved_moiety is not None + else [], "edges": [e.id for e in mi.kinetic_model.edges], "edges1": [e.id for e in mi.kinetic_model.edges], "unbalanced_mics": [ diff --git a/maud/running_stan.py b/maud/running_stan.py index e6128510..cfdef759 100644 --- a/maud/running_stan.py +++ b/maud/running_stan.py @@ -311,13 +311,13 @@ def predict( "output_dir": output_dir, "iter_warmup": 0, "iter_sampling": 1, + "adapt_engaged": False, "fixed_param": True, - "show_progress": False, } if mi.config.cmdstanpy_config_predict is not None: sample_args = { - **sample_args, **mi.config.cmdstanpy_config_predict, + **sample_args, } sample_args["chains"] = 1 mcmc_draw = model.sample(**sample_args) diff --git a/maud/stan/out_of_sample_model.stan b/maud/stan/out_of_sample_model.stan index 68cfc820..8da0cb74 100644 --- a/maud/stan/out_of_sample_model.stan +++ b/maud/stan/out_of_sample_model.stan @@ -2,9 +2,12 @@ data { // network properties int N_mic; + int N_pool; int N_edge_sub; int N_edge_prod; int N_unbalanced; + int N_independent; + int N_dependent; int N_metabolite; int N_km; int N_sub_km; @@ -19,7 +22,9 @@ data { int N_pme; // phosphorylation modifying enzyme int N_competitive_inhibition; matrix[N_mic, N_edge] S; - array[N_mic - N_unbalanced] int balanced_mic_ix; + matrix[N_pool, N_independent] left_nullspace_independent; + array[N_independent] int independent_bal_ix; + array[N_dependent] int dependent_bal_ix; array[N_unbalanced] int unbalanced_mic_ix; array[N_competitive_inhibition] int ci_mic_ix; array[N_edge] int edge_type; // 1 = reversible modular rate law, 2 = drain @@ -58,22 +63,21 @@ data { array[N_pme_knockout_test] int pme_knockout_test_long; array[N_experiment_test, 2] int pme_knockout_test_bounds; vector[N_experiment_test] temperature_test; - // hardcoded priors + // hardcoded priors (boundary conditions) array[2, N_experiment_test] vector[N_pme] priors_conc_phos_test; array[2, N_experiment_test] vector[N_unbalanced] priors_conc_unbalanced_test; + array[2, N_experiment_test] vector[N_pool] priors_conc_moiety_pool_test; array[2, N_experiment_test] vector[N_enzyme] priors_conc_enzyme_test; array[2, N_experiment_test] vector[N_drain] priors_drain_test; // configuration - array[N_experiment_test] vector[N_mic - N_unbalanced] conc_init; + real timepoint; + array[N_experiment_test] vector[N_independent] conc_init; real rel_tol_ode; real abs_tol_ode; int max_num_steps_ode; real rel_tol_alg; real abs_tol_alg; int max_num_steps_alg; - real steady_state_threshold_abs; - real steady_state_threshold_rel; - real steady_state_penalty_rel; int likelihood; // set to 0 for priors-only mode real drain_small_conc_corrector; int penalize_non_steady; @@ -95,6 +99,7 @@ generated quantities { array[N_experiment_test] vector[N_reaction] flux_test; array[N_experiment_test] vector[N_pme] conc_pme_test; array[N_experiment_test] vector[N_unbalanced] conc_unbalanced_test; + array[N_experiment_test] vector[N_pool] conc_moiety_pool_test; array[N_experiment_test] vector[N_enzyme] conc_enzyme_test; array[N_experiment_test] vector[N_drain] drain_test; array[N_experiment_test] vector[N_edge] free_enzyme_ratio_test; @@ -110,6 +115,8 @@ generated quantities { priors_conc_phos_test[2, e])); conc_unbalanced_test[e] = to_vector(lognormal_rng(priors_conc_unbalanced_test[1, e], priors_conc_unbalanced_test[2, e])); + conc_moiety_pool_test[e] = to_vector(lognormal_rng(priors_conc_moiety_pool_test[1, e], + priors_conc_moiety_pool_test[2, e])); conc_enzyme_test[e] = to_vector(lognormal_rng(priors_conc_enzyme_test[1, e], priors_conc_enzyme_test[2, e])); } @@ -118,7 +125,7 @@ generated quantities { flux_test[e] = rep_vector(0, N_reaction); vector[N_enzyme] conc_enzyme_experiment = conc_enzyme_test[e]; vector[N_pme] conc_pme_experiment = conc_pme_test[e]; - vector[N_mic - N_unbalanced] conc_balanced_experiment; + vector[N_independent] conc_independent_balanced_experiment; int N_eko_experiment = measure_ragged(enzyme_knockout_test_bounds, e); int N_pko_experiment = measure_ragged(pme_knockout_test_bounds, e); if (N_eko_experiment > 0) { @@ -134,42 +141,42 @@ generated quantities { e); conc_pme_experiment[pko_experiment] = rep_vector(0, N_pko_experiment); } - conc_balanced_experiment = solve_newton_tol(maud_ae_system, - conc_init[e], - rel_tol_alg, abs_tol_alg, - max_num_steps_alg, - rel_tol_ode, abs_tol_ode, - max_num_steps_ode, - conc_unbalanced_test[e], - balanced_mic_ix, - unbalanced_mic_ix, - conc_enzyme_experiment, - dgr_test[e], kcat, km, ki, - transfer_constant, - dissociation_constant, kcat_pme, - conc_pme_experiment, - drain_test[e], - temperature_test[e], - drain_small_conc_corrector, S, - subunits, edge_type, - edge_to_enzyme, edge_to_drain, - ci_mic_ix, sub_km_ix_by_edge_long, - sub_km_ix_by_edge_bounds, - prod_km_ix_by_edge_long, - prod_km_ix_by_edge_bounds, - sub_by_edge_long, - sub_by_edge_bounds, - prod_by_edge_long, - prod_by_edge_bounds, ci_ix_long, - ci_ix_bounds, allostery_ix_long, - allostery_ix_bounds, - allostery_type, allostery_mic, - edge_to_tc, - phosphorylation_ix_long, - phosphorylation_ix_bounds, - phosphorylation_type, - phosphorylation_pme); - conc_test[e, balanced_mic_ix] = conc_balanced_experiment; + conc_independent_balanced_experiment = ode_bdf_tol(dbalanced_dt, conc_init[e], 0, {timepoint}, + rel_tol_ode, abs_tol_ode, max_num_steps_ode, + conc_unbalanced_test[e], + conc_moiety_pool_test[e], + independent_bal_ix, + dependent_bal_ix, + unbalanced_mic_ix, + conc_enzyme_experiment, + dgr_test[e], kcat, km, ki, + transfer_constant, + dissociation_constant, kcat_pme, + conc_pme_experiment, + drain_test[e], + temperature_test[e], + drain_small_conc_corrector, S, + left_nullspace_independent, + subunits, edge_type, + edge_to_enzyme, edge_to_drain, + ci_mic_ix, sub_km_ix_by_edge_long, + sub_km_ix_by_edge_bounds, + prod_km_ix_by_edge_long, + prod_km_ix_by_edge_bounds, + sub_by_edge_long, + sub_by_edge_bounds, + prod_by_edge_long, + prod_by_edge_bounds, ci_ix_long, + ci_ix_bounds, allostery_ix_long, + allostery_ix_bounds, + allostery_type, allostery_mic, + edge_to_tc, + phosphorylation_ix_long, + phosphorylation_ix_bounds, + phosphorylation_type, + phosphorylation_pme)[1]; + conc_test[e, independent_bal_ix] = conc_independent_balanced_experiment; + conc_test[e, dependent_bal_ix] = conc_moiety_pool_test[e]- left_nullspace_independent * conc_independent_balanced_experiment; conc_test[e, unbalanced_mic_ix] = conc_unbalanced_test[e]; vector[N_edge] edge_flux = get_edge_flux(conc_test[e], conc_enzyme_experiment, diff --git a/pyproject.toml b/pyproject.toml index 2496fa95..7d42aeac 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -24,7 +24,7 @@ readme = "README.rst" requires-python = ">=3.10" dependencies = [ "pip >= 20", - "arviz >= 0.18.0", + "arviz>=0.18.0", "importlib_resources >= 3.2", "numpy", "scipy", @@ -32,7 +32,7 @@ dependencies = [ "pandas", "matplotlib", "toml", - "cmdstanpy >= 1.2.1", + "cmdstanpy>=1.2.3", "click", "depinfo == 1.7.0", "pydantic >= 2.0", @@ -118,6 +118,8 @@ markers = ["raises"] # Enable the following rules: # pycodestyle (`E`), Pyflakes (`F`), pycodestyle (`W`), flake8-bugbear (`B`) lint.select = ["E", "F", "B", "W", "D"] +# Enable the following rules: +# pycodestyle (`E`), Pyflakes (`F`), pycodestyle (`W`), flake8-bugbear (`B`) lint.ignore = [ "B905", # Use zip() without a `strict` parameter. "D107", # Init methods can be undocumented. @@ -125,7 +127,11 @@ lint.ignore = [ "D213", # https://beta.ruff.rs/docs/rules/multi-line-summary-second-line/ "D104", # __init__.py can be empty. ] +# Enable the following rules: +# pycodestyle (`E`), Pyflakes (`F`), pycodestyle (`W`), flake8-bugbear (`B`) lint.fixable = ["ALL"] +# Enable the following rules: +# pycodestyle (`E`), Pyflakes (`F`), pycodestyle (`W`), flake8-bugbear (`B`) lint.unfixable = [] line-length = 80 target-version = "py310"