cellnopt tutorial and added data ploting

corneto/methods/signal/cellnopt_ilp.py

@@ -352,20 +352,51 @@ def plot_solution_network_active_edges(G, P, iexp):
     )
 
 
-def plot_fitness(G, exp_list, P, measured_only=False):
+def plot_fitness(G, exp_list, P, measured_only=False, **kwargs):
     """Plot the fitness of the model simulation vs measurements.
 
     PARAMETERS:
     - G: corneto.Graph object
     - exp_list: dictionary of experiments
     - P: solution of the ILP model
     - measured_only: if True, plot only the measured nodes, otherwise plot all nodes
+    - **kwargs arguments are passed to subplots and figures of matplotlib
 
     TODO: there are some assumptions, like the first experiment is the reference experiment and it is called 'exp0'
     """
     import matplotlib.pyplot as plt
 
     N_exps = len(exp_list)
+
+    # Ensure that all experiments have the input and output variables
+    for exp in exp_list.values():
+        if "input" not in exp:
+            raise ValueError("Input not found in experiment")
+        if "output" not in exp:
+            raise ValueError("Output not found in experiment")
+
+    # Collect the input and output variables
+    input_matrix, input_vars = collect_field_into_matrix(exp_list, "input")
+    output_matrix, output_vars = collect_field_into_matrix(exp_list, "output")
+
+    # Check if inhibition is present in any of the experiments
+    inhibition_present = any("inhibition" in exp for exp in exp_list.values())
+    if inhibition_present:
+        inhibition_matrix, inhibition_vars = collect_field_into_matrix(
+            exp_list, "inhibition"
+        )
+        perturbation_matrix = np.hstack((input_matrix, inhibition_matrix))
+        perturbation_vars = input_vars + inhibition_vars
+    else:
+        perturbation_matrix = input_matrix
+        perturbation_vars = input_vars
+
+    # Create the figure
+    # Set colors: input colors are blue, inhibition colors are red
+    perturbation_colors = ["blue"] * len(input_vars)
+    if inhibition_present:
+        perturbation_colors += ["red"] * len(inhibition_vars)
+
     N_nodes = len(G.V)
     output_names = list(
         {key for exp in exp_list.values() for key in exp["output"].keys()}
@@ -374,13 +405,16 @@ def plot_fitness(G, exp_list, P, measured_only=False):
     # depending on the flag measured_only, we can plot only the measured nodes or all nodes
     if measured_only:
         fig, axs = plt.subplots(
-            N_exps - 1, len(output_names), figsize=(20, 5), squeeze=False
+            N_exps - 1, len(output_names) + 1, squeeze=False, **kwargs
         )
     else:
-        fig, axs = plt.subplots(N_exps - 1, N_nodes, figsize=(20, 5), squeeze=False)
+        fig, axs = plt.subplots(N_exps - 1, N_nodes + 1, squeeze=False, **kwargs)
 
     fig.tight_layout(pad=0.0)
 
+    # Adjust the space between subplots
+    plt.subplots_adjust(wspace=0.1, hspace=0.1)
+
     for exp, iexp in zip(exp_list, range(N_exps)):
         if iexp == 0:
             continue
@@ -396,9 +430,8 @@ def plot_fitness(G, exp_list, P, measured_only=False):
                         P.expr.vertex_value.value[imarker_inG, 0],
                         min(P.expr.vertex_value.value[imarker_inG, iexp], 1),
                     ],
+                    "bo-",
                     label=G.V[imarker_inG],
-                    color="blue",
-                    linestyle="-",
                 )
 
                 if G.V[imarker_inG] in exp_list[exp]["output"].keys():
@@ -410,11 +443,15 @@ def plot_fitness(G, exp_list, P, measured_only=False):
                         ],
                         "ro-",
                     )
-                axs[iexp - 1, imarker].set_ylim([-0.01, 1.1])
+                axs[iexp - 1, imarker].set_ylim([-0.05, 1.1])
                 if iexp == 1:
-                    axs[iexp - 1, imarker].set_title(G.V[imarker_inG])
+                    axs[iexp - 1, imarker].set_title(output_names[imarker])
+                if iexp != N_exps - 1:
+                    axs[iexp - 1, imarker].set_xticks([])
                 if imarker == 0:
-                    axs[iexp - 1, imarker].set_ylabel(f"Experiment {iexp}")
+                    axs[iexp - 1, imarker].set_ylabel(f"Exp. {iexp}")
+                else:
+                    axs[iexp - 1, imarker].set_yticks([])
         else:
             for imarker in range(N_nodes):
                 axs[iexp - 1, imarker].plot(
@@ -423,9 +460,10 @@ def plot_fitness(G, exp_list, P, measured_only=False):
                         P.expr.vertex_value.value[imarker, 0],
                         min(P.expr.vertex_value.value[imarker, iexp], 1),
                     ],
+                    "bo-",
                     label=G.V[imarker],
-                    color="blue",
-                    linestyle="-",
+                    # color="blue",
+                    # linestyle="o-",
                 )
 
                 if G.V[imarker] in exp_list[exp]["output"].keys():
@@ -437,10 +475,167 @@ def plot_fitness(G, exp_list, P, measured_only=False):
                         ],
                         "ro-",
                     )
-                axs[iexp - 1, imarker].set_ylim([-0.01, 1.1])
+                axs[iexp - 1, imarker].set_ylim([-0.05, 1.1])
                 if iexp == 1:
                     axs[iexp - 1, imarker].set_title(G.V[imarker])
+                if iexp != N_exps - 1:
+                    axs[iexp - 1, imarker].set_xticks([])
                 if imarker == 0:
