adjoint and algebra solver

README.rst

@@ -42,17 +42,17 @@ this command:
 
     pip install maud-metabolic-models
 
-Cmdstanpy depends on `cmdstan <https://github.com/stan-dev/cmdstan>`_, 
-which in turn requires a c++ toolchain. On some computers you will have to 
-install these in order to use Maud. You will hit an error at the next step if 
+Cmdstanpy depends on `cmdstan <https://github.com/stan-dev/cmdstan>`_,
+which in turn requires a c++ toolchain. On some computers you will have to
+install these in order to use Maud. You will hit an error at the next step if
 this applies to your computer. Luckily cmdstanpy comes with commands that
-can do the necessary installing for you. On windows the c++ toolchain can be installed with 
+can do the necessary installing for you. On windows the c++ toolchain can be installed with
 the following powershell commands:
 
 
 Usage
 =====
-Maud is used from the command line. To see all the available commands try running 
+Maud is used from the command line. To see all the available commands try running
 
 .. code-block:: console
 

docs/inputting.md

@@ -46,7 +46,27 @@ The following optional fields can also be specified:
 * `molecule_unit` A unit for counting molecules, like 'mol' or 'mmol'
 * `volume_unit` A unit for measuring volume, like 'L'
 * `energy_unit` A unit for measuring energy, like 'J' or 'kJ'
-
+* `algebra_solve` A boolean argument that specifies if the powell algebra solver will be used to determine the steady state
+* `adjoint_solve` A boolean argument that specifies if the adjoint ode solver will be used in place of the bdf solver.
+
+The following fields can be specified by the ODE solver. The variables and definitions can be found both in the [`Stan documentation`](`https://mc-stan.org/docs/functions-reference/functions-ode-solver.html`) with more detail in the [`Sundials documentation`](`https://sundials.readthedocs.io/en/latest/cvodes/Usage/SIM.html#c.CVodeSStolerances`):
+
+* `rel_tol` Relative tolerance of the BDF solver
+* `abs_tol` Absolute tolerance of the BDF solver -- applied to every state variable.
+* `max_num_steps` Maximum number of steps taken by both the algebra and ODE solvers.
+* `timepoint` The integration time solved for either ODE solver.
+* `rel_tol_forward` The relative tolerance used during the forward solve
+* `abs_tol_forward` The absolute tolerance used during the backwards solve -- can be converted into a vector
+* `rel_tol_backward` The relative tolerance used during the backwards solve
+* `abs_tol_backward` The absolute tolerance used during the backwards solve -- can be converted into a vector
+* `rel_tol_quadrature` The relative tolerance used in the backwards quadrature problem
+* `abs_tol_quadrature` The absolute tolerance used in the backwards quadrature problem
+* `num_steps_between_checkpoints` number of steps between checkpoints in forward solution
+* `interpolation_polynomial` 1 for hermite and 2 for polynomial
+* `solver_forward` Solver used for forward ODE problem: 1=Adams (non-stiff), 2=BDF (stiff)
+* `solver_backward` Solver used for backward ODE problem: 1=Adams (non-stiff), 2=BDF (stiff)
+* `algebra_func_tol` After convergence of the solver, the proposed solution is plugged into the algebraic system and its norm is compared to the function tolerance. If the norm is below the function tolerance, the solution is deemed acceptable
+* `algebra_rel_tol` The relative tolerance is the estimated relative error of the solver and serves to test if a satisfactory solution has been found
 
 Here is an example configuration file:
 
@@ -71,6 +91,18 @@ Here is an example configuration file:
     rel_tol = 1e-4
     max_num_steps = 1e6
     timepoint = 1e3
+    rel_tol_forward = 1e-9
+    abs_tol_forward = 1e-9
+    rel_tol_backward = 1e-9
+    abs_tol_backward = 1e-9
+    rel_tol_quadrature = 1e-9
+    abs_tol_quadrature = 1e-9
+    num_steps_between_checkpoints = int(100000)
+    interpolation_polynomial = 2
+    solver_forward = 2
+    solver_backward = 2
+    algebra_func_tol = 1e-6
+    algebra_rel_tol = 1e-10
 
 This file tells Maud that a file representing a kinetic model can be found at
 the relative path `kinetic_model.toml`, and that priors and experimental

docs/theory.md

@@ -44,7 +44,7 @@ want. Questions about the network can be formulated in terms of integrals over
 the posterior distribution. For example, suppose that our thermokinetic model
 says that the flux through a reaction is given by the function $flux:
 \Theta\rightarrow\mathbb{R}$; our model's probability, given measurements $y$,
