From b71fc0d7aec2ec272d0b0208a347e9332af0b771 Mon Sep 17 00:00:00 2001 From: Alon Grinberg Dana Date: Tue, 28 Nov 2023 09:20:50 +0200 Subject: [PATCH 1/6] When generating a flux d, Don't crash if the model cannot be simulated e.g., if the concentrations were fed incorrectly --- t3/utils/flux.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/t3/utils/flux.py b/t3/utils/flux.py index b84aedd0..1bff73c7 100644 --- a/t3/utils/flux.py +++ b/t3/utils/flux.py @@ -568,6 +568,10 @@ def create_digraph(flux_graph: dict, for rop_list in rxn_dict.values(): species_to_consider.update(rop_list[0]) xs = [v for k, v in profile['X'].items() if k in species_to_consider] + if not len(xs): + print(f'Could not create a flux diagram for observables {observables} at {time} s. ' + f'Could not simulate the system.') + return x_max, x_min = max(xs), min(xs) abs_rops = [abs(values[1]) for inner_dict in flux_graph.values() for values in inner_dict.values()] rop_min, rop_max = min(abs_rops), max(abs_rops) @@ -821,7 +825,7 @@ def get_flux_graph(profile: dict, graph[node][rxn][1] += rop else: graph[node][rxn] = [opposite_rxn_species, rop] - min_rop = min_rop * max_rop + min_rop = min_rop * max_rop if min_rop is not None else 0 return graph, nodes_to_explore, min_rop, max_rop From f86fa1c9a8d3018cb5fc04c24cffff260ec0143c Mon Sep 17 00:00:00 2001 From: Alon Grinberg Dana Date: Tue, 28 Nov 2023 09:21:36 +0200 Subject: [PATCH 2/6] Return width 1 as fallback in flux.get_width() --- t3/utils/flux.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/t3/utils/flux.py b/t3/utils/flux.py index 1bff73c7..34d932bc 100644 --- a/t3/utils/flux.py +++ b/t3/utils/flux.py @@ -694,6 +694,8 @@ def get_width(x: float, max_width, min_width = 4, 0.2 x, x_min, x_max = abs(x), abs(x_min), abs(x_max) if not log_scale: + if x == x_min == x_max: + return 1 return min_width + (x - x_min) * (max_width - min_width) / (x_max - x_min) return get_width(x=-np.log10(x_min) - np.log10(x_max) + np.log10(x), x_min=-np.log10(x_max), From 6e44261c6e64e88ad92361c911f6ef86cb79bd08 Mon Sep 17 00:00:00 2001 From: Alon Grinberg Dana Date: Tue, 28 Nov 2023 09:47:28 +0200 Subject: [PATCH 3/6] Added flux.count_species_in_well() To properly count species in cases where a species label is contained within another species label, e.g., 'NO' and 'NO2'. This is avoided when using RMG labels like 'S(1243)', but is not avoided when generating flux diagrams for non-RMG generated Cantera models --- t3/utils/flux.py | 25 ++++++++++++++++++++++++- 1 file changed, 24 insertions(+), 1 deletion(-) diff --git a/t3/utils/flux.py b/t3/utils/flux.py index 34d932bc..99981470 100644 --- a/t3/utils/flux.py +++ b/t3/utils/flux.py @@ -894,7 +894,7 @@ def get_opposite_rxn_species(rxn: str, spc: str) -> List[str]: """ arrow = ' <=> ' if ' <=> ' in rxn else ' => ' wells = rxn.split(arrow) - counts = wells[0].count(spc), wells[1].count(spc) + counts = count_species_in_well(well=wells[0], spc=spc), count_species_in_well(well=wells[1], spc=spc) i = int(counts[0] > counts[1]) species = wells[i].split(' + ') for token in [' + M', ' (+M)', 'M', ' + ']: @@ -902,6 +902,29 @@ def get_opposite_rxn_species(rxn: str, spc: str) -> List[str]: return [s for s in species if s != ''] +def count_species_in_well(well: str, + spc: str, + ) -> int: + """ + Count the number of times a species appears in a well. + + Args: + well (str): The well string. + spc (str): The species label. + + Returns: + int: The number of times a species appears in the well. + """ + count = 0 + for token in [' + M', ' (+M)', 'M']: + well = well.replace(token, '') + splits = well.split(' + ') + for s in splits: + if s == spc: + count += 1 + return count + + def unpack_stoichiometry(labels: List[str]) -> Tuple[List[str], List[int]]: """ Unpack stoichiometry. From 424311daf7447e8218aaaa391eb776f13e3426ff Mon Sep 17 00:00:00 2001 From: Alon Grinberg Dana Date: Thu, 30 Nov 2023 17:26:31 +0200 Subject: [PATCH 4/6] Tests: test_count_species_in_well --- tests/test_flux.py | 18 +++++++++++++++++- 1 file changed, 17 insertions(+), 1 deletion(-) diff --git a/tests/test_flux.py b/tests/test_flux.py index 023be4c0..acd0dfc3 100644 --- a/tests/test_flux.py +++ b/tests/test_flux.py @@ -361,7 +361,8 @@ def test_get_flux_graph(): assert len(profile['ROPs']) == 36 flux_graph, nodes_to_explore, min_rop, max_rop = flux.get_flux_graph(profile=profile, observables=['H4N2(1)']) if i == 0: - assert nodes_to_explore == {'H2(4)', 'H3N2(6)', 'ammonia(9)', 'H(3)', 'H2N2(7)', '2 NH2(5)', 'HN2(10)', 'N2(2)', 'NH2(5)'} + assert nodes_to_explore == {'H2(4)', '2 H3N2(6)', 'ammonia(9)', 'H(3)', 'H2N2(7)', '2 NH2(5)', + 'HN2(10)', 'N2(2)', 'NH2(5)', 'H3N2(6)'} assert almost_equal(min_rop, 3.27e-21, ratio=100) assert almost_equal(max_rop, 20.3659, places=3) assert list(flux_graph.keys()) == ['H4N2(1)', 'NH2(5)', 'HN2(10)', 'H(3)', 'H2N2(7)', 'H2(4)', 'ammonia(9)', 'H3N2(6)'] @@ -391,6 +392,21 @@ def test_get_opposite_rxn_species(): assert flux.get_opposite_rxn_species(rxn=rxn, spc='H(3)') == ['H2O2(18)'] assert flux.get_opposite_rxn_species(rxn=rxn, spc='H2O2(18)') == ['H(3)', 'HO2(6)'] + rxn = 'HO2 + NO <=> NO2 + OH' + assert flux.get_opposite_rxn_species(rxn=rxn, spc='HO2') == ['NO2', 'OH'] + assert flux.get_opposite_rxn_species(rxn=rxn, spc='NO') == ['NO2', 'OH'] + assert flux.get_opposite_rxn_species(rxn=rxn, spc='NO2') == ['HO2', 'NO'] + assert flux.get_opposite_rxn_species(rxn=rxn, spc='OH') == ['HO2', 'NO'] + + +def test_count_species_in_well(): + """Test counting a species in a well""" + assert flux.count_species_in_well(well='H2O + H + H', spc='H2O') == 1 + assert flux.count_species_in_well(well='H2O + H + H', spc='H') == 2 + assert flux.count_species_in_well(well='H2O + H + H', spc='H2') == 0 + assert flux.count_species_in_well(well='HO2 + NO', spc='NO') == 1 + assert flux.count_species_in_well(well='NO2 + OH', spc='NO') == 0 + def test_get_other_reactants_and_products(): """Test getting the reactants and products other than a given species.""" From a3fa93fa31fac96bf7e3d4ec6641ed8b739ce1e3 Mon Sep 17 00:00:00 2001 From: Alon Grinberg Dana Date: Wed, 29 Nov 2023 08:04:03 +0200 Subject: [PATCH 5/6] Added allowed_nodes to flux diagram To create edges only between a list of selected nodes if desired --- t3/utils/flux.py | 50 ++++++++++++++++++++++++++++++++++-------------- 1 file changed, 36 insertions(+), 14 deletions(-) diff --git a/t3/utils/flux.py b/t3/utils/flux.py index 99981470..e82b8e4a 100644 --- a/t3/utils/flux.py +++ b/t3/utils/flux.py @@ -34,6 +34,7 @@ def generate_flux(model_path: str, display_r_n_p: bool = True, scaling: Optional[float] = None, fix_cantera_model: bool = True, + allowed_nodes: Optional[List[str]] = None, ): """ Generate a flux diagram for a given model and composition. @@ -66,6 +67,8 @@ def generate_flux(model_path: str, display_r_n_p (bool, optional): Whether to display the other reactants and products on each arrow. scaling (Optional[float], optional): The scaling of the final image, 100 means no scaling. fix_cantera_model (bool, optional): Whether to fix the Cantera model before running the simulation. + allowed_nodes (Optional[List[str]], optional): A list of nodes to consider. + any node outside this list will not appear in the flux diagram. Structures: profiles: {