model

qemcmc.model

EnergyModel

EnergyModel(
    n,
    couplings=[],
    name=None,
    cost_function_signs=[-1, -1],
    model_type="ising",
)

A base class for energy models for a classical energy function over n spins, defined by arbitrary-order coupling tensors.

Parameters:

Name	Type	Description	Default
n	int	Number of spins	required
couplings	List[ndarray]	List of numpy arrays representing interaction tensors. A 1D array encodes linear terms, a 2D array encodes pairwise terms (expected symmetric), and higher-rank tensors encode higher-order interactions.	[]
name	str	Optional label for the model (used in plotting / logging).	None
cost_function_signs	list	Sign convention(s) used by downstream components (e.g. proposal/acceptance conventions). For example, [-1, -1] for linear and quadratic terms.	[-1, -1]
model_type	str	Type of model, either 'ising' or 'binary'. This determines how the binary states are interpreted and how the energy is calculated. 'ising' models use spin values {-1, +1}, while 'binary' models use binary values {0, 1}.	'ising'

Notes

• Energies are computed by mapping binary states {0,1} to spin values {-1,+1} internally.
• Brute-force methods such as get_all_energies scale as O(2^n) and are intended only for small systems.

Source code in src/qemcmc/model/energy_model.py

def __init__(self, n: int, couplings: List[np.ndarray] = [], name: str = None, cost_function_signs: list = [-1, -1], model_type: str = "ising"):
    self.n = n
    self.n_spins = n
    self.couplings = couplings
    self.name = name
    self.alpha = self.calculate_alpha(n, couplings)
    self.lowest_energy = None

    self.normalised_couplings = [self.couplings[i] * self.alpha for i in range(len(self.couplings))]
    self.cost_function_signs = cost_function_signs

    if model_type not in ["ising", "binary"]:
        raise ValueError(f"Invalid model_type '{model_type}'. Expected 'ising' or 'binary'.")
    else:
        self.model_type = model_type

    self.initial_state = self.get_initial_states(num_initial_states=100)

get_initial_states

get_initial_states(num_initial_states)

Generates a list of random initial states.

Parameters:

Name	Type	Description	Default
num_initial_states	int	The number of initial states to generate.	required

Returns:

Name	Type	Description
list		A list of random initial states.

Source code in src/qemcmc/model/energy_model.py

def get_initial_states(self, num_initial_states: int):
    """
    Generates a list of random initial states.

    Parameters
    ----------
    num_initial_states:
        The number of initial states to generate.

    Returns
    -------
    list:
        A list of random initial states.
    """

    init_states = []
    while len(init_states) < num_initial_states:
        state = get_random_state(self.n_spins)
        init_states.append(state)
    return init_states

get_ground_state

get_ground_state(num_reads=100, num_batches=10)

Finds an approximate ground state using Simulated Annealing.

Source code in src/qemcmc/model/energy_model.py