-that the flux is greater than $x$ is then $\int_x^\inftyp(\theta\mid
+that the flux is greater than $x$ is then $\int_x^\infty(\theta\mid
 y)d\theta$.
 
 The rest of this page sets out how Maud fleshes out the three sub-models and
@@ -128,7 +128,7 @@ metabolite's 'formation energy'.
 * $n$ is a vector representing the number of charges transported by each
 reaction.
 * $F$ is the Faraday constant (a number)
-* $\psi$ is a vector representin each reaction's membrane potential (these
+* $\psi$ is a vector representing each reaction's membrane potential (these
 numbers only matter for reactions that transport non-zero charge)
 
 Note that, for reactions with zero transported charge, the thermodynamic effect
@@ -277,7 +277,7 @@ model for these measurements is imperfect in at least these ways:
   absolute concentrations.
 
 These problems also apply to flux measurements, but there is another issue:
-flux meausrements are derived from k
+flux measurements are derived from k
 
 (sec-prior-model)=
 ### Prior model

maud/data/example_inputs/example_ode/config.toml

@@ -1,14 +1,30 @@
-name = "test-ode"
+name = "test-ode-adjoint"
 kinetic_model_file = "kinetic_model.toml"
 priors_file = "priors.toml"
 experiments_file = "experiments.toml"
 likelihood = true
+adjoint_solve = true
+algebra_solve = false
 
 [cmdstanpy_config]
 iter_warmup = 0
 iter_sampling = 1
 chains = 1
 
 [ode_config]
+rel_tol = 1e-9
+abs_tol = 1e-9
 max_num_steps = 1e9
-timepoint = 1e3
+timepoint = 500
+rel_tol_forward = 1e-12
+abs_tol_forward = 1e-12
+rel_tol_backward = 1e-12
+abs_tol_backward = 1e-12
+rel_tol_quadrature = 1e-12
+abs_tol_quadrature = 1e-12
+num_steps_between_checkpoints = 100000
+interpolation_polynomial = 2
+solver_forward = 2
+solver_backward = 2
+algebra_func_tol = 1e-6
+algebra_rel_tol = 1e-10

maud/data/example_inputs/example_ode/input_data_train_adjoint.json

