Daemontus/feat/set samplers (#50)

* feat: Add working samplers for VertexSet and ColorSet.

* feat: Add samplers for SpaceSet, ColoredVertexSet, and ColoredSpaceSet.

* feat: Add samplers for PerturbationSet and ColoredPerturbationSet.

* feat: Add example script for uniform random sampling.

Author: daemontus

Cargo.toml

@@ -28,13 +28,13 @@ serde = { version = "1.0", features = ["derive"] }
 serde_json = "1.0"
 regex = "1.12"
 pyo3 = { version = "0.27", features = ["abi3-py37", "extension-module", "num-bigint", "py-clone"] }
-biodivine-lib-param-bn = { version="0.7.0", features=["solver-z3", "serde"] }
+biodivine-lib-param-bn = { version="0.7.1", features=["solver-z3", "serde"] }
 biodivine-lib-bdd = { version = "0.6.2", features = ["serde"] }
 #biodivine-pbn-control = "0.3.1"
 biodivine-pbn-control = { git = "https://github.com/sybila/biodivine-pbn-control", rev = "9e31bab9d000266ea25c35cc752c69077973d43c" }
-biodivine-hctl-model-checker = "0.3.4"
+biodivine-hctl-model-checker = "0.3.5"
 biodivine-lib-io-bma = "0.1.2"
-biodivine-algo-bdd-scc = "0.3.1"
+biodivine-algo-bdd-scc = "0.3.2"
 # Cannot update `rand` because lib-bdd uses it to generate random big-int values.
 rand = "0.8"
 cancel-this = { version = "0.4.0", features = ["pyo3"] }

biodivine_aeon/__init__.pyi

@@ -1163,6 +1163,7 @@ class ColorSet:
     def extend_with_spaces(self, set: SpaceSet) -> ColoredSpaceSet: ...
     def extend_with_perturbations(self, set: PerturbationSet) -> ColoredPerturbationSet: ...
     def items(self, retained: Optional[Sequence[Union[VariableId, ParameterId, str]]] = None) -> Iterator[ColorModel]: ...
+    def sample_items(self, retained: Optional[Sequence[Union[VariableId, ParameterId, str]]] = None, seed: Optional[int] = None) -> Iterator[ColorModel]: ...
 
 class ColorModel:
     def __ctx__(self) -> SymbolicContext: ...
@@ -1218,6 +1219,7 @@ class VertexSet:
     def to_bdd(self) -> Bdd: ...
     def extend_with_colors(self, set: ColorSet) -> ColoredVertexSet: ...
     def items(self, retained: Optional[Sequence[Union[str, VariableId]]] = None) -> Iterator[VertexModel]: ...
+    def sample_items(self, retained: Optional[Sequence[Union[str, VariableId]]] = None, seed: Optional[int] = None) -> Iterator[VertexModel]: ...
     def enclosing_subspace(self) -> dict[VariableId, bool] | None: ...
     def enclosing_named_subspace(self) -> dict[str, bool] | None: ...
     def enclosed_subspace(self) -> dict[VariableId, bool] | None: ...
@@ -1269,6 +1271,7 @@ class SpaceSet:
     def to_bdd(self) -> Bdd: ...
     def extend_with_colors(self, set: ColorSet) -> ColoredSpaceSet: ...
     def items(self, retained: Optional[Sequence[Union[str, VariableId]]] = None) -> Iterator[SpaceModel]: ...
+    def sample_items(self, retained: Optional[Sequence[Union[str, VariableId]]] = None, seed: Optional[int] = None) -> Iterator[SpaceModel]: ...
     def with_all_sub_spaces(self) -> SpaceSet: ...
     def with_all_super_spaces(self) -> SpaceSet: ...
 
@@ -1329,6 +1332,11 @@ class ColoredVertexSet:
               retained_variables: Optional[Sequence[VariableIdType]] = None,
               retained_functions: Optional[Sequence[Union[VariableId, ParameterId, str]]] = None
               ) -> Iterator[tuple[ColorModel, VertexModel]]: ...
+    def sample_items(self,
+                     retained_variables: Optional[Sequence[VariableIdType]] = None,
+                     retained_functions: Optional[Sequence[Union[VariableId, ParameterId, str]]] = None,
+                     seed: Optional[int] = None
+                     ) -> Iterator[tuple[ColorModel, VertexModel]]: ...
 
 class ColoredSpaceSet:
     def __init__(self, ctx: SymbolicContext, bdd: Bdd):
