Skip to content

Commit

Permalink
fixing out_of_sample_model
Browse files Browse the repository at this point in the history
  • Loading branch information
@NicholasCowie
NicholasCowie committed Jun 7, 2024
1 parent 84dca7c commit 1e7f392
Show file tree
Hide file tree
Showing 5 changed files with 74 additions and 56 deletions.
15 changes: 8 additions & 7 deletions maud/cli.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"]
Expand Down
8 changes: 6 additions & 2 deletions maud/getting_idatas.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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],
Expand All @@ -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": [
Expand Down
4 changes: 2 additions & 2 deletions maud/running_stan.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
93 changes: 50 additions & 43 deletions maud/stan/out_of_sample_model.stan
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,9 +2,12 @@
data {
// network properties
int<lower=1> N_mic;
int<lower=0> N_pool;
int<lower=1> N_edge_sub;
int<lower=1> N_edge_prod;
int<lower=1> N_unbalanced;
int<lower=0> N_independent;
int<lower=0> N_dependent;
int<lower=1> N_metabolite;
int<lower=1> N_km;
int<lower=1> N_sub_km;
Expand All @@ -19,7 +22,9 @@ data {
int<lower=0> N_pme; // phosphorylation modifying enzyme
int<lower=0> N_competitive_inhibition;
matrix[N_mic, N_edge] S;
array[N_mic - N_unbalanced] int<lower=1, upper=N_mic> balanced_mic_ix;
matrix[N_pool, N_independent] left_nullspace_independent;
array[N_independent] int<lower=1, upper=N_mic> independent_bal_ix;
array[N_dependent] int<lower=1, upper=N_mic> dependent_bal_ix;
array[N_unbalanced] int<lower=1, upper=N_mic> unbalanced_mic_ix;
array[N_competitive_inhibition] int<lower=1, upper=N_mic> ci_mic_ix;
array[N_edge] int<lower=1, upper=3> edge_type; // 1 = reversible modular rate law, 2 = drain
Expand Down Expand Up @@ -58,22 +63,21 @@ data {
array[N_pme_knockout_test] int<lower=0, upper=N_pme> 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<lower=0>[N_mic - N_unbalanced] conc_init;
real timepoint;
array[N_experiment_test] vector<lower=0>[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<lower=0, upper=1> likelihood; // set to 0 for priors-only mode
real drain_small_conc_corrector;
int<lower=0, upper=1> penalize_non_steady;
Expand All @@ -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;
Expand All @@ -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]));
}
Expand All @@ -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) {
Expand All @@ -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,
Expand Down
10 changes: 8 additions & 2 deletions pyproject.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,15 +24,15 @@ readme = "README.rst"
requires-python = ">=3.10"
dependencies = [
"pip >= 20",
"arviz >= 0.18.0",
"arviz>=0.18.0",
"importlib_resources >= 3.2",
"numpy",
"scipy",
"sympy",
"pandas",
"matplotlib",
"toml",
"cmdstanpy >= 1.2.1",
"cmdstanpy>=1.2.3",
"click",
"depinfo == 1.7.0",
"pydantic >= 2.0",
Expand Down Expand Up @@ -118,14 +118,20 @@ 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.
"D203", # Class docstrints don't need blank lines before them.
"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"
Expand Down

0 comments on commit 1e7f392

Please sign in to comment.