@@ -0,0 +1 @@
+{"N_mic":4,"N_edge_sub":5,"N_edge_prod":5,"N_edge":5,"N_unbalanced":2,"N_enzyme":5,"N_phosphorylation":0,"N_pme":0,"N_competitive_inhibition":1,"N_allostery":2,"N_allosteric_enzyme":2,"N_drain":0,"N_km":10,"N_sub_km":5,"N_prod_km":5,"S":[[-1,-1,-1,0,0],[1,0,0,-1,0],[0,1,1,0,-1],[0,0,0,1,1]],"N_reaction":4,"N_metabolite":4,"balanced_mic_ix":[2,3],"unbalanced_mic_ix":[1,4],"ci_mic_ix":[4],"edge_type":[1,1,1,1,1],"edge_to_enzyme":[1,2,3,4,5],"edge_to_tc":[0,1,2,0,0],"edge_to_drain":[0,0,0,0,0],"edge_to_reaction":[1,2,2,3,4],"water_stoichiometry":[0.0,0.0,0.0,0.0,0.0],"transported_charge":[0.0,0.0,0.0,0.0,0.0],"mic_to_met":[1,2,3,4],"subunits":[1,1,1,1,1],"sub_by_edge_long":[1,1,1,2,3],"sub_by_edge_bounds":[[1,1],[2,2],[3,3],[4,4],[5,5]],"prod_by_edge_long":[2,3,3,4,4],"prod_by_edge_bounds":[[1,1],[2,2],[3,3],[4,4],[5,5]],"sub_km_ix_by_edge_long":[1,3,5,7,9],"sub_km_ix_by_edge_bounds":[[1,1],[2,2],[3,3],[4,4],[5,5]],"prod_km_ix_by_edge_long":[2,4,6,8,10],"prod_km_ix_by_edge_bounds":[[1,1],[2,2],[3,3],[4,4],[5,5]],"ci_ix_long":[1],"ci_ix_bounds":[[1,1],[2,1],[2,1],[2,1],[2,1]],"allostery_ix_long":[1,2],"allostery_ix_bounds":[[1,0],[1,1],[2,2],[3,2],[3,2]],"allostery_type":[1,2],"allostery_mic":[3,3],"phosphorylation_ix_long":[],"phosphorylation_ix_bounds":[[1,0],[1,0],[1,0],[1,0],[1,0]],"phosphorylation_type":[],"phosphorylation_pme":[],"priors_km":[[-0.6931471805599453,0.0,0.6931471805599453,-0.6931471805599453,0.0,0.6931471805599453,-0.6931471805599453,0.6931471805599453,0.0,1.0986122886681098],[0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2]],"priors_ki":[[0.0],[0.2]],"priors_kcat":[[-0.6931471805599453,0.6931471805599453,0.0,0.6931471805599453,0.0],[0.2,0.2,0.2,0.2,0.2]],"priors_dissociation_constant":[[-1.2039728043259361,-0.10536051565782628],[0.2,0.2]],"priors_transfer_constant":[[0.0,0.0],[0.2,0.2]],"priors_kcat_pme":[[],[]],"priors_drain_train":[[[]],[[]]],"priors_conc_enzyme_train":[[[0.0,0.6931471805599453,1.0986122886681098,0.6931471805599453,1.0986122886681098]],[[0.2,0.2,0.2,0.2,0.2]]],"priors_conc_unbalanced_train":[[[1.6094379124341003,-0.6931471805599453]],[[0.2,0.2]]],"priors_conc_pme_train":[[[]],[[]]],"priors_psi_train":[[0.0],[2.0]],"prior_loc_dgf":[10.0,2.0,5.0,0.0],"prior_cov_dgf":[[0.2,0.0,0.0,0.0],[0.0,0.2,0.0,0.0],[0.0,0.0,0.2,0.0],[0.0,0.0,0.0,0.2]],"N_experiment_train":1,"N_flux_measurement_train":2,"N_enzyme_measurement_train":0,"N_conc_measurement_train":2,"N_enzyme_knockout_train":0,"N_pme_knockout_train":0,"temperature_train":[298.15],"enzyme_knockout_train_long":[],"enzyme_knockout_train_bounds":[[1,0]],"pme_knockout_train_long":[],"pme_knockout_train_bounds":[[1,0]],"yconc_train":[0.323117,3.02187],"sigma_yconc_train":[0.1,0.1],"experiment_yconc_train":[1,1],"mic_ix_yconc_train":[2,3],"yflux_train":[0.421816,2.11674],"sigma_yflux_train":[0.01,0.01],"experiment_yflux_train":[1,1],"reaction_yflux_train":[3,4],"yenz_train":[],"sigma_yenz_train":[],"experiment_yenz_train":[],"enzyme_yenz_train":[],"likelihood":1,"drain_small_conc_corrector":1e-6,"reject_non_steady":1,"steady_state_threshold_abs":1e-8,"steady_state_threshold_rel":0.001,"rel_tol":1e-9,"abs_tol":1e-9,"timepoint":500.0,"max_num_steps":1000000000,"adjoint_solve":1,"algebra_solve":0,"rel_tol_forward":1e-12,"abs_tol_forward":1e-12,"rel_tol_backward":1e-12,"abs_tol_backward":1e-12,"rel_tol_quadrature":1e-12,"abs_tol_quadrature":1e-12,"num_steps_between_checkpoints":100000,"interpolation_polynomial":2,"solver_forward":2,"solver_backward":2,"algebra_func_tol":1e-6,"algebra_rel_tol":1e-10,"conc_init":[[0.323117,3.02187]]}

maud/data/example_inputs/example_ode/input_data_train_adjoint_algebra.json

