Merge pull request #28 from DenisTitovLab/prevent-param-snooping-in-q…

…ssa-equations-and-catch-wrong-param-range-entry

Prevent param snooping in qssa equations and catch wrong param range entry

src/data_driven_rate_equation_selection.jl

@@ -73,8 +73,14 @@ function data_driven_rate_equation_selection(
         all_param_removal_codes,
         param_names,
         param_removal_code_names,
+        metab_names,
         num_alpha_params,
     )
+    if isempty(starting_param_removal_codes)
+        @error "Equations for this enzymes with $(num_param_range[1]) parameters do not exist. One reason could be that some parameters are unidentifiable based on the data so upper bound of range_number_params should be reduced. Check number of practically_unidentifiable_params with find_practically_unidentifiable_params(data, param_names)."
+        return
+    end
+
 
     nt_param_removal_codes = starting_param_removal_codes
     nt_previous_param_removal_codes = similar(nt_param_removal_codes)
@@ -316,6 +322,7 @@ function calculate_all_parameter_removal_codes_w_num_params(
     all_param_removal_codes,
     param_names::Tuple{Symbol,Vararg{Symbol}},
     param_removal_code_names::Tuple{Symbol,Vararg{Symbol}},
+    metab_names::Tuple{Symbol,Vararg{Symbol}},
     num_alpha_params::Int,
 )
     codes_with_num_params = Tuple[]
@@ -338,7 +345,34 @@ function calculate_all_parameter_removal_codes_w_num_params(
     end
     nt_param_removal_codes =
         [NamedTuple{param_removal_code_names}(x) for x in unique(codes_with_num_params)]
-    return nt_param_removal_codes
+
+    # ensure that if K_S1 = Inf then all K_S1_S2 and all other K containing S1 in qssa cannot be 2, which stands for (K_S1_S2)^2 = K_S1 * K_S2
+    if any([occursin("allo", string(key)) for key in keys(nt_param_removal_codes[1])])
+        filtered_nt_param_removal_codes = nt_param_removal_codes
+    else
+        filtered_nt_param_removal_codes = NamedTuple[]
+        for nt_param_removal_code in nt_param_removal_codes
+            if all(
+                nt_param_removal_code[Symbol("K_" * string(metab))] != 1 for
+                metab in metab_names
+            )
+                push!(filtered_nt_param_removal_codes, nt_param_removal_code)
+            else
+                one_metab_codes = metab_names[findall(
+                    nt_param_removal_code[Symbol("K_" * string(metab))] == 1 for
+                    metab in metab_names
+                )]
+                if all(
+                    nt_param_removal_code[param_name] != 2 for
+                    param_name in keys(nt_param_removal_code) if
+                    any(occursin.(string.(one_metab_codes), string(param_name)))
+                )
+                    push!(filtered_nt_param_removal_codes, nt_param_removal_code)
+                end
+            end
+        end
+    end
+    return filtered_nt_param_removal_codes
 end
 
 """

test/tests_for_optimal_rate_eq_selection.jl

@@ -165,6 +165,7 @@ end
 
 #test calculate_all_parameter_removal_codes_w_num_params
 num_metabolites = rand(4:8)
+metab_names = Tuple(Symbol("Metabolite$(i)") for i = 1:num_metabolites)
 n_alphas = rand(1:4)
 param_names = (
     :L,
@@ -193,6 +194,7 @@ nt_param_subset_codes_w_num_params =
         all_param_removal_codes,
         param_names,
         param_removal_code_names,
+        metab_names,
         n_alphas,
     )
 #ensure that funct_output_param_subset_codes have the correct number of parameters