From c66e224345bb69461120227a7396f722193db6cf Mon Sep 17 00:00:00 2001
From: Morten Bormann Nielsen <mon@teknologisk.dk>
Date: Wed, 17 Dec 2025 10:27:44 +0100
Subject: [PATCH 01/19] Adds the Dirichlet Rescale Constraints (DRSC) algorithm
 for generating constrained points that sum to a constant value

---
 ProcessOptimizer/space/DRSC.py | 418 +++++++++++++++++++++++++++++++++
 1 file changed, 418 insertions(+)
 create mode 100644 ProcessOptimizer/space/DRSC.py

diff --git a/ProcessOptimizer/space/DRSC.py b/ProcessOptimizer/space/DRSC.py
new file mode 100644
index 00000000..49fc2f42
--- /dev/null
+++ b/ProcessOptimizer/space/DRSC.py
@@ -0,0 +1,418 @@
+from typing import List, Tuple, Callable, Optional
+
+import numpy as np
+
+from scipy.optimize import linprog
+
+
+class DRSCGenerator:
+    """
+    Dirichlet-Rescale-Constraints (DRSC) Algorithm
+    Generates n-dimensional vectors with fixed sum satisfying linear and nonlinear constraints.
+    
+    This implementation is based on the paper "Generating Random Vectors 
+    satisfying Linear and Nonlinear Constraints" by Rick S. H. Willemsen, 
+    Wilco van den Heuvel and Michel van de Velden. It can be found here:
+    https://arxiv.org/abs/2501.16936
+    """
+    
+    def __init__(
+        self, 
+        n: int,
+        bounds: np.ndarray,
+        linear_constraints: Optional[List[Tuple[np.ndarray, float]]] = None,
+        nonlinear_constraints: Optional[List[Callable]] = None,
+        epsilon_ineq: float = 1e-3,
+        epsilon_eq: float = 1e-2,
+        max_restarts: int = 10000
+    ):
+        """
+        Initialize the DRSC generator.
+        
+        Parameters
+        ----------
+        * `n` [int]: 
+            Dimension of vectors to be generated
+            
+        * `bounds` [np.ndarray]: 
+            Array of shape (n, 2) with [lower, upper] bounds for each dimension
+            (normalized to sum = 1 simplex)
+        
+        * `linear_constraints`: 
+            List of (a, b) tuples representing a^T x <= b
+            
+        * `nonlinear_constraints`: 
+            List of functions f(x) that should be <= 0
+            
+        * `epsilon_ineq`: 
+            Feasibility tolerance for inequality constraints
+            
+        * `epsilon_eq`:
+            Feasibility tolerance for equality constraints
+        
+        * `max_restarts`:
+            Maximum number of restarts before giving up
+        """
+        self.n = n
+        self.original_bounds = np.array(bounds)
+        self.linear_constraints = linear_constraints or []
+        self.nonlinear_constraints = nonlinear_constraints or []
+        self.epsilon_ineq = epsilon_ineq
+        self.epsilon_eq = epsilon_eq
+        self.max_restarts = max_restarts
+        
+        # Compute optimal dimension ordering based on constrained dimension width
+        self._sort_order, self._inverse_order = self._compute_optimal_order()
+        
+        # Reorder bounds and constraints according to sort order
+        self._sorted_bounds = self.original_bounds[self._sort_order]
+        self._sorted_linear_constraints = self._reorder_linear_constraints()
+        self._sorted_nonlinear_constraints = self._reorder_nonlinear_constraints()
+        
+        # Find induced simplices using sorted dimensions (Algorithm 3, line 1)
+        self._thetas = self._find_induced_simplices()
+
+    
+    def _compute_optimal_order(self) -> Tuple[np.ndarray, np.ndarray]:
+        """
+        Compute optimal dimension ordering based on width.
+        
+        Dimensions with smaller widths (less flexible) are processed first.
+        This maximizes sampling efficiency by leaving flexible dimensions
+        to absorb variation from random sampling.
+        
+        Returns:
+            sort_order: Indices to sort dimensions by width (ascending)
+            inverse_order: Indices to restore original dimension order
+        """
+        widths = self.original_bounds[:, 1] - self.original_bounds[:, 0]
+        sort_order = np.argsort(widths)
+        inverse_order = np.argsort(sort_order)
+        
+        return sort_order, inverse_order
+    
+    
+    def _reorder_linear_constraints(self) -> List[Tuple[np.ndarray, float]]:
+        """Reorder linear constraint coefficients according to sort order."""
+        reordered = []
+        for a, b in self.linear_constraints:
+            a_reordered = a[self._sort_order]
+            reordered.append((a_reordered, b))
+        return reordered
+    
+    
+    def _reorder_nonlinear_constraints(self) -> List[Callable]:
+        """Wrap nonlinear constraints to work with sorted dimensions."""
+        wrapped = []
+        for f in self.nonlinear_constraints:
+            def wrapped_f(x_sorted, func=f):
+                # Convert from sorted order to original order before evaluating
+                x_original = x_sorted[self._inverse_order]
+                return func(x_original)
+            wrapped.append(wrapped_f)
+        return wrapped
+    
+    
+    def _find_induced_simplices(self) -> np.ndarray:
+        """
+        Find induced simplices Si for each dimension using LP formulation (1)-(7).
+        Uses sorted dimension order for optimal efficiency.
+        
+        For each dimension i, find theta_i such that the induced simplex
+        S_i = {x in S : x_i >= theta_i} covers part of the infeasible region
+        while not overlapping with previous induced simplices.
+        
+        The key insight from the paper is that induced simplices should cover
+        the INFEASIBLE region, so that when a point lands there, we can 
+        transform it back to the standard simplex.
+        
+        Returns
+        -------
+        Array of theta values for each dimension
+        """
+        thetas = np.zeros(self.n)
+        
+        for i in range(self.n):
+            theta_i = self._solve_lp_for_dimension(i, thetas)
+            thetas[i] = theta_i
+        
+        return thetas
+    
+    
+    def _solve_lp_for_dimension(
+            self,
+            dim: int,
+            previous_thetas: np.ndarray
+        ) -> float:
+        """
+        Solve LP formulation (1)-(7) for dimension dim (in sorted order).
+        
+        The goal is to find the largest theta_i such that the induced simplex
+        S_i = {x : x_i >= theta_i, sum(x) = 1, x >= 0} lies within the 
+        infeasible region AND doesn't overlap with previous induced simplices.
+        
+        According to the paper, we maximize x_i subject to:
+        - sum(x) = 1 (on the simplex)
+        - 0 <= x <= 1 (standard simplex bounds)
+        - All linear constraints in J_lin are satisfied
+        - x_k < theta_k for all k < i (non-overlapping with previous simplices)
+        
+        Parameters
+        ----------
+            dim : int
+                Current dimension index (in sorted order)
+            previous_thetas : np.ndarray
+                Array of theta values from previous dimensions
+            
+        Returns
+        -------
+            theta_i: float
+                theta_i value for the induced simplex
+        """
+        # Objective: maximize x_i (we'll negate for minimization)
+        c = np.zeros(self.n)
+        c[dim] = -1  # Negative because linprog minimizes
+        
+        # Equality constraint: sum(x) = 1
+        A_eq = np.ones((1, self.n))
+        b_eq = np.array([1.0])
+        
+        # Inequality constraints
+        A_ub = []
+        b_ub = []
+        
+        # Box constraints: l <= x <= u (from the sorted bounds)
+        for i in range(self.n):
+            # Lower bound: x[i] >= lower  →  -x[i] <= -lower
+            a_lower = np.zeros(self.n)
+            a_lower[i] = -1.0
+            A_ub.append(a_lower)
+            b_ub.append(-self._sorted_bounds[i, 0] + self.epsilon_ineq)
+            
+            # Upper bound: x[i] <= upper
+            a_upper = np.zeros(self.n)
+            a_upper[i] = 1.0
+            A_ub.append(a_upper)
+            b_ub.append(self._sorted_bounds[i, 1] + self.epsilon_ineq)
+        
+        # Additional linear constraints (already reordered)
+        for a, b in self._sorted_linear_constraints:
+            A_ub.append(a)
+            b_ub.append(b + self.epsilon_ineq)
+        
+        # Non-overlapping constraints: x_k < theta_k for k < dim
+        # This ensures the new induced simplex doesn't overlap with previous ones
+        # We enforce x[k] <= theta_k - epsilon (strictly less than)
+        for k in range(dim):
+            if previous_thetas[k] > 0:
+                constraint = np.zeros(self.n)
+                constraint[k] = 1.0
+                A_ub.append(constraint)
+                # x[k] <= theta_k - small_epsilon to ensure non-overlap
+                b_ub.append(previous_thetas[k] - 1e-9)
+        
+        # Convert to arrays
+        A_ub = np.array(A_ub) if A_ub else None
+        b_ub = np.array(b_ub) if b_ub else None
+        
+        # Bounds: 0 <= x_i <= 1
+        bounds = [(0, 1) for _ in range(self.n)]
+        
+        # Solve LP
+        result = linprog(
+            c, 
+            A_ub=A_ub, 
+            b_ub=b_ub, 
+            A_eq=A_eq, 
+            b_eq=b_eq, 
+            bounds=bounds, 
+            method='highs',
+        )
+        
+        if result.success:
+            return result.x[dim]
+        else:
+            return 0.0
+    
+    
+    def _get_simplex_vertices(self, dim: int, theta: float) -> np.ndarray:
+        """
+        Get vertices of the induced simplex for dimension dim with threshold theta.
+        
+        The induced simplex is defined by x_dim >= theta and sum(x) = 1, 0 <= x <= 1.
+        """        
+        induced_vertices = []
+        
+        # Vertex where x_dim = 1, all others = 0
+        v = np.zeros(self.n)
+        v[dim] = 1.0
+        induced_vertices.append(v)
+        
+        # Vertices where x_dim = theta, distributed among other dimensions
+        for j in range(self.n):
+            if j != dim:
+                v = np.zeros(self.n)
+                v[dim] = theta
+                v[j] = 1.0 - theta
+                induced_vertices.append(v)
+        
+        return np.array(induced_vertices)
+    
+    
+    def _sample_flat_dirichlet(self) -> np.ndarray:
+        """
+        Sample from flat Dirichlet distribution (uniform on standard simplex).
+        """
+        y = np.random.exponential(1.0, self.n)
+        return y / np.sum(y)
+    
+    
+    def _check_all_constraints(self, x: np.ndarray) -> bool:
+        """
+        Check if x satisfies ALL constraints.
+        
+        For the DRSC algorithm to work correctly, we must check:
+        1. Box bounds (essential - feasible region is subset of standard simplex)
+        2. Additional linear constraints (if any)
+        3. Nonlinear constraints (if any)
+        
+        Parameters
+        ----------
+        x : np.ndarray
+            Point to check (in sorted dimension order)
+            
+        Returns
+        -------
+        bool
+            True if all constraints are satisfied
+        """
+        # Check box bounds - ESSENTIAL
+        # The standard simplex is [0,1]^n with sum=1, but our feasible region
+        # is defined by tighter bounds
+        for i in range(self.n):
+            if x[i] < self._sorted_bounds[i, 0] - self.epsilon_ineq:
+                return False
+            if x[i] > self._sorted_bounds[i, 1] + self.epsilon_ineq:
+                return False
+        
+        # Check additional linear constraints (if any)
+        for a, b in self._sorted_linear_constraints:
+            if np.dot(a, x) > b + self.epsilon_ineq:
+                return False
+        
+        # Check nonlinear constraints (if any)
+        for f in self._sorted_nonlinear_constraints:
+            if f(x) > self.epsilon_ineq:
+                return False
+        
+        return True
+    
+    
+    def _is_in_induced_simplex(self, x: np.ndarray) -> Optional[int]:
+        """
+        Check if x is in any induced simplex.
+        
+        Returns:
+            Index of the induced simplex, or None if not in any
+        """
+        for idx, (theta, vertices, dim) in enumerate(self.induced_simplices):
+            if x[dim] >= theta:
+                return idx
+        return None
+    
+    
+    def _find_containing_induced_simplex(self, x: np.ndarray) -> Optional[int]:
+        """
+        Find which induced simplex contains point x.
+        
+        An induced simplex S_i is defined as {x in S : x_i >= theta_i}.
+        Due to the non-overlapping construction, a point can be in at most
+        one induced simplex.
+        
+        Returns
+        -------
+            Index of the induced simplex containing x, or None if not in any
+        """
+        for i in range(self.n):
+            if self._thetas[i] > 0 and x[i] >= self._thetas[i]:
+                return i
+        return None
+    
+    
+    def _affine_transform_to_standard_simplex(self, x: np.ndarray, dim: int) -> np.ndarray:
+        """
+        Apply affine transformation from induced simplex S_dim to standard simplex.
+        
+        Parameters
+        ----------
+            x: Point in induced simplex (sorted order)
+            dim: Index of the induced simplex dimension
+            
+        Returns
+        -------
+            Transformed point in standard simplex (sorted order)
+        """
+        theta = self._thetas[dim]
+        
+        # Create the transformed point
+        x_new = x.copy()
+        x_new[dim] = x[dim] - theta
+        
+        # Rescale to sum to 1
+        # The sum is now 1 - theta, so we divide by (1 - theta)
+        scale = 1.0 - theta
+        if scale > 1e-10:
+            x_new = x_new / scale
+        
+        return x_new
+    
+    
+    def generate(self) -> np.ndarray:
+        """
+        Generate a single vector using the DRSC algorithm (Algorithm 3).
+        
+        Returns
+        -------
+            n-dimensional vector with sum=1 satisfying all constraints.
+            Vector is returned in ORIGINAL dimension order
+        """
+        restarts = 0
+        
+        while restarts < self.max_restarts:
+            # Line 2: Sample from flat Dirichlet (in sorted dimension space)
+            x = self._sample_flat_dirichlet()
+
+            while True:
+                # Line 3: Check if all constraints satisfied
+                if self._check_all_constraints(x):
+                    return x[self._inverse_order]
+                
+                # Line 5: Check if x is inside an induced simplex
+                simplex_idx = self._find_containing_induced_simplex(x)
+                
+                if simplex_idx is not None:
+                    # Line 6: Apply affine transformation
+                    x = self._affine_transform_to_standard_simplex(x, simplex_idx)
+                else:
+                    # Line 8: Restart - violated constraint but not in any induced simplex
+                    restarts += 1
+                    break
+        
+        raise RuntimeError(
+            f"Failed to generate valid vector after {self.max_restarts} restarts. "
+            f"The constraint region may be very small or empty."
+        )
+    
+    def generate_sample(self, size: int) -> np.ndarray:
+        """
+        Generate multiple vectors.
+        
+        Parameters
+        ----------
+            size: Number of vectors to generate
+            
+        Returns
+        -------
+            Array of shape (size, n) with vectors in ORIGINAL dimension order
+        """
+        return np.array([self.generate() for _ in range(size)])
\ No newline at end of file

