Repository: MachineLearningLifeScience/poli (language: python)

Issue:
Add Ehrlich functions
[We could add the latest work by Stanton et al. as a black box](https://arxiv.org/abs/2407.0023)


Discussion / hints:

Fix patch:
diff --git a/src/poli/objective_repository/__init__.py b/src/poli/objective_repository/__init__.py
index e6ba21c8..66debbe2 100644
--- a/src/poli/objective_repository/__init__.py
+++ b/src/poli/objective_repository/__init__.py
@@ -11,6 +11,9 @@
     ToyContinuousBlackBox,
 )
 
+# Discrete toy examples
+from .ehrlich.register import EhrlichBlackBox, EhrlichProblemFactory
+
 # Chemistry black boxes
 from .dockstring.register import DockstringProblemFactory, DockstringBlackBox
 from .drd3_docking.register import DRD3ProblemFactory, DRD3BlackBox
@@ -142,6 +145,7 @@
 
 AVAILABLE_PROBLEM_FACTORIES = {
     "aloha": AlohaProblemFactory,
+    "ehrlich": EhrlichProblemFactory,
     "dockstring": DockstringProblemFactory,
     "drd3_docking": DRD3ProblemFactory,
     "foldx_rfp_lambo": FoldXRFPLamboProblemFactory,
@@ -185,6 +189,7 @@
 
 AVAILABLE_BLACK_BOXES = {
     "aloha": AlohaBlackBox,
+    "ehrlich": EhrlichBlackBox,
     "dockstring": DockstringBlackBox,
     "drd3_docking": DRD3BlackBox,
     "foldx_rfp_lambo": FoldXRFPLamboBlackBox,
diff --git a/src/poli/objective_repository/ehrlich/__init__.py b/src/poli/objective_repository/ehrlich/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/src/poli/objective_repository/ehrlich/_construct_feasibility_matrix.py b/src/poli/objective_repository/ehrlich/_construct_feasibility_matrix.py
new file mode 100644
index 00000000..e6600948
--- /dev/null
+++ b/src/poli/objective_repository/ehrlich/_construct_feasibility_matrix.py
@@ -0,0 +1,63 @@
+from __future__ import annotations
+
+import numpy as np
+
+
+def _construct_banded_matrix(size: int) -> np.ndarray:
+    """
+    Constructs a matrix of zeroes and ones, where
+    the ones are bands that can loop around.
+    """
+    matrix = np.zeros((size, size), dtype=int)
+    band_index = 0
+    band_length = size - 1
+    for row_i in range(size):
+        indices_for_positions_that_will_be_1 = list(
+            range(band_index, band_index + band_length)
+        )
+
+        # Looping the ones that go over the limit
+        for i in range(len(indices_for_positions_that_will_be_1)):
+            if indices_for_positions_that_will_be_1[i] >= size:
+                indices_for_positions_that_will_be_1[i] -= size
+
+        matrix[row_i, indices_for_positions_that_will_be_1] = 1
+
+        band_index += 1
+
+    return matrix
+
+
+def _construct_binary_mask(size: int) -> np.ndarray:
+    banded_matrix = _construct_banded_matrix(size)
+
+    # Shuffle its rows
+    random_indices_for_rows = np.random.permutation(size)
+    binary_mask_matrix = banded_matrix[random_indices_for_rows]
+
+    # Making sure that the diagonal is full
+    # of ones
+    binary_mask_matrix[np.diag_indices(size)] = 1
+
+    return binary_mask_matrix
+
+
+def _construct_transition_matrix(size: int, seed: int | None = None) -> np.ndarray:
+    binary_mask_matrix = _construct_binary_mask(size)
+
+    # Creating a random state and matrix
+    random_state = np.random.RandomState(seed)
+    random_matrix = random_state.randn(size, size)
+
+    # Softmax it
+    transition_matrix = np.exp(random_matrix) / np.sum(np.exp(random_matrix), axis=0)
+
+    # Mask it
+    masked_transition_matrix = transition_matrix * binary_mask_matrix
+
+    # Normalize it
+    normalized_transition_matrix = masked_transition_matrix / np.sum(
+        masked_transition_matrix, axis=1, keepdims=True
+    )
+
+    return normalized_transition_matrix
diff --git a/src/poli/objective_repository/ehrlich/information.py b/src/poli/objective_repository/ehrlich/information.py
new file mode 100644
index 00000000..45e4032a
--- /dev/null
+++ b/src/poli/objective_repository/ehrlich/information.py
@@ -0,0 +1,15 @@
+import numpy as np
+
+from poli.core.black_box_information import BlackBoxInformation
+
+ehrlich_info = BlackBoxInformation(
+    name="ehrlich",
+    max_sequence_length=np.inf,
+    aligned=True,
+    fixed_length=True,
+    deterministic=True,  # ?
+    alphabet=None,  # TODO: add alphabet once we settle for one for SMLIES/SELFIES.
+    log_transform_recommended=False,
+    discrete=True,
+    padding_token="",
+)
diff --git a/src/poli/objective_repository/ehrlich/register.py b/src/poli/objective_repository/ehrlich/register.py
new file mode 100644
index 00000000..70e89bb2
--- /dev/null
+++ b/src/poli/objective_repository/ehrlich/register.py
@@ -0,0 +1,451 @@
+"""
+This module implements Ehrlich functions as black boxes in poli.
+
+Ehrlich functions were proposed by Stanton et al. [1] as a quick-and-easy
+alternative for testing discrete sequence optimizers (with protein
+optimization in mind). They are deviced to 
+
+(i) be easy to query,
+(ii) have feasible and unfeasible sequences,
+(iii) have uninformative random samples (i.e. randomly sampling
+     and evaluating should not be competitive, as many of these
+     should be unfeasible).
+(iv) be maximized when certain motifs are present in the sequence.
+     These motifs can be long-range within the sequence, and are
+     meant to be non-additive.
+
+Check the references for details on the implementation.
+
+References
+----------
+[1] Stanton, S., Alberstein, R., Frey, N., Watkins, A., & Cho, K. (2024).
+    Closed-Form Test Functions for Biophysical Sequence Optimization Algorithms.
+    arXiv preprint arXiv:2407.00236. https://arxiv.org/abs/2407.00236
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+from poli.core.abstract_black_box import AbstractBlackBox
+from poli.core.abstract_problem_factory import AbstractProblemFactory
+
+from poli.core.black_box_information import BlackBoxInformation
+from poli.core.problem import Problem
+from poli.objective_repository.ehrlich._construct_feasibility_matrix import (
+    _construct_transition_matrix,
+)
+from poli.objective_repository.ehrlich.information import (
+    ehrlich_info,
+)
+
+from poli.core.util.seeding import seed_python_numpy_and_torch
+
+from poli.core.util.proteins.defaults import AMINO_ACIDS
+
+
+class EhrlichBlackBox(AbstractBlackBox):
+    """
+    Ehrlich functions were proposed by Stanton et al. [1] as a quick-and-easy
+    alternative for testing discrete sequence optimizers (with protein
+    optimization in mind). They are deviced to
+
+    (i) be easy to query,
+    (ii) have feasible and unfeasible sequences,
+    (iii) have uninformative random samples (i.e. randomly sampling
+        and evaluating should not be competitive, as many of these
+        should be unfeasible).
+    (iv) be maximized when certain motifs are present in the sequence.
+        These motifs can be long-range within the sequence, and are
+        meant to be non-additive.
+
+    Check the references for details on the implementation.
+
+    Parameters
+    ----------
+    sequence_length : int
+        The length of the sequence to be optimized. This length is fixed, and
+        _only_ sequences of this length are considered.
+    motif_length : int
+        The length of the motifs.
+    n_motifs : int
+        The number of motifs.
+    quantization : int, optional
+        The quantization parameter. This parameter must be between 1 and the
+        motif length, and the motif length must be divisible by the quantization.
+        By default, it is None (which corresponds to the motif length).
+    seed : int, optional
+        The seed for the random number generator. By default, it is None
+        (i.e. no seed is set).
+    alphabet : list of str, optional
+        The alphabet to be used for the sequences. By default, it is the
+        of 20 amino acids.
+    batch_size : int, optional
+        The batch size for the black box. By default, it is None (i.e. all
+        sequences are evaluated in a vectorized way).
+    parallelize : bool, optional
+        Whether to parallelize the evaluation of the black box. By default,
+        it is False.
+    num_workers : int, optional
+        The number of processors used in parallelization.
+    evaluation_budget : int, optional
+        The evaluation budget for the black box. By default, it is infinite.
+
+    References
+    ----------
+    [1] Stanton, S., Alberstein, R., Frey, N., Watkins, A., & Cho, K. (2024).
+        Closed-Form Test Functions for Biophysical Sequence Optimization Algorithms.
+        arXiv preprint arXiv:2407.00236. https://arxiv.org/abs/2407.00236
+
+    """
+
+    def __init__(
+        self,
+        sequence_length: int,
+        motif_length: int,
+        n_motifs: int,
+        quantization: int | None = None,
+        seed: int = None,
+        alphabet: list[str] = AMINO_ACIDS,
+        batch_size: int = None,
+        parallelize: bool = False,
+        num_workers: int = None,
+        evaluation_budget: int = float("inf"),
+    ):
+        super().__init__(batch_size, parallelize, num_workers, evaluation_budget)
+        self.alphabet = alphabet
+        self.sequence_length = sequence_length
+
+        if seed is not None:
+            seed_python_numpy_and_torch(seed)
+
+        if motif_length * n_motifs > sequence_length:
+            raise ValueError(
+                "The total length of the motifs is greater than the sequence length."
+            )
+
+        if quantization is None:
+            quantization = motif_length
+
+        if not (1 <= quantization <= motif_length) or motif_length % quantization != 0:
+            raise ValueError(
+                "The quantization parameter must be between 1 and the motif length, "
+                "and the motif length must be divisible by the quantization."
+            )
+
+        self.motif_length = motif_length
+        self.n_motifs = n_motifs
+        self.quantization = quantization
+
+        self.transition_matrix = _construct_transition_matrix(
+            size=len(alphabet),
+            seed=seed,
+        )
+
+        self.motifs = self.construct_random_motifs(
+            motif_length=motif_length,
+            n_motifs=n_motifs,
+            seed=seed,
+        )
+        self.offsets = self.construct_random_offsets(
+            motif_length=motif_length,
+            n_motifs=n_motifs,
+            seed=seed,
+        )
+
+    def _sample_random_sequence(
+        self,
+        length: int | None = None,
+        random_state: int | np.random.RandomState | None = None,
+        repeating_allowed: bool = True,
+    ) -> str:
+        """
+        Uses the sparse transition matrix to generate a random sequence
+        of a given length.
+        """
+        if length is None:
+            length = self.sequence_length
+
+        if random_state is None:
+            random_state = np.random.RandomState()
+
+        if isinstance(random_state, int):
+            random_state = np.random.RandomState(random_state)
+        elif isinstance(random_state, np.random.RandomState):
+            pass
+        else:
+            raise ValueError(
+                "The random_state parameter must be an integer or an instance of "
+                "np.random.RandomState."
+            )
+
+        sequence = self.alphabet[random_state.randint(len(self.alphabet))]
+        current_state = self.alphabet.index(sequence)
+
+        for _ in range(length - 1):
+            next_state = random_state.choice(
+                len(self.alphabet), p=self.transition_matrix[current_state]
+            )
+            if not repeating_allowed:
+                while next_state == current_state:
+                    next_state = random_state.choice(
+                        len(self.alphabet),
+                        p=self.transition_matrix[current_state],
+                    )
+            sequence += self.alphabet[next_state]
+            current_state = next_state
+
+        return sequence
+
+    def _is_feasible(self, sequence: str | np.ndarray) -> bool:
+        """
+        Checks whether a sequence (str or array of one sequence) is feasible
+        under the transition matrix. This is done by looping through
+        the sequence and determining whether the transition probabilities
+        are non-zero.
+        """
+        if isinstance(sequence, np.ndarray):
+            assert sequence.ndim == 1 or sequence.shape[0] == 1
+            sequence = "".join(sequence.flatten())
+
+        current_state = self.alphabet.index(sequence[0])
+        for i in range(1, len(sequence)):
+            next_state = self.alphabet.index(sequence[i])
+
+            if np.isclose(self.transition_matrix[current_state, next_state], 0.0):
+                return False
+            current_state = next_state
+
+        return True
+
+    def construct_random_motifs(
+        self, motif_length: int, n_motifs: int, seed: int = None
+    ) -> np.ndarray:
+        """
+        Creates a given number of random motifs of a certain length.
+        """
+        assert motif_length * n_motifs <= self.sequence_length
+
+        random_state = np.random.RandomState(seed)
+
+        # Sampling a sequence of length motif_length * n_motifs
+        sequence = self._sample_random_sequence(
+            length=motif_length * n_motifs,
+            random_state=random_state,
+            repeating_allowed=False,
+        )
+
+        # Chunking it into n_motifs
+        motifs = np.array(
+            [
+                list(sequence[i * motif_length : (i + 1) * motif_length])
+                for i in range(0, n_motifs)
+            ]
+        )
+
+        return motifs
+
+    def construct_random_offsets(
+        self,
+        motif_length: int,
+        n_motifs: int,
+        seed: int = None,
+    ) -> np.ndarray:
+        """
+        Creates a given number of random offsets for the motifs.
+        """
+        all_motifs_length = motif_length * n_motifs
+
+        # For each motif, we sample weights in the simplex
+        # from a uniform dirichlet
+        random_state = np.random.RandomState(seed)
+
+        offsets = []
+        for _ in range(n_motifs):
+            weights = random_state.dirichlet(np.ones(motif_length - 1))
+            _offset_for_motif = [0]
+            for weight in weights:
+                _offset_for_motif.append(
+                    1
+                    + np.floor(
+                        weight * (self.sequence_length - all_motifs_length) // n_motifs
+                    )
+                )
+
+            offset_for_motif = np.cumsum(np.array(_offset_for_motif, dtype=int))
+
+            offsets.append(offset_for_motif)
+
+        return np.array(offsets)
+
+    def construct_optimal_solution(
+        self, motifs: np.ndarray | None = None, offsets: np.ndarray | None = None
+    ) -> np.ndarray:
+        """
+        Constructs an optimal solution for a given set of motifs and offsets.
+
+        If None are provided, then the motifs and offsets of the black box
+        are used.
+        """
+        if motifs is None:
+            motifs = self.motifs
+
+        if offsets is None:
+            offsets = self.offsets
+
+        # de-cumsum the offsets
+        offsets = np.diff(offsets, prepend=0)
+        optimal_sequence = []
+
+        for motif, offset in zip(motifs, offsets):
+            # We write first l - 1 characters according to the offsets,
+            # and wait to write the last character at the end
+            for character, next_offset_value in zip(motif, offset[1:]):
+                # Put the current character in the current position all the way through just before the offset
+                optimal_sequence += [character] * next_offset_value
+
+            # Write the last character
+            optimal_sequence += [motif[-1]]
+
+        # We pad until the sequence length with the last character
+        # of the last motif
+        optimal_sequence += [motifs[-1][-1]] * (
+            self.sequence_length - len(optimal_sequence)
+        )
+
+        return np.array(optimal_sequence).reshape(1, -1)
+
+    def _maximal_motif_matches(
+        self, sequence: np.ndarray, motif: np.ndarray, offset: np.ndarray
+    ) -> int:
+        """
+        Counts the maximal motif match.
+        """
+        assert sequence.ndim == 1 or sequence.shape[0] == 1
+        sequence = "".join(sequence.flatten())
+        maximal_match = 0
+        for seq_idx in range(len(sequence) - max(offset)):
+            matches = 0
+            sequence_at_offset = np.array(
+                [sequence[seq_idx + offset_value] for offset_value in offset]
+            )
+            matches = sum(sequence_at_offset == motif)
+
+            maximal_match = max(maximal_match, matches)
+
+        return maximal_match
+
+    def _black_box(self, x: np.ndarray, context=None) -> np.ndarray:
+        """
+        Evaluates the sequences in x by checking maximal matches and multiplying.
+        """
+        values = []
+        for sequence in x:
+            if not self._is_feasible(sequence):
+                values.append(-np.inf)
+                continue
+
+            value = 1.0
+            for motif, offset in zip(self.motifs, self.offsets):
+                maximal_matches = self._maximal_motif_matches(sequence, motif, offset)
+                value *= (
+                    maximal_matches // (self.motif_length / self.quantization)
+                ) / self.quantization
+
+            values.append(value)
+
+        return np.array(values).reshape(-1, 1)
+
+    @staticmethod
+    def get_black_box_info() -> BlackBoxInformation:
+        return ehrlich_info
+
+
+class EhrlichProblemFactory(AbstractProblemFactory):
+    """
+    A factory for creating Ehrlich functions and initial conditions.
+
+    References
+    ----------
+    [1] Stanton, S., Alberstein, R., Frey, N., Watkins, A., & Cho, K. (2024).
+        Closed-Form Test Functions for Biophysical Sequence Optimization Algorithms.
+        arXiv preprint arXiv:2407.00236. https://arxiv.org/abs/2407.00236
+    """
+
+    def __init__(self) -> None:
+        super().__init__()
+
+    def get_setup_information(self) -> BlackBoxInformation:
+        return ehrlich_info
+
+    def create(
+        self,
+        sequence_length: int,
+        motif_length: int,
+        n_motifs: int,
+        quantization: int | None = None,
+        seed: int = None,
+        alphabet: list[str] = AMINO_ACIDS,
+        batch_size: int = None,
+        parallelize: bool = False,
+        num_workers: int = None,
+        evaluation_budget: int = float("inf"),
+    ) -> Problem:
+        """
+        Creates an Ehrlich function problem (containing an Ehrlich black box and
+        an initial condition).
+
+        Parameters
+        ----------
+        sequence_length : int
+            The length of the sequence to be optimized. This length is fixed, and
+            _only_ sequences of this length are considered.
+        motif_length : int
+            The length of the motifs.
+        n_motifs : int
+            The number of motifs.
+        quantization : int, optional
+            The quantization parameter. This parameter must be between 1 and the
+            motif length, and the motif length must be divisible by the quantization.
+            By default, it is None (which corresponds to the motif length).
+        seed : int, optional
+            The seed for the random number generator. By default, it is None
+            (i.e. no seed is set).
+        alphabet : list of str, optional
+            The alphabet to be used for the sequences. By default, it is the
+            of 20 amino acids.
+        batch_size : int, optional
+            The batch size for the black box. By default, it is None (i.e. all
+            sequences are evaluated in a vectorized way).
+        parallelize : bool, optional
+            Whether to parallelize the evaluation of the black box. By default,
+            it is False.
+        num_workers : int, optional
+            The number of processors used in parallelization.
+        evaluation_budget : int, optional
+            The evaluation budget for the black box. By default, it is infinite.
+
+        References
+        ----------
+        [1] Stanton, S., Alberstein, R., Frey, N., Watkins, A., & Cho, K. (2024).
+            Closed-Form Test Functions for Biophysical Sequence Optimization Algorithms.
+            arXiv preprint arXiv:2407.00236. https://arxiv.org/abs/2407.00236
+        """
+        if seed is not None:
+            seed_python_numpy_and_torch(seed)
+
+        f = EhrlichBlackBox(
+            sequence_length=sequence_length,
+            motif_length=motif_length,
+            n_motifs=n_motifs,
+            quantization=quantization,
+            seed=seed,
+            alphabet=alphabet,
+            batch_size=batch_size,
+            parallelize=parallelize,
+            num_workers=num_workers,
+            evaluation_budget=evaluation_budget,
+        )
+        x0 = np.array([list(f._sample_random_sequence())])
+
+        return Problem(f, x0)


Test patch:
diff --git a/src/poli/tests/registry/toy_discrete_problems/__init__.py b/src/poli/tests/registry/toy_discrete_problems/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/src/poli/tests/registry/toy_discrete_problems/test_ehrlich.py b/src/poli/tests/registry/toy_discrete_problems/test_ehrlich.py
new file mode 100644
index 00000000..f21415c0
--- /dev/null
+++ b/src/poli/tests/registry/toy_discrete_problems/test_ehrlich.py
@@ -0,0 +1,166 @@
+"""
+This is a suite of tests for the Ehrlich functions
+proposed by Stanton et al. (2024).
+"""
+
+import numpy as np
+
+import pytest
+
+from poli.objective_repository.ehrlich._construct_feasibility_matrix import (
+    _construct_transition_matrix,
+)
+from poli.repository import EhrlichBlackBox, EhrlichProblemFactory
+
+
+@pytest.mark.parametrize("size", [3, 5, 8, 10])
+@pytest.mark.parametrize("seed", [1, 2, 3, 4])
+def test_sparse_matrix_construction_is_ergodic_and_aperiodic(size: int, seed: int):
+    sparse_transition_matrix = _construct_transition_matrix(size, seed=seed)
+
+    # Checking with the Perron-Frobenius condition
+    assert (
+        np.linalg.matrix_power(sparse_transition_matrix, (size - 1) ** 2 + 1) > 0.0
+    ).all()
+
+
+@pytest.mark.parametrize("seed", [1, 4])
+@pytest.mark.parametrize("sequence_length", [10, 20, 50, 100])
+@pytest.mark.parametrize("motif_length", [2, 3, 4, 5])
+@pytest.mark.parametrize("n_motifs", [1, 2, 3, 4])
+def test_ehrlich_motifs_and_offsets_are_deterministic(
+    seed: int, sequence_length: int, motif_length: int, n_motifs: int
+):
+    if n_motifs * motif_length > sequence_length:
+        pytest.skip(
+            "The total length of the motifs is greater than the sequence length."
+        )
+    ehrlich_1 = EhrlichBlackBox(
+        sequence_length=sequence_length,
+        motif_length=motif_length,
+        n_motifs=n_motifs,
+        seed=seed,
+    )
+    ehrlich_2 = EhrlichBlackBox(
+        sequence_length=sequence_length,
+        motif_length=motif_length,
+        n_motifs=n_motifs,
+        seed=seed,
+    )
+
+    assert np.all(ehrlich_1.motifs == ehrlich_2.motifs)
+    assert np.all(ehrlich_1.offsets == ehrlich_2.offsets)
+
+
+@pytest.mark.parametrize("sequence_length", [10, 20, 50, 100])
+@pytest.mark.parametrize("motif_length", [2, 3, 4, 5])
+@pytest.mark.parametrize("n_motifs", [1, 2, 3, 4])
+def test_ehrlich_gives_different_motifs_for_different_seeds(
+    sequence_length, motif_length, n_motifs
+):
+    if n_motifs * motif_length > sequence_length:
+        pytest.skip(
+            "The total length of the motifs is greater than the sequence length."
+        )
+
+    ehrlich_1 = EhrlichBlackBox(
+        sequence_length=sequence_length,
+        motif_length=motif_length,
+        n_motifs=n_motifs,
+        seed=1,
+    )
+    ehrlich_2 = EhrlichBlackBox(
+        sequence_length=sequence_length,
+        motif_length=motif_length,
+        n_motifs=n_motifs,
+        seed=2,
+    )
+
+    assert not np.all(ehrlich_1.motifs == ehrlich_2.motifs)
+
+
+@pytest.mark.parametrize("sequence_length", [10, 20, 50, 100])
+@pytest.mark.parametrize("motif_length", [2, 3, 4, 5])
+@pytest.mark.parametrize("n_motifs", [1, 2, 3, 4])
+def test_ehrlich_function_produces_optimal_sequences(
+    sequence_length: int, motif_length: int, n_motifs: int
+):
+    if n_motifs * motif_length > sequence_length:
+        pytest.skip(
+            "The total length of the motifs is greater than the sequence length."
+        )
+
+    ehrlich = EhrlichBlackBox(
+        sequence_length=sequence_length,
+        motif_length=motif_length,
+        n_motifs=n_motifs,
+        seed=1,
+    )
+
+    optimal_sequence = ehrlich.construct_optimal_solution()
+    assert ehrlich._is_feasible(optimal_sequence)
+    assert ehrlich(optimal_sequence) == 1.0
+
+
+def test_consistency_of_ehrlich_function_motif_matching():
+    ehrlich = EhrlichBlackBox(
+        sequence_length=10,
+        motif_length=3,
+        n_motifs=2,
+        quantization=3,
+        seed=1,
+    )
+
+    one_sequence = np.array(["E"] * 10).reshape(1, 10)
+    motif_matches = ehrlich._maximal_motif_matches(
+        one_sequence, np.array(["E", "V", "D"]), np.array([0, 1, 3])
+    )
+    assert motif_matches == 1
+
+    another_sequence = "EVEEEEEEEE"
+    another_sequence = np.array(list(another_sequence)).reshape(1, 10)
+    assert (
+        ehrlich._maximal_motif_matches(
+            another_sequence, np.array(["E", "V", "D"]), np.array([0, 1, 3])
+        )
+        == 2
+    )
+
+    yet_another_sequence = "EEEEVEDEEE"
+    yet_another_sequence = np.array(list(yet_another_sequence)).reshape(1, 10)
+    assert (
+        ehrlich._maximal_motif_matches(
+            yet_another_sequence, np.array(["E", "V", "D"]), np.array([0, 1, 3])
+        )
+        == 3
+    )
+
+
+def test_creating_a_problem_with_a_factory():
+    problem_factory = EhrlichProblemFactory()
+
+    problem = problem_factory.create(
+        sequence_length=10,
+        motif_length=3,
+        n_motifs=2,
+        quantization=3,
+        seed=1,
+    )
+
+    f, x0 = problem.black_box, problem.x0
+    y0 = f(x0)
+
+
+def test_creating_with_create():
+    from poli import create
+
+    problem = create(
+        name="ehrlich",
+        sequence_length=10,
+        motif_length=3,
+        n_motifs=2,
+        quantization=3,
+        seed=1,
+    )
+    f, x0 = problem.black_box, problem.x0
+    y0 = f(x0)