def get_ground_state(self, num_reads=100, num_batches=10):
    """
    Finds an approximate ground state using Simulated Annealing.
    """

    h, J = self.couplings

    h_dict = {i: h[i] for i in range(self.n_spins)}
    J_dict = {(i, j): J[i, j] for i in range(self.n_spins) for j in range(i + 1, self.n_spins)}

    bqm = dimod.BinaryQuadraticModel.from_ising(h_dict, J_dict)

    sampler = dimod.SimulatedAnnealingSampler()
    best_energy = float("inf")

    reads_per_batch = max(1, num_reads // num_batches)

    for _ in tqdm(range(num_batches), desc=f"Annealing ({num_reads} total reads)"):
        response = sampler.sample(bqm, num_reads=reads_per_batch)

        current_best = response.first
        if current_best.energy < best_energy:
            best_energy = current_best.energy

    print("\n--- Simulated Annealing Results ---")
    print(f"Lowest Energy Found: {best_energy:.4f}")

    return best_energy

calculate_energy

calculate_energy(state, couplings, cost_function_signs)

Calculate the energy of a given state for an arbitrary-order Ising/binary model.

Parameters:

Name	Type	Description	Default
state	array - like(str, list, tuple, array)	State configuration. Can be: - Binary: "011", [0,1,1], (0,1,1), etc. - Spin: [-1,1,1], (-1,1,1), etc.	required
couplings	list of numpy arrays	List of coupling tensors where: - 1D arrays represent linear terms (h_i) - 2D arrays represent quadratic terms (J_ij) - 3D arrays represent cubic terms, etc.	required

Returns:

Name	Type	Description
float	Total energy of the state

Source code in src/qemcmc/model/energy_model.py

def calculate_energy(self, state, couplings, cost_function_signs):
    """
    Calculate the energy of a given state for an arbitrary-order Ising/binary model.

    Parameters
    -----------
    state : array-like (str, list, tuple, np.array)

        State configuration. Can be:
        - Binary: "011", [0,1,1], (0,1,1), etc.
        - Spin: [-1,1,1], (-1,1,1), etc.

    couplings : list of numpy arrays

        List of coupling tensors where:
        - 1D arrays represent linear terms (h_i)
        - 2D arrays represent quadratic terms (J_ij)
        - 3D arrays represent cubic terms, etc.

    Returns
    -------
    float : Total energy of the state
    """

    if isinstance(state, str):
        state = np.array([int(bit) for bit in state])
    if isinstance(state, (list, tuple, np.ndarray)):
        state = np.array(state)
        if state.dtype == np.bool_:
            state = state.astype(int)
    else:
        raise TypeError(f"State must be a string, list, tuple, or numpy array, but got {type(state)}")
    if self.model_type == "binary":
        if not np.all(np.isin(state, [0, 1])):
            # If not already in binary format, try interpreting as spin values and converting
            if not np.all(np.isin(state, [-1, 1])):
                # not in spin format either?
                raise ValueError("Spin configuration must be in spin (-1/+1) or binary (0/1) format, but got values outside these sets.")
            else:
                # Convert to binary
                state = np.array([(1-spin) // 2 for spin in state])

    elif self.model_type == "ising":
        if not np.all(np.isin(state, [-1, 1])):
            # If not already in spin format, try interpreting as binary and converting
            if not np.all(np.isin(state, [0, 1])):
                # not in binary format either?
                raise ValueError("Spin configuration must be in binary (0/1) or spin (-1/+1) format, but got values outside these sets.")
            else:
                # Convert to spin
                state = np.array([2 * int(bit) - 1 for bit in state])
                #state = np.array([1-2 * int(bit)  for bit in state])

    total_energy = 0.0
    for term_index, coupling in enumerate(couplings):
        coupling = np.array(coupling)
        order = coupling.ndim
        if order == 0:
            total_energy += cost_function_signs[term_index] * coupling
        elif order == 1:
            total_energy += cost_function_signs[term_index] * np.dot(coupling, state)
        elif order == 2:
            total_energy += cost_function_signs[term_index] * 0.5 * np.einsum("ij, i, j->", coupling, state, state)
        else:
            # General case for any order >=3 (cubic, quartic etc.)
            indices = "".join(chr(97 + i) for i in range(order))  # 'abc...', 'ijkl...'
            einsum_str = f"{indices}," + ",".join([indices[i] for i in range(order)]) + "->"
            coefficient = 1.0 / math.factorial(order)
            total_energy += cost_function_signs[term_index] * coefficient * np.einsum(einsum_str, coupling, *([state] * order))

    return total_energy

get_subgroup_couplings

get_subgroup_couplings(
    subgroup, current_state, coupling_weights=None
)

Calculates local couplings for a subgroup of spins.

Spins outside the specified subgroup are treated as frozen constants, and their values are absorbed into the couplings of the subgroup. This is useful for calculating local effective Hamiltonians.

Parameters:

Name	Type	Description	Default
subgroup	List[int]	A list of indices of the spins in the subgroup.	required
current_state	str	The current state of the full system, as a binary string.	required
coupling_weights	List[float]	Optional weights for the couplings. This can be used to effectively remove certain couplings from the proposal (e.g., constraints) by setting their weight to 0. It is important not to normalize the couplings after applying these weights.	None

Returns:

Type	Description
List[ndarray]	A new list of coupling tensors for the subgroup.

Notes

It might seem off to do the reweighting at this point, but here are reasons why I [SF] did it!
Basically, when you get the subgroup couplings, the terms jumble up, so the returned local
coupling list necessarily does not respect the order etc. of the different terms in the original.

Source code in src/qemcmc/model/energy_model.py

def get_subgroup_couplings(self, subgroup: List[int], current_state: str, coupling_weights: List[float] = None):
    """
    Calculates local couplings for a subgroup of spins.

    Spins outside the specified subgroup are treated as frozen constants, and their
    values are absorbed into the couplings of the subgroup. This is useful for
    calculating local effective Hamiltonians.

    Parameters
    ----------
    subgroup : List[int]
        A list of indices of the spins in the subgroup.
    current_state : str
        The current state of the full system, as a binary string.
    coupling_weights : List[float], optional
        Optional weights for the couplings. This can be used to effectively
        remove certain couplings from the proposal (e.g., constraints) by setting
        their weight to 0. It is important not to normalize the couplings after
        applying these weights.

    Returns
    -------
    List[np.ndarray]
        A new list of coupling tensors for the subgroup.

    Notes
    -----
        It might seem off to do the reweighting at this point, but here are reasons why I [SF] did it!
        Basically, when you get the subgroup couplings, the terms jumble up, so the returned local
        coupling list necessarily does not respect the order etc. of the different terms in the original.
    """

    if coupling_weights is None:
        coupling_weights = [1.0] * len(self.normalised_couplings)

    n_sub = len(subgroup)
    subgroup_set = set(subgroup)
    g_to_l = {g_idx: l_idx for l_idx, g_idx in enumerate(subgroup)}

    # Map bitstring '0'/'1' to spin values -1/+1
    state_vals = np.array([1 if b == "1" else -1 for b in current_state])
    max_order = max(c.ndim for c in self.normalised_couplings)
    new_couplings = [np.zeros((n_sub,) * d) for d in range(1, max_order + 1)]

    for couplings_index, coupling in enumerate(self.normalised_couplings):
        coupling = np.array(coupling) * coupling_weights[couplings_index]
        # only loop over elements that actually exist (non-zero)
        non_zero_indices = np.transpose(np.nonzero(coupling))

        for indices in non_zero_indices:
            indices = tuple(indices)
            coeff = coupling[indices]

            if len(set(indices)) != len(indices):
                continue

            in_group = [i for i in indices if i in subgroup_set]
            out_group = [i for i in indices if i not in subgroup_set]

            # Multiply coefficient by values of fixed spins outside the subgroup
            multiplier = np.prod(state_vals[out_group])
            effective_coeff = coeff * multiplier

            if in_group:
                new_order = len(in_group)
                local_indices = tuple(g_to_l[i] for i in in_group)
                new_couplings[new_order - 1][local_indices] += effective_coeff

    return new_couplings

calculate_alpha

calculate_alpha(n, couplings=None, eps=1e-15)

Compute alpha = sqrt(n) / sqrt(sum of squares of UNIQUE coupling coefficients), assuming coupling tensors are symmetric representations.

Any non-symmetric 2-body input raises ValueError.

Parameters:

Name	Type	Description	Default
n	int	Number of spins.	required
couplings	list[ndarray]	Couplings to use. Defaults to self.couplings if not provided.	None
eps	float	Small threshold to avoid division by zero.	1e-15

Returns:

Type	Description
np.ndarray :	An array of normalising factors for each term in the couplings list.

Source code in src/qemcmc/model/energy_model.py

def calculate_alpha(self, n: int, couplings: list = None, eps: float = 1e-15) -> np.ndarray:
    """
    Compute alpha = sqrt(n) / sqrt(sum of squares of UNIQUE coupling coefficients),
    assuming coupling tensors are symmetric representations.

    Any non-symmetric 2-body input raises ValueError.

    Parameters
    ----------
    n : int
        Number of spins.
    couplings : list[np.ndarray], optional
        Couplings to use. Defaults to self.couplings if not provided.
    eps : float
        Small threshold to avoid division by zero.

    Returns
    -------
    np.ndarray :
        An array of normalising factors for each term in the couplings list.
    """

    if couplings is None:
        couplings = self.couplings

    norm_sq_arr = np.zeros(len(couplings))

    for T_ind, T in enumerate(couplings):
        norm_sq = 0.0
        T = np.asarray(T)
        order = T.ndim

        if order == 0:
            pass

        # 1-body: h_i
        elif order == 1:
            if T.shape != (n,):
                raise ValueError(f"1-body tensor has shape {T.shape}, expected ({n},)")
            for i in range(n):
                c = float(T[i])
                norm_sq += c * c

        # 2-body: symmetric J
        elif order == 2:
            if T.shape != (n, n):
                raise ValueError(f"2-body tensor has shape {T.shape}, expected ({n},{n})")

            # Enforce symmetry (reject pure upper/lower triangular)
            if not np.allclose(T, T.T):
                raise ValueError("Non-symmetric J provided. This alpha function only accepts symmetric J.")

            # Count each interaction once: i<j
            for i in range(n):
                for j in range(i + 1, n):
                    c = float(T[i, j])
                    if c != 0.0:
                        norm_sq += c * c

        # Order >= 3: count each unordered interaction once using i1<i2<...<ik
        else:
            if T.shape != (n,) * order:
                raise ValueError(f"{order}-body tensor has shape {T.shape}, expected {(n,) * order}")

            for comb in itertools.combinations(range(n), order):
                c = float(T[comb])
                if c != 0.0:
                    norm_sq += c * c

        norm_sq_arr[T_ind] = norm_sq

    norm_sq_tot = np.sum(norm_sq_arr)
    if norm_sq_tot < eps:
        print("Warning: Couplings are all zero (or very close to zero). Returning alpha=1 to avoid division by zero.")
        print("Input of a zero-coupling may be intended to sample Uniform distributions, however if this is not your intention, please check your couplings.")
        return np.ones(len(couplings))

    return np.sqrt(n) / np.sqrt(norm_sq_tot)

get_energy

get_energy(state)

Returns the energy of a given state

Source code in src/qemcmc/model/energy_model.py

def get_energy(self, state: str) -> float:
    """
    Returns the energy of a given state
    """

    if not isinstance(state, str):
        raise TypeError(f"State must be a string, but got {type(state)}")
    energy = self.calculate_energy(state, self.couplings, self.cost_function_signs)
    return energy

get_all_energies

get_all_energies()

Calculate the energies for all possible spin states. This method generates all possible spin states for the system and calculates the energy for each state.

Returns:

Type	Description
np.ndarray :	An array containing the energies of all possible spin states.

Source code in src/qemcmc/model/energy_model.py

def get_all_energies(self) -> np.ndarray:
    """
    Calculate the energies for all possible spin states.
    This method generates all possible spin states for the
    system and calculates the energy for each state.

    Returns
    -------
    np.ndarray :
        An array containing the energies of all possible spin states.
    """

    self.S = ["".join(i) for i in itertools.product("01", repeat=self.n)]
    all_energies = np.zeros(len(self.S))
    for state in self.S:
        all_energies[int(state, 2)] = self.calculate_energy(state, self.couplings, self.cost_function_signs)
    return all_energies

get_lowest_energies

get_lowest_energies(
    num_states, return_configurations=False
)

Retrieve the lowest energy states and their degeneracies. This method computes all possible energies and then finds the specified number of lowest energy states along with their degeneracies.

Parameters:

Name	Type	Description	Default
num_states	int		required
return_configurations	bool		False

Returns:

Type	Description
Tuple[ndarray, ndarray]	The first array contains the lowest energy values. The second array contains the degeneracies of the corresponding energy values.

Notes

This method is intended for small instances due to its brute-force nature, which is extremely memory intensive and slow.

Source code in src/qemcmc/model/energy_model.py

def get_lowest_energies(self, num_states: int, return_configurations: bool = False) -> typing.Tuple[np.ndarray, np.ndarray]:
    """
    Retrieve the lowest energy states and their degeneracies.
    This method computes all possible energies and then finds the specified number
    of lowest energy states along with their degeneracies.

    Parameters
    ----------
    num_states (int): The number of lowest energy states to retrieve.
    return_configurations (bool): Whether to also return the corresponding configurations of the lowest energy states. Defaults to False.

    Returns
    -------
    Tuple[np.ndarray, np.ndarray]
        - The first array contains the lowest energy values.
        - The second array contains the degeneracies of the corresponding energy values.

    Notes
    -----
        This method is intended for small instances due to its brute-force nature, which is extremely memory intensive and slow.
    """

    # only to be used for small instances, it is just brute force so extremely memory intensive and slow
    all_energies = self.get_all_energies()

    # very slow (sorts whole array)
    lowest_energies, lowest_energy_degeneracy = self.find_lowest_values(all_energies, num_values=num_states)
    if return_configurations:
        lowest_configs = []
        for energy in lowest_energies:
            configs = [self.S[i] for i, e in enumerate(all_energies) if e == energy]
            lowest_configs.append(configs)
        self.lowest_energy = lowest_energies[0]
        return lowest_energies, lowest_energy_degeneracy, lowest_configs
    else:
        self.lowest_energy = lowest_energies[0]
        return lowest_energies, lowest_energy_degeneracy

find_lowest_values

find_lowest_values(arr, num_values=5)

Find the lowest unique values in an array and their degeneracies.

Parameters:

Name	Type	Description	Default
arr	ndarray		required
num_values	int		5

Returns:

Type	Description
Tuple[lowest_values(ndarray), degeneracy(ndarray)]	lowest_values (np.ndarray): The lowest unique values in the array. degeneracy (np.ndarray): The counts of each of the lowest unique values.

Source code in src/qemcmc/model/energy_model.py

def find_lowest_values(self, arr: np.ndarray, num_values: int = 5):
    """
    Find the lowest unique values in an array and their degeneracies.

    Parameters
    ----------
    arr (np.ndarray): The input array from which to find the lowest values.
    num_values (int, optional): The number of lowest unique values to find. Defaults to 5.

    Returns
    -------
    Tuple[lowest_values (np.ndarray), degeneracy (np.ndarray)]
        - lowest_values (np.ndarray): The lowest unique values in the array.
        - degeneracy (np.ndarray): The counts of each of the lowest unique values.
    """

    # Count the occurrences of each value
    unique_values, counts = np.unique(arr, return_counts=True)
    # Sort the unique values and counts by value
    sorted_indices = np.argsort(unique_values)
    unique_values_sorted = unique_values[sorted_indices]
    counts_sorted = counts[sorted_indices]
    # Find the first num_values
    lowest_values = unique_values_sorted[:num_values]
    degeneracy = counts_sorted[:num_values]
    return lowest_values, degeneracy

get_lowest_energy

get_lowest_energy()

Calculate and return the lowest energy from all possible energies. This method uses a brute force approach to find the lowest energy, making it extremely memory intensive and slow. It is recommended to use this method only for small instances.

Returns:

Name	Type	Description
float	The lowest energy value.

Notes

If the lowest energy has already been calculated and stored in self.lowest_energy, it will return that value directly to save computation time.

Source code in src/qemcmc/model/energy_model.py

def get_lowest_energy(self):
    """
    Calculate and return the lowest energy from all possible energies.
    This method uses a brute force approach to find the lowest energy,
    making it extremely memory intensive and slow. It is recommended
    to use this method only for small instances.

    Returns
    -------
    float: The lowest energy value.

    Notes
    -----
    If the lowest energy has already been calculated and stored
    in `self.lowest_energy`, it will return that value directly
    to save computation time.
    """

    # Only to be used for small instances, it is just brute force so extremely memory intensive and slow
    if self.lowest_energy is not None:
        return self.lowest_energy
    else:
        all_energies = self.get_all_energies()

    lowest_energy = np.min(all_energies)
    self.lowest_energy = lowest_energy
    return lowest_energy

get_boltzmann_factor

get_boltzmann_factor(state, beta=1.0)

Get un-normalised boltzmann probability of a given state

Parameters:

Name	Type	Description	Default
state	str	configuration of spins for which probability is to be calculated	required
beta	float	inverse temperature (1/T) at which the probability is to be calculated.	1.0

Returns:

Type	Description
float corresponding to the un-normalised boltzmann probability of the given state.

Source code in src/qemcmc/model/energy_model.py

def get_boltzmann_factor(self, state: str, beta: float = 1.0) -> float:
    """
    Get un-normalised boltzmann probability of a given state

    Parameters
    ----------
    state : str
        configuration of spins for which probability is to be calculated
    beta : float
        inverse temperature (1/T) at which the probability is to be calculated.

    Returns
    -------
        float corresponding to the un-normalised boltzmann probability of the given state.
    """

    E = self.get_energy(state)
    boltzmann_factor = np.exp(-1 * beta * E, dtype=np.longdouble)

    return boltzmann_factor

get_boltzmann_factor_from_energy

get_boltzmann_factor_from_energy(E, beta=1.0)

Get un-normalized Boltzmann probability for a given energy.

Parameters:

Name	Type	Description	Default
E	float	Energy for which the Boltzmann factor is to be calculated.	required
beta	float	Inverse temperature (1/T) at which the probability is to be calculated.	1.0

Returns:

Name	Type	Description
float	float	The un-normalized Boltzmann probability for the given energy.

Source code in src/qemcmc/model/energy_model.py

def get_boltzmann_factor_from_energy(self, E, beta: float = 1.0) -> float:
    """
    Get un-normalized Boltzmann probability for a given energy.

    Parameters
    ----------
    E : float
        Energy for which the Boltzmann factor is to be calculated.
    beta : float
        Inverse temperature (1/T) at which the probability is to be calculated.

    Returns
    -------
    float:
        The un-normalized Boltzmann probability for the given energy.
    """

    return np.exp(-1 * beta * E, dtype=np.longdouble)

ModelMaker

ModelMaker(
    n_spins, model_type, name, cost_function_signs=[-1, -1]
)

Utility class for constructing standard Ising energy models used in simulations and experiments.

This class constructs predefined random energy models used for testing and experimentation. Depending on the chosen model_type, it generates coupling tensors and initialises an :class:EnergyModel instance.

Source code in src/qemcmc/model/model_maker.py

def __init__(self, n_spins: int, model_type: str, name: str, cost_function_signs: list = [-1, -1]):
    self.name = name
    self.n_spins = n_spins
    self.cost_function_signs = cost_function_signs or [-1, -1]

    if not isinstance(model_type, str):
        raise TypeError("model_type must be a string")

    if model_type == "Fully Connected Ising":
        self.make_fully_connected_ising()
    elif model_type == "Fully Connected QUBO":
        self.make_fully_connected_binary()
    else:
        raise ValueError(f"Unknown model_type: {model_type}")

make_fully_connected_binary

make_fully_connected_binary()

Transforms the existing Ising couplings into an mathematically equivalent QUBO model via s = 2x - 1.

Source code in src/qemcmc/model/model_maker.py

def make_fully_connected_binary(self):
    """
    Transforms the existing Ising couplings into an mathematically 
    equivalent QUBO model via s = 2x - 1.
    """
    couplings = self.make_fully_connected_ising(return_couplings=True)
    # all_energies_ising = self.model.get_all_energies()


    h,J = couplings

    #constant = 2*np.sum(J)-np.sum(h)
    Q_binary = 4 * J
    q_binary = 2 * h - 2 * np.sum(J, axis=1)
    binary_couplings = [np.round(q_binary, 4), np.round(Q_binary, 4)]

    self.model = EnergyModel(
        n=self.n_spins, 
        couplings=binary_couplings, 
        name=self.name, 
        cost_function_signs = [-1,-1],
        model_type="binary"
    )

ConstraintModel

ConstraintModel(
    n,
    constraint_couplings,
    constraint_signs,
    couplings,
    constraint_func,
    **kwargs
)

Bases: EnergyModel

A subclass of EnergyModel that incorporates a constraint function to define valid configurations.

The constraint function takes a state as input and returns True if the state is valid (satisfies the constraint) and False otherwise. The energy of invalid states is set to infinity, effectively excluding them from the Boltzmann distribution.

Parameters:

Name	Type	Description	Default
n	int	Number of spins in the model.	required
constraint_couplings	list[ndarray]	List of coupling tensors (numpy arrays) defining the constraint.	required
constraint_signs	list[float]	Sign convention(s) for the constraint couplings.	required
couplings	list[ndarray]	List of coupling tensors (numpy arrays) defining the energy function.	required
constraint_func	Callable[[str], bool]	A function that takes a state (string representation of spin configuration) and returns True if the state satisfies the constraint, and False otherwise.	required
name	str	Optional label for the model (used in plotting / logging).	required
cost_function_signs	list[float]	Sign convention(s) used by downstream components (e.g. proposal/acceptance conventions).	required
model_type	str	Type of model, either 'ising' or 'binary'. This determines how the binary states are interpreted and how the energy is calculated. 'ising' models use spin values {-1, +1}, while 'binary' models use binary values {0, 1}.	required

Notes

• The energy of any state that does not satisfy the constraint is set to infinity, which means such states will have zero probability in the Boltzmann distribution.
• This class can be used to model systems with hard constraints on the configurations, such as certain combinatorial optimization problems or physical systems with forbidden states.

Source code in src/qemcmc/model/constraint_model.py

def __init__(self, n: int, constraint_couplings: list[np.ndarray], constraint_signs: list[float], couplings: list[np.ndarray], constraint_func: typing.Callable[[str], bool], **kwargs):
    if not callable(constraint_func):
        raise ValueError("constraint_func must be a callable function that takes a state as input and returns True/False.")
    self.constraint_func = constraint_func

    # This model requires a special way to generate initial states that respect the constraint.
    self.get_initial_states = self.get_initial_states_constraint

    if couplings is not None:
        super().__init__(n=n, couplings=couplings, **kwargs)
    else:
        super().__init__(n=n, couplings=[], **kwargs)

    self.constraint_couplings = constraint_couplings
    self.constraint_signs = constraint_signs

    # Calculate normalization factors for constraint couplings
    self.constraint_coupling_alpha = self.calculate_alpha(n, constraint_couplings)

    # These are the couplings used in quantum proposals
    # Combine and normalize the energy and constraint couplings
    self.normalised_couplings = (
        [self.couplings[i] * self.alpha for i in range(len(self.couplings))]
        +
        # [self.constraint_couplings[i] for i in range(len(self.constraint_couplings))]
        [self.constraint_couplings[i] * self.constraint_coupling_alpha for i in range(len(self.constraint_couplings))]
    )

    # Store the un-normalized total couplings
    self.total_couplings = self.couplings + self.constraint_couplings

get_initial_states_constraint

get_initial_states_constraint(num_initial_states)

Generates a list of random initial states that satisfy the constraint function.

Parameters:

Name	Type	Description	Default
num_initial_states	int	The number of initial states to generate.	required

Returns:

Name	Type	Description
list		A list of random initial states that are valid according to the constraint.

Source code in src/qemcmc/model/constraint_model.py

def get_initial_states_constraint(self, num_initial_states: int):
    """
    Generates a list of random initial states that satisfy the constraint function.

    Parameters
    ----------
    num_initial_states:
        The number of initial states to generate.

    Returns
    -------
    list:
        A list of random initial states that are valid according to the constraint.
    """
    init_states = []
    counter = 0
    max_attempts = 1000
    while len(init_states) < num_initial_states and counter < max_attempts:
        state = get_random_state(self.n_spins)
        if self.constraint_func(state):
            init_states.append(state)
        counter += 1

    if len(init_states) < num_initial_states:
        warnings.warn(
            f"Could not find enough valid initial states satisfying the constraint. "
            f"Found {len(init_states)} valid states after {max_attempts} attempts. "
            f"You may want to provide them manually if more are needed."
        )
    return init_states

get_constraint_energy

get_constraint_energy(state)

Calculate the energy contribution from the constraint couplings for a given state.

Parameters:

Name	Type	Description	Default
state	str	The state for which to calculate the constraint energy.	required

Returns:

Name	Type	Description
float	float	The energy contribution from the constraint couplings for the given state.

Source code in src/qemcmc/model/constraint_model.py

def get_constraint_energy(self, state: str) -> float:
    """
    Calculate the energy contribution from the constraint couplings for a given state.

    Parameters
    ----------
    state : str
        The state for which to calculate the constraint energy.

    Returns
    -------
    float:
        The energy contribution from the constraint couplings for the given state.
    """
    return self.calculate_energy(state, self.constraint_couplings, self.constraint_signs)

get_total_energy

get_total_energy(state)

Calculate the total energy of a given state, including both the regular energy and the constraint energy.

Parameters:

Name	Type	Description	Default
state	str	The state for which to calculate the total energy.	required

Returns:

Name	Type	Description
float	float	The total energy of the given state, including contributions from both the regular couplings and the constraint couplings.

Source code in src/qemcmc/model/constraint_model.py

def get_total_energy(self, state: str) -> float:
    """
    Calculate the total energy of a given state, including both the regular energy and the constraint energy.

    Parameters
    ----------
    state : str
        The state for which to calculate the total energy.

    Returns
    -------
    float:
        The total energy of the given state, including contributions from both the regular couplings and the constraint couplings.
    """
    return self.get_energy(state) + self.get_constraint_energy(state)