Implementing an alternative network sampler#
In this tutorial, we demonstrate how to use CORNETO to implement a simple sampler to explore alternative solutions for network inference problems. By perturbing the objective function, we can explore the solution space and identify robust as well as variable components of the network.
We will illustrate this process using a simple example of a network inference problem with CARNIVAL, using the data from the CARNIVAL transcriptomics tutorial.
NOTE: This notebook uses
gurobias the solver andcvxpyas the backend to accelerate the generation of alternative solutions using thewarm_startoption from cvxpy
import corneto as cn
cn.info()
|
|
from corneto.data import GraphData
# Load dataset from tutorial
dataset = GraphData.load("data/carnival_transcriptomics_dataset.zip")
from corneto.methods import CarnivalFlow
m = CarnivalFlow(lambda_reg=0.1)
P = m.build(dataset.graph, dataset.data)
P.objectives
Unreachable vertices for sample: 0
[error_sample1_0: Expression(AFFINE, UNKNOWN, ()),
regularization_edge_has_signal: Expression(AFFINE, NONNEGATIVE, ())]
# We will sample solutions by perturbing the edge_has_signal variable
P.expr
{'_flow': _flow: Variable((3251,), _flow),
'edge_inhibits': edge_inhibits: Variable((3251, 1), edge_inhibits, boolean=True),
'edge_activates': edge_activates: Variable((3251, 1), edge_activates, boolean=True),
'const0x103326ca2e6db1fe': const0x103326ca2e6db1fe: Constant(CONSTANT, NONNEGATIVE, (937, 3251)),
'_dag_layer': _dag_layer: Variable((937, 1), _dag_layer),
'const0x18ac965210dd36c7': const0x18ac965210dd36c7: Constant(CONSTANT, NONNEGATIVE, (937, 3251)),
'flow': _flow: Variable((3251,), _flow),
'vertex_value': Expression(AFFINE, UNKNOWN, (937, 1)),
'vertex_activated': Expression(AFFINE, NONNEGATIVE, (937, 1)),
'vertex_inhibited': Expression(AFFINE, NONNEGATIVE, (937, 1)),
'edge_value': Expression(AFFINE, UNKNOWN, (3251, 1)),
'edge_has_signal': Expression(AFFINE, NONNEGATIVE, (3251, 1)),
'vertex_max_depth': _dag_layer: Variable((937, 1), _dag_layer)}
Sampling Solutions Using Random Perturbation#
We’ll implement a generic approach to sample alternative solutions by:
Introducing small random perturbations to the objective function
Re-optimizing the model with these perturbations
Collecting solutions that remain close to optimal (within a tolerance threshold)
This technique allows us to explore the space of near-optimal solutions and assess the robustness of different model components.
import numpy as np
import pandas as pd
def sample_alternative_solutions(
problem,
target_variable,
percentage=0.1,
scale=0.03,
rel_opt_tol=0.05,
max_samples=30,
time_limit=60,
):
"""Sample alternative solutions by perturbing the objective function.
Parameters:
-----------
problem : CORNETO optimization problem
The optimization problem to sample
target_variable : variable expression
The variable to perturb (e.g., P.expr.edge_has_signal or P.expr.vertex_value)
percentage : float
Percentage of elements to perturb
scale : float
Scale parameter for normal distribution perturbation
rel_opt_tol : float
Relative optimality tolerance
max_samples : int
Maximum number of samples to generate
time_limit : int
Time limit for solver
Returns:
--------
dict : A dictionary containing the results
"""
# Save the original problem state
original_problem = problem.solve(solver="gurobi", verbosity=1)
# Store original objective values
orig_objectives = []
for o in problem.objectives:
print(o.name, o.value)
orig_objectives.append(o.value)
# Prepare for perturbation
shape = target_variable.shape
num_elements = shape[0]
n_elements = int(num_elements * percentage)
vec = np.zeros(shape=num_elements)
pert = problem.backend.Parameter("perturbation", shape=vec.shape, value=vec)
# Add perturbation objective
problem.add_objective(target_variable.T @ pert, name="perturbation")
# Sample solutions
selected_edges = []
edge_values = []
vertex_values = []
for i in range(max_samples):
# Generate random perturbation
c = np.random.normal(scale=scale, size=n_elements)
random_indices = np.random.choice(vec.shape[0], size=n_elements, replace=False)
vec[:] = 0 # Reset vector
vec[random_indices] = c
# Update parameter and solve
pert.value = vec
solved_problem = problem.solve(
solver="gurobi",
warm_start=True,
ignore_dpp=True,
TimeLimit=time_limit,
verbosity=0,
)
print(f"Sample {i + 1}/{max_samples}: Objective value = {solved_problem.value:.4f}")
# Check if solution is within tolerance
accept = True
for o_orig, o in zip(orig_objectives, problem.objectives):
if o.name == "perturbation":
continue
relative_error = np.abs(o.value - o_orig) / np.abs(o_orig + 1e-10)
print(f" {o.name}: {o.value:.4f}, rel error: {relative_error:.4f}")
if relative_error > rel_opt_tol:
print(f" > Rejected! ({relative_error:.4f} >= {rel_opt_tol})")
accept = False
break
if accept:
selected_edges.append(problem.expr.edge_has_signal.value > 0.5)
edge_values.append(problem.expr.edge_value.value)
vertex_values.append(problem.expr.vertex_value.value)
return {
"selected_edges": selected_edges,
"edge_values": edge_values,
"vertex_values": vertex_values,
}
Edge-Based Perturbation#
First, we’ll sample solutions by perturbing the edge variables. This allows us to explore alternative network topologies where different edges might be selected.
# Create a fresh model
m = CarnivalFlow(lambda_reg=0.1)
P = m.build(dataset.graph, dataset.data)
# Sample solutions by perturbing edge_has_signal
edge_results = sample_alternative_solutions(
problem=P,
target_variable=P.expr.edge_has_signal,
percentage=0.10,
scale=0.03,
rel_opt_tol=0.05,
max_samples=30,
)
# Analyze edge-based results
df_sols = pd.DataFrame(np.concatenate(edge_results["selected_edges"], axis=1), index=m.processed_graph.E).astype(int)
print(f"Generated {df_sols.shape[1]} alternative solutions")
df_sols.head()
Unreachable vertices for sample: 0
===============================================================================
CVXPY
v1.7.1
===============================================================================
-------------------------------------------------------------------------------
Compilation
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Numerical solver
-------------------------------------------------------------------------------
Set parameter Username
Set parameter LicenseID to value 2775543
Academic license - for non-commercial use only - expires 2027-02-05
Set parameter OutputFlag to value 1
Set parameter QCPDual to value 1
Gurobi Optimizer version 12.0.3 build v12.0.3rc0 (mac64[arm] - Darwin 25.5.0 25F84)
CPU model: Apple M4 Pro
Thread count: 12 physical cores, 12 logical processors, using up to 12 threads
Non-default parameters:
QCPDual 1
Optimize a model with 30607 rows, 10690 columns and 108040 nonzeros
Model fingerprint: 0xf48e83e6
Variable types: 4188 continuous, 6502 integer (6502 binary)
Coefficient statistics:
Matrix range [1e+00, 9e+02]
Objective range [1e-01, 7e+00]
Bounds range [1e+00, 1e+00]
RHS range [1e+00, 1e+03]
Found heuristic solution: objective 0.1000000
Presolve removed 17840 rows and 1466 columns
Presolve time: 0.24s
Presolved: 12767 rows, 9224 columns, 81311 nonzeros
Variable types: 3504 continuous, 5720 integer (5714 binary)
Root relaxation: objective -1.186277e+02, 1171 iterations, 0.05 seconds (0.07 work units)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 -118.62771 0 53 0.10000 -118.62771 - - 0s
0 0 -118.51700 0 62 0.10000 -118.51700 - - 0s
H 0 0 -3.7564221 -118.51700 3055% - 0s
H 0 0 -76.6465356 -118.51700 54.6% - 0s
H 0 0 -76.7465356 -118.51700 54.4% - 0s
0 0 -118.21224 0 60 -76.74654 -118.21224 54.0% - 0s
H 0 0 -79.0760121 -118.21224 49.5% - 0s
H 0 0 -113.3642782 -118.12057 4.20% - 0s
0 0 -118.12057 0 55 -113.36428 -118.12057 4.20% - 0s
0 0 -118.07057 0 41 -113.36428 -118.07057 4.15% - 0s
H 0 0 -113.4642782 -118.07057 4.06% - 0s
0 0 -118.07057 0 36 -113.46428 -118.07057 4.06% - 0s
0 0 -118.07057 0 33 -113.46428 -118.07057 4.06% - 0s
0 0 -118.07057 0 18 -113.46428 -118.07057 4.06% - 0s
0 0 -118.07057 0 18 -113.46428 -118.07057 4.06% - 0s
H 0 0 -117.5705690 -118.07057 0.43% - 0s
0 0 -118.07057 0 18 -117.57057 -118.07057 0.43% - 0s
0 0 -118.07057 0 13 -117.57057 -118.07057 0.43% - 1s
0 0 -118.07057 0 13 -117.57057 -118.07057 0.43% - 1s
0 0 -118.07057 0 2 -117.57057 -118.07057 0.43% - 1s
0 0 -118.07057 0 37 -117.57057 -118.07057 0.43% - 1s
0 0 -118.07057 0 25 -117.57057 -118.07057 0.43% - 1s
0 0 -118.07057 0 20 -117.57057 -118.07057 0.43% - 1s
0 0 -118.07057 0 2 -117.57057 -118.07057 0.43% - 1s
H 0 0 -117.8705690 -118.07057 0.17% - 1s
0 0 -118.07057 0 11 -117.87057 -118.07057 0.17% - 1s
0 0 -118.07057 0 11 -117.87057 -118.07057 0.17% - 1s
0 0 -118.07057 0 11 -117.87057 -118.07057 0.17% - 1s
0 0 -118.07057 0 44 -117.87057 -118.07057 0.17% - 2s
H 0 0 -117.9705690 -118.07057 0.08% - 2s
0 0 -118.07057 0 38 -117.97057 -118.07057 0.08% - 2s
0 0 -118.07057 0 31 -117.97057 -118.07057 0.08% - 2s
0 0 -118.07057 0 23 -117.97057 -118.07057 0.08% - 2s
0 0 -118.07057 0 22 -117.97057 -118.07057 0.08% - 2s
0 0 -118.07057 0 22 -117.97057 -118.07057 0.08% - 2s
0 0 -118.07057 0 22 -117.97057 -118.07057 0.08% - 2s
0 0 -118.07057 0 13 -117.97057 -118.07057 0.08% - 2s
0 0 -118.07057 0 13 -117.97057 -118.07057 0.08% - 2s
Cutting planes:
Implied bound: 17
MIR: 8
Inf proof: 1
Zero half: 1
RLT: 12
Relax-and-lift: 4
Explored 1 nodes (11239 simplex iterations) in 2.68 seconds (2.40 work units)
Thread count was 12 (of 12 available processors)
Solution count 10: -117.971 -117.971 -117.871 ... -76.7465
No other solutions better than -117.971
Optimal solution found (tolerance 1.00e-04)
Best objective -1.179705689649e+02, best bound -1.179705689649e+02, gap 0.0000%
-------------------------------------------------------------------------------
Summary
-------------------------------------------------------------------------------
error_sample1_0 0.0
regularization_edge_has_signal 54.0
Sample 1/30: Objective value = 5.2910
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 2/30: Objective value = 5.3342
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 3/30: Objective value = 5.2666
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 4/30: Objective value = 5.2464
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 5/30: Objective value = 5.3448
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 6/30: Objective value = 5.1948
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 7/30: Objective value = 5.3077
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 8/30: Objective value = 5.2081
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 9/30: Objective value = 5.3122
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 10/30: Objective value = 5.2652
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 11/30: Objective value = 5.2137
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 12/30: Objective value = 5.3676
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 13/30: Objective value = 5.2529
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 14/30: Objective value = 5.4194
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 15/30: Objective value = 5.2722
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 16/30: Objective value = 5.1633
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 17/30: Objective value = 5.4053
error_sample1_0: -0.0000, rel error: 315.2032
> Rejected! (315.2032 >= 0.05)
Sample 18/30: Objective value = 5.2273
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 19/30: Objective value = 5.3283
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 20/30: Objective value = 5.1877
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 21/30: Objective value = 5.2407
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 22/30: Objective value = 5.3266
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 23/30: Objective value = 5.3726
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 24/30: Objective value = 5.3392
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 25/30: Objective value = 5.2540
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 26/30: Objective value = 5.2754
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 27/30: Objective value = 5.2294
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 28/30: Objective value = 5.2531
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 29/30: Objective value = 5.3124
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 30/30: Objective value = 5.3054
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 55.0000, rel error: 0.0185
Generated 29 alternative solutions
| 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | ... | 19 | 20 | 21 | 22 | 23 | 24 | 25 | 26 | 27 | 28 | ||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| (SMAD3) | (MYOD1) | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| (GRK2) | (BDKRB2) | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| (MAPK14) | (MAPKAPK2) | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| (DEPTOR_EEF1A1_MLST8_MTOR_PRR5_RICTOR) | (FBXW8) | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| (SLK) | (MAP3K5) | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
5 rows × 29 columns
# Analyze variability in edge selection
df_var = pd.concat([df_sols.mean(axis=1), df_sols.std(axis=1)], axis=1)
df_var.columns = ["mean", "std"]
print("Edges with highest variability across solutions:")
df_var.sort_values(by="std", ascending=False).head(50).sort_values(by="mean")
Edges with highest variability across solutions:
| mean | std | ||
|---|---|---|---|
| (TBK1) | (AKT1) | 0.172414 | 0.384426 |
| (PPP1CA) | (AKT1) | 0.172414 | 0.384426 |
| (MAPK14) | (SMAD3) | 0.172414 | 0.384426 |
| (CASP3) | 0.206897 | 0.412251 | |
| (PRKCA) | (TP53) | 0.206897 | 0.412251 |
| (PPP1CA) | (TP53) | 0.206897 | 0.412251 |
| (TBK1) | (RELA) | 0.206897 | 0.412251 |
| (TAB1) | (ROR2) | 0.206897 | 0.412251 |
| (MAPK8) | (TP53) | 0.206897 | 0.412251 |
| (PRKCD) | (MAPK8) | 0.206897 | 0.412251 |
| (MAP3K7) | (MAP2K3) | 0.206897 | 0.412251 |
| (RELA) | (JUN) | 0.206897 | 0.412251 |
| (MAP2K3) | (MAPK14) | 0.206897 | 0.412251 |
| (TRAF6) | (TAB1) | 0.206897 | 0.412251 |
| (MAPK8) | (MYC) | 0.241379 | 0.435494 |
| (TRAF2) | (TICAM1) | 0.241379 | 0.435494 |
| (PPP2CA) | (TRAF2) | 0.241379 | 0.435494 |
| (GSK3B) | (RELA) | 0.241379 | 0.435494 |
| (AIMP1) | (SMURF2) | 0.241379 | 0.435494 |
| (MAPK8) | (AKT1) | 0.241379 | 0.435494 |
| (AIMP1) | 0.241379 | 0.435494 | |
| (SMURF2) | (SMAD5) | 0.241379 | 0.435494 |
| (ROR2) | (MAPK8) | 0.241379 | 0.435494 |
| (MAP2K6) | (MAPK14) | 0.241379 | 0.435494 |
| (MAP3K7) | (MAP2K6) | 0.241379 | 0.435494 |
| (TBK1) | (NFKB1_RELA) | 0.275862 | 0.454859 |
| (MAPK14) | (CDKN1A) | 0.275862 | 0.454859 |
| (CDK1) | (TP53) | 0.310345 | 0.470824 |
| (MAPK1) | (JUN) | 0.379310 | 0.493804 |
| (MAP3K7) | (MAP2K4) | 0.413793 | 0.501230 |
| (MAP2K4) | (MAPK14) | 0.413793 | 0.501230 |
| (CDK1) | (AKAP12) | 0.448276 | 0.506120 |
| (PHLPP1) | (PRKCA) | 0.448276 | 0.506120 |
| (GSK3B) | (PHLPP1) | 0.448276 | 0.506120 |
| (SMAD3) | (SMAD4) | 0.448276 | 0.506120 |
| (AKAP12) | (PRKCA) | 0.448276 | 0.506120 |
| (RELA) | (EGR1) | 0.482759 | 0.508548 |
| (NFKB1_RELA) | (EGR1) | 0.517241 | 0.508548 |
| (SMAD2) | (SMAD4) | 0.551724 | 0.506120 |
| (GSK3B) | (CDKN1A) | 0.689655 | 0.470824 |
| (PPP1CA) | (MAPK1) | 0.724138 | 0.454859 |
| (CDK1) | (PPP1CA) | 0.724138 | 0.454859 |
| (JUN) | (SPI1) | 0.758621 | 0.435494 |
| (GSK3B) | (PRKCE) | 0.758621 | 0.435494 |
| (PRKCE) | (TICAM2) | 0.758621 | 0.435494 |
| (MAPK8) | (STAT1) | 0.758621 | 0.435494 |
| (MAPK1) | (MYC) | 0.758621 | 0.435494 |
| (AKT1) | (SKI) | 0.758621 | 0.435494 |
| (SKI) | (SMAD5) | 0.758621 | 0.435494 |
| (TICAM2) | (TICAM1) | 0.758621 | 0.435494 |
# Analyze edge values across solutions
df_edge_val = pd.DataFrame(np.concatenate(edge_results["edge_values"], axis=1), index=m.processed_graph.E).astype(int)
df_edge_val.head()
| 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | ... | 19 | 20 | 21 | 22 | 23 | 24 | 25 | 26 | 27 | 28 | ||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| (SMAD3) | (MYOD1) | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| (GRK2) | (BDKRB2) | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| (MAPK14) | (MAPKAPK2) | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| (DEPTOR_EEF1A1_MLST8_MTOR_PRR5_RICTOR) | (FBXW8) | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| (SLK) | (MAP3K5) | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
5 rows × 29 columns
Vertex-Based Perturbation#
Next, we’ll sample solutions by perturbing vertex variables. This approach can highlight alternative active node sets that are consistent with the data.
# Create a fresh model
m = CarnivalFlow(lambda_reg=0.1)
P = m.build(dataset.graph, dataset.data)
# Sample solutions by perturbing vertex_value
vertex_results = sample_alternative_solutions(
problem=P,
target_variable=P.expr.vertex_value,
percentage=0.10,
scale=0.03,
rel_opt_tol=0.05,
max_samples=30,
)
# Analyze vertex-based results
df_vertex_val = pd.DataFrame(np.concatenate(vertex_results["vertex_values"], axis=1), index=m.processed_graph.V).astype(
int
)
print(f"Generated {df_vertex_val.shape[1]} alternative solutions")
df_vertex_val.head()
Unreachable vertices for sample: 0
===============================================================================
CVXPY
v1.7.1
===============================================================================
-------------------------------------------------------------------------------
Compilation
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Numerical solver
-------------------------------------------------------------------------------
Set parameter OutputFlag to value 1
Set parameter QCPDual to value 1
Gurobi Optimizer version 12.0.3 build v12.0.3rc0 (mac64[arm] - Darwin 25.5.0 25F84)
CPU model: Apple M4 Pro
Thread count: 12 physical cores, 12 logical processors, using up to 12 threads
Non-default parameters:
QCPDual 1
Optimize a model with 30607 rows, 10690 columns and 108040 nonzeros
Model fingerprint: 0xf48e83e6
Variable types: 4188 continuous, 6502 integer (6502 binary)
Coefficient statistics:
Matrix range [1e+00, 9e+02]
Objective range [1e-01, 7e+00]
Bounds range [1e+00, 1e+00]
RHS range [1e+00, 1e+03]
Found heuristic solution: objective 0.1000000
Presolve removed 17840 rows and 1466 columns
Presolve time: 0.20s
Presolved: 12767 rows, 9224 columns, 81311 nonzeros
Variable types: 3504 continuous, 5720 integer (5714 binary)
Root relaxation: objective -1.186277e+02, 1171 iterations, 0.04 seconds (0.07 work units)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 -118.62771 0 53 0.10000 -118.62771 - - 0s
0 0 -118.51700 0 62 0.10000 -118.51700 - - 0s
H 0 0 -3.7564221 -118.51700 3055% - 0s
H 0 0 -76.6465356 -118.51700 54.6% - 0s
H 0 0 -76.7465356 -118.51700 54.4% - 0s
0 0 -118.21224 0 60 -76.74654 -118.21224 54.0% - 0s
H 0 0 -79.0760121 -118.21224 49.5% - 0s
H 0 0 -113.3642782 -118.12057 4.20% - 0s
0 0 -118.12057 0 55 -113.36428 -118.12057 4.20% - 0s
0 0 -118.07057 0 41 -113.36428 -118.07057 4.15% - 0s
H 0 0 -113.4642782 -118.07057 4.06% - 0s
0 0 -118.07057 0 36 -113.46428 -118.07057 4.06% - 0s
0 0 -118.07057 0 33 -113.46428 -118.07057 4.06% - 0s
0 0 -118.07057 0 18 -113.46428 -118.07057 4.06% - 0s
0 0 -118.07057 0 18 -113.46428 -118.07057 4.06% - 0s
H 0 0 -117.5705690 -118.07057 0.43% - 0s
0 0 -118.07057 0 18 -117.57057 -118.07057 0.43% - 0s
0 0 -118.07057 0 13 -117.57057 -118.07057 0.43% - 0s
0 0 -118.07057 0 13 -117.57057 -118.07057 0.43% - 0s
0 0 -118.07057 0 2 -117.57057 -118.07057 0.43% - 1s
0 0 -118.07057 0 37 -117.57057 -118.07057 0.43% - 1s
0 0 -118.07057 0 25 -117.57057 -118.07057 0.43% - 1s
0 0 -118.07057 0 20 -117.57057 -118.07057 0.43% - 1s
0 0 -118.07057 0 2 -117.57057 -118.07057 0.43% - 1s
H 0 0 -117.8705690 -118.07057 0.17% - 1s
0 0 -118.07057 0 11 -117.87057 -118.07057 0.17% - 1s
0 0 -118.07057 0 11 -117.87057 -118.07057 0.17% - 1s
0 0 -118.07057 0 11 -117.87057 -118.07057 0.17% - 1s
0 0 -118.07057 0 44 -117.87057 -118.07057 0.17% - 1s
H 0 0 -117.9705690 -118.07057 0.08% - 1s
0 0 -118.07057 0 38 -117.97057 -118.07057 0.08% - 1s
0 0 -118.07057 0 31 -117.97057 -118.07057 0.08% - 1s
0 0 -118.07057 0 23 -117.97057 -118.07057 0.08% - 1s
0 0 -118.07057 0 22 -117.97057 -118.07057 0.08% - 1s
0 0 -118.07057 0 22 -117.97057 -118.07057 0.08% - 1s
0 0 -118.07057 0 22 -117.97057 -118.07057 0.08% - 1s
0 0 -118.07057 0 13 -117.97057 -118.07057 0.08% - 1s
0 0 -118.07057 0 13 -117.97057 -118.07057 0.08% - 1s
Cutting planes:
Implied bound: 17
MIR: 8
Inf proof: 1
Zero half: 1
RLT: 12
Relax-and-lift: 4
Explored 1 nodes (11239 simplex iterations) in 1.86 seconds (2.40 work units)
Thread count was 12 (of 12 available processors)
Solution count 10: -117.971 -117.971 -117.871 ... -76.7465
No other solutions better than -117.971
Optimal solution found (tolerance 1.00e-04)
Best objective -1.179705689649e+02, best bound -1.179705689649e+02, gap 0.0000%
-------------------------------------------------------------------------------
Summary
-------------------------------------------------------------------------------
error_sample1_0 0.0
regularization_edge_has_signal 54.0
Sample 1/30: Objective value = 5.3506
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 2/30: Objective value = 5.3926
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 3/30: Objective value = 5.4280
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 4/30: Objective value = 5.3630
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 5/30: Objective value = 5.4156
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 6/30: Objective value = 5.3873
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 7/30: Objective value = 5.3842
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 8/30: Objective value = 5.3906
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 9/30: Objective value = 5.3064
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 10/30: Objective value = 5.4230
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 11/30: Objective value = 5.3525
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 12/30: Objective value = 5.3631
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 13/30: Objective value = 5.3528
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 14/30: Objective value = 5.3961
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 15/30: Objective value = 5.3628
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 16/30: Objective value = 5.3982
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 17/30: Objective value = 5.3655
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 55.0000, rel error: 0.0185
Sample 18/30: Objective value = 5.2853
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 19/30: Objective value = 5.3232
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 20/30: Objective value = 5.3768
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 21/30: Objective value = 5.2673
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 22/30: Objective value = 5.3205
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 23/30: Objective value = 5.3125
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 55.0000, rel error: 0.0185
Sample 24/30: Objective value = 5.3489
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 25/30: Objective value = 5.4452
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 26/30: Objective value = 5.3848
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 27/30: Objective value = 5.2272
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 28/30: Objective value = 5.3814
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 29/30: Objective value = 5.3650
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Sample 30/30: Objective value = 5.2891
error_sample1_0: 0.0000, rel error: 0.0000
regularization_edge_has_signal: 54.0000, rel error: 0.0000
Generated 30 alternative solutions
| 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | ... | 20 | 21 | 22 | 23 | 24 | 25 | 26 | 27 | 28 | 29 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| SMAD2_SMAD4 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| NFKBIB | -1 | -1 | -1 | -1 | -1 | -1 | -1 | -1 | -1 | -1 | ... | -1 | -1 | -1 | -1 | -1 | -1 | -1 | -1 | -1 | -1 |
| MDC1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| BRSK1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| EREG | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
5 rows × 30 columns
# Focus on signaling proteins with no measurement data
df_sig_prot_pred = df_vertex_val.loc[df_vertex_val.index.difference(dataset.data.query.pluck_features())]
df_sig_prot_pred.head()
| 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | ... | 20 | 21 | 22 | 23 | 24 | 25 | 26 | 27 | 28 | 29 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| AAK1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| ABI1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| ABL1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| ABRAXAS1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| ACAT1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
5 rows × 30 columns
# Visualize the variability in predicted signaling protein activities
import matplotlib.pyplot as plt
plt.figure(figsize=(12, 6))
df_sig_prot_pred.std(axis=1).sort_values().tail(30).plot.bar(title="Variability in Signaling Protein Activity")
plt.xlabel("Protein")
plt.ylabel("Standard Deviation")
plt.tight_layout()
plt.show()
Summary#
In this notebook, we:
Created a generic sampling method for exploring alternative network solutions
Applied the method to both edge and vertex perturbations
Analyzed the variability in solutions to identify:
Edges with high uncertainty (appearing in some but not all solutions)
Signaling proteins with variable activity across solutions
This approach helps assess the robustness of network inference results and can guide further experimental validation efforts by highlighting the most uncertain parts of the network.