From b8e4b12512058b23794ee5a13b3f0ab6554ee02e Mon Sep 17 00:00:00 2001
From: Morten Bormann Nielsen <mon@teknologisk.dk>
Date: Wed, 17 Dec 2025 10:29:17 +0100
Subject: [PATCH 02/19] Expands the SumEquals constraint type to use the new
 DRSC algorithm by default

---
 ProcessOptimizer/space/constraints.py | 203 +++++++++++++++++++++++++-
 1 file changed, 197 insertions(+), 6 deletions(-)

diff --git a/ProcessOptimizer/space/constraints.py b/ProcessOptimizer/space/constraints.py
index 760bea09..7c8ddf4a 100644
--- a/ProcessOptimizer/space/constraints.py
+++ b/ProcessOptimizer/space/constraints.py
@@ -1,12 +1,17 @@
-from sklearn.utils import check_random_state
-from .space import Real, Integer, Categorical, Space, Task
+from typing import Union, List, Optional, Callable, Tuple
+
 import numpy as np
+
+from sklearn.utils import check_random_state
 from scipy import linalg
-from typing import Union, List
+
+from .space import Real, Integer, Categorical, Space, Task
+from .DRSC import DRSCGenerator
+
 
 class Constraints:
     def __init__(self, constraints_list, space):
-        """Constraints used when sampling for the aqcuisiiton function
+        """Constraints used when sampling for the aqcuisition function
 
         Parameters
         ----------
