updating plot look for direct simulation

crnt4sbml/bistability_analysis.py

@@ -765,8 +765,8 @@ def __parent_initialize_direct_simulation(self):
         global numpy
         global warnings
         global pandas
-        global ggplot, aes, geom_line, ylim, scale_color_distiller
-        global facet_wrap, theme_bw, geom_path, geom_point, labs, annotate
+        global ggplot, aes, geom_path, arrow, geom_point, scale_alpha_manual, theme, element_blank
+        global element_text, element_rect, element_line, guides, guide_legend
         global itg
         global sympy
         global time
@@ -778,8 +778,8 @@ def __parent_initialize_direct_simulation(self):
         import numpy
         import warnings
         import pandas
-        from plotnine import ggplot, aes, geom_line, ylim, scale_color_distiller, facet_wrap, theme_bw, geom_path, \
-            geom_point, labs, annotate
+        from plotnine import ggplot, aes, geom_path, arrow, geom_point, scale_alpha_manual, theme, element_blank, \
+            element_text, element_rect, element_line, guides, guide_legend
         import scipy.integrate as itg
         import time
 
@@ -1197,25 +1197,49 @@ def jac_f(t, cs, ks, ode_lambda_functions, jacobian):
     def __plot_direct_simulation(self, pcp_scan, forward_scan, reverse_scan, path, index):
 
         out = pandas.DataFrame(columns=['dir', 'signal'] + [self.__response])
+        fwd_scan = []
         for i in range(len(forward_scan)):
             out_i = pandas.DataFrame([forward_scan[i]], columns=[out.columns[2]])
+            fwd_scan.append(forward_scan[i])
             out_i['signal'] = pcp_scan[i]
             out_i['dir'] = 'Forward scan'
             out = pandas.concat([out, out_i[out.columns]])
+        rvrs_scan = []
         for i in range(len(reverse_scan)):
             out_i = pandas.DataFrame([reverse_scan[i]], columns=[out.columns[2]])
+            rvrs_scan.append(reverse_scan[i])
             out_i['signal'] = pcp_scan[i]
             out_i['dir'] = 'Reverse scan'
             out = pandas.concat([out, out_i[out.columns]])
 
-        g = (ggplot(out)
-             + aes(x='signal', y=self.__response, color='dir')
-             + labs(x=f"{self.__signal} total", y=f"[{self.__response}]", color="")
-             + geom_path(size=2, alpha=0.5)
-             + geom_point(color="black")
-             + theme_bw()
-             + geom_point(color="black"))
+        def mask(df, key, value):
+            return df[df[key] == value]
+
+        pandas.DataFrame.mask = mask
 
-        g.save(filename=path + f"/sim_bif_diag_{index}.png", format="png", width=6, height=4, units='in', verbose=False)
+        # flip one of the directions
+        xf = out.mask('dir', 'Forward scan')
+        xr = out.mask('dir', 'Reverse scan')
+        xr = xr.iloc[::-1]
 
-        return g
+        out = pandas.concat([xf, xr])
+
+        g = (ggplot(out)
+             + aes(x='signal', y='s2', color='dir', alpha='dir', group='dir')
+             + geom_path(size=4, arrow=arrow(), show_legend=False)
+             + geom_point(size=2, shape='s', stroke=0.0)
+             + geom_point(color="black", size=2, alpha=1.0)
+             + scale_alpha_manual(values=[0.85, 0.35], guide=False)
+             # + scale_color_manual(values=["red", "blue"], labels=["d", "k"])
+             + guides(color=guide_legend(override_aes={'size': 6, 'alpha': [0.85, 0.35]}))
+             + theme(legend_title=element_blank(),
+                     axis_text=element_text(size=0.8 * 20), legend_key=element_rect(color='None', fill='None'),
+                     panel_background=element_rect(fill='None'),
+                     panel_border=element_rect(fill='None', color='#7f7f7f'),
+                     panel_grid_major=element_line(color='#E5E5E5', size=0.8),
+                     panel_grid_minor=element_line(color='#FAFAFA', size=1),
+                     strip_background=element_rect(fill='#CCCCCC', color='#7F7F7F', size=1)))
+
+        g.save(filename=path + f"/sim_bif_diag_{index}.png", format="png", width=8, height=5, units='in', verbose=False)
+
+        return g