@@ -0,0 +1 @@
+{"N_mic":4,"N_edge_sub":5,"N_edge_prod":5,"N_edge":5,"N_unbalanced":2,"N_enzyme":5,"N_phosphorylation":0,"N_pme":0,"N_competitive_inhibition":1,"N_allostery":2,"N_allosteric_enzyme":2,"N_drain":0,"N_km":10,"N_sub_km":5,"N_prod_km":5,"S":[[-1,-1,-1,0,0],[1,0,0,-1,0],[0,1,1,0,-1],[0,0,0,1,1]],"N_reaction":4,"N_metabolite":4,"balanced_mic_ix":[2,3],"unbalanced_mic_ix":[1,4],"ci_mic_ix":[4],"edge_type":[1,1,1,1,1],"edge_to_enzyme":[1,2,3,4,5],"edge_to_tc":[0,1,2,0,0],"edge_to_drain":[0,0,0,0,0],"edge_to_reaction":[1,2,2,3,4],"water_stoichiometry":[0.0,0.0,0.0,0.0,0.0],"transported_charge":[0.0,0.0,0.0,0.0,0.0],"mic_to_met":[1,2,3,4],"subunits":[1,1,1,1,1],"sub_by_edge_long":[1,1,1,2,3],"sub_by_edge_bounds":[[1,1],[2,2],[3,3],[4,4],[5,5]],"prod_by_edge_long":[2,3,3,4,4],"prod_by_edge_bounds":[[1,1],[2,2],[3,3],[4,4],[5,5]],"sub_km_ix_by_edge_long":[1,3,5,7,9],"sub_km_ix_by_edge_bounds":[[1,1],[2,2],[3,3],[4,4],[5,5]],"prod_km_ix_by_edge_long":[2,4,6,8,10],"prod_km_ix_by_edge_bounds":[[1,1],[2,2],[3,3],[4,4],[5,5]],"ci_ix_long":[1],"ci_ix_bounds":[[1,1],[2,1],[2,1],[2,1],[2,1]],"allostery_ix_long":[1,2],"allostery_ix_bounds":[[1,0],[1,1],[2,2],[3,2],[3,2]],"allostery_type":[1,2],"allostery_mic":[3,3],"phosphorylation_ix_long":[],"phosphorylation_ix_bounds":[[1,0],[1,0],[1,0],[1,0],[1,0]],"phosphorylation_type":[],"phosphorylation_pme":[],"priors_km":[[-0.6931471805599453,0.0,0.6931471805599453,-0.6931471805599453,0.0,0.6931471805599453,-0.6931471805599453,0.6931471805599453,0.0,1.0986122886681098],[0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2]],"priors_ki":[[0.0],[0.2]],"priors_kcat":[[-0.6931471805599453,0.6931471805599453,0.0,0.6931471805599453,0.0],[0.2,0.2,0.2,0.2,0.2]],"priors_dissociation_constant":[[-1.2039728043259361,-0.10536051565782628],[0.2,0.2]],"priors_transfer_constant":[[0.0,0.0],[0.2,0.2]],"priors_kcat_pme":[[],[]],"priors_drain_train":[[[]],[[]]],"priors_conc_enzyme_train":[[[0.0,0.6931471805599453,1.0986122886681098,0.6931471805599453,1.0986122886681098]],[[0.2,0.2,0.2,0.2,0.2]]],"priors_conc_unbalanced_train":[[[1.6094379124341003,-0.6931471805599453]],[[0.2,0.2]]],"priors_conc_pme_train":[[[]],[[]]],"priors_psi_train":[[0.0],[2.0]],"prior_loc_dgf":[10.0,2.0,5.0,0.0],"prior_cov_dgf":[[0.2,0.0,0.0,0.0],[0.0,0.2,0.0,0.0],[0.0,0.0,0.2,0.0],[0.0,0.0,0.0,0.2]],"N_experiment_train":1,"N_flux_measurement_train":2,"N_enzyme_measurement_train":0,"N_conc_measurement_train":2,"N_enzyme_knockout_train":0,"N_pme_knockout_train":0,"temperature_train":[298.15],"enzyme_knockout_train_long":[],"enzyme_knockout_train_bounds":[[1,0]],"pme_knockout_train_long":[],"pme_knockout_train_bounds":[[1,0]],"yconc_train":[0.323117,3.02187],"sigma_yconc_train":[0.1,0.1],"experiment_yconc_train":[1,1],"mic_ix_yconc_train":[2,3],"yflux_train":[0.421816,2.11674],"sigma_yflux_train":[0.01,0.01],"experiment_yflux_train":[1,1],"reaction_yflux_train":[3,4],"yenz_train":[],"sigma_yenz_train":[],"experiment_yenz_train":[],"enzyme_yenz_train":[],"likelihood":1,"drain_small_conc_corrector":1e-6,"reject_non_steady":1,"steady_state_threshold_abs":1e-8,"steady_state_threshold_rel":0.001,"rel_tol":1e-9,"abs_tol":1e-9,"timepoint":500.0,"max_num_steps":1000000000,"adjoint_solve":1,"algebra_solve":1,"rel_tol_forward":1e-12,"abs_tol_forward":1e-12,"rel_tol_backward":1e-12,"abs_tol_backward":1e-12,"rel_tol_quadrature":1e-12,"abs_tol_quadrature":1e-12,"num_steps_between_checkpoints":100000,"interpolation_polynomial":2,"solver_forward":2,"solver_backward":2,"algebra_func_tol":1e-6,"algebra_rel_tol":1e-10,"conc_init":[[0.323117,3.02187]]}