@@ -1372,6 +1380,11 @@ class ColoredSpaceSet:
               retained_variables: Optional[Sequence[VariableIdType]] = None,
               retained_functions: Optional[Sequence[Union[VariableId, ParameterId, str]]] = None
               ) -> Iterator[tuple[ColorModel, SpaceModel]]: ...
+    def sample_items(self,
+                     retained_variables: Optional[Sequence[VariableIdType]] = None,
+                     retained_functions: Optional[Sequence[Union[VariableId, ParameterId, str]]] = None,
+                     seed: Optional[int] = None
+                     ) -> Iterator[tuple[ColorModel, SpaceModel]]: ...
     def with_all_sub_spaces(self) -> ColoredSpaceSet: ...
     def with_all_super_spaces(self) -> ColoredSpaceSet: ...
 
@@ -1740,6 +1753,10 @@ class PerturbationSet:
     def to_internal(self) -> ColoredVertexSet: ...
     def extend_with_colors(self, set: ColorSet) -> ColoredPerturbationSet: ...
     def items(self, retained: Optional[Sequence[VariableIdType]] = None) -> Iterator[PerturbationModel]: ...
+    def sample_items(self,
+                     retained: Optional[Sequence[VariableIdType]] = None,
+                     seed: Optional[int] = None
+                     ) -> Iterator[PerturbationModel]: ...
 
 class PerturbationModel:
     def __ctx__(self) -> AsynchronousPerturbationGraph: ...
@@ -1805,6 +1822,11 @@ class ColoredPerturbationSet:
               retained_variables: Optional[Sequence[VariableIdType]] = None,
               retained_functions: Optional[Sequence[Union[VariableId, ParameterId, str]]] = None
               ) -> Iterator[tuple[ColorModel, PerturbationModel]]: ...
+    def sample_items(self,
+                     retained_variables: Optional[Sequence[VariableIdType]] = None,
+                     retained_functions: Optional[Sequence[Union[VariableId, ParameterId, str]]] = None,
+                     seed: Optional[int] = None
+                     ) -> Iterator[tuple[ColorModel, PerturbationModel]]: ...
 
 class AsynchronousPerturbationGraph(AsynchronousGraph):
     def __init__(self, network: BooleanNetwork, perturb: Optional[Sequence[VariableIdType]] = None):

example/script/uniform_random_sampling.py

@@ -0,0 +1,45 @@
+from biodivine_aeon import *
+
+# This script showcases the ability to randomly uniformly sample
+# from symbolic sets. This property then allows building more complex
+# statistics workflows on top of it.
+
+bn = BooleanNetwork.from_file('../case-study/butanol-production/butanol-pathway-f1-f2-f3.aeon')
+graph = AsynchronousGraph(bn)
+
+print("Total number of interpretations:", graph.mk_unit_colors().cardinality())
+
+# The network has three unknown functions: f_1(2), f_2(2) and f_3(3).
+# First, let us find all colors that exhibit the `butanol` phenotype:
+
+butanol = graph.mk_subspace({ 'butanol': 1 })
+print("Colors where butanol=1: ", butanol.colors())
+butanol_phenotype = Reachability.forward_subset(graph, butanol)
+print("Colors where exists trap set with butanol=1: ", butanol_phenotype.colors())
+
+# Now, we want to test a hypothesis that f_1 and f_2 are independent of each other.
+# In particular, both f_1 and f_2 can be either conjunction or disjunction. We thus
+# want to test whether the choice in f_1 influences the choice in f_2.
+
+# (Note that we could also test this particular hypothesis exactly through
+# careful manipulation of the symbolic set, but it is a useful motivating example)
+
+color_sampler = butanol_phenotype.colors().sample_items(retained=['f_1', 'f_2'], seed=1)
+samples = 0
+same = 0
+different = 0
+while samples < 10000:
+    samples += 1
+    sample = next(color_sampler)
+    f_1_expr = sample['f_1']
+    f_2_expr = sample['f_2']
+    if f_1_expr.is_and() == f_2_expr.is_and():
+        same += 1
+    else:
+        different += 1
+    if samples % 1000 == 0:
+        print(f"{same/samples:.2f} / {different/samples:.2f}")
+
+# If all went according to plan, we should get a roughly 50/50 distribution in the end.
+
+# Note that similar random sampling functions are available on all symbolic sets.

src/bindings/lib_param_bn/symbolic/set_color.rs

@@ -2,7 +2,17 @@ use std::collections::hash_map::DefaultHasher;
 use std::hash::{Hash, Hasher};
 use std::ops::Not;
 