crnt4sbml/c_graph.py

@@ -523,6 +523,9 @@ def __simplify_nullspace(self, the_null_space):
         return temp_null
 
     def __create_non_negative_b_matrix(self, the_null_space, sizes):
+        '''
+        Idea inspired by https://r.789695.n4.nabble.com/convex-nonnegative-basis-vectors-in-nullspace-of-matrix-td4548822.html.
+        '''
 
         a_null = numpy.zeros((len(self.__species), sizes))
         for i in range(len(self.__species)):
@@ -532,10 +535,11 @@ def __create_non_negative_b_matrix(self, the_null_space, sizes):
             a_null[i, :] = temp_vec
 
         a_eq = numpy.array([numpy.sum(a_null, axis=0)])
+
         # must multiply by negative one because in optimization we have the inequality <= 0.0
         a_ub = -1.0*a_null
         b_ub = numpy.zeros(len(self.__species))
-        b_eq = numpy.array([1.0])        
+        b_eq = numpy.array([1.0])
 
         # defining the number of solutions to simulate
         num_sols = ((a_ub.shape[0]+1)*(a_ub.shape[1]+1))*10
@@ -561,7 +565,7 @@ def __create_non_negative_b_matrix(self, the_null_space, sizes):
         num_rows = numpy.size(unique_vecs, 0)
 
         # taking the smallest nonzero entry of the unique vectors
-        # and dividing by it this will hopfully create nice
+        # and dividing by it this will hopefully create nice
         # looking vectors that are whole numbers
         for i in range(num_rows):
             minval = numpy.min(unique_vecs[i, numpy.nonzero(unique_vecs[i, :])])

crnt4sbml/general_approach.py

