From dcff7b7d634c08f4970072e83935cbe48623ab1a Mon Sep 17 00:00:00 2001 From: thalassemia <67928790+thalassemia@users.noreply.github.com> Date: Thu, 18 Jan 2024 16:45:13 -0800 Subject: [PATCH 01/11] Cache expensive sparse array operation --- reconstruction/ecoli/dataclasses/process/transcription.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/reconstruction/ecoli/dataclasses/process/transcription.py b/reconstruction/ecoli/dataclasses/process/transcription.py index fcb3cd10a..3b987b3cc 100755 --- a/reconstruction/ecoli/dataclasses/process/transcription.py +++ b/reconstruction/ecoli/dataclasses/process/transcription.py @@ -5,6 +5,7 @@ TODO: handle ppGpp and DksA-ppGpp regulation separately """ +from functools import cache from typing import cast import numpy as np @@ -908,6 +909,7 @@ def cistron_id_to_rna_indexes(self, cistron_id): return self.cistron_tu_mapping_matrix.getrow( self._cistron_id_to_index[cistron_id]).nonzero()[1] + @cache def rna_id_to_cistron_indexes(self, rna_id): """ Returns the indexes of cistrons that constitute the given transcription From 7accb6adc1df57022751477a90c9e6652f4664b3 Mon Sep 17 00:00:00 2001 From: thalassemia <67928790+thalassemia@users.noreply.github.com> Date: Thu, 18 Jan 2024 16:47:25 -0800 Subject: [PATCH 02/11] Remove need for expensive array intersect --- wholecell/utils/polymerize.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/wholecell/utils/polymerize.py b/wholecell/utils/polymerize.py index 40134010e..b3ab1d78f 100644 --- a/wholecell/utils/polymerize.py +++ b/wholecell/utils/polymerize.py @@ -155,6 +155,8 @@ def _prepare_running_values(self): self._progress = np.zeros(self._nSequences, np.int64) self._activeSequencesIndexes = self._activeSequencesIndexes[ self._sequenceReactions[:, self._currentStep]] + self.elongation_rates = self.elongation_rates[ + self._sequenceReactions[:, self._currentStep]] self._update_elongation_resource_demands() self._monomerHistory = np.empty((self._maxElongation, self._nMonomers), np.int64) @@ -221,8 +223,8 @@ def _elongate_to_limit(self): step = self._currentStep + projectionIndex level = self.elongation_rates * (step + 1) last_unit = (level - np.floor(level)) - elongating = np.where(self.elongation_rates > last_unit)[0] - active = np.intersect1d(self._activeSequencesIndexes, elongating) + active = self._activeSequencesIndexes[ + self.elongation_rates > last_unit] else: active = self._activeSequencesIndexes @@ -345,6 +347,7 @@ def _finalize_resource_limited_elongations(self): # Cull sequences self._activeSequencesIndexes = self._activeSequencesIndexes[~sequencesToCull] + self.elongation_rates = self.elongation_rates[~sequencesToCull] def _update_elongation_resource_demands(self): ''' From 1da4703a535a4e996e91fa89b535c893f7047b2a Mon Sep 17 00:00:00 2001 From: thalassemia <67928790+thalassemia@users.noreply.github.com> Date: Thu, 18 Jan 2024 17:10:52 -0800 Subject: [PATCH 03/11] Use Numba JIT for tRNA charging IVP --- ecoli/library/initial_conditions.py | 2 +- ecoli/models/polypeptide_elongation_models.py | 64 ++++++--- .../ecoli/dataclasses/process/metabolism.py | 132 ++++++++++++------ 3 files changed, 134 insertions(+), 64 deletions(-) diff --git a/ecoli/library/initial_conditions.py b/ecoli/library/initial_conditions.py index 1e8fc65e8..61bdb17b3 100644 --- a/ecoli/library/initial_conditions.py +++ b/ecoli/library/initial_conditions.py @@ -1629,7 +1629,7 @@ def initialize_trna_charging(bulk_state, unique_molecules, sim_data, variable_el 'unit_conversion': metabolism.get_amino_acid_conc_conversion(MICROMOLAR_UNITS) } fraction_charged, *_ = calculate_trna_charging(synthetase_conc, uncharged_trna_conc, - charged_trna_conc, aa_conc, ribosome_conc, f, charging_params) + charged_trna_conc, aa_conc, ribosome_conc[0], f, charging_params) # Update counts of tRNA to match charging total_trna_counts = uncharged_trna + charged_trna diff --git a/ecoli/models/polypeptide_elongation_models.py b/ecoli/models/polypeptide_elongation_models.py index 9e32b816a..09487049e 100644 --- a/ecoli/models/polypeptide_elongation_models.py +++ b/ecoli/models/polypeptide_elongation_models.py @@ -2,6 +2,7 @@ from typing import Callable, Optional, Tuple +from numba import njit import numpy as np from scipy.integrate import solve_ivp from vivarium.library.units import units as vivunits @@ -565,7 +566,8 @@ def isTimeStepShortEnough(self, inputTimeStep, timeStepSafetyFraction): short_enough = False return short_enough - + + def ppgpp_metabolite_changes(uncharged_trna_conc, charged_trna_conc, ribosome_conc, f, rela_conc, spot_conc, ppgpp_conc, counts_to_molar, v_rib, charging_params, ppgpp_params, time_step, @@ -700,6 +702,7 @@ def ppgpp_metabolite_changes(uncharged_trna_conc, charged_trna_conc, return delta_metabolites, n_syn_reactions, n_deg_reactions, v_rela_syn, v_spot_syn, v_deg, v_deg_inhibited + def calculate_trna_charging(synthetase_conc, uncharged_trna_conc, charged_trna_conc, aa_conc, ribosome_conc, f, params, supply=None, time_limit=1000, limit_v_rib=False, use_disabled_aas=False): ''' @@ -762,34 +765,17 @@ def dcdt(t, c): ndarray[float]: dc/dt for tRNA concentrations dims: 2 * number of amino acids (uncharged tRNA come first, then charged) ''' + v_charging, dtrna, daa = dcdt_jit(t, c, n_aas_masked, n_aas, mask, + params['kS'], synthetase_conc, params['KMaa'], params['KMtf'], + f, params['krta'], params['krtf'], params['max_elong_rate'], + ribosome_conc, limit_v_rib, aa_rate_limit, v_rib_max) - uncharged_trna_conc = c[:n_aas_masked] - charged_trna_conc = c[n_aas_masked:2*n_aas_masked] - aa_conc = c[2*n_aas_masked:2*n_aas_masked+n_aas] - masked_aa_conc = aa_conc[mask] - - v_charging = (params['kS'] * synthetase_conc * uncharged_trna_conc * masked_aa_conc / (params['KMaa'][mask] * params['KMtf'][mask]) - / (1 + uncharged_trna_conc/params['KMtf'][mask] + masked_aa_conc/params['KMaa'][mask] + uncharged_trna_conc*masked_aa_conc/params['KMtf'][mask]/params['KMaa'][mask])) - with np.errstate(divide='ignore'): - numerator_ribosome = 1 + np.sum(f * (params['krta'] / charged_trna_conc + uncharged_trna_conc / charged_trna_conc * params['krta'] / params['krtf'])) - v_rib = params['max_elong_rate'] * ribosome_conc / numerator_ribosome - - # Handle case when f is 0 and charged_trna_conc is 0 - if not np.isfinite(v_rib): - v_rib = 0 - - # Limit v_rib and v_charging to the amount of available amino acids - if limit_v_rib: - v_charging = np.fmin(v_charging, aa_rate_limit) - v_rib = min(v_rib, v_rib_max) - - dtrna = v_charging - v_rib*f - daa = np.zeros(n_aas) if supply is None: v_synthesis = np.zeros(n_aas) v_import = np.zeros(n_aas) v_export = np.zeros(n_aas) else: + aa_conc = c[2*n_aas_masked:2*n_aas_masked+n_aas] v_synthesis, v_import, v_export = supply(unit_conversion * aa_conc) v_supply = v_synthesis + v_import - v_export daa[mask] = v_supply[mask] - v_charging @@ -854,6 +840,38 @@ def dcdt(t, c): return new_fraction_charged, v_rib, total_synthesis, total_import, total_export + +@njit(error_model='numpy') +def dcdt_jit(t, c, n_aas_masked, n_aas, mask, + kS, synthetase_conc, KMaa, KMtf, + f, krta, krtf, max_elong_rate, + ribosome_conc, limit_v_rib, aa_rate_limit, v_rib_max +): + uncharged_trna_conc = c[:n_aas_masked] + charged_trna_conc = c[n_aas_masked:2*n_aas_masked] + aa_conc = c[2*n_aas_masked:2*n_aas_masked+n_aas] + masked_aa_conc = aa_conc[mask] + + v_charging = (kS * synthetase_conc * uncharged_trna_conc * masked_aa_conc / (KMaa[mask] * KMtf[mask]) + / (1 + uncharged_trna_conc/KMtf[mask] + masked_aa_conc/KMaa[mask] + uncharged_trna_conc*masked_aa_conc/KMtf[mask]/KMaa[mask])) + numerator_ribosome = 1 + np.sum(f * (krta / charged_trna_conc + uncharged_trna_conc / charged_trna_conc * krta / krtf)) + v_rib = max_elong_rate * ribosome_conc / numerator_ribosome + + # Handle case when f is 0 and charged_trna_conc is 0 + if not np.isfinite(v_rib): + v_rib = 0 + + # Limit v_rib and v_charging to the amount of available amino acids + if limit_v_rib: + v_charging = np.fmin(v_charging, aa_rate_limit) + v_rib = min(v_rib, v_rib_max) + + dtrna = v_charging - v_rib*f + daa = np.zeros(n_aas) + + return v_charging, dtrna, daa + + def get_charging_supply_function( supply_in_charging: bool, mechanistic_supply: bool, diff --git a/reconstruction/ecoli/dataclasses/process/metabolism.py b/reconstruction/ecoli/dataclasses/process/metabolism.py index 4b747cdf2..77fead761 100755 --- a/reconstruction/ecoli/dataclasses/process/metabolism.py +++ b/reconstruction/ecoli/dataclasses/process/metabolism.py @@ -13,6 +13,7 @@ import re from typing import Any, cast, Dict, Iterable, List, Optional, Set, Tuple, Union +from numba import njit import numpy as np import sympy as sp from sympy.parsing.sympy_parser import parse_expr @@ -1317,23 +1318,11 @@ def amino_acid_synthesis(self, counts_per_aa_fwd: np.ndarray, counts_per_aa_rev: # Convert to appropriate arrays if units.hasUnit(aa_conc): aa_conc = aa_conc.asNumber(METABOLITE_CONCENTRATION_UNITS) - - km_saturation = np.product(1 / (1 + self.aa_upstream_kms / aa_conc), axis=1) - - # Determine saturation fraction for reactions - forward_fraction = 1 / (1 + aa_conc / self.aa_kis) * km_saturation - reverse_fraction = 1 / (1 + self.aa_reverse_kms / aa_conc) - deg_fraction = 1 / (1 + self.aa_degradation_kms / aa_conc) - loss_fraction = reverse_fraction + deg_fraction - - # Calculate synthesis rate - synthesis = ( - self.aa_forward_stoich @ (self.aa_kcats_fwd * counts_per_aa_fwd * forward_fraction) - - self.aa_reverse_stoich @ (self.aa_kcats_rev * counts_per_aa_rev * reverse_fraction) - - self.aa_kcats_rev * counts_per_aa_rev * deg_fraction - ) - - return synthesis, forward_fraction, loss_fraction + + return amino_acid_synthesis_jit(counts_per_aa_fwd, counts_per_aa_rev, aa_conc, + self.aa_upstream_kms, self.aa_kis, self.aa_reverse_kms, + self.aa_degradation_kms, self.aa_forward_stoich, self.aa_kcats_fwd, + self.aa_reverse_stoich, self.aa_kcats_rev) def amino_acid_export(self, aa_transporters_counts: np.ndarray, aa_conc: Union[units.Unum, np.ndarray], mechanistic_uptake: bool): @@ -1349,20 +1338,12 @@ def amino_acid_export(self, aa_transporters_counts: np.ndarray, rate of export for each amino acid. array is unitless but represents counts of amino acid per second """ - - if mechanistic_uptake: - # Export based on mechanistic model - if units.hasUnit(aa_conc): - aa_conc = aa_conc.asNumber(METABOLITE_CONCENTRATION_UNITS) - trans_counts_per_aa = self.aa_to_exporters_matrix @ aa_transporters_counts - export_rates = self.export_kcats_per_aa * trans_counts_per_aa / (1 + self.aa_export_kms / aa_conc) - else: - # Export is lumped with specific uptake rates in amino_acid_import - # and not dependent on internal amino acid concentrations or - # explicitly considered here - export_rates = np.zeros(len(aa_conc)) - - return export_rates + if units.hasUnit(aa_conc): + aa_conc = aa_conc.asNumber(METABOLITE_CONCENTRATION_UNITS) + + return amino_acid_export_jit(aa_transporters_counts, aa_conc, + mechanistic_uptake, self.aa_to_exporters_matrix, + self.export_kcats_per_aa, self.aa_export_kms) def amino_acid_import(self, aa_in_media: np.ndarray, dry_mass: units.Unum, internal_aa_conc: Union[units.Unum, np.ndarray], aa_transporters_counts: np.ndarray, @@ -1385,15 +1366,11 @@ def amino_acid_import(self, aa_in_media: np.ndarray, dry_mass: units.Unum, if units.hasUnit(internal_aa_conc): internal_aa_conc = internal_aa_conc.asNumber(METABOLITE_CONCENTRATION_UNITS) - saturation = 1 / (1 + internal_aa_conc / self.aa_import_kis) - if mechanisitic_uptake: - # Uptake based on mechanistic model - counts_per_aa = self.aa_to_importers_matrix @ aa_transporters_counts - import_rates = self.import_kcats_per_aa * counts_per_aa - else: - import_rates = self.max_specific_import_rates * dry_mass.asNumber(DRY_MASS_UNITS) - - return import_rates * saturation * aa_in_media + dry_mass = dry_mass.asNumber(DRY_MASS_UNITS) + return amino_acid_import_jit(aa_in_media, dry_mass, internal_aa_conc, + aa_transporters_counts, mechanisitic_uptake, self.aa_import_kis, + self.aa_to_importers_matrix, self.import_kcats_per_aa, + self.max_specific_import_rates) def get_amino_acid_conc_conversion(self, conc_units): return units.strip_empty_units(conc_units / METABOLITE_CONCENTRATION_UNITS) @@ -2224,6 +2201,81 @@ def _is_transport_rxn(self, stoich): return is_transport_rxn +@njit +def np_apply_along_axis(func1d, axis, arr): + assert arr.ndim == 2 + assert axis in [0, 1] + if axis == 0: + result = np.empty(arr.shape[1]) + for i in range(len(result)): + result[i] = func1d(arr[:, i]) + else: + result = np.empty(arr.shape[0]) + for i in range(len(result)): + result[i] = func1d(arr[i, :]) + return result + + +@njit +def np_prod(array, axis): + return np_apply_along_axis(np.prod, axis, array) + + +@njit +def amino_acid_synthesis_jit(counts_per_aa_fwd, counts_per_aa_rev, aa_conc, + aa_upstream_kms, aa_kis, aa_reverse_kms, + aa_degradation_kms, aa_forward_stoich, aa_kcats_fwd, + aa_reverse_stoich, aa_kcats_rev): + km_saturation = np_prod(1 / (1 + aa_upstream_kms / aa_conc), axis=1) + + # Determine saturation fraction for reactions + forward_fraction = 1 / (1 + aa_conc / aa_kis) * km_saturation + reverse_fraction = 1 / (1 + aa_reverse_kms / aa_conc) + deg_fraction = 1 / (1 + aa_degradation_kms / aa_conc) + loss_fraction = reverse_fraction + deg_fraction + + # Calculate synthesis rate + synthesis = ( + aa_forward_stoich.astype(np.float64) @ (aa_kcats_fwd * counts_per_aa_fwd * forward_fraction).astype(np.float64) + - aa_reverse_stoich.astype(np.float64) @ (aa_kcats_rev * counts_per_aa_rev * reverse_fraction).astype(np.float64) + - aa_kcats_rev * counts_per_aa_rev * deg_fraction + ) + + return synthesis, forward_fraction, loss_fraction + + +@njit +def amino_acid_export_jit(aa_transporters_counts, aa_conc, + mechanistic_uptake, aa_to_exporters_matrix, + export_kcats_per_aa, aa_export_kms): + if mechanistic_uptake: + # Export based on mechanistic model + trans_counts_per_aa = aa_to_exporters_matrix.astype(np.float64) @ aa_transporters_counts.astype(np.float64) + export_rates = export_kcats_per_aa * trans_counts_per_aa / (1 + aa_export_kms / aa_conc) + else: + # Export is lumped with specific uptake rates in amino_acid_import + # and not dependent on internal amino acid concentrations or + # explicitly considered here + export_rates = np.zeros(len(aa_conc)) + return export_rates + + +@njit +def amino_acid_import_jit(aa_in_media, dry_mass, internal_aa_conc, + aa_transporters_counts, mechanisitic_uptake, aa_import_kis, + aa_to_importers_matrix, import_kcats_per_aa, + max_specific_import_rates): + saturation = 1 / (1 + internal_aa_conc / aa_import_kis) + if mechanisitic_uptake: + # Uptake based on mechanistic model + counts_per_aa = aa_to_importers_matrix.astype(np.float64) @ aa_transporters_counts.astype(np.float64) + import_rates = import_kcats_per_aa * counts_per_aa + else: + import_rates = max_specific_import_rates * dry_mass + + return import_rates * saturation * aa_in_media + + # Class used to update metabolite concentrations based on the current nutrient conditions class ConcentrationUpdates(object): def __init__(self, concDict, relative_changes, equilibriumReactions, exchange_data_dict, all_metabolite_ids, linked_metabolites): From 5fbabbe7687a93f21dd62d9155d9b610bb394430 Mon Sep 17 00:00:00 2001 From: thalassemia <67928790+thalassemia@users.noreply.github.com> Date: Thu, 18 Jan 2024 17:11:02 -0800 Subject: [PATCH 04/11] Update all requirements --- requirements.txt | 274 +++++++++++++++++++++++++---------------------- 1 file changed, 143 insertions(+), 131 deletions(-) diff --git a/requirements.txt b/requirements.txt index e245161a3..52d5bdba6 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,193 +1,205 @@ -absl-py==1.4.0 +absl-py==2.1.0 access==1.1.9 -aesara==2.9.0 +aesara==2.9.3 affine==2.4.0 -anyio==3.7.0 -argon2-cffi==21.3.0 +anyio==4.2.0 +argon2-cffi==23.1.0 argon2-cffi-bindings==21.2.0 -arrow==1.2.3 -asttokens==2.2.1 -attrs==23.1.0 +arrow==1.3.0 +asttokens==2.4.1 +async-lru==2.0.4 +attrs==23.2.0 +Babel==2.14.0 backcall==0.2.0 -beautifulsoup4==4.12.2 -biopython==1.81 -bleach==6.0.0 -certifi==2023.5.7 -cffi==1.15.1 -charset-normalizer==3.1.0 -click==8.1.3 +beautifulsoup4==4.12.3 +biopython==1.83 +bleach==6.1.0 +certifi==2023.11.17 +cffi==1.16.0 +charset-normalizer==3.3.2 +clarabel==0.6.0 +click==8.1.7 click-plugins==1.1.1 cligj==0.7.2 -comm==0.1.3 +comm==0.2.1 cons==0.4.6 -contourpy==1.1.0 -coverage==7.2.7 +contourpy==1.2.0 +coverage==7.4.0 cvxpy==1.4.1 -cycler==0.11.0 -Cython==0.29.35 -debugpy==1.6.7 +cycler==0.12.1 +Cython==3.0.8 +debugpy==1.8.0 decorator==5.1.1 defusedxml==0.7.1 deprecation==2.1.0 -dill==0.3.6 -dnspython==2.3.0 +dill==0.3.7 +dnspython==2.4.2 ecos==2.0.12 -esda==2.4.3 +esda==2.5.1 ete3==3.1.3 etuples==0.3.9 -executing==1.2.0 -fastjsonschema==2.17.1 -filelock==3.12.2 -Fiona==1.9.4.post1 -fonttools==4.40.0 +executing==2.0.1 +fastjsonschema==2.19.1 +filelock==3.13.1 +fiona==1.9.5 +fonttools==4.47.2 fqdn==1.5.1 -geopandas==0.13.2 -giddy==2.3.4 -idna==3.4 -imageio==2.31.1 -inequality==1.0.0 +geopandas==0.14.2 +giddy==2.3.5 +idna==3.6 +imageio==2.33.1 +inequality==1.0.1 iniconfig==2.0.0 ipdb==0.13.13 -ipykernel==6.23.3 -ipython==8.14.0 +ipykernel==6.29.0 +ipython==8.20.0 ipython-genutils==0.2.0 -ipywidgets==8.0.6 +ipywidgets==8.1.1 isoduration==20.11.0 iteround==1.0.4 -jedi==0.18.2 -Jinja2==3.1.2 -joblib==1.2.0 +jedi==0.19.1 +Jinja2==3.1.3 +joblib==1.3.2 +json5==0.9.14 jsonpointer==2.4 -jsonschema==4.17.3 +jsonschema==4.21.0 +jsonschema-specifications==2023.12.1 jupyter==1.0.0 jupyter-console==6.6.3 -jupyter-events==0.6.3 -jupyter_client==8.3.0 -jupyter_core==5.3.1 -jupyter_server==2.7.0 -jupyter_server_terminals==0.4.4 -jupyterlab-pygments==0.2.2 -jupyterlab-widgets==3.0.7 -kiwisolver==1.4.4 -lazy_loader==0.2 -libpysal==4.7.0 -llvmlite==0.40.1 +jupyter-events==0.9.0 +jupyter-lsp==2.2.2 +jupyter_client==8.6.0 +jupyter_core==5.7.1 +jupyter_server==2.12.5 +jupyter_server_terminals==0.5.1 +jupyterlab==4.0.10 +jupyterlab-widgets==3.0.9 +jupyterlab_pygments==0.3.0 +jupyterlab_server==2.25.2 +kiwisolver==1.4.5 +lazy_loader==0.3 +libpysal==4.9.2 +llvmlite==0.41.1 logical-unification==0.4.6 -mapclassify==2.5.0 +mapclassify==2.6.1 MarkupSafe==2.1.3 -matplotlib==3.7.1 +matplotlib==3.8.2 matplotlib-inline==0.1.6 -mgwr==2.1.2 +mgwr==2.2.1 miniKanren==1.0.3 -mistune==3.0.1 -momepy==0.6.0 +mistune==3.0.2 +momepy==0.7.0 mpmath==1.3.0 multipledispatch==1.0.0 nbclassic==1.0.0 -nbclient==0.8.0 -nbconvert==7.6.0 -nbformat==5.9.0 -nest-asyncio==1.5.6 -networkx==3.1 -notebook==6.5.4 +nbclient==0.9.0 +nbconvert==7.14.2 +nbformat==5.9.2 +nest-asyncio==1.5.9 +networkx==3.2.1 +notebook==7.0.6 notebook_shim==0.2.3 -numba==0.57.1 -numpy==1.24.4 -opencv-python==4.8.0.76 -orjson==3.9.1 -ortools==9.7.2996 +numba==0.58.1 +numpy==1.26.3 +opencv-python==4.9.0.80 +orjson==3.9.12 +ortools==9.8.3296 osqp==0.6.3 -overrides==7.3.1 -packaging==23.1 -pandas==2.0.2 -pandocfilters==1.5.0 +overrides==7.4.0 +packaging==23.2 +pandas==2.1.4 +pandocfilters==1.5.1 parso==0.8.3 -patsy==0.5.3 -pexpect==4.8.0 +patsy==0.5.6 +pexpect==4.9.0 pickleshare==0.7.5 -Pillow==9.5.0 -Pint==0.22 -platformdirs==3.8.0 -pluggy==1.2.0 -pointpats==2.3.0 -prometheus-client==0.17.0 -prompt-toolkit==3.0.38 -protobuf==4.23.3 -psutil==5.9.5 +pillow==10.2.0 +Pint==0.23 +platformdirs==4.1.0 +pluggy==1.3.0 +pointpats==2.4.0 +prometheus-client==0.19.0 +prompt-toolkit==3.0.43 +protobuf==4.25.2 +psutil==5.9.7 ptyprocess==0.7.0 -PuLP==2.7.0 +PuLP==2.8.0 pure-eval==0.2.2 +pybind11==2.11.1 pycparser==2.21 -Pygments==2.15.1 -pymongo==4.4.0 -pymunk==6.5.1 -pyparsing==3.1.0 -pyproj==3.6.0 -PyQt5==5.15.9 +Pygments==2.17.2 +pymongo==4.6.1 +pymunk==6.6.0 +pyparsing==3.1.1 +pyproj==3.6.1 +PyQt5==5.15.10 PyQt5-Qt5==5.15.2 -PyQt5-sip==12.12.1 -pyrsistent==0.19.3 -pysal==23.1 -pytest==7.4.0 +PyQt5-sip==12.13.0 +pyrsistent==0.20.0 +pysal==23.7 +pytest==7.4.4 pytest-cov==4.1.0 python-dateutil==2.8.2 python-json-logger==2.0.7 -pytz==2023.3 -PyWavelets==1.4.1 -PyYAML==6.0 -pyzmq==25.1.0 -qdldl==0.1.7 -qtconsole==5.4.3 -QtPy==2.3.1 +pytz==2023.3.post1 +PyWavelets==1.5.0 +PyYAML==6.0.1 +pyzmq==25.1.2 +qdldl==0.1.7.post0 +qtconsole==5.5.1 +QtPy==2.4.1 quantecon==0.7.1 -rasterio==1.3.7 +rasterio==1.3.9 rasterstats==0.19.0 +referencing==0.32.1 requests==2.31.0 rfc3339-validator==0.1.4 rfc3986-validator==0.1.1 -Rtree==1.0.1 -scikit-image==0.21.0 -scikit-learn==1.2.2 -scipy==1.11.0 -scs==3.2.3 -seaborn==0.12.2 -segregation==2.4.2 +rpds-py==0.17.1 +Rtree==1.1.0 +scikit-image==0.22.0 +scikit-learn==1.4.0 +scipy==1.11.4 +scs==3.2.4.post1 +seaborn==0.13.1 +segregation==2.5 Send2Trash==1.8.2 -shapely==2.0.1 -simplejson==3.19.1 +shapely==2.0.2 +simplejson==3.19.2 six==1.16.0 sniffio==1.3.0 snuggs==1.4.7 -soupsieve==2.4.1 +soupsieve==2.5 spaghetti==1.7.4 -spglm==1.0.8 +spglm==1.1.0 spint==1.0.7 splot==1.1.5.post1 -spopt==0.5.0 -spreg==1.3.2 +spopt==0.6.0 +spreg==1.4.2 spvcm==0.3.0 -stack-data==0.6.2 -statsmodels==0.14.0 +stack-data==0.6.3 +statsmodels==0.14.1 stochastic-arrow==1.0.0 -swiglpk==5.0.8 +swiglpk==5.0.10 sympy==1.12 -terminado==0.17.1 -threadpoolctl==3.1.0 -tifffile==2023.4.12 +terminado==0.18.0 +threadpoolctl==3.2.0 +tifffile==2023.12.9 tinycss2==1.2.1 -tobler==0.10 +tobler==0.11.2 toolz==0.12.0 -tornado==6.3.2 -tqdm==4.65.0 -traitlets==5.9.0 -typing_extensions==4.6.3 -tzdata==2023.3 +tornado==6.4 +tqdm==4.66.1 +traitlets==5.14.1 +types-python-dateutil==2.8.19.20240106 +typing_extensions==4.9.0 +tzdata==2023.4 Unum==4.2.1 uri-template==1.3.0 -urllib3==2.0.3 -vivarium-core==1.6.0 -wcwidth==0.2.6 +urllib3==2.1.0 +vivarium-core==1.6.1 +wcwidth==0.2.13 webcolors==1.13 webencodings==0.5.1 -websocket-client==1.6.1 -widgetsnbextension==4.0.7 +websocket-client==1.7.0 +widgetsnbextension==4.0.9 From b9f75a740acea3b791f2c49631ab0c99459a453d Mon Sep 17 00:00:00 2001 From: thalassemia <67928790+thalassemia@users.noreply.github.com> Date: Thu, 18 Jan 2024 17:15:20 -0800 Subject: [PATCH 05/11] Do not allocate dense array for FBA --- wholecell/utils/_netflow/nf_glpk.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/wholecell/utils/_netflow/nf_glpk.py b/wholecell/utils/_netflow/nf_glpk.py index fb6b746df..8af3db53a 100644 --- a/wholecell/utils/_netflow/nf_glpk.py +++ b/wholecell/utils/_netflow/nf_glpk.py @@ -403,15 +403,19 @@ def buildEqConst(self): n_flows = len(self._flows) self._add_rows(n_coeffs) - A = np.zeros((n_coeffs, n_flows)) + A_row = [] + A_col = [] + A_data = [] # avoid creating duplicate constraints self._materialIdxLookup = {} for materialIdx, (material, pairs) in enumerate(sorted(self._materialCoeffs.items())): self._materialIdxLookup[material] = materialIdx for pair in pairs: - A[materialIdx, pair[1]] = pair[0] + A_row.append(materialIdx) + A_col.append(pair[1]) + A_data.append(pair[0]) - A_coo = coo_matrix(A) + A_coo = coo_matrix((A_data, (A_row, A_col)), shape=(n_coeffs, n_flows)) rowIdxs = _toIndexArray(A_coo.row) colIdxs = _toIndexArray(A_coo.col) data = _toDoubleArray(A_coo.data) From a88478fa691294e891d3f385268155af2b9f57e8 Mon Sep 17 00:00:00 2001 From: thalassemia <67928790+thalassemia@users.noreply.github.com> Date: Thu, 18 Jan 2024 17:16:54 -0800 Subject: [PATCH 06/11] Ribosome conc is scalar --- ecoli/library/initial_conditions.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ecoli/library/initial_conditions.py b/ecoli/library/initial_conditions.py index 61bdb17b3..1e8fc65e8 100644 --- a/ecoli/library/initial_conditions.py +++ b/ecoli/library/initial_conditions.py @@ -1629,7 +1629,7 @@ def initialize_trna_charging(bulk_state, unique_molecules, sim_data, variable_el 'unit_conversion': metabolism.get_amino_acid_conc_conversion(MICROMOLAR_UNITS) } fraction_charged, *_ = calculate_trna_charging(synthetase_conc, uncharged_trna_conc, - charged_trna_conc, aa_conc, ribosome_conc[0], f, charging_params) + charged_trna_conc, aa_conc, ribosome_conc, f, charging_params) # Update counts of tRNA to match charging total_trna_counts = uncharged_trna + charged_trna From 409e5f88052276f522b6c21c24dec00ee1af7081 Mon Sep 17 00:00:00 2001 From: thalassemia <67928790+thalassemia@users.noreply.github.com> Date: Thu, 18 Jan 2024 17:26:00 -0800 Subject: [PATCH 07/11] Downgrade ortools for cvxpy --- requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index 52d5bdba6..fb6314de9 100644 --- a/requirements.txt +++ b/requirements.txt @@ -103,7 +103,7 @@ numba==0.58.1 numpy==1.26.3 opencv-python==4.9.0.80 orjson==3.9.12 -ortools==9.8.3296 +ortools==9.7.2996 osqp==0.6.3 overrides==7.4.0 packaging==23.2 From d15dcc17803a92879c7d979fd4ae264c634b46a5 Mon Sep 17 00:00:00 2001 From: thalassemia <67928790+thalassemia@users.noreply.github.com> Date: Thu, 18 Jan 2024 17:31:22 -0800 Subject: [PATCH 08/11] Cleanup whitespace --- ecoli/models/polypeptide_elongation_models.py | 4 +-- .../ecoli/dataclasses/process/metabolism.py | 32 +++++++++---------- 2 files changed, 18 insertions(+), 18 deletions(-) diff --git a/ecoli/models/polypeptide_elongation_models.py b/ecoli/models/polypeptide_elongation_models.py index 09487049e..d9bd973cf 100644 --- a/ecoli/models/polypeptide_elongation_models.py +++ b/ecoli/models/polypeptide_elongation_models.py @@ -765,7 +765,7 @@ def dcdt(t, c): ndarray[float]: dc/dt for tRNA concentrations dims: 2 * number of amino acids (uncharged tRNA come first, then charged) ''' - v_charging, dtrna, daa = dcdt_jit(t, c, n_aas_masked, n_aas, mask, + v_charging, dtrna, daa = dcdt_jit(t, c, n_aas_masked, n_aas, mask, params['kS'], synthetase_conc, params['KMaa'], params['KMtf'], f, params['krta'], params['krtf'], params['max_elong_rate'], ribosome_conc, limit_v_rib, aa_rate_limit, v_rib_max) @@ -842,7 +842,7 @@ def dcdt(t, c): @njit(error_model='numpy') -def dcdt_jit(t, c, n_aas_masked, n_aas, mask, +def dcdt_jit(t, c, n_aas_masked, n_aas, mask, kS, synthetase_conc, KMaa, KMtf, f, krta, krtf, max_elong_rate, ribosome_conc, limit_v_rib, aa_rate_limit, v_rib_max diff --git a/reconstruction/ecoli/dataclasses/process/metabolism.py b/reconstruction/ecoli/dataclasses/process/metabolism.py index 77fead761..adcc9a95c 100755 --- a/reconstruction/ecoli/dataclasses/process/metabolism.py +++ b/reconstruction/ecoli/dataclasses/process/metabolism.py @@ -1318,9 +1318,9 @@ def amino_acid_synthesis(self, counts_per_aa_fwd: np.ndarray, counts_per_aa_rev: # Convert to appropriate arrays if units.hasUnit(aa_conc): aa_conc = aa_conc.asNumber(METABOLITE_CONCENTRATION_UNITS) - - return amino_acid_synthesis_jit(counts_per_aa_fwd, counts_per_aa_rev, aa_conc, - self.aa_upstream_kms, self.aa_kis, self.aa_reverse_kms, + + return amino_acid_synthesis_jit(counts_per_aa_fwd, counts_per_aa_rev, aa_conc, + self.aa_upstream_kms, self.aa_kis, self.aa_reverse_kms, self.aa_degradation_kms, self.aa_forward_stoich, self.aa_kcats_fwd, self.aa_reverse_stoich, self.aa_kcats_rev) @@ -1340,9 +1340,9 @@ def amino_acid_export(self, aa_transporters_counts: np.ndarray, """ if units.hasUnit(aa_conc): aa_conc = aa_conc.asNumber(METABOLITE_CONCENTRATION_UNITS) - - return amino_acid_export_jit(aa_transporters_counts, aa_conc, - mechanistic_uptake, self.aa_to_exporters_matrix, + + return amino_acid_export_jit(aa_transporters_counts, aa_conc, + mechanistic_uptake, self.aa_to_exporters_matrix, self.export_kcats_per_aa, self.aa_export_kms) def amino_acid_import(self, aa_in_media: np.ndarray, dry_mass: units.Unum, @@ -1367,9 +1367,9 @@ def amino_acid_import(self, aa_in_media: np.ndarray, dry_mass: units.Unum, internal_aa_conc = internal_aa_conc.asNumber(METABOLITE_CONCENTRATION_UNITS) dry_mass = dry_mass.asNumber(DRY_MASS_UNITS) - return amino_acid_import_jit(aa_in_media, dry_mass, internal_aa_conc, - aa_transporters_counts, mechanisitic_uptake, self.aa_import_kis, - self.aa_to_importers_matrix, self.import_kcats_per_aa, + return amino_acid_import_jit(aa_in_media, dry_mass, internal_aa_conc, + aa_transporters_counts, mechanisitic_uptake, self.aa_import_kis, + self.aa_to_importers_matrix, self.import_kcats_per_aa, self.max_specific_import_rates) def get_amino_acid_conc_conversion(self, conc_units): @@ -2222,8 +2222,8 @@ def np_prod(array, axis): @njit -def amino_acid_synthesis_jit(counts_per_aa_fwd, counts_per_aa_rev, aa_conc, - aa_upstream_kms, aa_kis, aa_reverse_kms, +def amino_acid_synthesis_jit(counts_per_aa_fwd, counts_per_aa_rev, aa_conc, + aa_upstream_kms, aa_kis, aa_reverse_kms, aa_degradation_kms, aa_forward_stoich, aa_kcats_fwd, aa_reverse_stoich, aa_kcats_rev): km_saturation = np_prod(1 / (1 + aa_upstream_kms / aa_conc), axis=1) @@ -2245,8 +2245,8 @@ def amino_acid_synthesis_jit(counts_per_aa_fwd, counts_per_aa_rev, aa_conc, @njit -def amino_acid_export_jit(aa_transporters_counts, aa_conc, - mechanistic_uptake, aa_to_exporters_matrix, +def amino_acid_export_jit(aa_transporters_counts, aa_conc, + mechanistic_uptake, aa_to_exporters_matrix, export_kcats_per_aa, aa_export_kms): if mechanistic_uptake: # Export based on mechanistic model @@ -2261,9 +2261,9 @@ def amino_acid_export_jit(aa_transporters_counts, aa_conc, @njit -def amino_acid_import_jit(aa_in_media, dry_mass, internal_aa_conc, - aa_transporters_counts, mechanisitic_uptake, aa_import_kis, - aa_to_importers_matrix, import_kcats_per_aa, +def amino_acid_import_jit(aa_in_media, dry_mass, internal_aa_conc, + aa_transporters_counts, mechanisitic_uptake, aa_import_kis, + aa_to_importers_matrix, import_kcats_per_aa, max_specific_import_rates): saturation = 1 / (1 + internal_aa_conc / aa_import_kis) if mechanisitic_uptake: From d88c43e7257c07e7482d7eb2e4c0e8cd3fef5109 Mon Sep 17 00:00:00 2001 From: thalassemia Date: Thu, 18 Jan 2024 18:10:45 -0800 Subject: [PATCH 09/11] Upgrade PyQt5-Qt5 --- requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index fb6314de9..008b47325 100644 --- a/requirements.txt +++ b/requirements.txt @@ -133,7 +133,7 @@ pymunk==6.6.0 pyparsing==3.1.1 pyproj==3.6.1 PyQt5==5.15.10 -PyQt5-Qt5==5.15.2 +PyQt5-Qt5==5.15.12 PyQt5-sip==12.13.0 pyrsistent==0.20.0 pysal==23.7 From 509760f82710d610169fb8c13a9aa7e2995648fb Mon Sep 17 00:00:00 2001 From: thalassemia Date: Fri, 19 Jan 2024 11:39:43 -0800 Subject: [PATCH 10/11] PyQt5-Qt5 5.15.12 only on macOS --- requirements.txt | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index 008b47325..8a8983f36 100644 --- a/requirements.txt +++ b/requirements.txt @@ -133,7 +133,8 @@ pymunk==6.6.0 pyparsing==3.1.1 pyproj==3.6.1 PyQt5==5.15.10 -PyQt5-Qt5==5.15.12 +PyQt5-Qt5==5.15.12; sys_platform == 'darwin' +PyQt5-Qt5==5.15.2; sys_platform != 'darwin' PyQt5-sip==12.13.0 pyrsistent==0.20.0 pysal==23.7 From e0a99b10ee17e2505369088e05896a31aa2a681d Mon Sep 17 00:00:00 2001 From: thalassemia <67928790+thalassemia@users.noreply.github.com> Date: Fri, 19 Jan 2024 11:59:40 -0800 Subject: [PATCH 11/11] Remove asserts in JIT functions (breaks with pytest) --- .../ecoli/dataclasses/process/metabolism.py | 26 ++++++++++--------- 1 file changed, 14 insertions(+), 12 deletions(-) diff --git a/reconstruction/ecoli/dataclasses/process/metabolism.py b/reconstruction/ecoli/dataclasses/process/metabolism.py index adcc9a95c..77880286d 100755 --- a/reconstruction/ecoli/dataclasses/process/metabolism.py +++ b/reconstruction/ecoli/dataclasses/process/metabolism.py @@ -2203,22 +2203,24 @@ def _is_transport_rxn(self, stoich): @njit def np_apply_along_axis(func1d, axis, arr): - assert arr.ndim == 2 - assert axis in [0, 1] - if axis == 0: - result = np.empty(arr.shape[1]) - for i in range(len(result)): - result[i] = func1d(arr[:, i]) - else: - result = np.empty(arr.shape[0]) - for i in range(len(result)): - result[i] = func1d(arr[i, :]) - return result + if arr.ndim != 2: + raise RuntimeError("Array must have 2 dimensions.") + if axis not in [0, 1]: + raise RuntimeError("Axis must be 0 or 1.") + if axis == 0: + result = np.empty(arr.shape[1]) + for i in range(len(result)): + result[i] = func1d(arr[:, i]) + else: + result = np.empty(arr.shape[0]) + for i in range(len(result)): + result[i] = func1d(arr[i, :]) + return result @njit def np_prod(array, axis): - return np_apply_along_axis(np.prod, axis, array) + return np_apply_along_axis(np.prod, axis, array) @njit