-use biodivine_lib_bdd::Bdd as RsBdd;
+use crate::AsNative;
+use crate::bindings::lib_bdd::bdd::Bdd;
+use crate::bindings::lib_param_bn::symbolic::model_color::ColorModel;
+use crate::bindings::lib_param_bn::symbolic::set_colored_space::ColoredSpaceSet;
+use crate::bindings::lib_param_bn::symbolic::set_colored_vertex::ColoredVertexSet;
+use crate::bindings::lib_param_bn::symbolic::set_spaces::SpaceSet;
+use crate::bindings::lib_param_bn::symbolic::set_vertex::VertexSet;
+use crate::bindings::lib_param_bn::symbolic::symbolic_context::SymbolicContext;
+use crate::bindings::pbn_control::{ColoredPerturbationSet, PerturbationSet};
+use biodivine_lib_bdd::random_sampling::UniformValuationSampler;
+use biodivine_lib_bdd::{Bdd as RsBdd, BddPartialValuation};
 use biodivine_lib_param_bn::biodivine_std::traits::Set;
 use biodivine_lib_param_bn::symbolic_async_graph::projected_iteration::{
     OwnedRawSymbolicIterator, RawProjection,
@@ -15,16 +25,8 @@ use pyo3::IntoPyObjectExt;
 use pyo3::basic::CompareOp;
 use pyo3::prelude::*;
 use pyo3::types::PyList;
-
-use crate::AsNative;
-use crate::bindings::lib_bdd::bdd::Bdd;
-use crate::bindings::lib_param_bn::symbolic::model_color::ColorModel;
-use crate::bindings::lib_param_bn::symbolic::set_colored_space::ColoredSpaceSet;
-use crate::bindings::lib_param_bn::symbolic::set_colored_vertex::ColoredVertexSet;
-use crate::bindings::lib_param_bn::symbolic::set_spaces::SpaceSet;
-use crate::bindings::lib_param_bn::symbolic::set_vertex::VertexSet;
-use crate::bindings::lib_param_bn::symbolic::symbolic_context::SymbolicContext;
-use crate::bindings::pbn_control::{ColoredPerturbationSet, PerturbationSet};
+use rand::SeedableRng;
+use rand::prelude::StdRng;
 
 /// A symbolic representation of a set of "colors", i.e., interpretations of explicit and
 /// implicit parameters within a particular `BooleanNetwork`.
@@ -44,6 +46,20 @@ pub struct _ColorModelIterator {
     retained_implicit: Vec<biodivine_lib_param_bn::VariableId>,
 }
 
+/// An internal class used for random uniform sampling of `ColorModel` instances from a `ColorSet`.
+/// On the Python side, it just looks like an infinite iterator.
+///
+/// Similar to [`_ColorModelIterator`], but uses a sampler to get valuations from the projection
+/// BDD instead of iterating it fully.
+#[pyclass(module = "biodivine_aeon")]
+pub struct _ColorModelSampler {
+    ctx: Py<SymbolicContext>,
+    projection: RawProjection,
+    sampler: UniformValuationSampler<StdRng>,
+    retained_explicit: Vec<biodivine_lib_param_bn::ParameterId>,
+    retained_implicit: Vec<biodivine_lib_param_bn::VariableId>,
+}
+
 impl AsNative<GraphColors> for ColorSet {
     fn as_native(&self) -> &GraphColors {
         &self.native
@@ -238,6 +254,75 @@ impl ColorSet {
     #[pyo3(signature = (retained = None))]
     pub fn items(&self, retained: Option<&Bound<'_, PyList>>) -> PyResult<_ColorModelIterator> {
         let ctx = self.ctx.get();
+        let (retained, implicit, explicit) = Self::read_retained_functions(ctx, retained)?;
+
+        let projection = RawProjection::new(retained, self.as_native().as_bdd());
+        Ok(_ColorModelIterator {
+            ctx: self.ctx.clone(),
+            native: projection.into_iter(),
+            retained_implicit: implicit,
+            retained_explicit: explicit,
+        })
+    }
+
+    /// Returns a sampler for random uniform sampling of colors from this `ColorSet` with an
+    /// optional projection to a subset of uninterpreted functions. **If a projection is specified,
+    /// the sampling is uniform with respect to the projected set.**
+    ///
+    /// See also the `items` method regarding the `retained` projection set.
+    ///
+    /// You can specify an optional seed to make the sampling random but deterministic.
+    #[pyo3(signature = (retained = None, seed = None))]
+    pub fn sample_items(
+        &self,
+        retained: Option<&Bound<'_, PyList>>,
+        seed: Option<u64>,
+    ) -> PyResult<_ColorModelSampler> {
+        let ctx = self.ctx.get();
+        let (retained, implicit, explicit) = Self::read_retained_functions(ctx, retained)?;
+        let projection = RawProjection::new(retained, self.as_native().as_bdd());
+        let rng = StdRng::seed_from_u64(seed.unwrap_or_default());
+        let sampler = projection.bdd().mk_uniform_valuation_sampler(rng);
+        Ok(_ColorModelSampler {
+            ctx: self.ctx.clone(),
+            projection,
+            sampler,
+            retained_explicit: explicit,
+            retained_implicit: implicit,
+        })
+    }
+}
+
+impl ColorSet {
+    pub fn mk_native(ctx: Py<SymbolicContext>, native: GraphColors) -> Self {
+        Self { ctx, native }
+    }
+
+    pub fn mk_derived(&self, native: GraphColors) -> ColorSet {
+        ColorSet {
+            ctx: self.ctx.clone(),
+            native,
+        }
+    }
+
+    pub fn semantic_eq(a: &ColorSet, b: &ColorSet) -> bool {
+        let a = a.as_native().as_bdd();
+        let b = b.as_native().as_bdd();
+        if a.num_vars() != b.num_vars() {
+            return false;
+        }
+
+        RsBdd::binary_op_with_limit(1, a, b, biodivine_lib_bdd::op_function::xor).is_some()
+    }
+
+    pub fn read_retained_functions(
+        ctx: &SymbolicContext,
+        retained: Option<&Bound<'_, PyList>>,
+    ) -> PyResult<(
+        Vec<biodivine_lib_bdd::BddVariable>,
+        Vec<biodivine_lib_param_bn::VariableId>,
+        Vec<biodivine_lib_param_bn::ParameterId>,
+    )> {
         let mut retained_explicit = Vec::new();
         let mut retained_implicit = Vec::new();
         let retained = if let Some(retained) = retained {
@@ -266,41 +351,12 @@ impl ColorSet {
         } else {
             retained_explicit.append(&mut ctx.as_native().network_parameters().collect::<Vec<_>>());
             retained_implicit.append(&mut ctx.as_native().network_implicit_parameters());
-            self.ctx.get().as_native().parameter_variables().clone()
+            ctx.as_native().parameter_variables().clone()
         };
         retained_explicit.sort();
         retained_implicit.sort();
 
-        let projection = RawProjection::new(retained, self.as_native().as_bdd());
-        Ok(_ColorModelIterator {
-            ctx: self.ctx.clone(),
-            native: projection.into_iter(),
-            retained_implicit,
-            retained_explicit,
-        })
-    }
-}
-
-impl ColorSet {
-    pub fn mk_native(ctx: Py<SymbolicContext>, native: GraphColors) -> Self {
-        Self { ctx, native }
-    }
-
-    pub fn mk_derived(&self, native: GraphColors) -> ColorSet {
-        ColorSet {
-            ctx: self.ctx.clone(),
-            native,
-        }
-    }
-
-    pub fn semantic_eq(a: &ColorSet, b: &ColorSet) -> bool {
-        let a = a.as_native().as_bdd();
-        let b = b.as_native().as_bdd();
-        if a.num_vars() != b.num_vars() {
-            return false;
-        }
-
-        RsBdd::binary_op_with_limit(1, a, b, biodivine_lib_bdd::op_function::xor).is_some()
+        Ok((retained, retained_implicit, retained_explicit))
     }
 }
 
@@ -326,3 +382,34 @@ impl _ColorModelIterator {
         self.__next__()
     }
 }
+
+#[pymethods]
+impl _ColorModelSampler {
+    fn __iter__(self_: Py<Self>) -> Py<Self> {
+        self_
+    }
+
+    fn __next__(&mut self) -> Option<ColorModel> {
+        self.projection
+            .bdd()
+            .random_valuation_sample(&mut self.sampler)
+            .map(|it| {
+                let mut retained = BddPartialValuation::empty();
+                for x in self.projection.retained_variables() {
+                    retained.set_value(*x, it[*x]);
+                }
+
+                ColorModel::new_native(
+                    self.ctx.clone(),
+                    retained,
+                    self.retained_implicit.clone(),
+                    self.retained_explicit.clone(),
+                )
+            })
+    }
+
+    #[allow(clippy::should_implement_trait)]
+    pub fn next(&mut self) -> Option<ColorModel> {
+        self.__next__()
+    }
+}