@@ -731,6 +731,8 @@ def run_direct_simulation(self, params_for_global_min=None, dir_path="./dir_sim_
 
         Note: This routine is more expensive than the numerical continuation routines, but can provide solutions
         when the Jacobian of the ODE system is always singular. A parallel version of this routine is available.
+        The routine automatically produces plots of the direct simulation runs and puts them in the user specified
+        dir_path.
 
         Parameters
         ------------
@@ -755,6 +757,9 @@ def run_direct_simulation(self, params_for_global_min=None, dir_path="./dir_sim_
             right_multiplier: float
                 A float value that determines the percentage of the signal that will be searched to the right of the signal
                 value. For example, the upperbound for the signal range will be signal_value + signal_value*right_multiplier.
+        Returns
+        ---------
+            list_of_ggplots: list of ggplots produced by plotnine
 
         Example
         ---------

crnt4sbml/mass_conservation_approach.py

@@ -1233,6 +1233,8 @@ def run_direct_simulation(self, response=None, signal=None, params_for_global_mi
 
         Note: This routine is more expensive than the numerical continuation routines, but can provide solutions
         when the Jacobian of the ODE system is always singular. A parallel version of this routine is available.
+        The routine automatically produces plots of the direct simulation runs and puts them in the user specified
+        dir_path.
 
         Parameters
         ------------
@@ -1261,6 +1263,9 @@ def run_direct_simulation(self, response=None, signal=None, params_for_global_mi
             right_multiplier: float
                 A float value that determines the percentage of the signal that will be searched to the right of the signal
                 value. For example, the upperbound for the signal range will be signal_value + signal_value*right_multiplier.
+        Returns
+        ---------
+            list_of_ggplots: list of ggplots produced by plotnine
 
         Example
         ---------

example_scripts/Fig1C_GB.py

@@ -14,12 +14,13 @@
 ds16 = re3*s6*s7 + s16*(-re3r - re4)
 ds15 = re5*s1*s6 + s15*(-re5r - re6)
 
-replacements = [[s1, (C3 - 2.0*s15 - s16 - s3 - s6)]]
+# replacements = [[s1, (C3 - 2.0*s15 - s16 - s3 - s6)]]
+
+# ds2 = sympy.simplify(ds2.subs(s1, (C3 - 2.0*s15 - s16 - s3 - s6)))
+# ds3 = sympy.simplify(ds3.subs(s1, (C3 - 2.0*s15 - s16 - s3 - s6)))
+# ds6 = sympy.simplify(ds6.subs(s1, (C3 - 2.0*s15 - s16 - s3 - s6)))
+# ds15 = sympy.simplify(ds15.subs(s1, (C3 - 2.0*s15 - s16 - s3 - s6)))
 
-ds2 = sympy.simplify(ds2.subs(s1, (C3 - 2.0*s15 - s16 - s3 - s6)))
-ds3 = sympy.simplify(ds3.subs(s1, (C3 - 2.0*s15 - s16 - s3 - s6)))
-ds6 = sympy.simplify(ds6.subs(s1, (C3 - 2.0*s15 - s16 - s3 - s6)))
-ds15 = sympy.simplify(ds15.subs(s1, (C3 - 2.0*s15 - s16 - s3 - s6)))
 
 
 # k1, k2, k3, k4, s, y, x, v0 = symbols('k1 k2 k3 k4 s y x v0', real=True)
@@ -36,11 +37,21 @@
 # ds16 = re3*s6*s7 + s16*(-re3r - re4)
 # ds15 = re5*s1*s6 + s15*(-re5r - re6)
 
-# ds15 = sympy.simplify(sympy.expand(re5*s6*(C3 - 2.0*s15 - s16 - s3 - s6) + s15*(-re5r - re6)))
-# ds3 = sympy.simplify(sympy.expand(re1*(C2 - s3)*(C3 - 2.0*s15 - s16 - s3 - s6) + s3*(-re1r - re2)))
-# ds6 = sympy.simplify(sympy.expand(re2*s3 - re3*s6*(C1 - s16) + re3r*s16 - re5*s6*(C3 - 2.0*s15 - s16 - s3 - s6) + s15*(re5r + 2.0*re6)))
-# ds16 = sympy.simplify(sympy.expand(re3*s6*(C1 - s16) + s16*(-re3r - re4)))
+ds15 = sympy.simplify(sympy.expand(re5*s6*(C3 - 2.0*s15 - s16 - s3 - s6) + s15*(-re5r - re6)))
+ds3 = sympy.simplify(sympy.expand(re1*(C2 - s3)*(C3 - 2.0*s15 - s16 - s3 - s6) + s3*(-re1r - re2)))
+ds6 = sympy.simplify(sympy.expand(re2*s3 - re3*s6*(C1 - s16) + re3r*s16 - re5*s6*(C3 - 2.0*s15 - s16 - s3 - s6) + s15*(re5r + 2.0*re6)))
+ds16 = sympy.simplify(sympy.expand(re3*s6*(C1 - s16) + s16*(-re3r - re4)))
+
+print(ds3)
+print("")
+print(ds6)
+print("")
+print(ds16)
+print("")
+print(ds15)
+print("")
 
+sys.exit()
 
 # in CoCoLib:
 # use QQab ::= QQ[re1, re1r, re2, re3, re3r, re4, re5, re5r, re6, C1, C2, C3];
@@ -50,6 +61,16 @@
 # SetVerbosityLevel(100);
 # ReducedGBasis(I);
 
+# nondimensional Fig1Ci
+# in CoCoLib:
+# use QQab ::= QQ[a,b,c,d,e,f,g,h,k,l];
+# K := NewFractionField(QQab);
+# use K[B, AE1, BE2, AB], lex;
+# I := ideal(AE1 + b*BE2 + AB*d + 2*AB*e + 2*AB*B*f + AE1*B*f + f*B^2 + B*BE2*f - B*h + B*BE2*h - B*f*l, g*k*2*AB + g*k*AE1 + g*k*B + g*k*BE2 - g*k*l + AE1 + AE1*a - AE1*2*AB*g - g*AE1^2 - AE1*B*g - AE1*BE2*g + AE1*g*l, -b*BE2 - BE2*c + B*h - B*BE2*h, AB*d + AB*e + 2*AB*B*f + AE1*B*f + f*B^2 + B*BE2*f - B*f*l);
+# SetVerbosityLevel(100);
+# ReducedGBasis(I);
+
+
 # SageMath Fig1Ci:
 # A.<re1, re1r, re2, re3, re3r, re4, re5, re5r, re6, C1, C2, C3> = PolynomialRing(QQ)
 # F = A.fraction_field()
@@ -65,6 +86,14 @@
 # J=ideal(C2*C3*re1 - 2*C2*re1*s15 - C2*re1*s16 - C2*re1*s3 - C2*re1*s6 - C3*re1*s3 + 2*re1*s15*s3 + re1*s16*s3 + re1*s3^2 + re1*s3*s6 - re1r*s3 - re2*s3, -C1*re3*s6 - C3*re5*s6 + re2*s3 + re3*s16*s6 + re3r*s16 + 2*re5*s15*s6 + re5*s16*s6 + re5*s3*s6 + re5*s6^2 + re5r*s15 + 2*re6*s15, C1*re3*s6 - re3*s16*s6 - re3r*s16 - re4*s16, C3*re5*s6 - 2*re5*s15*s6 - re5*s16*s6 - re5*s3*s6 - re5*s6^2 - re5r*s15 - re6*s15)
 # G = ideal groebnerBasis(J, Strategy=>"F4")
 
+# Macaulay Fig1Ci full non dimensional:
+# K = frac(QQ[a,b,c,d,e,f,g,h,k,l]);
+# R = K[BE2, E2, AE1, E1, A, B, AB, MonomialOrder => Lex];
+# gbTrace = 3
+# J=ideal(-a*AE1 - BE2*c - AB*d + A*B*f + A*E1*g, -AE1 - b*BE2 - AB*d - 2*AB*e + A*B*f + B*E2*h, -AE1 - a*AE1 + A*E1*g, -b*BE2 - BE2*c + B*E2*h, -AE1 - a*AE1 + A*E1*g, -b*BE2 - BE2*c + B*E2*h, -AB*d - AB*e + A*B*f, -1 - BE2 + E2, -AE1 + E1 - k, A - 2*AB - AE1 - B - BE2 - l)
+# G = ideal groebnerBasis(J, Strategy=>"F4")
+
+
 # SageMath simple example:
 # A.<k1, k2, k3, k4> = PolynomialRing(QQ)
 # F = A.fraction_field()
@@ -172,7 +201,7 @@
 species = [s15, s3, s6, s16]
 # # equations = [ds1, ds2, ds3, ds6, ds7, ds15, ds16]
 equations = [ds3, ds6, ds16, ds15]
-
+print(equations)
 start = time.time()
 gb = groebner(equations, species, order='lex', method='f5b')
 end = time.time()

example_scripts/run_general_quickstart.py

@@ -3,6 +3,7 @@
 network = crnt4sbml.CRNT("../sbml_files/Fig1Ci.xml")
 
 approach = network.get_general_approach()
+
 bnds = approach.get_optimization_bounds()
 
 approach.initialize_general_approach(signal="C3", response="s15", fix_reactions=True)

example_scripts/run_mass_conservation_quickstart.py

@@ -1,3 +1,5 @@
+import sys
+sys.path.insert(0, "..")
 import crnt4sbml
 
 network = crnt4sbml.CRNT("../sbml_files/Fig1Ci.xml")