@@ -138,6 +143,8 @@ def sumequal_sampling(
 
         The samples are in the original space. They need to be transformed
         before being passed to a model or minimizer by `space.transform()`.
+        
+        Uses DRSC algorithm if available, falls back to null-space method.
 
         Parameters
         ----------
@@ -154,6 +161,76 @@ def sumequal_sampling(
            Points sampled from the space.
         """
         
+        # Check if we should use DRSC
+        if (len(self.sum_equals) == 1 and 
+            self.sum_equals[0].sampler == 'DRSC'):
+            
+            return self._drsc_sampling(n_samples, random_state)
+        else:
+            # Use existing null-space method
+            return self._nullspace_sampling(n_samples, random_state)
+        
+    def _drsc_sampling(self, n_samples, random_state):
+        """Sample using DRSC algorithm."""
+        rng = check_random_state(random_state)
+        constraint = self.sum_equals[0]
+        constrained_dims = constraint.dimensions
+        d = len(constrained_dims)
+        # Get bounds for constrained dimensions
+        bounds = np.array([self.space.bounds[dim] for dim in constrained_dims])
+        
+        # Get or create DRSC generator (handles sorting internally)
+        drsc_gen = constraint.get_drsc_generator(self.space)
+        
+        # Set numpy random seed for DRSC
+        np.random.seed(rng.randint(0, 2**31))
+        
+        samples = []
+        # Tolerance for clipping
+        eps = 1e-10
+        
+        for _ in range(n_samples):
+            # Generate sample on simplex
+            x_simplex = drsc_gen.generate()
+            
+            # Convert to original space
+            x_constrained = x_simplex * constraint.value
+            
+            # Clip to bounds (handles numerical tolerance)
+            x_clipped = np.clip(x_constrained, bounds[:, 0] - eps, bounds[:, 1] + eps)
+            
+            # Ensure exact bounds (remove the epsilon slack)
+            x_clipped = np.clip(x_clipped, bounds[:, 0], bounds[:, 1])
+            
+            # Verify sum constraint is still approximately satisfied
+            current_sum = np.sum(x_clipped)
+            if not np.isclose(current_sum, constraint.value, rtol=1e-6):
+                # Rescale proportionally to restore sum
+                x_clipped = x_clipped * (constraint.value / current_sum)
+                
+            # Final clip to ensure we are in bounds
+            x_constrained = np.clip(x_clipped, bounds[:, 0], bounds[:, 1])        
+                        
+            # Build full sample
+            if d < self.space.n_dims:
+                # Generate random settings across all factors
+                full_sample = self.space.rvs(n_samples=1, random_state=rng)[0]
+                
+                # Overwrite the random setting values for the constrained factors
+                for i, dim_idx in enumerate(constrained_dims):
+                    full_sample[dim_idx] = x_constrained[i]
+                    
+                samples.append(full_sample)
+            else:
+                samples.append(x_constrained.tolist())
+        
+        return samples
+        
+    def _nullspace_sampling(self, n_samples, random_state):
+        """
+        Original null-space based sampling method. Consider deprecating
+        immediately before the full pull request.
+        """
         def null_space(A, rcond=None) -> np.ndarray:
             """Helper function to calculate the null space of a matrix
 
@@ -663,7 +740,14 @@ def __eq__(self, other):
             return False
 
 class SumEquals():
-    def __init__(self, dimensions: List[int], value: Union[float, int]):
+    def __init__(
+        self, 
+        dimensions: List[int], 
+        value: Union[float, int],
+        sampler: str = "DRSC",
+        linear_constraints: Optional[List[Tuple[np.ndarray, float]]] = None,
+        nonlinear_constraints: Optional[List[Callable]] = None,
+    ):
         """Constraint class of type SumEquals.
 
         This constraint enforces that the sum of all values drawn for the
@@ -678,6 +762,17 @@ def __init__(self, dimensions: List[int], value: Union[float, int]):
 
         * `value` [float or int]:
             The value for which the sum should be equal to.
+        
+        * `sampler` [str, default='DRSC']:
+            Sampling algorithm to use. Options:
+            - 'DRSC': Dirichlet-Rescale-Constraints (fast, handles additional constraints)
+            - 'nullspace': Original null-space based sampler
+            
+        * `linear_constraints` [list of (a, b) tuples, optional]:
+            Additional linear constraints a^T x <= b (only used with DRSC)
+            
+        * `nonlinear_constraints` [list of callables, optional]:
+            Additional nonlinear constraints f(x) <= 0 (only used with DRSC)
 
         """
 
@@ -694,9 +789,105 @@ def __init__(self, dimensions: List[int], value: Union[float, int]):
 
         self.dimensions = tuple(dimensions)
         self.value = value
+        self.sampler = sampler
+        self.linear_constraints = linear_constraints or []
+        self.nonlinear_constraints = nonlinear_constraints or []
+        # Lazy initialization of DRSC generator
+        self._drsc_generator = None
 
         self.validate_sample = self._validate_sample
-
+        
+    def get_drsc_generator(self, space) -> DRSCGenerator:
+        """
+        Get or create DRSC generator with space information.
+        
+        Parameters
+        ----------
+            space: Space object containing dimension bounds
+            
+        Returns
+        -------
+            DRSCGenerator instance configured for this constraint
+        """
+        if self._drsc_generator is None and self.sampler == 'DRSC':
+            # Extract bounds for constrained dimensions
+            bounds = [space.bounds[dim] for dim in self.dimensions]
+            
+            # Normalize bounds to simplex (sum=1)
+            normalized_bounds = [(low / self.value, high / self.value)  for low, high in bounds]
+            
+            # Create DRSC generator
+            self._drsc_generator = DRSCGenerator(
+                n=len(self.dimensions),
+                bounds=normalized_bounds,
+                linear_constraints=self._convert_linear_constraints(bounds),
+                nonlinear_constraints=self._wrap_nonlinear_constraints(bounds),
+            )
+        return self._drsc_generator
+    
+    def _convert_linear_constraints(self, bounds):
+        """
+        Convert user-provided linear constraints to normalized form for use 
+        in the DRSC format.
+        """
+        # DRSC works on unit simplex, need to normalize constraints
+        converted = []
+        
+        # Add box constraints from dimension bounds
+        for i, (low, high) in enumerate(bounds):
+            # x_i >= low becomes constraint on simplex
+            # After scaling: x_simplex[i] * value >= low
+            # In simplex coords: -x[i] <= -low/value
+            a = np.zeros(len(self.dimensions))
+            a[i] = -1.0
+            b = -low / self.value
+            converted.append((a, b))
+            
+            # x_i <= high becomes: x[i] <= high/value
+            a = np.zeros(len(self.dimensions))
+            a[i] = 1.0
+            b = high / self.value
+            converted.append((a, b))
+        
+        # Add user-provided linear constraints (if any)
+        for a, b in self.linear_constraints:
+            # Extract only the constrained dimensions
+            a_constrained = np.array([a[dim] for dim in self.dimensions])
+            # Normalize by the sum value
+            a_normalized = a_constrained / self.value
+            b_normalized = b / self.value
+            converted.append((a_normalized, b_normalized))
+        
+        return converted
+    
+    def _wrap_nonlinear_constraints(self, bounds):
+        """Wrap nonlinear constraints to work on simplex."""
+        if not self.nonlinear_constraints:
+            return []
+        
+        wrapped = []
+        for f in self.nonlinear_constraints:
+            def wrapped_f(x_simplex):
+                # Convert from simplex to original space
+                x_original = self._simplex_to_original(x_simplex, bounds)
+                return f(x_original)
+            wrapped.append(wrapped_f)
+        return wrapped
+    
+    def _simplex_to_original(self, x_simplex, bounds):
+        """Convert from unit simplex coordinates to original space."""
+        x_original = np.zeros(len(self.dimensions))
+        for i, (low, high) in enumerate(bounds):
+            # Scale from simplex (sums to 1) to original space (sums to value)
+            scaled_value = x_simplex[i] * self.value
+            # Clip to dimension bounds to handle numerical tolerance and avoid errors
+            x_original[i] = np.clip(scaled_value, low, high)
+        return x_original
+    
+    def _original_to_simplex(self, x_original):
+        """Convert from original space to unit simplex."""
+        return np.array(x_original) / self.value
+    
     def _validate_sample(self, sample: List[int]) -> bool:
         # Returns True if sample does not violate the constraints aside from 
         # floating point errors

From adff4b29974e6a752423fa41f97136e2f1fb6864 Mon Sep 17 00:00:00 2001
From: Morten Bormann Nielsen <mon@teknologisk.dk>
Date: Wed, 17 Dec 2025 11:26:05 +0100
Subject: [PATCH 03/19] Adjusts tolerance in the test of factor sum constraint
 check to match the default tolerance of the DRSC algorithm

---
 ProcessOptimizer/tests/test_constraints.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/ProcessOptimizer/tests/test_constraints.py b/ProcessOptimizer/tests/test_constraints.py
index 5124d96f..ac8815d0 100644
--- a/ProcessOptimizer/tests/test_constraints.py
+++ b/ProcessOptimizer/tests/test_constraints.py
@@ -228,7 +228,8 @@ def test_SumEquals():
     samples = cons.sumequal_sampling(n_samples=1000)
     for sample in samples:
         factor_sum = np.sum(sample[0] + sample[1])
-        assert np.isclose(factor_sum, cons.sum_equals[0].value)
+        # The default tightness of the SumEquals sampling should be within 1e-4
+        assert np.isclose(factor_sum, cons.sum_equals[0].value, rtol=1e-4)
     # Check that the samples also have settings for the other dimensions
     assert len(samples[0]) == 4 and len(samples[999]) == 4
     # Check that the other dimensions have correct type

From cb40d5ef55cb0614445e65874de7807231644b03 Mon Sep 17 00:00:00 2001
From: Morten Bormann Nielsen <mon@teknologisk.dk>
Date: Wed, 17 Dec 2025 11:26:53 +0100
Subject: [PATCH 04/19] Updates function documentation format to be consistent
 with function descriptions across the package

---
 ProcessOptimizer/space/DRSC.py | 48 ++++++++++++++++++++--------------
 1 file changed, 29 insertions(+), 19 deletions(-)

diff --git a/ProcessOptimizer/space/DRSC.py b/ProcessOptimizer/space/DRSC.py
index 49fc2f42..1a641eca 100644
--- a/ProcessOptimizer/space/DRSC.py
+++ b/ProcessOptimizer/space/DRSC.py
@@ -81,9 +81,12 @@ def _compute_optimal_order(self) -> Tuple[np.ndarray, np.ndarray]:
         This maximizes sampling efficiency by leaving flexible dimensions
         to absorb variation from random sampling.
         
-        Returns:
-            sort_order: Indices to sort dimensions by width (ascending)
-            inverse_order: Indices to restore original dimension order
+        Returns
+        -------
+            `sort_order` [np.ndarray]: 
+                Indices to sort dimensions by width (ascending)
+            `inverse_order`[np.ndarray]: 
+                Indices to restore original dimension order
         """
         widths = self.original_bounds[:, 1] - self.original_bounds[:, 0]
         sort_order = np.argsort(widths)
@@ -128,7 +131,8 @@ def _find_induced_simplices(self) -> np.ndarray:
         
         Returns
         -------
-        Array of theta values for each dimension
+            `thetas` [np.ndarray]:
+                Array of theta values for each dimension
         """
         thetas = np.zeros(self.n)
         
@@ -159,14 +163,14 @@ def _solve_lp_for_dimension(
         
         Parameters
         ----------
-            dim : int
+            `dim` [int]:
                 Current dimension index (in sorted order)
-            previous_thetas : np.ndarray
+            `previous_thetas` [np.ndarray]:
                 Array of theta values from previous dimensions
             
         Returns
         -------
-            theta_i: float
+            `theta_i` [float]:
                 theta_i value for the induced simplex
         """
         # Objective: maximize x_i (we'll negate for minimization)
@@ -278,12 +282,12 @@ def _check_all_constraints(self, x: np.ndarray) -> bool:
         
         Parameters
         ----------
-        x : np.ndarray
+        * `x` [np.ndarray]:
             Point to check (in sorted dimension order)
             
         Returns
         -------
-        bool
+        * bool
             True if all constraints are satisfied
         """
         # Check box bounds - ESSENTIAL
@@ -312,8 +316,9 @@ def _is_in_induced_simplex(self, x: np.ndarray) -> Optional[int]:
         """
         Check if x is in any induced simplex.
         
-        Returns:
-            Index of the induced simplex, or None if not in any
+        Returns
+        -------
+            * Index of the induced simplex, or None if not in any
         """
         for idx, (theta, vertices, dim) in enumerate(self.induced_simplices):
             if x[dim] >= theta:
@@ -331,7 +336,7 @@ def _find_containing_induced_simplex(self, x: np.ndarray) -> Optional[int]:
         
         Returns
         -------
-            Index of the induced simplex containing x, or None if not in any
+            * Index of the induced simplex containing x, or None if not in any
         """
         for i in range(self.n):
             if self._thetas[i] > 0 and x[i] >= self._thetas[i]:
@@ -345,12 +350,15 @@ def _affine_transform_to_standard_simplex(self, x: np.ndarray, dim: int) -> np.n
         
         Parameters
         ----------
-            x: Point in induced simplex (sorted order)
-            dim: Index of the induced simplex dimension
+            * `x` [np.ndarray]: 
+                Point in induced simplex (sorted order)
+            * `dim` [int]: 
+                Index of the induced simplex dimension
             
         Returns
         -------
-            Transformed point in standard simplex (sorted order)
+            * `x_new` [np.ndarray]:
+                Transformed point in standard simplex (sorted order)
         """
         theta = self._thetas[dim]
         
@@ -373,8 +381,9 @@ def generate(self) -> np.ndarray:
         
         Returns
         -------
-            n-dimensional vector with sum=1 satisfying all constraints.
-            Vector is returned in ORIGINAL dimension order
+            * `x` [np.ndarray]:
+                n-dimensional vector with sum=1 satisfying all constraints. 
+                Vector is returned in ORIGINAL dimension order.
         """
         restarts = 0
         
@@ -409,10 +418,11 @@ def generate_sample(self, size: int) -> np.ndarray:
         
         Parameters
         ----------
-            size: Number of vectors to generate
+            * `size`[int]: 
+                Number of vectors to generate
             
         Returns
         -------
-            Array of shape (size, n) with vectors in ORIGINAL dimension order
+            * Array of shape (size, n) with vectors in ORIGINAL dimension order
         """
         return np.array([self.generate() for _ in range(size)])
\ No newline at end of file

From 69ebc8628fb342d3d9d33595dfb1439588a87d27 Mon Sep 17 00:00:00 2001
From: Morten Bormann Nielsen <mon@teknologisk.dk>
Date: Wed, 17 Dec 2025 15:05:52 +0100
Subject: [PATCH 05/19] Adds tests for asking for multiple new suggestions with
 SumEquals constraints, and tests of seeding of SumEquals sampling

---
 ProcessOptimizer/tests/test_constraints.py | 43 ++++++++++++++++++++++
 1 file changed, 43 insertions(+)

diff --git a/ProcessOptimizer/tests/test_constraints.py b/ProcessOptimizer/tests/test_constraints.py
index ac8815d0..8a338b52 100644
--- a/ProcessOptimizer/tests/test_constraints.py
+++ b/ProcessOptimizer/tests/test_constraints.py
@@ -277,6 +277,49 @@ def test_SumEquals():
     # The values suggested for the constrained dimensions should be the same 
     # irrespective of their order in the dimension list
     assert x1[:3] == x2[1:]
+    
+    # Test that we can ask for all initial points as long as we use strategy=cl_min
+    assert opt.ask(5, strategy="cl_min")
+    # We can also ask for points beyond the initial ones without data
+    assert opt.ask(6, strategy="cl_min")
+    # Give the optimizer some data and test that we can ask for multiple points
+    # after we have data
+    y = 0
+    for _ in range(10):
+        x = opt.ask()
+        opt.tell(x, y)
+        y -= 0.1
+    assert opt.ask(3)
+    
+    # Test that different seeds of the Optimizer provides different initial points
+    space = [
+        (50., 100.),
+        (1., 5.),
+        (3., 6.),
+        (1., 35.),
+        (1., 5.),
+        ("A", "B", "C"),
+        (1, 10),
+    ]
+    cons = [SumEquals(dimensions=[0, 1, 2, 3, 4], value=100.0, sampler="DRSC")]
+    # Build optimziers with different seeds
+    opt1 = Optimizer(space, lhs=False, n_initial_points=10, random_state=1)
+    opt1.set_constraints(cons)
+    opt2 = Optimizer(space, lhs=False, n_initial_points=10, random_state=2)
+    opt2.set_constraints(cons)
+    # Ask for the initial points
+    x1 = opt1.ask(10, strategy="cl_min")
+    x2 = opt2.ask(10, strategy="cl_min")
+    assert not x1 == x2
+    
+    # Check that seeding the generator directly inside each optimizer produces
+    # consistent points, despite the optimizers having different random states
+    constraint1 = opt1.get_constraints()
+    constraint2 = opt2.get_constraints()
+    # Ask for lots of points
+    x1 = constraint1.sumequal_sampling(n_samples=1000, random_state=31031988)
+    x2 = constraint2.sumequal_sampling(n_samples=1000, random_state=31031988)
+    assert x1 == x2
 
 @pytest.mark.fast_test
 def test_Conditional():

From e34a6b7522bfa1b6058546641d2acd332ed54396 Mon Sep 17 00:00:00 2001
From: Morten Bormann Nielsen <mon@teknologisk.dk>
Date: Thu, 18 Dec 2025 11:08:14 +0100
Subject: [PATCH 06/19] Remove unneeded function from first try at implementing
 the algorithm

---
 ProcessOptimizer/space/DRSC.py | 24 ------------------------
 1 file changed, 24 deletions(-)

diff --git a/ProcessOptimizer/space/DRSC.py b/ProcessOptimizer/space/DRSC.py
index 1a641eca..2ddca11b 100644
--- a/ProcessOptimizer/space/DRSC.py
+++ b/ProcessOptimizer/space/DRSC.py
@@ -239,30 +239,6 @@ def _solve_lp_for_dimension(
             return 0.0
     
     
-    def _get_simplex_vertices(self, dim: int, theta: float) -> np.ndarray:
-        """
-        Get vertices of the induced simplex for dimension dim with threshold theta.
-        
-        The induced simplex is defined by x_dim >= theta and sum(x) = 1, 0 <= x <= 1.
-        """        
-        induced_vertices = []
-        
-        # Vertex where x_dim = 1, all others = 0
-        v = np.zeros(self.n)
-        v[dim] = 1.0
-        induced_vertices.append(v)
-        
-        # Vertices where x_dim = theta, distributed among other dimensions
-        for j in range(self.n):
-            if j != dim:
-                v = np.zeros(self.n)
-                v[dim] = theta
-                v[j] = 1.0 - theta
-                induced_vertices.append(v)
-        
-        return np.array(induced_vertices)
-    
-    
     def _sample_flat_dirichlet(self) -> np.ndarray:
         """
         Sample from flat Dirichlet distribution (uniform on standard simplex).

From 13bf424b83bcedd1588bca0e0f91f3b6b60596e2 Mon Sep 17 00:00:00 2001
From: Morten Bormann Nielsen <mon@teknologisk.dk>
Date: Thu, 18 Dec 2025 11:09:11 +0100
Subject: [PATCH 07/19] Add a guardrail to the main algorithm to handle rare
 cases where the transformations cause the x-values to blow up outside the
 simplex

---
 ProcessOptimizer/space/DRSC.py | 12 ++++++++++--
 1 file changed, 10 insertions(+), 2 deletions(-)

diff --git a/ProcessOptimizer/space/DRSC.py b/ProcessOptimizer/space/DRSC.py
index 2ddca11b..adc2b893 100644
--- a/ProcessOptimizer/space/DRSC.py
+++ b/ProcessOptimizer/space/DRSC.py
@@ -341,7 +341,7 @@ def _affine_transform_to_standard_simplex(self, x: np.ndarray, dim: int) -> np.n
         # Create the transformed point
         x_new = x.copy()
         x_new[dim] = x[dim] - theta
-        
+    
         # Rescale to sum to 1
         # The sum is now 1 - theta, so we divide by (1 - theta)
         scale = 1.0 - theta
@@ -362,12 +362,14 @@ def generate(self) -> np.ndarray:
                 Vector is returned in ORIGINAL dimension order.
         """
         restarts = 0
+        max_transforms_per_sample = 1000
         
         while restarts < self.max_restarts:
             # Line 2: Sample from flat Dirichlet (in sorted dimension space)
             x = self._sample_flat_dirichlet()
+            transforms = 0
 
-            while True:
+            while transforms < max_transforms_per_sample:
                 # Line 3: Check if all constraints satisfied
                 if self._check_all_constraints(x):
                     return x[self._inverse_order]
@@ -378,6 +380,12 @@ def generate(self) -> np.ndarray:
                 if simplex_idx is not None:
                     # Line 6: Apply affine transformation
                     x = self._affine_transform_to_standard_simplex(x, simplex_idx)
+                    transforms += 1
+                    
+                    # Safety check: if values are getting too large, restart
+                    if np.any(np.abs(x) > 1e6):
+                        restarts += 1
+                        break  # Will trigger restart
                 else:
                     # Line 8: Restart - violated constraint but not in any induced simplex
                     restarts += 1

From 6cf95265d6b8ec2af4700feff044b76c168b5e66 Mon Sep 17 00:00:00 2001
From: Morten Bormann Nielsen <mon@teknologisk.dk>
Date: Thu, 18 Dec 2025 11:10:13 +0100
Subject: [PATCH 08/19] Adjusts how clipping is done to ensure that we both
 preserve the desired sum, but also that we don't violate our dimension bounds
 when we fix the sum. This can happen if the space is very unbalanced.

---
 ProcessOptimizer/space/constraints.py | 84 +++++++++++++++++++--------
 1 file changed, 61 insertions(+), 23 deletions(-)

diff --git a/ProcessOptimizer/space/constraints.py b/ProcessOptimizer/space/constraints.py
index 7c8ddf4a..a27cc40b 100644
--- a/ProcessOptimizer/space/constraints.py
+++ b/ProcessOptimizer/space/constraints.py
@@ -186,46 +186,84 @@ def _drsc_sampling(self, n_samples, random_state):
         np.random.seed(rng.randint(0, 2**31))
         
         samples = []
-        # Tolerance for clipping
-        eps = 1e-10
         
         for _ in range(n_samples):
-            # Generate sample on simplex
+            # Generate sample on simplex, then convert to original space
             x_simplex = drsc_gen.generate()
+            x_constrained = x_simplex * constraint.value              
             
-            # Convert to original space
-            x_constrained = x_simplex * constraint.value
-            
-            # Clip to bounds (handles numerical tolerance)
-            x_clipped = np.clip(x_constrained, bounds[:, 0] - eps, bounds[:, 1] + eps)
-            
-            # Ensure exact bounds (remove the epsilon slack)
-            x_clipped = np.clip(x_clipped, bounds[:, 0], bounds[:, 1])
-            
-            # Verify sum constraint is still approximately satisfied
-            current_sum = np.sum(x_clipped)
-            if not np.isclose(current_sum, constraint.value, rtol=1e-6):
-                # Rescale proportionally to restore sum
-                x_clipped = x_clipped * (constraint.value / current_sum)
-                
-            # Final clip to ensure we are in bounds
-            x_constrained = np.clip(x_clipped, bounds[:, 0], bounds[:, 1])        
+            # If the dimensions in the space have very different widths we risk 
+            # violating the space bounds due to the numerical tolerance of the
+            # DRSC algorithm. Fix this, while preserving the sum
+            x_constrained = self._fix_bounds_preserving_sum(
+                x_constrained, 
+                bounds, 
+                constraint.value
+            )
                         
             # Build full sample
             if d < self.space.n_dims:
                 # Generate random settings across all factors
                 full_sample = self.space.rvs(n_samples=1, random_state=rng)[0]
-                
                 # Overwrite the random setting values for the constrained factors
                 for i, dim_idx in enumerate(constrained_dims):
-                    full_sample[dim_idx] = x_constrained[i]
-                    
+                    full_sample[dim_idx] = x_constrained[i]   
                 samples.append(full_sample)
             else:
                 samples.append(x_constrained.tolist())
         
         return samples
+    
+   
+    def _fix_bounds_preserving_sum(self, x, bounds, target_sum, max_iterations=100):
+        """
+        Adjust x to satisfy bounds while preserving the sum constraint.
+        
+        When a value exceeds its bounds, the excess is redistributed to other
+        dimensions that have room to absorb it.
+        """
+        x = x.copy()
+        # Rescale proportionally to restore sum of x to target value
+        x = x * (target_sum / np.sum(x))
         
+        for _ in range(max_iterations):
+            all_satisfied = True
+            
+            for i in range(len(x)):
+                if x[i] < bounds[i, 0]:
+                    deficit = bounds[i, 0] - x[i]
+                    x[i] = bounds[i, 0]
+                    
+                    # Subtract from other dimensions that have room
+                    for j in range(len(x)):
+                        if j != i and deficit > 0:
+                            room = x[j] - bounds[j, 0]
+                            transfer = min(room, deficit)
+                            x[j] -= transfer
+                            deficit -= transfer
+                    
+                    all_satisfied = False
+                    
+                elif x[i] > bounds[i, 1]:
+                    excess = x[i] - bounds[i, 1]
+                    x[i] = bounds[i, 1]
+                    
+                    # Add to other dimensions that have room
+                    for j in range(len(x)):
+                        if j != i and excess > 0:
+                            room = bounds[j, 1] - x[j]
+                            transfer = min(room, excess)
+                            x[j] += transfer
+                            excess -= transfer
+                    
+                    all_satisfied = False
+            
+            if all_satisfied:
+                break
+        
+        return x
+
+    
     def _nullspace_sampling(self, n_samples, random_state):
         """
         Original null-space based sampling method. Consider deprecating

From 14e11c7647abb2cd4d42482dc2bd1ece9e23214e Mon Sep 17 00:00:00 2001
From: Morten Bormann Nielsen <mon@teknologisk.dk>
Date: Thu, 18 Dec 2025 11:12:14 +0100
Subject: [PATCH 09/19] Increased the number of samples to check during testing
 to catch rarer violations

---
 ProcessOptimizer/tests/test_constraints.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/ProcessOptimizer/tests/test_constraints.py b/ProcessOptimizer/tests/test_constraints.py
index 8a338b52..1fff742b 100644
--- a/ProcessOptimizer/tests/test_constraints.py
+++ b/ProcessOptimizer/tests/test_constraints.py
@@ -216,7 +216,7 @@ def test_SumEquals():
     assert not cons.validate_sample([3.00006, 2.0, "A"])    
     
     # Check again that only valid samples are drawn
-    samples = cons.sumequal_sampling(n_samples=1000)
+    samples = cons.sumequal_sampling(n_samples=10000)
     for sample in samples:
         factor_sum = np.sum(sample[0] + sample[1])
         assert np.isclose(factor_sum, cons.sum_equals[0].value)
@@ -225,7 +225,7 @@ def test_SumEquals():
     space = Space([[10.0, 20.0], [1000.0, 2000.0], [1, 5], ["A", "B"]])
     cons = Constraints([SumEquals((0, 1), 1234)], space)
     # Check again that only valid samples are drawn
-    samples = cons.sumequal_sampling(n_samples=1000)
+    samples = cons.sumequal_sampling(n_samples=10000)
     for sample in samples:
         factor_sum = np.sum(sample[0] + sample[1])
         # The default tightness of the SumEquals sampling should be within 1e-4
@@ -317,8 +317,8 @@ def test_SumEquals():
     constraint1 = opt1.get_constraints()
     constraint2 = opt2.get_constraints()
     # Ask for lots of points
-    x1 = constraint1.sumequal_sampling(n_samples=1000, random_state=31031988)
-    x2 = constraint2.sumequal_sampling(n_samples=1000, random_state=31031988)
+    x1 = constraint1.sumequal_sampling(n_samples=10000, random_state=31031988)
+    x2 = constraint2.sumequal_sampling(n_samples=10000, random_state=31031988)
     assert x1 == x2
 
 @pytest.mark.fast_test

From 87eb9d1907fc4ac4b6b3ba3edf59e52fbc294f89 Mon Sep 17 00:00:00 2001
From: Morten Bormann Nielsen <mon@teknologisk.dk>
Date: Thu, 18 Dec 2025 11:24:55 +0100
Subject: [PATCH 10/19] Fix small omission in test call of opt.ask(3) during
 SumEquals testing

---
 ProcessOptimizer/tests/test_constraints.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/ProcessOptimizer/tests/test_constraints.py b/ProcessOptimizer/tests/test_constraints.py
index 1fff742b..d459bb1f 100644
--- a/ProcessOptimizer/tests/test_constraints.py
+++ b/ProcessOptimizer/tests/test_constraints.py
@@ -289,7 +289,7 @@ def test_SumEquals():
         x = opt.ask()
         opt.tell(x, y)
         y -= 0.1
-    assert opt.ask(3)
+    assert opt.ask(3, strategy="cl_min")
     
     # Test that different seeds of the Optimizer provides different initial points
     space = [

From 3c5748db0a5ba02e7dfc90bd9a89b8ff3a131659 Mon Sep 17 00:00:00 2001
From: Morten Bormann Nielsen <mon@teknologisk.dk>
Date: Thu, 18 Dec 2025 12:51:55 +0100
Subject: [PATCH 11/19] Fixed a bunch of spelling mistakes through the file

---
 ProcessOptimizer/space/constraints.py | 31 ++++++++++++++-------------
 1 file changed, 16 insertions(+), 15 deletions(-)

diff --git a/ProcessOptimizer/space/constraints.py b/ProcessOptimizer/space/constraints.py
index a27cc40b..2ef1d4f1 100644
--- a/ProcessOptimizer/space/constraints.py
+++ b/ProcessOptimizer/space/constraints.py
@@ -11,7 +11,7 @@
 
 class Constraints:
     def __init__(self, constraints_list, space):
-        """Constraints used when sampling for the aqcuisition function
+        """Constraints used when sampling for the acquisition function
 
         Parameters
         ----------
@@ -123,7 +123,7 @@ def rvs(
 
             n_samples_candidates += n_samples
             if n_samples_candidates > 100000 and len(rows) < 100:
-                # If we have less than a 1/10.000 succes rate on the sampling
+                # If we have less than a 1/10.000 success rate on the sampling
                 # we throw an error
                 raise RuntimeError(
                     '''Could not find valid samples in constrained space.
@@ -435,7 +435,7 @@ def validate_sample(self, sample: List) -> bool:
         """
 
         # We iterate through all the dimensions and check the the type of
-        # constriants that are applied to a single dimensions, i.e Single,
+        # constraints that are applied to a single dimensions, i.e Single,
         # Exclusive and Inclusive
         #
         # We iterate through all samples which corresponds to number of
@@ -449,7 +449,7 @@ def validate_sample(self, sample: List) -> bool:
             # Inclusive constraints.
             # Check if there is a least one inclusive constraint:
             if self.inclusive[dim]:
-                # We go through all inlcusive constraints for this dimension
+                # We go through all inclusive constraints for this dimension
                 # and if the value is not found to be included in any of the
                 # bounds of the inclusive constraints we return false.
                 value_is_valid = False
@@ -467,11 +467,11 @@ def validate_sample(self, sample: List) -> bool:
                 if not constraint.validate_constraint(sample[dim]):
                     return False
 
-        # We iterate through sum constriants
+        # We iterate through sum constraints
         for constraint in self.sum:
             if not constraint.validate_sample(sample):
                 return False
-        # We iterate through sum_equals constriants
+        # We iterate through sum_equals constraints
         for constraint in self.sum_equals:
             if not constraint.validate_sample(sample):
                 return False
@@ -479,8 +479,7 @@ def validate_sample(self, sample: List) -> bool:
         for constraint in self.conditional:
             if not constraint.validate_sample(sample):
                 return False
-        # If we we did not find any violaiton of the constraints we return
-        # True.
+        # If we did not find any violation of the constraints we return True
         return True
 
     def __repr__(self):
@@ -611,7 +610,7 @@ def __init__(self, dimension, bounds, dimension_type):
             consist of floats.
 
             For 'categorical' dimensions the tuple must be of length < number
-            of dimensions and lenght > 1.
+            of dimensions and length > 1.
 
             The tuple can contain any combination of str, int or float
 
@@ -678,7 +677,7 @@ def __init__(self, dimension, bounds, dimension_type):
             consist of floats.
 
             For 'categorical' dimensions the tuple must be of length < number
-            of dimensions and lenght > 1.
+            of dimensions and length > 1.
 
             The tuple can contain any combination of str, int or float
 
@@ -731,7 +730,7 @@ def __init__(self, dimensions, value, less_than=True):
         Parameters
         ----------
         * `dimensions` [list of ints]:
-            A list of integers coresponding to the index of the dimensions that
+            A list of integers corresponding to the index of the dimensions that
             should be summed
 
         * `value` [float or int]:
@@ -1068,7 +1067,7 @@ def check_constraints(space, constraints):
                 raise IndexError('Dimension index exceeds number of dimensions')
             for ind_dim in constraint.dimensions:
                 if isinstance(space.dimensions[ind_dim], Categorical):
-                    raise ValueError('Sum constraint can not be applid to categorical dimension: {}'.format(space.dimensions[ind_dim]))
+                    raise ValueError('Sum constraint can not be applied to categorical dimension: {}'.format(space.dimensions[ind_dim]))
         elif isinstance(constraint, SumEquals):
             # Check that there is only one SumEquals constraint being applied
             if len(constraints) > 1:
@@ -1103,7 +1102,9 @@ def check_constraints(space, constraints):
                 for constraint_if_false in constraint.if_false:
                     check_constraints(space, [constraint_if_false])
         else:
-            raise TypeError('Constraints must be of type "Single", "Exlusive", "Inclusive", "Sum", "SumEquals" or "Conditional". Got {}'.format(type(constraint)))
+            raise TypeError(
+                'Constraints must be of type "Single", "Exclusive", "Inclusive", "Sum", "SumEquals" or "Conditional". Got {}'.format(type(constraint))
+            )
 
 
 def check_dim_and_space(space, constraint):
@@ -1170,7 +1171,7 @@ def check_value(dim, value):
             raise ValueError('Value {} exceeds bounds of space {}'.format(value, [dim.low, dim.high]))
     else:  # Categorical dimension.
         if value not in dim.categories:
-            raise ValueError('Categorical value {} is not in space with categoreis {}'.format(value, dim.categories))
+            raise ValueError('Categorical value {} is not in space with categories {}'.format(value, dim.categories))
 
 
 def check_is_constraint(constraint):
@@ -1180,4 +1181,4 @@ def check_is_constraint(constraint):
     if not (isinstance(constraint, Single) or isinstance(constraint, Inclusive)
             or isinstance(constraint, Exclusive) or isinstance(constraint, Sum)
             or isinstance(constraint, Conditional)):
-        raise TypeError('Constraint must be of type Inclusive, Exlusive, Single, Sum or Conditional. Got {}'.format(type(constraint)))
+        raise TypeError('Constraint must be of type Inclusive, Exclusive, Single, Sum or Conditional. Got {}'.format(type(constraint)))

From 181e244a7e1904cc6c06a9b2659476a8a33eaadd Mon Sep 17 00:00:00 2001
From: Morten Bormann Nielsen <mon@teknologisk.dk>
Date: Thu, 18 Dec 2025 15:24:07 +0100
Subject: [PATCH 12/19] Updates changelog for DRSC algorithm addition

---
 CHANGELOG.md | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 3997cde5..1b41f4d2 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,15 @@
 # Release history
 
+## Version 1.1.3 [unpublished]
+
+### Changes
+
+- Implemented the DRSC algorithm for SumEquals constraints. A side-effect of this implementation is that you can now
+  ask for multiple points with opt.ask() while using this type of constraint (previously you could only ask for one).
+
+### Bugfixes
+
+
 ## Version 1.1.2 [unpublished]
 
 ### Changes

From 0785a72e715af8b394b883e05641babe233e1e11 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=B8ren=20Furbo?= <srfu@novonordisk.com>
Date: Thu, 5 Feb 2026 13:04:27 +0100
Subject: [PATCH 13/19] fix: Swtiched evaluation order

---
 ProcessOptimizer/optimizer/optimizer.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/ProcessOptimizer/optimizer/optimizer.py b/ProcessOptimizer/optimizer/optimizer.py
index 9dd8e5ce..39c7f371 100644
--- a/ProcessOptimizer/optimizer/optimizer.py
+++ b/ProcessOptimizer/optimizer/optimizer.py
@@ -269,7 +269,7 @@ def __init__(
             )
 
         # check if regressor
-        if not is_regressor(base_estimator) and base_estimator is not None:
+        if base_estimator is not None and not is_regressor(base_estimator):
             raise ValueError("%s has to be a regressor." % base_estimator)
 
         # treat per second acqusition function specially
@@ -1309,12 +1309,12 @@ def __MinimalDistance(X, Y):
     # This function calls NSGAII to estimate the Pareto Front
     def NSGAII(self, MU=40):
         from ._NSGA2 import NSGAII
-        
+
         if MU % 4 != 0:
             raise ValueError(
                 "Number of simulated points must be divisible by 4 for the NSGAII algorithm"
             )
-            
+
         pop, logbook, front = NSGAII(
             self.n_objectives,
             self.__ObjectiveGP,

From 7ec5caac6f17e61a955d8f29ac0215fc05656eaf Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=B8ren=20Furbo?= <srfu@novonordisk.com>
Date: Thu, 5 Feb 2026 13:19:45 +0100
Subject: [PATCH 14/19] fix: Catching error and raising the correct type

---
 ProcessOptimizer/optimizer/optimizer.py | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/ProcessOptimizer/optimizer/optimizer.py b/ProcessOptimizer/optimizer/optimizer.py
index 39c7f371..229e1bf9 100644
--- a/ProcessOptimizer/optimizer/optimizer.py
+++ b/ProcessOptimizer/optimizer/optimizer.py
@@ -269,9 +269,11 @@ def __init__(
             )
 
         # check if regressor
-        if base_estimator is not None and not is_regressor(base_estimator):
-            raise ValueError("%s has to be a regressor." % base_estimator)
-
+        try:
+            if base_estimator is not None and not is_regressor(base_estimator):
+                raise ValueError("%s has to be a regressor." % base_estimator)
+        except AttributeError as me:
+            raise ValueError("%s has to be a regressor." % base_estimator) from me
         # treat per second acqusition function specially
         is_multi_regressor = isinstance(base_estimator, MultiOutputRegressor)
         if "ps" in self.acq_func and not is_multi_regressor:

From bd3e3ce891089fac6369652b9c4837bd3628f9f8 Mon Sep 17 00:00:00 2001
From: Morten Bormann Nielsen <mon@teknologisk.dk>
Date: Mon, 9 Feb 2026 14:21:47 +0100
Subject: [PATCH 15/19] Updates changelog and adds a demonstration script for
 the DRSC algorithm

---
 CHANGELOG.md                                  |   1 +
 ProcessOptimizer/plots.py                     |   3 +-
 .../DRSC SumEqual implementation.py           | 275 ++++++++++++++++++
 3 files changed, 278 insertions(+), 1 deletion(-)
 create mode 100644 examples/design_qualifications/DRSC SumEqual implementation.py

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 1b41f4d2..32adbe21 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -9,6 +9,7 @@
 
 ### Bugfixes
 
+- Fixed TypeError when combining categorical and numerical dimensions in get_Brownie_Bee_Pareto
 
 ## Version 1.1.2 [unpublished]
 
diff --git a/ProcessOptimizer/plots.py b/ProcessOptimizer/plots.py
index 20a8d053..417df47e 100644
--- a/ProcessOptimizer/plots.py
+++ b/ProcessOptimizer/plots.py
@@ -2719,7 +2719,8 @@ def get_Brownie_Bee_Pareto(optimizer, n_points=100):
     front_x = np.asarray(
         optimizer.space.inverse_transform(
             front_x.reshape(len(front_x), optimizer.space.transformed_n_dims)
-        )
+        ),
+        dtype=object,
     )
     
     # Sort the points in ascending order on objective 1
diff --git a/examples/design_qualifications/DRSC SumEqual implementation.py b/examples/design_qualifications/DRSC SumEqual implementation.py
new file mode 100644
index 00000000..3fa3242a
--- /dev/null
+++ b/examples/design_qualifications/DRSC SumEqual implementation.py	
@@ -0,0 +1,275 @@
+# -*- coding: utf-8 -*-
+"""
+A script demonstrating the features of the DRSC algorithm for use with the
+SumEquals constraint in ProcessOptimizer.
+
+Author: Morten Bormann Nielsen
+December 2025
+"""
+
+
+import time
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib as mpl
+
+from ProcessOptimizer import Optimizer
+from ProcessOptimizer.space import Real
+from ProcessOptimizer.space.constraints import Constraints, SumEquals
+
+from ProcessOptimizer.space.DRSC import DRSCGenerator
+
+from ProcessOptimizer.model_systems import hart3
+from ProcessOptimizer.plots import plot_objective_1d
+
+#%% Demonstrate speed of generator
+space = [
+    (52.9, 96.3),
+    (0.1, 0.5),
+    (3.6, 6.6),
+    (1.0, 5.0),
+    (1.0, 35.0),        
+]
+cons = [SumEquals(dimensions=[0, 1, 2, 3, 4], value=100.0, sampler="DRSC")]
+opt = Optimizer(space, lhs=False, n_initial_points=10)
+opt.set_constraints(cons)
+
+constraint = opt.get_constraints().sum_equals[0]
+
+drsc_gen = constraint.get_drsc_generator(opt.space)
+# Demonstrate speed of algorithm by generating 10,000 points. Takes about 5 seconds.
+start_time = time.time()
+x_simplex = drsc_gen.generate_sample(10000)
+print("Generated 10,000 points in the space in %s seconds" % (time.time() - start_time))
+
+#%% Use the new sampling method to generate a set of initial samples
+
+space = [
+    (52.9, 96.3),
+    (0.1, 0.5),
+    (3.6, 6.6),
+    (1.0, 5.0),
+    (1.0, 35.0),        
+]
+
+cons = [SumEquals(dimensions=[0, 1, 2, 3, 4], value=100.0, sampler="DRSC")]
+opt = Optimizer(space, lhs=False, n_initial_points=10, random_state=31031988)
+opt.set_constraints(cons)
+
+print(time.strftime("Starting opt.ask calculation with new method at:") + " " + time.strftime("%H:%M:%S"))
+x = opt.ask(10, strategy="cl_min")
+print(time.strftime("Finished opt.ask calculation with new method at:") + " " + time.strftime("%H:%M:%S"))
+
+#%% Generate a new optimizer with the same seed and constraint
+
+
+opt1 = Optimizer(space, lhs=False, n_initial_points=10, random_state=1)
+opt1.set_constraints(cons)
+opt2 = Optimizer(space, lhs=False, n_initial_points=10, random_state=2)
+opt2.set_constraints(cons)
+
+x1 = opt1.ask(10, strategy="cl_min")
+x2 = opt2.ask(10, strategy="cl_min")
+
+x!=x2
+
+#%% Test MANY more points
+
+constraint1 = opt1.get_constraints()
+constraint2 = opt2.get_constraints()
+
+x1 = constraint1.sumequal_sampling(n_samples=1000, random_state=31031988)
+x2 = constraint2.sumequal_sampling(n_samples=1000, random_state=31031988)
+
+# You can now test for x1 == x2
+
+#%% Test with categorical dimensions and integer dimensions present
+space = [
+    ("A", "B"),
+    (52.9, 96.3),
+    (0.1, 0.5),
+    (3.6, 6.6),
+    (1.0, 5.0),
+    (1.0, 35.0),
+    ("C", "D", "E"),
+    (1, 5),
+]
+cons = [SumEquals(dimensions=[1, 2, 3, 4, 5], value=100.0, sampler="DRSC")]
+opt = Optimizer(space, lhs=False, n_initial_points=10)
+opt.set_constraints(cons)
+
+print(time.strftime("Starting opt.ask calculation with new method at:") + " " + time.strftime("%H:%M:%S"))
+x = opt.ask(10, strategy="cl_min")
+print(time.strftime("Finished opt.ask calculation with new method at:") + " " + time.strftime("%H:%M:%S"))
+
+#%% Use the new sampling method with simulated data
+
+# Build Hartmann 3D ModelSystem object
+hart3_model = hart3.create_hart3(noise=True)
+# Define space
+space = [
+    Real(0., 1., name='x0'),
+    Real(0., 1., name='x1'),
+    Real(0., 1., name='x2'),
+]
+
+cons = [SumEquals(dimensions=[0, 1, 2], value=1.0, sampler="DRSC")]
+opt = Optimizer(space, lhs=False, n_initial_points=20)
+opt.set_constraints(cons)
+
+# Run optimization, first through the initial points, then using EI  
+for _ in range(30):
+    x = opt.ask(5, strategy="cl_min")
+    y = [hart3_model.get_score(x) for x in x]
+    res = opt.tell(x, y)
+
+# res = opt.get_result()
+plot_objective_1d(res, pars="expected_minimum")
+
+#%% Plot points
+fig = plt.figure(figsize=(8,6))
+ax = fig.add_subplot(projection='3d')
+cmap = mpl.colormaps["viridis"].resampled(len(opt.Xi))
+
+idx = np.arange(len(opt.Xi))
+x = np.array(opt.Xi)
+y = np.array(opt.yi)
+
+ax.scatter(x[:,0], x[:, 1], x[:, 2], c=idx, cmap=cmap, alpha=1)
+
+ax.set_xlim(0, 1)
+ax.set_ylim(0, 1)
+ax.set_zlim(0, 1)
+#% Manual check for speed of generating valid samples outside of debug mode
+
+# constraint = opt.get_constraints().sum_equals[0]
+
+# drsc_gen = constraint.get_drsc_generator(opt.space)
+# print(time.strftime("Generation 1000 points in the space. Starting calculation at:") + " " + time.strftime("%H:%M:%S"))
+# x_simplex = drsc_gen.generate_sample(1000)
+# print(time.strftime("Finished calculation at:") + " " + time.strftime("%H:%M:%S"))
+
+
+#%% Testing of a non-linear constraint 
+
+# Define space
+space = [
+    Real(0., 1., name='x0'),
+    Real(0., 1., name='x1'),
+    Real(0., 1., name='x2'),
+]
+
+sum_value = 1.0
+
+# Nonlinear constraint: x0 * x1 < 0.1
+# Convention: function returns value that should be <= 0 when SATISFIED
+def product_constraint(x):
+    """x[0] * x[1] <= 0.1"""
+    return x[0] * x[1] - 0.1
+
+# Create SumEquals constraint with nonlinear constraint
+cons = [
+    SumEquals(
+        dimensions=[0, 1, 2],
+        value=sum_value,
+        nonlinear_constraints=[product_constraint]
+    )
+]
+
+# Create optimizer
+opt = Optimizer(space, lhs=False, n_initial_points=10)
+opt.set_constraints(cons)
+
+# Test sampling
+constraint = opt.get_constraints().sum_equals[0]
+drsc_gen = constraint.get_drsc_generator(opt.space)
+
+# Generate samples
+samples = drsc_gen.generate_sample(1000)
+
+# Convert to original space and verify constraints
+samples_original = samples * sum_value
+
+x = samples_original
+fig = plt.figure(figsize=(8,6))
+ax = fig.add_subplot(projection='3d')
+ax.scatter(x[:, 0], x[:, 1], x[:, 2],)
+ax.set_xlim(0, 1)
+ax.set_ylim(0, 1)
+ax.set_zlim(0, 1)
+
+
+#%% Asking for an initial set of points with the old method
+
+# cons = [SumEquals(dimensions=[0, 1, 2, 3, 4], value=100.0, sampler="nullspace")]
+# opt = Optimizer(space, lhs=False, n_initial_points=10)
+# opt.set_constraints(cons)
+
+# # This takes about 25 seconds on MON's laptop
+# print(time.strftime("Starting opt.ask calculation with old method at:") + " " + time.strftime("%H:%M:%S"))
+# x = opt.ask(5, strategy="cl_min")
+# print(time.strftime("Finished opt.ask calculation with old method at:") + " " + time.strftime("%H:%M:%S"))
+
+
+#%% Create an example with three constrained dimensions
+
+# Sample 1000 points and plot them to show it works
+space = [
+    (0., 10.),
+    (0., 10.),
+    (0., 10.),
+    ("A", "B", "C"),
+    ]
+
+cons = [SumEquals(dimensions=[0, 1, 2], value=15.0, sampler="DRSC")]
+opt = Optimizer(space, lhs=False, n_initial_points=100)
+opt.set_constraints(cons)
+
+constraint = opt.get_constraints().sum_equals[0]
+
+x = opt.ask(100, strategy="cl_min")
+colors = []
+for xx in x:
+    if xx[3] == "A":
+        colors.append("r")
+    elif xx[3] == "B":
+        colors.append("g")
+    else:
+        colors.append("b")
+    
+#%%
+
+fig = plt.figure(figsize=(8,6))
+ax = fig.add_subplot(projection='3d')
+[ax.scatter(x[0], x[1], x[2], c=colors[i]) for i, x in enumerate(x)]
+ax.set_xlim(0, 10)
+ax.set_ylim(0, 10)
+ax.set_zlim(0, 10)
+
+
+#%%
+
+space = [
+    (0., 7.0),
+    (0., 7.0),
+    (0., 7.0),
+    (0., 7.0),
+    (0., 7.0),
+    (0., 7.0),
+    (0., 7.0),
+    (5.0, 20.0),
+    (4.0, 10.0),
+]
+
+cons = [SumEquals(dimensions=[0, 1, 2, 3, 4, 5, 6], value=7.0)]
+
+opt = Optimizer(space, lhs=False, n_initial_points=10,)
+opt.set_constraints(cons)
+start_time = time.time()
+x = opt.ask(10, strategy="cl_min")
+print("Generated 10 points in the space in %s seconds" % (time.time() - start_time))
+
+x = np.array(x)
+fig = plt.figure(figsize=(8,6))
+for i in range(7):
+    plt.scatter(np.arange(10), x[:, i])
\ No newline at end of file

From 3697a424b7330973c09c40fe7dae3569aa617edc Mon Sep 17 00:00:00 2001
From: Morten Bormann Nielsen <mon@teknologisk.dk>
Date: Mon, 9 Feb 2026 14:40:05 +0100
Subject: [PATCH 16/19] Removes the old null_space sampler option from
 SumEquals sampling

---
 ProcessOptimizer/space/constraints.py      | 179 +--------------------
 ProcessOptimizer/tests/test_constraints.py |   2 +-
 2 files changed, 5 insertions(+), 176 deletions(-)

diff --git a/ProcessOptimizer/space/constraints.py b/ProcessOptimizer/space/constraints.py
index 2ef1d4f1..1c2e4846 100644
--- a/ProcessOptimizer/space/constraints.py
+++ b/ProcessOptimizer/space/constraints.py
@@ -3,7 +3,6 @@
 import numpy as np
 
 from sklearn.utils import check_random_state
-from scipy import linalg
 
 from .space import Real, Integer, Categorical, Space, Task
 from .DRSC import DRSCGenerator
@@ -139,12 +138,11 @@ def sumequal_sampling(
             n_samples: int = 1,
             random_state: Union[int, np.random.RandomState, None] = None,
         ) -> List:
-        """Draw samples that respect SumEquals constraints.
+        """Draw samples that respect SumEquals constraints using the DRSC
+        algorithm.
 
         The samples are in the original space. They need to be transformed
         before being passed to a model or minimizer by `space.transform()`.
-        
-        Uses DRSC algorithm if available, falls back to null-space method.
 
         Parameters
         ----------
@@ -160,15 +158,8 @@ def sumequal_sampling(
         * `points`: [list of lists, shape=(n_points, n_dims)]
            Points sampled from the space.
         """
-        
-        # Check if we should use DRSC
-        if (len(self.sum_equals) == 1 and 
-            self.sum_equals[0].sampler == 'DRSC'):
-            
-            return self._drsc_sampling(n_samples, random_state)
-        else:
-            # Use existing null-space method
-            return self._nullspace_sampling(n_samples, random_state)
+        return self._drsc_sampling(n_samples, random_state)
+
         
     def _drsc_sampling(self, n_samples, random_state):
         """Sample using DRSC algorithm."""
@@ -263,161 +254,6 @@ def _fix_bounds_preserving_sum(self, x, bounds, target_sum, max_iterations=100):
         
         return x
 
-    
-    def _nullspace_sampling(self, n_samples, random_state):
-        """
-        Original null-space based sampling method. Consider deprecating
-        immediately before the full pull request.
-        """
-        def null_space(A, rcond=None) -> np.ndarray:
-            """Helper function to calculate the null space of a matrix
-
-            Parameters
-            ----------
-            A : numpy array
-                The matrix to calculate the null space of
-            rcond : float, optional
-                The tolerance for determining the effective rank of A.
-
-            Returns
-            -------
-            Q : numpy array
-                The null space of A, as a matrix with orthonormal columns.
-            """
-            
-            u, s, vh = np.linalg.svd(A, full_matrices=True) # A = u*s*vh
-            M, N = u.shape[0], vh.shape[1] # M = number of rows, N = number of columns
-            if rcond is None: # default value of rcond
-                rcond = np.finfo(s.dtype).eps * max(M, N) # machine precision times max dimension
-            tol = np.amax(s) * rcond # tolerance for singular values
-            num = np.sum(s > tol, dtype=int) # number of singular values greater than tol
-            Q = vh[num:,:].T.conj() # columns of vh corresponding to singular values greater than tol
-            return Q # return the null space of A
-
-        rng = check_random_state(random_state)
-                
-        # Find a point on the plane defined by the SumEquals constraint where
-        # A + B + ... = value. We do this by asking where the diagonal between
-        # the origin and A_max, B_max, ... intersects the plane.
-        d = len(self.sum_equals[0].dimensions)
-        origin = np.array(
-            [self.space.bounds[dim][0] 
-             for dim in self.sum_equals[0].dimensions]
-        )
-        delta = np.array(
-            [self.space.bounds[dim][1] - self.space.bounds[dim][0] 
-             for dim in self.sum_equals[0].dimensions]
-        )
-        
-        A = np.zeros((d,d))
-        B = np.zeros(d)
-        # Row representing the sum constraint
-        A[0,:] = 1
-        B[0] = self.sum_equals[0].value
-        # Rows that define the linear equation for the diagonal along the 
-        # constrained dimensions
-        for i in range(1,d):
-            A[i, 0] = -delta[i]/delta[0]
-            A[i, i] = 1
-            B[i] = origin[i]
-        # Identify the point that lies on the constraint plane and on the diagonal
-        point = np.linalg.solve(A, B)
-        # Use the fact that the vector [1, 1, ...] (a 1 for each constrained 
-        # dimension) is normal to the plane defined by A + B + ... to build
-        # basis-vectors inside the plane, using the null_space function
-        N = np.array(np.ones(d))
-        ns = null_space(N[np.newaxis, :])
-        # We only need to simulate points up to a distance of half the diagonal
-        # from the origin to A_max, B_max, etc.
-        sim_distance = np.sqrt(np.sum(delta**2)) / 2
-        
-        # To avoid "clustering" of points in the constrained plane, we will 
-        # create samples using low discrepancy quasirandom sequences, see:
-        # http://extremelearning.com.au/unreasonable-effectiveness-of-quasirandom-sequences/
-        # for background on this method
-        
-        # Helper function for calculating the generalized golden ratio
-        def phi(d):
-            x = 2.0
-            for i in range(10): 
-                x = pow(1+x,1/(d+1)) 
-            return x
-        # Golden ratio for our present dimensionality (the constrained space)
-        g = phi(d-1)
-        alpha = np.zeros(d-1)
-        for j in range(d-1):
-            alpha[j] = pow(1/g, j+1) %1
-        vec_comp = np.zeros((1, d-1))
-        # Choose seed (starting location) in the normalized space
-        seed = 0.5
-        
-        # Build a list of samples
-        samples = []
-        j = 0
-        while len(samples) < n_samples:
-            # Generate next step in the sequence
-            vec_comp = (seed + alpha*(j+1)) %1
-            j += 1
-            # Center these components on zero
-            vec_comp = vec_comp - 0.5
-            # Simulate lengths of each null_space vector to add to our point
-            vec_comp = vec_comp * sim_distance * 2
-            sample_candidate = point[:, None].T + ns @ vec_comp.T
-            # Generate the correct shape 
-            sample_candidate = sample_candidate[0]
-            # Check that the candidate is inside the original parameter space
-            inspace = [
-                (sample_candidate[i] >= self.space.bounds[dim][0]) and
-                (sample_candidate[i] <= self.space.bounds[dim][1])
-                for i, dim in enumerate(self.sum_equals[0].dimensions)
-            ]
-            # Only accept the candidate if it is in our space
-            if all(inspace):
-                samples.append(sample_candidate)      
-        # Convert the list of arrays to a list of lists
-        samples = [arr.tolist() for arr in samples]
-        
-        # Create settings for the dimensions that are not part of the constraint
-        if d < self.space.n_dims:
-            remaining_dimensions = [
-                i for i in range(self.space.n_dims) 
-                if i not in self.sum_equals[0].dimensions
-            ]
-            # Convert our list of samples to an array
-            samples = np.array(samples)
-            # Expand the sample array to make space for settings of the 
-            # unconstrained dimensions
-            for i in remaining_dimensions:
-                samples = np.insert(
-                    samples,
-                    i, 
-                    np.zeros(len(samples)), 
-                    axis=1
-                )
-            # Convert back to list of lists
-            samples = samples.tolist()
-            
-            # Generate random settings across all factors
-            full_sample = self.space.rvs(n_samples=n_samples, random_state=rng)
-            # Sort the settings in each column, which will ensure that the
-            # unconstrained settings are distributed in a space-filling way too
-            transposed_sample = list(zip(*full_sample))
-            sort_trans_sample = [sorted(inner_tuple) for inner_tuple in transposed_sample]
-            # Place the settings back in a list of lists
-            full_sample = list(map(list, zip(*sort_trans_sample)))
-            
-            # Overwrite the random setting values for the constrained factors
-            for j in remaining_dimensions:
-                for i in range(len(samples)):
-                    samples[i][j] = full_sample[i][j]
-            
-            # Shuffle the order of the samples, otherwise the unconstrained
-            # settings will be returned in a sorted order. Use seeding to 
-            # provide consistent initial samples
-            rng2 = np.random.default_rng(seed=42)
-            rng2.shuffle(samples)
-                    
-        return samples
 
     def validate_sample(self, sample: List) -> bool:
         """ Validates a sample of parameter values in regards to the
@@ -781,7 +617,6 @@ def __init__(
         self, 
         dimensions: List[int], 
         value: Union[float, int],
-        sampler: str = "DRSC",
         linear_constraints: Optional[List[Tuple[np.ndarray, float]]] = None,
         nonlinear_constraints: Optional[List[Callable]] = None,
     ):
@@ -799,11 +634,6 @@ def __init__(
 
         * `value` [float or int]:
             The value for which the sum should be equal to.
-        
-        * `sampler` [str, default='DRSC']:
-            Sampling algorithm to use. Options:
-            - 'DRSC': Dirichlet-Rescale-Constraints (fast, handles additional constraints)
-            - 'nullspace': Original null-space based sampler
             
         * `linear_constraints` [list of (a, b) tuples, optional]:
             Additional linear constraints a^T x <= b (only used with DRSC)
@@ -826,7 +656,6 @@ def __init__(
 
         self.dimensions = tuple(dimensions)
         self.value = value
-        self.sampler = sampler
         self.linear_constraints = linear_constraints or []
         self.nonlinear_constraints = nonlinear_constraints or []
         # Lazy initialization of DRSC generator
diff --git a/ProcessOptimizer/tests/test_constraints.py b/ProcessOptimizer/tests/test_constraints.py
index d459bb1f..10e62dde 100644
--- a/ProcessOptimizer/tests/test_constraints.py
+++ b/ProcessOptimizer/tests/test_constraints.py
@@ -301,7 +301,7 @@ def test_SumEquals():
         ("A", "B", "C"),
         (1, 10),
     ]
-    cons = [SumEquals(dimensions=[0, 1, 2, 3, 4], value=100.0, sampler="DRSC")]
+    cons = [SumEquals(dimensions=[0, 1, 2, 3, 4], value=100.0)]
     # Build optimziers with different seeds
     opt1 = Optimizer(space, lhs=False, n_initial_points=10, random_state=1)
     opt1.set_constraints(cons)

From e6a47b9642c8541d119c4d77491b0d9ba2122a3f Mon Sep 17 00:00:00 2001
From: Morten Bormann Nielsen <mon@teknologisk.dk>
Date: Mon, 9 Feb 2026 14:41:46 +0100
Subject: [PATCH 17/19] Removes need to specify sampler when getting DRSC
 generator

---
 ProcessOptimizer/space/constraints.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/ProcessOptimizer/space/constraints.py b/ProcessOptimizer/space/constraints.py
index 1c2e4846..c6e33bd6 100644
--- a/ProcessOptimizer/space/constraints.py
+++ b/ProcessOptimizer/space/constraints.py
@@ -675,7 +675,7 @@ def get_drsc_generator(self, space) -> DRSCGenerator:
         -------
             DRSCGenerator instance configured for this constraint
         """
-        if self._drsc_generator is None and self.sampler == 'DRSC':
+        if self._drsc_generator is None:
             # Extract bounds for constrained dimensions
             bounds = [space.bounds[dim] for dim in self.dimensions]
             

From e77b5c110a009831c00c5f25b5732b1944be5648 Mon Sep 17 00:00:00 2001
From: Morten Bormann Nielsen <mon@teknologisk.dk>
Date: Mon, 9 Feb 2026 14:45:29 +0100
Subject: [PATCH 18/19] Adds an example script demonstrating the use of
 SumEquals sampling with the DRSC algorithm

---
 .../DRSC SumEqual implementation.py           | 113 ++++++------------
 1 file changed, 34 insertions(+), 79 deletions(-)

diff --git a/examples/design_qualifications/DRSC SumEqual implementation.py b/examples/design_qualifications/DRSC SumEqual implementation.py
index 3fa3242a..a22fcba3 100644
--- a/examples/design_qualifications/DRSC SumEqual implementation.py	
+++ b/examples/design_qualifications/DRSC SumEqual implementation.py	
@@ -3,7 +3,7 @@
 A script demonstrating the features of the DRSC algorithm for use with the
 SumEquals constraint in ProcessOptimizer.
 
-Author: Morten Bormann Nielsen
+Author: Morten Bormann Nielsen, Danish Technological Institute
 December 2025
 """
 
@@ -15,14 +15,15 @@
 
 from ProcessOptimizer import Optimizer
 from ProcessOptimizer.space import Real
-from ProcessOptimizer.space.constraints import Constraints, SumEquals
-
-from ProcessOptimizer.space.DRSC import DRSCGenerator
+from ProcessOptimizer.space.constraints import SumEquals
 
 from ProcessOptimizer.model_systems import hart3
 from ProcessOptimizer.plots import plot_objective_1d
 
-#%% Demonstrate speed of generator
+#%% Setup of SumEquals sampling for very oblong space
+
+# The DRSC generator for SumEquals sampling is able to generate valid points
+# extremely fast
 space = [
     (52.9, 96.3),
     (0.1, 0.5),
@@ -30,20 +31,23 @@
     (1.0, 5.0),
     (1.0, 35.0),        
 ]
-cons = [SumEquals(dimensions=[0, 1, 2, 3, 4], value=100.0, sampler="DRSC")]
+# Define the constraint as all factor settings adding up to 100
+cons = [SumEquals(dimensions=[0, 1, 2, 3, 4], value=100.0)]
 opt = Optimizer(space, lhs=False, n_initial_points=10)
 opt.set_constraints(cons)
 
+# Extract the generator itself
 constraint = opt.get_constraints().sum_equals[0]
-
 drsc_gen = constraint.get_drsc_generator(opt.space)
-# Demonstrate speed of algorithm by generating 10,000 points. Takes about 5 seconds.
+# Demonstrate speed of algorithm by generating 10,000 points. This call takes
+# about 12 seconds on an Intel i7-12850HX.
 start_time = time.time()
 x_simplex = drsc_gen.generate_sample(10000)
 print("Generated 10,000 points in the space in %s seconds" % (time.time() - start_time))
 
-#%% Use the new sampling method to generate a set of initial samples
+#%% Demonstrate typical use for a real user with fewer points
 
+# Define the space
 space = [
     (52.9, 96.3),
     (0.1, 0.5),
@@ -51,17 +55,16 @@
     (1.0, 5.0),
     (1.0, 35.0),        
 ]
-
-cons = [SumEquals(dimensions=[0, 1, 2, 3, 4], value=100.0, sampler="DRSC")]
+# Set up the constraint
+cons = [SumEquals(dimensions=[0, 1, 2, 3, 4], value=100.0)]
 opt = Optimizer(space, lhs=False, n_initial_points=10, random_state=31031988)
 opt.set_constraints(cons)
-
-print(time.strftime("Starting opt.ask calculation with new method at:") + " " + time.strftime("%H:%M:%S"))
+# Ask for settings for the initial experiments
+print(time.strftime("Starting opt.ask calculation with DRSC method at:") + " " + time.strftime("%H:%M:%S"))
 x = opt.ask(10, strategy="cl_min")
-print(time.strftime("Finished opt.ask calculation with new method at:") + " " + time.strftime("%H:%M:%S"))
-
-#%% Generate a new optimizer with the same seed and constraint
+print(time.strftime("Finished opt.ask calculation with DRSC method at:") + " " + time.strftime("%H:%M:%S"))
 
+#%% Demonstrate that different seeds lead to different points
 
 opt1 = Optimizer(space, lhs=False, n_initial_points=10, random_state=1)
 opt1.set_constraints(cons)
@@ -73,7 +76,7 @@
 
 x!=x2
 
-#%% Test MANY more points
+#%% Demonstrate that identical seeds lead to identical points
 
 constraint1 = opt1.get_constraints()
 constraint2 = opt2.get_constraints()
@@ -81,9 +84,9 @@
 x1 = constraint1.sumequal_sampling(n_samples=1000, random_state=31031988)
 x2 = constraint2.sumequal_sampling(n_samples=1000, random_state=31031988)
 
-# You can now test for x1 == x2
+x1 == x2
 
-#%% Test with categorical dimensions and integer dimensions present
+#%% Demonstrate how to use SumEqual constraints when categoricals are present
 space = [
     ("A", "B"),
     (52.9, 96.3),
@@ -94,13 +97,13 @@
     ("C", "D", "E"),
     (1, 5),
 ]
-cons = [SumEquals(dimensions=[1, 2, 3, 4, 5], value=100.0, sampler="DRSC")]
+cons = [SumEquals(dimensions=[1, 2, 3, 4, 5], value=100.0)]
 opt = Optimizer(space, lhs=False, n_initial_points=10)
 opt.set_constraints(cons)
 
-print(time.strftime("Starting opt.ask calculation with new method at:") + " " + time.strftime("%H:%M:%S"))
+print(time.strftime("Starting opt.ask calculation with DRSC method at:") + " " + time.strftime("%H:%M:%S"))
 x = opt.ask(10, strategy="cl_min")
-print(time.strftime("Finished opt.ask calculation with new method at:") + " " + time.strftime("%H:%M:%S"))
+print(time.strftime("Finished opt.ask calculation with DRSC method at:") + " " + time.strftime("%H:%M:%S"))
 
 #%% Use the new sampling method with simulated data
 
@@ -113,7 +116,7 @@
     Real(0., 1., name='x2'),
 ]
 
-cons = [SumEquals(dimensions=[0, 1, 2], value=1.0, sampler="DRSC")]
+cons = [SumEquals(dimensions=[0, 1, 2], value=1.0)]
 opt = Optimizer(space, lhs=False, n_initial_points=20)
 opt.set_constraints(cons)
 
@@ -123,10 +126,10 @@
     y = [hart3_model.get_score(x) for x in x]
     res = opt.tell(x, y)
 
-# res = opt.get_result()
+# Show the system
 plot_objective_1d(res, pars="expected_minimum")
 
-#%% Plot points
+# Show the location of the sampled points in this experiment
 fig = plt.figure(figsize=(8,6))
 ax = fig.add_subplot(projection='3d')
 cmap = mpl.colormaps["viridis"].resampled(len(opt.Xi))
@@ -140,17 +143,9 @@
 ax.set_xlim(0, 1)
 ax.set_ylim(0, 1)
 ax.set_zlim(0, 1)
-#% Manual check for speed of generating valid samples outside of debug mode
-
-# constraint = opt.get_constraints().sum_equals[0]
 
-# drsc_gen = constraint.get_drsc_generator(opt.space)
-# print(time.strftime("Generation 1000 points in the space. Starting calculation at:") + " " + time.strftime("%H:%M:%S"))
-# x_simplex = drsc_gen.generate_sample(1000)
-# print(time.strftime("Finished calculation at:") + " " + time.strftime("%H:%M:%S"))
 
-
-#%% Testing of a non-linear constraint 
+#%% Demonstrate the use of a non-linear SumEquals constraint 
 
 # Define space
 space = [
@@ -191,6 +186,7 @@ def product_constraint(x):
 samples_original = samples * sum_value
 
 x = samples_original
+# Create a plot showing where the samples lie in the space
 fig = plt.figure(figsize=(8,6))
 ax = fig.add_subplot(projection='3d')
 ax.scatter(x[:, 0], x[:, 1], x[:, 2],)
@@ -199,19 +195,7 @@ def product_constraint(x):
 ax.set_zlim(0, 1)
 
 
-#%% Asking for an initial set of points with the old method
-
-# cons = [SumEquals(dimensions=[0, 1, 2, 3, 4], value=100.0, sampler="nullspace")]
-# opt = Optimizer(space, lhs=False, n_initial_points=10)
-# opt.set_constraints(cons)
-
-# # This takes about 25 seconds on MON's laptop
-# print(time.strftime("Starting opt.ask calculation with old method at:") + " " + time.strftime("%H:%M:%S"))
-# x = opt.ask(5, strategy="cl_min")
-# print(time.strftime("Finished opt.ask calculation with old method at:") + " " + time.strftime("%H:%M:%S"))
-
-
-#%% Create an example with three constrained dimensions
+#%% Demonstrate the use of three constrained dimensions and a categorical
 
 # Sample 1000 points and plot them to show it works
 space = [
@@ -219,9 +203,9 @@ def product_constraint(x):
     (0., 10.),
     (0., 10.),
     ("A", "B", "C"),
-    ]
+]
 
-cons = [SumEquals(dimensions=[0, 1, 2], value=15.0, sampler="DRSC")]
+cons = [SumEquals(dimensions=[0, 1, 2], value=15.0)]
 opt = Optimizer(space, lhs=False, n_initial_points=100)
 opt.set_constraints(cons)
 
@@ -237,39 +221,10 @@ def product_constraint(x):
     else:
         colors.append("b")
     
-#%%
-
+# Create a plot showing the distribution of points in this space
 fig = plt.figure(figsize=(8,6))
 ax = fig.add_subplot(projection='3d')
 [ax.scatter(x[0], x[1], x[2], c=colors[i]) for i, x in enumerate(x)]
 ax.set_xlim(0, 10)
 ax.set_ylim(0, 10)
 ax.set_zlim(0, 10)
-
-
-#%%
-
-space = [
-    (0., 7.0),
-    (0., 7.0),
-    (0., 7.0),
-    (0., 7.0),
-    (0., 7.0),
-    (0., 7.0),
-    (0., 7.0),
-    (5.0, 20.0),
-    (4.0, 10.0),
-]
-
-cons = [SumEquals(dimensions=[0, 1, 2, 3, 4, 5, 6], value=7.0)]
-
-opt = Optimizer(space, lhs=False, n_initial_points=10,)
-opt.set_constraints(cons)
-start_time = time.time()
-x = opt.ask(10, strategy="cl_min")
-print("Generated 10 points in the space in %s seconds" % (time.time() - start_time))
-
-x = np.array(x)
-fig = plt.figure(figsize=(8,6))
-for i in range(7):
-    plt.scatter(np.arange(10), x[:, i])
\ No newline at end of file

From dd32362445c093b3ed8952bef4bf4f1ddfee406d Mon Sep 17 00:00:00 2001
From: Morten Bormann Nielsen <mon@teknologisk.dk>
Date: Mon, 9 Feb 2026 14:58:55 +0100
Subject: [PATCH 19/19] Updates changelog and pyproject.toml to make clear we
 no longer support Python 3.9

---
 CHANGELOG.md   | 1 +
 pyproject.toml | 1 -
 2 files changed, 1 insertion(+), 1 deletion(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 32adbe21..6cdf2ffc 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -6,6 +6,7 @@
 
 - Implemented the DRSC algorithm for SumEquals constraints. A side-effect of this implementation is that you can now
   ask for multiple points with opt.ask() while using this type of constraint (previously you could only ask for one).
+- Dropped support for Python 3.9
 
 ### Bugfixes
 
diff --git a/pyproject.toml b/pyproject.toml
index 2e94aa42..0730b9d1 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -19,7 +19,6 @@ classifiers = [
   "License :: OSI Approved :: BSD License",
   "Operating System :: OS Independent",
   "Programming Language :: Python :: 3",
-  "Programming Language :: Python :: 3.9",
   "Programming Language :: Python :: 3.10",
   "Programming Language :: Python :: 3.11",
   "Programming Language :: Python :: 3.12",