Introduce the DRSC algorithm for constraint handling

[dk-teknologisk-mon]

This pull request implements the Dirichlet Rescale Constraints (DRSC) algorithm from this paper. This algorithm generates random n-dimensional vectors with a fixed sum in a way that is many, many orders of magnitude faster than the present way we handle SumEquals constrains.

A bit of motivation, using a formulation example:

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)]
opt = Optimizer(space, lhs=False, n_initial_points=10)
opt.set_constraints(cons)
x = opt.ask(5, strategy="cl_min")

Running this code takes about 25 seconds on my laptop despite the internals only calling this function asking for 5 points:

self._constraints.sumequal_sampling(n_samples=self.n_points, random_state=self.rng)

Things become unsustainable if you get as far as providing data, because at that stage opt._tell() will ask for a set of candidate points for acquisition where it will call sumequal_sampling for 10,000 points. This would take about 8 hours to complete.

With the new algorithm however:

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)
x_simplex = drsc_gen.generate_sample(10000)

This completes in about 13 seconds!

=======================
Some added benefits of the way I've implemented this:

1. You can now ask for multiple points with SumEquals constraints to your hearts desire, as long as you pass strategy="cl_min" top opt.ask().
2. The implementation makes it feasible to introduce Steinerberger sampling for SumEquals constraints as well (but this should be its own pull request).
3. The change is entirely backwards compatible and we don't have to change anything at all inside the optimizer class.
4. You can use nonlinear inequality constraints on top of the SumEquals constraint, like the constraint x0*x1 <= 0.1, shown below with 1000 samples:

Nonlinear constraints require the user to set up a function that returns a value <= 0 when it is satisfied, like this example:

# 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):
    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]
    )
]

closes #325 #256 #254

[dk-teknologisk-mon]

Some questions I'd like input on:

1. Should we just kill the old implementation of SumEquals sampling entirely? Right now I've kept it as an option you can pass to the SumEquals constraint in the form of sampler="nullspace" but truth be told it is so inefficient that I can't think of any reason to keep it around?
2. What on earth is going on in python 3.11 and 3.12?

[SRFU-NN]

Really nice