-                    axs[iexp - 1, imarker].set_ylabel(f"Experiment {iexp}")
+                    axs[iexp - 1, imarker].set_ylabel(f"Exp. {iexp}")
+                else:
+                    axs[iexp - 1, imarker].set_yticks([])
+        # Plot perturbation
+        if measured_only:
+            plot_location = len(output_names)
+        else:
+            plot_location = N_nodes
+
+        axs[iexp - 1, plot_location].bar(
+            range(len(perturbation_vars)),
+            perturbation_matrix[iexp],
+            color=perturbation_colors,
+        )
+        if iexp == N_exps - 1:
+            axs[iexp - 1, plot_location].set_xticks(range(len(perturbation_vars)))
+            axs[iexp - 1, plot_location].set_xticklabels(
+                perturbation_vars, rotation=45
+            )
+        else:
+            # No xtick label
+            axs[iexp - 1, plot_location].set_xticks([])
+
+        axs[iexp - 1, plot_location].set_ylim([-0.01, 1.1])
+        axs[iexp - 1, plot_location].set_yticks([])
+        if iexp == 1:
+            axs[iexp - 1, plot_location].set_title("Pert.")
+
+    plt.show()
+
+
+def collect_field_into_matrix(experiments, field_name="input"):
+    """Collects the field_name values into matrix.
+
+    Collects the field_name values (input, inhibition etc) from a dictionary
+    of experiments and returns them as a numpy array.
+
+    PARAMETERS:
+    - experiments: dictionary of experiments containing input values
+    - field_name: name of the field to collect (default: 'input')
+
+    Returns:
+    - input_matrix: numpy array of input values
+    - input_vars: list of unique input variable names
+
+    """
+    # Collect all unique input names
+    input_vars = set()
+    for exp in experiments.values():
+        # ensure field_name exists
+        if field_name not in exp:
+            raise ValueError("Field name (" + field_name + ")  not found in experiment")
+        input_vars.update(exp[field_name].keys())
+
+    input_vars = sorted(input_vars)  # Sorting to keep a consistent order
+    data = []
+
+    # Collect input values for each experiment
+    for exp in experiments.values():
+        row = [exp[field_name].get(var, 0) for var in input_vars]
+        data.append(row)
+
+    # Convert the data to a numpy array
+    input_matrix = np.array(data)
+
+    return input_matrix, input_vars
+
+
+def plot_data(exp_list):
+    """Plot the data.
+
+    PARAMETERS:
+    - exp_list: dictionary of experiments
+
+    """
+    import matplotlib.pyplot as plt
+
+    N_exps = len(exp_list)
+
+    # Ensure that all experiments have the input and output variables
+    for exp in exp_list.values():
+        if "input" not in exp:
+            raise ValueError("Input not found in experiment")
+        if "output" not in exp:
+            raise ValueError("Output not found in experiment")
+
+    # Collect the input and output variables
+    input_matrix, input_vars = collect_field_into_matrix(exp_list, "input")
+    output_matrix, output_vars = collect_field_into_matrix(exp_list, "output")
+
+    # Check if inhibition is present in any of the experiments
+    inhibition_present = any("inhibition" in exp for exp in exp_list.values())
+    if inhibition_present:
+        inhibition_matrix, inhibition_vars = collect_field_into_matrix(
+            exp_list, "inhibition"
+        )
+        perturbation_matrix = np.hstack((input_matrix, inhibition_matrix))
+        perturbation_vars = input_vars + inhibition_vars
+    else:
+        perturbation_matrix = input_matrix
+        perturbation_vars = input_vars
+
+    # Create the figure
+    # Set colors: input colors are blue, inhibition colors are red
+    perturbation_colors = ["blue"] * len(input_vars)
+    if inhibition_present:
+        perturbation_colors += ["red"] * len(inhibition_vars)
+
+    fig, axs = plt.subplots(N_exps - 1, len(output_vars) + 1, squeeze=False)
+
+    fig.tight_layout(pad=0.0)
+    # Adjust the space between subplots
+    plt.subplots_adjust(wspace=0.1, hspace=0.1)
+
+    for exp, iexp in zip(exp_list, range(N_exps)):
+        if iexp == 0:
+            continue
+
+        for imarker in range(len(output_vars)):
+            # output_names[imarker] is the name of the output node, find the position in the graph
+            imarker_name = output_vars[imarker]
+
+            axs[iexp - 1, imarker].plot(
+                [0, 10],
+                [output_matrix[0, imarker], output_matrix[iexp, imarker]],
+                "ro-",
+            )
+            axs[iexp - 1, imarker].set_ylim([-0.01, 1.1])
+
+            if iexp == 1:
+                axs[iexp - 1, imarker].set_title(imarker_name)
+            if iexp != N_exps - 1:
+                axs[iexp - 1, imarker].set_xticks([])
+            if imarker == 0:
+                axs[iexp - 1, imarker].set_ylabel(f"Exp. {iexp}")
+            else:
+                axs[iexp - 1, imarker].set_yticks([])
+
+        # Plot perturbation
+        axs[iexp - 1, len(output_vars)].bar(
+            range(len(perturbation_vars)),
+            perturbation_matrix[iexp],
+            color=perturbation_colors,
+        )
+        if iexp == N_exps - 1:
+            axs[iexp - 1, len(output_vars)].set_xticks(range(len(perturbation_vars)))
+            axs[iexp - 1, len(output_vars)].set_xticklabels(
+                perturbation_vars, rotation=45
+            )
+        else:
+            # No xtick label
+            axs[iexp - 1, len(output_vars)].set_xticks([])
 
+        axs[iexp - 1, len(output_vars)].set_ylim([-0.01, 1.1])
+        axs[iexp - 1, len(output_vars)].set_yticks([])
+        if iexp == 1:
+            axs[iexp - 1, len(output_vars)].set_title("Pert.")
     plt.show()