Skip to content

qemcmc

qemcmc

Proposal 

Proposal(model)

Bases: ABC

Abstract base class for producing proposals for Markov Chain Monte Carlo algorithms.

Subclasses implement the proposal mechanism by defining an update(state) method that generates a candidate state from the current one (e.g. single-spin flips, block updates, or quantum proposals).

Parameters:

Name Type Description Default
model EnergyModel

Energy model defining the target distribution over spin configurations.

 required
Source code in src/qemcmc/sampler/proposal.py 
20
21
22
def __init__(self, model: EnergyModel):
    self.model = model
    self.n_spins = model.n_spins

update abstractmethod 

update(state)

Generate a candidate state from the current state using the proposal mechanism.

This method should be implemented by subclasses to define the specific proposal strategy (e.g., single-spin flips, block updates, or quantum proposals).

Parameters:

Name Type Description Default
state MCMCState

The current state of the Markov chain.

 required

Returns:

Type Description
MCMCState

A new candidate state.
Source code in src/qemcmc/sampler/proposal.py 
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
@abstractmethod
def update(self, state: MCMCState) -> MCMCState:
    """
    Generate a candidate state from the current state using the proposal mechanism.

    This method should be implemented by subclasses to define the specific proposal strategy
    (e.g., single-spin flips, block updates, or quantum proposals).

    Parameters
    ----------
    state : MCMCState
        The current state of the Markov chain.

    Returns
    -------
    MCMCState
        A new candidate state.
    """
    pass

ClassicalProposal 

ClassicalProposal(model, method='uniform')

Bases: Proposal

Classical Markov Chain Monte Carlo proposer.

This class implements purely classical proposal mechanisms for MCMC. New candidate states are generated either by sampling a completely random (uniform) configuration, or by performing a local single-spin or two-spin flip.

Parameters:

Name Type Description Default
model EnergyModel

Energy model defining the target Boltzmann distribution.

 required
method str

Proposal mechanism used to generate candidate states.

  • "uniform" : propose a completely random spin configuration.
  • "local" : flip a single randomly chosen spin.
  • "2-local" : flip two randomly chosen spins.

Default is "uniform".

 'uniform'
Source code in src/qemcmc/sampler/classical_proposal.py 
32
33
34
35
36
def __init__(self, model: EnergyModel, method: str = "uniform"):
    if method not in self.METHODS:
        raise ValueError(f"Method '{method}' is not supported. Choose from 'uniform', 'local', or '2-local'.")
    super().__init__(model)
    self.method = method

update_uniform 

update_uniform(current_state_bitstring)

Proposes a new state by generating a random bitstring.

Parameters:

Name Type Description Default
current_state_bitstring str

The current state represented as a bitstring (not used, but required by the API).

 required

Returns:

Type Description
str

A new random state bitstring of the same length.
Source code in src/qemcmc/sampler/classical_proposal.py 
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
def update_uniform(self, current_state_bitstring: str) -> str:
    """
    Proposes a new state by generating a random bitstring.

    Parameters
    ----------
    current_state_bitstring : str
        The current state represented as a bitstring (not used, but required by the API).

    Returns
    -------
    str
        A new random state bitstring of the same length.
    """

    return get_random_state(len(current_state_bitstring))

update_local 

update_local(current_state_bitstring)

Proposes a new state by flipping a single randomly chosen spin.

Parameters:

Name Type Description Default
current_state_bitstring str

The current state represented as a bitstring.

 required

Returns:

Type Description
str

The new state bitstring after flipping one spin.
Source code in src/qemcmc/sampler/classical_proposal.py 
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
def update_local(self, current_state_bitstring: str) -> str:
    """
    Proposes a new state by flipping a single randomly chosen spin.

    Parameters
    ----------
    current_state_bitstring : str
        The current state represented as a bitstring.

    Returns
    -------
    str
        The new state bitstring after flipping one spin.
    """

    # Randomly choose which spin to flip
    choice = np.random.randint(0, self.n_spins)

    # Flip the chosen spin
    c_s = list(current_state_bitstring)
    c_s[choice] = str(int(c_s[choice]) ^ 1)

    # Return the new state as a bitstring
    return "".join(c_s)

update_2local 

update_2local(current_state_bitstring)

Proposes a new state by flipping two randomly chosen spins.

Parameters:

Name Type Description Default
current_state_bitstring str

The current state represented as a bitstring.

 required

Returns:

Type Description
str

The new state bitstring after flipping two spins.
Source code in src/qemcmc/sampler/classical_proposal.py 
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
def update_2local(self, current_state_bitstring: str) -> str:
    """
    Proposes a new state by flipping two randomly chosen spins.

    Parameters
    ----------
    current_state_bitstring : str
        The current state represented as a bitstring.

    Returns
    -------
    str
        The new state bitstring after flipping two spins.
    """
    # Randomly choose which two spins to flip
    choices = np.random.choice(self.n_spins, size=2, replace=False)

    # Flip the chosen spins
    c_s = list(current_state_bitstring)
    for choice in choices:
        c_s[choice] = str(int(c_s[choice]) ^ 1)

    # Return the new state as a bitstring
    return "".join(c_s)

QeProposal 

QeProposal(
    model,
    gamma,
    time,
    r=None,
    delta_t=None,
    coarse_graining=None,
    coupling_weights=None,
    m=1,
)

Bases: Proposal

Quantum-enhanced Markov Chain Monte Carlo sampler.

This class implements the proposal mechanism of the quantum-enhanced MCMC algorithm. Candidate states are generated by simulating the quantum time evolution of a transverse-field Hamiltonian.

The quantum proposal circuit is constructed using :class:CircuitMaker and can optionally operate on coarse-grained subgroups of spins to improve scalability.

Parameters:

Name Type Description Default
model EnergyModel

Energy model defining the target Boltzmann distribution.

 required
gamma float | tuple[float, float]

Transverse field strength (Γ). If a tuple is provided, a value is sampled uniformly from the range [min, max] at each step.

 required
time float | tuple[float, float]

Total evolution time. If a tuple is provided, a value is sampled uniformly from the range [min, max] at each step.

 required
r int | None

Number of Trotter steps. Mutually exclusive with delta_t. If provided, the Trotter step size is derived as Δt = t / r and r stays fixed while Δt varies with t. A warning is issued at init if the resulting Δt range falls outside [0.1, 2.0]. By default None.

 None
delta_t float | None

Trotter step size. Mutually exclusive with r. If provided, the number of Trotter steps is derived as r = floor(t / Δt) per step, so r varies with t while Δt stays fixed. A warning is issued if Δt falls outside [0.1, 2.0]. By default None.

 None
coarse_graining CoarseGraining | None

A coarse-graining scheme to define spin subgroups. If None, no coarse-graining is used. By default None.

 None
coupling_weights list[float | tuple[float, float]] | None

Weights for the coupling tensors. If a tuple is provided, a weight is sampled uniformly from the range. This allows for dynamic adjustment of the influence of different terms in the Hamiltonian during the proposal step. By default None. Note that this is further adjusted by (1 - gamma) to balance the influence of the problem Hamiltonian with the mixer term. Divide by (1 - gamma) if needed. Also, note that the coupling weights should include the mixing term.

 None
m int

Number of subgroups to partition the spins into for sequential updates. By default 1.

 1
Notes

The proposal step simulates the time evolution under a transverse-field Hamiltonian defined by the energy model and measures the resulting state to produce a candidate configuration. This proposal is then accepted or rejected using the Metropolis criterion to ensure convergence to the target Boltzmann distribution.

Exactly one of r or delta_t should be passed in as an argument. If neither is given, a default Δt of 0.8 is used.

Source code in src/qemcmc/sampler/qe_proposal.py 
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
def __init__(
    self,
    model: EnergyModel,
    gamma: float | tuple[float, float],
    time: float | tuple[float, float],
    r: Optional[int] = None,
    delta_t: Optional[float] = None,
    coarse_graining: Optional[CoarseGraining] = None,
    coupling_weights: Optional[list[float | tuple[float, float]]] = None,
    m: int = 1,
):
    super().__init__(model)

    if coupling_weights is not None:
        if len(coupling_weights) != len(model.normalised_couplings):
            raise ValueError(
                f"Length of coupling_weights ({len(coupling_weights)}) must match the number of couplings in the model ({len(model.normalised_couplings)}), including constraint terms."
            )
        self.coupling_weights = coupling_weights
    else:
        self.coupling_weights = list(np.ones(len(model.normalised_couplings)))

    self.gamma = self._validate_gamma(gamma)
    self.time = self._validate_time(time)
    self.r = self._validate_r(r)
    self.delta_t = self._validate_delta_t(delta_t)
    self._validate_trotter_params()
    self.m = m
    self.method = "quantum"

    self.CM = CircuitMaker(self.model)

    self.cg = coarse_graining or CoarseGraining(model.n_spins)

update 

update(current_state)

Perform 'm' sequential quantum updates across non-overlapping subgroups.

Parameters:

Name Type Description Default
current_state str

Current state of the system as a bitstring.

 required

Returns:

Type Description
str

Updated state of the system as a bitstring.
Source code in src/qemcmc/sampler/qe_proposal.py 
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
def update(self, current_state: str) -> str:
    """
    Perform 'm' sequential quantum updates across non-overlapping subgroups.

    Parameters
    ----------
    current_state : str
        Current state of the system as a bitstring.

    Returns
    -------
    str
        Updated state of the system as a bitstring.
    """

    if not isinstance(current_state, str):
        raise TypeError(f"Bitstring must be of type str, got {type(current_state)}: {current_state!r}")

    # Generate m disjoint partitions of the spins
    partitions = self.cg.get_partitions(m=self.m)

    final_coupling_weights, g, t, r = self.sample_hyperparams()

    working_state = current_state
    for subgroup in partitions:
        # Recalculate local couplings as spins outside the subgroup may have flipped
        local_couplings = self.model.get_subgroup_couplings(subgroup=subgroup, current_state=working_state, coupling_weights=final_coupling_weights)
        working_state = self.CM.update(s=working_state, subgroup_choice=subgroup, local_couplings=local_couplings, gamma=g, time=t, r=r)
    return working_state

sample_hyperparams 

sample_hyperparams()

Sample hyperparameters (coupling weights, gamma, time, r) for the current proposal step.

Returns:

Type Description
tuple[list[float], float, float, int]

(final_coupling_weights, gamma, time, r) for this proposal step.
Source code in src/qemcmc/sampler/qe_proposal.py 
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
def sample_hyperparams(self) -> tuple[list[float], float, float, int]:
    """
    Sample hyperparameters (coupling weights, gamma, time, r) for the current proposal step.

    Returns
    -------
    tuple[list[float], float, float, int]
        (final_coupling_weights, gamma, time, r) for this proposal step.
    """

    # Sample gamma and time once per full proposal step
    g = self.gamma if not isinstance(self.gamma, tuple) else np.random.uniform(*self.gamma)
    t = self.time if not isinstance(self.time, tuple) else np.random.uniform(*self.time)

    # Determine r (number of Trotter steps) for this step
    if self.r is not None:
        # Fixed r mode: r stays constant, Δt varies with t
        r = self.r
        dt = t / r
    else:
        # Fixed Δt mode (explicit or default): r varies with t
        dt = self.delta_t if self.delta_t is not None else DEFAULT_DELTA_T
        r = max(1, int(np.floor(t / dt)))

    # Pre-sample weights for all unique coupling weight settings (tuples or floats)
    unique_items = sorted(list(set(self.coupling_weights)), key=lambda x: str(x))
    unique_index = [[i for i, x in enumerate(self.coupling_weights) if x == item] for item in unique_items]

    sampled_weights = np.ones(len(unique_items), dtype=float)
    for i, item in enumerate(unique_items):
        if isinstance(item, (int, float)):
            sampled_weights[i] = item
        elif isinstance(item, tuple):
            sampled_weights[i] = np.random.uniform(*item)

    # Adjust the sampled weights by (1 - gamma) to balance their influence with the mixer term in the Hamiltonian
    sampled_weights_adjusted_by_gamma = (1 - g) * sampled_weights  # Adjust weights by (1 - gamma) to balance with the mixer term

    # Construct the final list of coupling weights for this update step
    final_coupling_weights = np.ones(len(self.coupling_weights))
    for i, indices in enumerate(unique_index):
        for index in indices:
            final_coupling_weights[index] = sampled_weights_adjusted_by_gamma[i]

    return final_coupling_weights, g, t, r

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 
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
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 
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
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 
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
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 
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
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 
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
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 
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
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 
341
342
343
344
345
346
347
348
349
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 
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
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 
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
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 
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
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 
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
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 
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
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 
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
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)

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 
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
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 
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
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 
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
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 
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
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)

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 
19
20
21
22
23
24
25
26
27
28
29
30
31
32
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 
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
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"
    )

CircuitMaker 

CircuitMaker(model)

Constructs and simulates quantum circuits used to generate QeMCMC proposals.

This class builds the Hamiltonian corresponding to a given energy model and simulates its time evolution using PennyLane. Starting from a classical bitstring configuration, the circuit performs Trotterised quantum evolution and samples a new configuration from the resulting quantum state.

The generated sample serves as the proposal state in the quantum-enhanced MCMC algorithm.

Parameters:

Name Type Description Default
model EnergyModel

Energy model defining the problem Hamiltonian.

 required
Notes

The total Hamiltonian simulated by the circuit is

``H = γ H_mixer + (1 - γ) α H_problem``

where H_problem encodes the classical energy model and H_mixer corresponds to a transverse-field term. The evolution time t and number of Trotter steps r are passed per call, giving an effective Trotter step size of δt = t / r. The evolution is approximated using Trotterisation via qml.ApproxTimeEvolution.

Source code in src/qemcmc/circuits.py 
40
41
42
43
44
45
46
47
def __init__(self, model: EnergyModel):
    self.model = model
    self.n_qubits = model.n
    self.dev = qml.device("lightning.qubit", wires=self.n_qubits)
    # self.dev = qml.device("default.tensor", wires=self.n_qubits, method="mps", max_bond_dim=2, contract="auto-mps")

    self.model_type = model.model_type
    self.devices = {}  # cache devices for dynamic subgroup sizes if needed

get_problem_hamiltonian 

get_problem_hamiltonian(couplings, sign=1)

Construct the problem Hamiltonian from symmetric coupling tensors.

This method supports both 'ising' (-1/+1) and 'binary' (0/1) models.

Parameters:

Name Type Description Default
couplings List[ndarray]

A list of coupling tensors.

 required
sign int

A sign to apply to the Hamiltonian. Default is 1.

 1

Returns:

Type Description
Hamiltonian

The problem Hamiltonian.
Source code in src/qemcmc/circuits.py 
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
def get_problem_hamiltonian(self, couplings: List[np.ndarray], sign: int = 1) -> qml.Hamiltonian:
    """
    Construct the problem Hamiltonian from symmetric coupling tensors.

    This method supports both 'ising' (-1/+1) and 'binary' (0/1) models.

    Parameters
    ----------
    couplings : List[np.ndarray]
        A list of coupling tensors.
    sign : int, optional
        A sign to apply to the Hamiltonian. Default is ``1``.

    Returns
    -------
    qml.Hamiltonian
        The problem Hamiltonian.
    """
    total_hamiltonian = 0.0 * qml.Identity(0)
    for coupling_tensor in couplings:
        coupling_tensor = np.asarray(coupling_tensor)
        order = coupling_tensor.ndim
        if order == 0:
            continue

        spin_sign = (-1) ** order if self.model_type == "ising" else 1
        non_zero_indices = np.transpose(np.nonzero(coupling_tensor))
        for index_tuple in non_zero_indices:
            index_tuple = tuple(int(i) for i in index_tuple)

            if len(set(index_tuple)) != len(index_tuple):  # skip repeated indices
                continue
            if index_tuple != tuple(sorted(index_tuple)):  # keep only strictly increasing i1 < i2 < ... < ik
                continue

            coeff = float(coupling_tensor[index_tuple])
            if coeff == 0.0:
                continue

            if self.model_type == "ising":
                term = qml.PauliZ(index_tuple[0])
                for q in index_tuple[1:]:
                    term = term @ qml.PauliZ(q)
                total_hamiltonian += (sign * spin_sign * coeff) * term

            elif self.model_type == "binary":
                # 0.5 * (I - Z) for first variable
                #term = 0.5 * (qml.Identity(index_tuple[0]) + qml.PauliZ(index_tuple[0]))
                term = 0.5 * (qml.Identity(index_tuple[0]) - qml.PauliZ(index_tuple[0]))
                # multiply by 0.5 * (I - Z) for rest
                for q in index_tuple[1:]:
                    next_var = 0.5 * (qml.Identity(q) - qml.PauliZ(q))
                    #next_var = 0.5 * (qml.Identity(q) + qml.PauliZ(q))
                    term = term @ next_var

                total_hamiltonian += (sign * coeff) * term
    simplified_H = qml.simplify(total_hamiltonian)

    coeffs, ops = simplified_H.terms()
    return qml.Hamiltonian(coeffs, ops)

get_mixer_hamiltonian 

get_mixer_hamiltonian(num_wires=None)

Constructs the Mixer Hamiltonian: Σ X_i.

This can be for the full system or a subgroup.

Parameters:

Name Type Description Default
num_wires int

The number of wires (qubits) for the mixer. If None, uses the total number of qubits.

 None

Returns:

Type Description
Hamiltonian

The mixer Hamiltonian.
Source code in src/qemcmc/circuits.py 
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
def get_mixer_hamiltonian(self, num_wires: int = None) -> qml.Hamiltonian:
    """
    Constructs the Mixer Hamiltonian: Σ X_i.

    This can be for the full system or a subgroup.

    Parameters
    ----------
    num_wires : int, optional
        The number of wires (qubits) for the mixer. If None, uses the total number of qubits.

    Returns
    -------
    qml.Hamiltonian
        The mixer Hamiltonian.
    """
    if num_wires is None:
        num_wires = self.n_qubits
    return qml.Hamiltonian([1.0] * num_wires, [qml.PauliX(i) for i in range(num_wires)])

get_state_vector 

get_state_vector(s, weights, time, r, mix_weight)

Evolve the initial state and return the final state vector.

Parameters:

Name Type Description Default
s str

Input bitstring representing the initial state.

 required
weights list of float

Coefficients for the problem Hamiltonian terms.

 required
time float

Total evolution time.

 required
r int

Number of Trotter steps for the approximate time evolution.

 required
mix_weight float

Coefficient for the mixer Hamiltonian.

 required

Returns:

Type Description
ndarray

The final state vector.
Notes

The total Hamiltonian simulated by the circuit is a weighted sum of the problem Hamiltonian terms and the mixer Hamiltonian:

In Ferguson et al. (2025) [arXiv:2506.19538], we use gammas to weight the entire problem Hamiltonian vs the mixer vs the constraint Hamiltonian, such that the total Hamiltonian is:

H = g_p * H_p + g_m * H_m + g_c * H_c

but here we allow for separate weights for each coupling tensor term, as well as a separate gamma for the mixer. The total Hamiltonian is then:

H = (w_b1*H_b1 + w_b2*H_b2 + ... + w_b2*H_bm) + g_m * H_m

In other words, the constraint hamiltonian is absorbed in the coupling list, and weighted by the corresponding gamma in the gammas list. This allows for more flexible weighting of different terms, and also allows us to use the same code for both constrained and unconstrained problems (by simply including or excluding the constraint Hamiltonian in the coupling list and adjusting the gammas accordingly).

Note that it is assumed that each term is already normalised appropriately, so the gammas can be interpreted as the relative weights of each term in the total Hamiltonian.

Source code in src/qemcmc/circuits.py 
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
def get_state_vector(self, s: str, weights: List[float], time: float, r: int, mix_weight: float) -> np.ndarray:
    """
    Evolve the initial state and return the final state vector.

    Parameters
    ----------
    s : str
        Input bitstring representing the initial state.
    weights : list of float
        Coefficients for the problem Hamiltonian terms.
    time : float
        Total evolution time.
    r : int
        Number of Trotter steps for the approximate time evolution.
    mix_weight : float
        Coefficient for the mixer Hamiltonian.

    Returns
    -------
    np.ndarray
        The final state vector.


    Notes
    -----
    The total Hamiltonian simulated by the circuit is a weighted sum of the problem Hamiltonian terms and the mixer Hamiltonian:

    In Ferguson et al. (2025) [arXiv:2506.19538], we use gammas to weight the entire problem Hamiltonian vs the mixer vs the constraint Hamiltonian, such that the total Hamiltonian is:

    ``H = g_p * H_p + g_m * H_m + g_c * H_c``

    but here we allow for separate weights for each coupling tensor term, as well as a separate gamma for the mixer. The total Hamiltonian is then:

    ``H = (w_b1*H_b1 + w_b2*H_b2 + ... + w_b2*H_bm) + g_m * H_m``

    In other words, the constraint hamiltonian is absorbed in the coupling list, and weighted by the corresponding gamma in the gammas list.
    This allows for more flexible weighting of different terms, and also allows us to use the same code for both constrained and unconstrained problems (by simply including or excluding the constraint Hamiltonian in the coupling list and adjusting the gammas accordingly).

    Note that it is assumed that each term is already normalised appropriately, so the gammas can be interpreted as the relative weights of each term in the total Hamiltonian.

    """


    num_wires = len(s)
    dev = self._get_device(num_wires)


    #if mix_weight < 0 or mix_weight > 1:
    #    raise ValueError(f"mix_weight must be between 0 and 1. Got {mix_weight}")

    #if np.any(np.array(weights) < 0):
    #    raise ValueError(f"Weights must be non-negative. Got {weights}")

    coeff_mixer = mix_weight
    coeff_problem = weights

    H_total = qml.Hamiltonian(
        [coeff_mixer] + list(np.ones(len(self.model.normalised_couplings))),
        [self.get_mixer_hamiltonian(num_wires)]
        + [self.get_problem_hamiltonian(couplings=[self.model.normalised_couplings[i]], sign=coeff_problem[i]) for i in range(len(self.model.normalised_couplings))],
    )

    @qml.qnode(dev)
    def quantum_evolution(input_string):
        for i, bit in enumerate(input_string):
            if bit == "1":
                qml.PauliX(i)
        qml.ApproxTimeEvolution(H_total, time, r)
        return qml.state()

    state_vector = quantum_evolution(s)
    return state_vector

get_sample 

get_sample(
    s_cg, time, r, mix_weight, local_couplings, weights=None
)

Generate a single sample by evolving the system and measuring.

Parameters:

Name Type Description Default
s_cg str

Input bitstring for the subgroup.

 required
time float

Total evolution time.

 required
r int

Number of Trotter steps for the approximate time evolution.

 required
mix_weight float

Coefficient for the mixer Hamiltonian.

 required
local_couplings list

Coupling tensors for the subgroup.

 required
weights list of float

Coefficients for the problem Hamiltonian terms. Defaults to ones.

 None

Returns:

Type Description
str

A single bitstring sample.
Source code in src/qemcmc/circuits.py 
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
def get_sample(self, s_cg: str, time: float, r: int, mix_weight: float, local_couplings: list, weights: List[float] = None) -> str:
    """
    Generate a single sample by evolving the system and measuring.

    Parameters
    ----------
    s_cg : str
        Input bitstring for the subgroup.
    time : float
        Total evolution time.
    r : int
        Number of Trotter steps for the approximate time evolution.
    mix_weight : float
        Coefficient for the mixer Hamiltonian.
    local_couplings : list
        Coupling tensors for the subgroup.
    weights : list of float, optional
        Coefficients for the problem Hamiltonian terms. Defaults to ones.

    Returns
    -------
    str
        A single bitstring sample.
    """
    if weights is None:
        weights = np.ones(len(local_couplings))

    num_wires = len(s_cg)
    dev = self._get_device(num_wires)

    #if mix_weight < 0 or mix_weight > 1:
    #    raise ValueError(f"mix_weight must be between 0 and 1. Got {mix_weight}")

    #if np.any(np.array(weights) < 0):
    #    raise ValueError(f"Weights must be non-negative. Got {weights}")

    coeff_mixer = mix_weight
    coeff_problem = weights

    # coeff_problem = [coeff for i, coeff in enumerate(coeff_problem) if local_couplings[i].ndim > 0 and np.any(local_couplings[i] != 0)]
    # local_couplings = [coupling for coupling in local_couplings if coupling.ndim > 0 and np.any(coupling != 0)]
    if len(local_couplings) == 0:
        # If there are no local couplings, just return the input state
        print("no non-zero local couplings, skipping evolution and returning input state")
        return s_cg

    H_total = qml.Hamiltonian(
        [coeff_mixer] + list(np.ones(len(local_couplings))),
        [self.get_mixer_hamiltonian(num_wires)] + [self.get_problem_hamiltonian(couplings=[local_couplings[i]], sign=coeff_problem[i]) for i in range(len(local_couplings))],
    )

    # set qnode to use our device with dynamically chosen wires
    @qml.qnode(dev, shots=1)
    def quantum_evolution(input_string):
        for i, bit in enumerate(input_string):
            if bit == "1":
                qml.PauliX(i)
        qml.ApproxTimeEvolution(H_total, time, r)
        return qml.sample()

    """print(f"Simulating quantum circuit with {num_wires} qubits, time={time}, num_trotter_steps={num_trotter_steps}")
    compiled_circuit = qml.compile(quantum_evolution)

    specs_dict = qml.specs(compiled_circuit)(s_cg)
    specs_dict = specs_dict['resources']
    # 3. Print the relevant metrics
    print(f"Circuit Depth: {specs_dict['depth']}")
    print(f"Total Gates:   {specs_dict['num_gates']}")

    # To see the count of two-qubit gates specifically:
    two_qubit_gates = specs_dict['gate_sizes'].get(2, 0)
    print(f"Two-qubit Gates: {two_qubit_gates}")

    # To see the breakdown by gate type (e.g., CNOT, RZ, etc.)
    print("\nGate Breakdown:")
    for gate, count in specs_dict['gate_types'].items():
        print(f"- {gate}: {count}")"""
    # Get the first shot from the sample
    compiled_circuit = qml.compile(quantum_evolution)
    sample = compiled_circuit(s_cg)[0]
    # sample = quantum_evolution(s_cg)[0]
    bitstring = "".join(str(int(b)) for b in sample)
    return bitstring

update 

update(s, subgroup_choice, local_couplings, gamma, time, r)

Update a bitstring by evolving a subgroup.

This performs a time evolution on a coarse-grained Hamiltonian to get s' from s.

Parameters:

Name Type Description Default
s str

The initial bitstring.

 required
subgroup_choice list of int

Indices of the subgroup to evolve.

 required
local_couplings list

Coupling tensors for the subgroup.

 required
gamma float

Mixing parameter.

 required
time float

Evolution time.

 required
r int

Number of Trotter steps for the approximate time evolution.

 required

Returns:

Type Description
str

The updated bitstring s'.
Source code in src/qemcmc/circuits.py 
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
def update(self, s: str, subgroup_choice: List[int], local_couplings: list, gamma: float, time: float, r: int) -> str:
    """
    Update a bitstring by evolving a subgroup.

    This performs a time evolution on a coarse-grained Hamiltonian to get s' from s.

    Parameters
    ----------
    s : str
        The initial bitstring.
    subgroup_choice : list of int
        Indices of the subgroup to evolve.
    local_couplings : list
        Coupling tensors for the subgroup.
    gamma : float
        Mixing parameter.
    time : float
        Evolution time.
    r : int
        Number of Trotter steps for the approximate time evolution.

    Returns
    -------
    str
        The updated bitstring s'.
    """

    self._validate_bitstring(s)
    # Get s_cg' for the subgroup and reconstruct full s' using s and s_cg'
    s_cg = "".join([s[i] for i in subgroup_choice])
    s_cg_prime = self.get_sample(s_cg, time, r, gamma, local_couplings)
    s_list = list(s)
    for i, global_index in enumerate(subgroup_choice):
        s_list[global_index] = s_cg_prime[i]
    return "".join(s_list)

check_Hamiltonian 

check_Hamiltonian(
    basis_states_ints=None, couplings=None, plot=False
)

Utility function to print the Energy of the Hamiltonian to be simulated. Naturally, for a classical problem, the energy of the Hamiltonian wrt to some basis state b H|b> = E|b> should equal the energy of the corresponding state in the classical model. This is a useful sanity check to ensure that the Hamiltonian is being constructed correctly from the coupling tensors.

Source code in src/qemcmc/circuits.py 
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
def check_Hamiltonian(self, basis_states_ints: List[str] = None, couplings: List[np.ndarray] = None, plot =  False) -> None:
    """
    Utility function to print the Energy of the Hamiltonian to be simulated. Naturally, for a classical problem, the energy of the Hamiltonian wrt to some basis state b H|b> = E|b> should equal the energy of the corresponding state in the classical model. This is a useful sanity check to ensure that the Hamiltonian is being constructed correctly from the coupling tensors.


    """

    # Get hamiltonian (not normalised)
    Hamiltonian = self.get_problem_hamiltonian(couplings)
    energies = []
    quantum_energies = []
    # Evaluate the Hamiltonian on some randomly chosen bitstring states, and compare to the classical energy of those states in the model
    for state in basis_states_ints:
        print("state int: ", state)
        state =bin(state)[2:].zfill(self.n_qubits)
        classical_energy = self.model.get_energy(state)

        dev = qml.device("lightning.qubit", wires=self.n_qubits)
        @qml.qnode(dev)
        def evaluate_energy(basis_state):
            # Prepares the state |1, 0> if basis_state=[1, 0]
            #print(basis_state)
            #print(len(bin(state)[2:].zfill(self.n_qubits)))
            #qml.BasisState(, wires=range(self.n_qubits))

            for i, bit in enumerate(basis_state):
                if bit == "1":
                    qml.PauliX(i)
            return qml.expval(Hamiltonian)

        quantum_energy = evaluate_energy(state)

        if type(self.model) is ConstraintModel:
            constraint_energy = self.model.get_constraint_energy(state)
            print(f"State: {state}, Classical Energy: {np.round(classical_energy,2)}, Quantum Hamiltonian Energy: {np.round(quantum_energy,2)}, Constraint energu {np.round(constraint_energy,2)},Constraint Satisfaction: {self.model.constraint_func(state)}")
        else:
            print(f"State: {state}, Classical Energy: {np.round(classical_energy,2)}, Quantum Hamiltonian Energy: {np.round(quantum_energy,2)}")
        energies.append(classical_energy)
        quantum_energies.append(quantum_energy)
    if plot:
        from matplotlib import pyplot as plt
        plt.scatter(energies, quantum_energies)
        plt.xlabel("Classical Energy")
        plt.ylabel("Quantum Energy")
        plt.show()
        plt.plot(np.arange(0,2**self.n_qubits), np.array(energies), label = "classical")
        plt.plot(np.arange(0,2**self.n_qubits), np.array(quantum_energies), label = "quantum")
        plt.show()


    return

MCMCState dataclass 

MCMCState(bitstring, accepted, energy=None, position=None)

Represents a single step in an MCMC trajectory.

Stores the proposed configuration, whether it was accepted by the Metropolis rule, its energy, and the position of the step in the chain.

MCMCChain dataclass 

MCMCChain(states=None, name='MCMC')

Container for the sequence of states produced during an MCMC run.

This class records all proposed states, tracks accepted configurations, and provides helper methods for extracting trajectories, energies, and empirical distributions from the Markov chain.

Source code in src/qemcmc/utils/helpers.py 
63
64
65
66
67
68
69
70
71
72
73
74
75
76
def __init__(self, states: Optional[List[MCMCState]] = None, name: Optional[str] = "MCMC"):
    self.name = name

    if states is None:
        self._states: List[MCMCState] = []
        self._current_state: MCMCState = None
        self._states_accepted: List[MCMCState] = []
        self.markov_chain: List[str] = []

    else:
        self._states = states
        self._current_state: MCMCState = next((s for s in self._states[::-1] if s.accepted), None)
        self._states_accepted: List[MCMCState] = [state for state in states if state.accepted]
        self.markov_chain: List[str] = self.get_list_markov_chain()

SpectralGap 

SpectralGap(proposal, model, temp=1.0)

Bases: Runner

Class that finds the spectral gap, and the acceptance and proposal matrices for a given mcmc.

Source code in src/qemcmc/spectralgap.py 
15
16
17
18
def __init__(self, proposal: Proposal, model: EnergyModel, temp: float = 1.0):
    self.proposal = proposal
    self.model = model
    self.temp = temp

find_acceptance_matrix 

find_acceptance_matrix()

Function to find the acceptance matrix for a given model instance.

Returns: A (np.ndarray): The acceptance matrix for the mcmc

Source code in src/qemcmc/spectralgap.py 
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
def find_acceptance_matrix(self):
    """
    Function to find the acceptance matrix for a given model instance.

    Returns:
        A (np.ndarray): The acceptance matrix for the mcmc

    """

    num_states = 2**self.proposal.n_spins

    A = np.zeros((num_states, num_states))

    energies = self.model.get_all_energies()

    for i in range(num_states):
        for j in range(num_states):
            if i != j:
                A[i][j] = self.get_acceptance_probability(energies[i], energies[j], temperature=self.temp)
            else:
                A[i][j] = 0

    return A

find_proposal_matrix_local 

find_proposal_matrix_local()

Function to find the proposal matrix for a given local chain.

Returns: Q (np.ndarray): The Q matrix for local proposal

Source code in src/qemcmc/spectralgap.py 
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
def find_proposal_matrix_local(self):
    """
    Function to find the proposal matrix for a given local chain.

    Returns:
        Q (np.ndarray): The Q matrix for local proposal
    """

    possible_states = self.model.S
    # TODO: define S separately, not inside the model since S is the state space and will be pretty inefficient otherwise

    Q = np.zeros((2**self.model.n_spins, 2**self.model.n_spins))

    # loop throguh and find the difference in bitstrings.
    # When the ith bitstring is different (by a valua of 1) from the jth bitstring add 1 to Q[i,j]
    for i in range(2**self.model.n_spins):
        for j in range(2**self.model.n_spins):
            # TODO: what's sm? maybe rename to something like diff_count or something descriptive
            sm = 0
            for k in range(self.model.n_spins):
                sm += abs(int(possible_states[i][k]) - int(possible_states[j][k]))

            # ie if the number of different strings is the size of the cluster (= 1 for local)
            if sm == 1:
                Q[i, j] = 1

    row_sums = Q.sum(axis=1)
    Q = Q / row_sums[:, np.newaxis]

    return Q

find_proposal_matrix_uniform 

find_proposal_matrix_uniform()

Function to find the proposal matrix for a given uniform chain.

Returns: Q (np.ndarray): The Q matrix for uniform proposal

Source code in src/qemcmc/spectralgap.py 
75
76
77
78
79
80
81
82
83
84
85
86
87
def find_proposal_matrix_uniform(self):
    """
    Function to find the proposal matrix for a given uniform chain.

    Returns:
        Q (np.ndarray): The Q matrix for uniform proposal
    """

    Q = np.ones((2**self.model.n_spins, 2**self.model.n_spins)) / (self.proposal.n_spins**2 - 1)
    row_sums = Q.sum(axis=1)
    Q = Q / row_sums[:, np.newaxis]

    return Q

find_proposal_matrix_quantum 

find_proposal_matrix_quantum()

Function to find the proposal matrix for a given QeMCMCChain object.

Returns: Q (np.ndarray): The Q matrix for quantum proposal

Source code in src/qemcmc/spectralgap.py 
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
def find_proposal_matrix_quantum(self):
    """
    Function to find the proposal matrix for a given QeMCMCChain object.

    Returns:
        Q (np.ndarray): The Q matrix for quantum proposal
    """
    if not isinstance(self.proposal.gamma, (int, float)):
        raise ValueError("gamma must be a float to find the proposal matrix for quantum proposal, not a tuple. Got: ", self.proposal.gamma, "of type: ", type(self.proposal.gamma))
    if not isinstance(self.proposal.time, (int, float)):
        raise ValueError("time must be a number to find the proposal matrix for quantum proposal, not a tuple. Got: ", self.proposal.time, "of type: ", type(self.proposal.time))

    # Determine r to compute spectral gap
    if self.proposal.r is not None:
        r = self.proposal.r
    else:
        dt = self.proposal.delta_t if self.proposal.delta_t is not None else DEFAULT_DELTA_T
        r = max(1, int(np.floor(self.proposal.time / dt)))

    Q = np.zeros((2**self.model.n_spins, 2**self.model.n_spins))

    for i in range(2**self.model.n_spins):
        Q[i, :] += abs(self.proposal.CM.get_state_vector(self.model.S[i], self.proposal.coupling_weights, self.proposal.time, r, self.proposal.gamma)) ** 2
        # get_output_statevector(self.proposal.model.S[i])
    Q = Q

    return Q

find_proposal_matrix 

find_proposal_matrix()

Function to find the proposal matrix for a given mcmc. This is not done by brute force

Source code in src/qemcmc/spectralgap.py 
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
def find_proposal_matrix(self):
    """
    Function to find the proposal matrix for a given mcmc.
    This is not done by brute force
    """

    if self.proposal.method == "local":
        Q = self.find_proposal_matrix_local()
    elif self.proposal.method == "uniform":
        Q = self.find_proposal_matrix_uniform()
    elif self.proposal.method == "quantum":
        Q = self.find_proposal_matrix_quantum()
    else:
        raise ValueError("Method not recognised. Only 'local', 'uniform' or 'quantum' proposal methods are implimented in find_proposal_method.")

    return Q

find_spectral_gap 

find_spectral_gap(A=None, Q=None)

Function to find the spectral gap of a given mcmc.

Parameters:

Name Type Description Default
A

The acceptance matrix for the mcmc (optional, if not given, will be calculated)

 None
Q

The proposal matrix for the mcmc (optional, if not given, will be calculated)

 None
Source code in src/qemcmc/spectralgap.py 
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
def find_spectral_gap(self, A=None, Q=None):
    """
    Function to find the spectral gap of a given mcmc.

    Parameters
    ----------
    A (np.ndarray):
        The acceptance matrix for the mcmc (optional, if not given, will be calculated)
    Q (np.ndarray):
        The proposal matrix for the mcmc (optional, if not given, will be calculated)
    """

    if A is None:
        A = self.find_acceptance_matrix()
    if Q is None:
        Q = self.find_proposal_matrix()
        # Q = self.find_proposal_matrix_brute_force(multiple = 10*2**self.mcmc.n_spins)

    P = np.multiply(Q, A)

    # account for rejected swaps to add to s = s' matrix element
    for i in range(P.shape[0]):
        s = np.sum(P[i, :]) - P[i, i]
        P[i, i] = 1 - s

    # find eigenvalues
    e_vals, e_vecs = sp.linalg.eig(P)
    e_vals = np.flip(np.sort(abs(e_vals)))

    # find spectral gap
    delta = e_vals[1]
    delta = 1 - delta

    return delta

CoarseGraining 

CoarseGraining(
    n, subgroups=None, subgroup_probs=None, repeated=True
)

CoarseGraining class to generate partitions of spins for quantum proposals.

Parameters:

Name Type Description Default
n int

Number of spins in the system.

 required
subgroups list[list[int]]

A list of subgroups, where each subgroup is a list of spin indices. If None, the entire set of spins is treated as one subgroup.

 None
subgroup_probs list[float]

A list of probabilities corresponding to each subgroup, used for weighted random selection. Must sum to 1. If None, subgroups are selected uniformly at random.

 None
repeated bool

If True, then multiple subgroups are run on the quantum computer in serial, if not then only one subgroup is selected at random and run on the quantum computer. Default is True.

 True
Source code in src/qemcmc/coarse_grain.py 
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
def __init__(self, n, subgroups=None, subgroup_probs=None, repeated=True):

    self._user_specified = subgroups is not None

    if subgroups is None:
        subgroups = [list(range(n))]
        subgroup_probs = [1.0]
    else:
        if subgroup_probs is None:
            subgroup_probs = [1.0 / len(subgroups)] * len(subgroups)

    validate_subgroups(subgroups=subgroups, subgroup_probs=subgroup_probs, n_spins=n)

    self.n = n
    self.subgroups = subgroups
    self.subgroup_probs = subgroup_probs
    self.repeated = repeated

    if self.repeated and self._user_specified:
        if not np.all([len(s) == len(subgroups[0]) for s in subgroups]):
            raise NotImplementedError("If providing own subgroups, repeated must be set to False as it is not obvious how one performs repeated coarse-graining without knowledge of the subgroups. If length, l, of each subrgoup is equal, we continue assuming l = n//m.")
        else:
            print("Provided subgroups are all of same length, so repeated coarse-graining will be performed accross a random choice of subgroups. This means that the same variables may be chosen multiple times.")

sample 

sample()

Randomly samples a subgroup according to the specified probability distribution.

Source code in src/qemcmc/coarse_grain.py 
45
46
47
48
49
50
def sample(self) -> list[int]:
    """
    Randomly samples a subgroup according to the specified probability distribution.
    """
    idx = np.random.choice(len(self.subgroups), p=self.subgroup_probs)
    return self.subgroups[idx]

get_partitions 

get_partitions(m)

Returns partitions of spins for sequential quantum updates.

If the user provided explicit subgroups at initialisation, those are sampled accordingly. returned directly (all if repeated=True, otherwise just the first). If no subgroups were specified, random disjoint partitions of approximate size n/m are generated.

Parameters:

Name Type Description Default
m int

The number of partitions to generate. Ignored if user-specified subgroups are provided.

 required
Source code in src/qemcmc/coarse_grain.py 
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
def get_partitions(self, m: int) -> list[list[int]]:
    """
    Returns partitions of spins for sequential quantum updates.

    If the user provided explicit subgroups at initialisation, those are sampled accordingly.
    returned directly (all if ``repeated=True``, otherwise just the first).
    If no subgroups were specified, random disjoint partitions of
    approximate size n/m are generated.

    Parameters
    ----------
    m : int
        The number of partitions to generate. Ignored if user-specified subgroups are provided.
    """

    if self._user_specified:

        if self.repeated:
            # choose m subgroups from self.subgroups
            if len(self.subgroups) < m:
                raise ValueError("Number of partitions, m, must be less than or equal to the number of provided subgroups.")
            choices = np.random.choice(np.arange(0,len(self.subgroups)), size = m, p = self.subgroup_probs, replace=False)

            chosen_subgroups = []
            for choice in choices:
                chosen_subgroups.append(self.subgroups[choice])
            return chosen_subgroups
        else:
            chosen_subgroup = np.random.choice(self.subgroups, p=self.subgroup_probs)

            return [self.subgroups[chosen_subgroup]]

    spins = np.arange(self.n)
    np.random.shuffle(spins)
    # array_split handles uneven divisions gracefully
    chunks = [list(chunk) for chunk in np.array_split(spins, m)]
    if self.repeated:
        return chunks
    else:
        return [chunks[0]]

get_random_state 

get_random_state(num_spins)

Generate a random state for a given number of spins.

Parameters: num_spins (int): The number of spins in the system. Returns: str: A bitstring representing the random state.

Source code in src/qemcmc/utils/helpers.py 
188
189
190
191
192
193
194
195
196
197
198
199
200
201
def get_random_state(num_spins: int) -> str:
    """
    Generate a random state for a given number of spins.

    Parameters:
        num_spins (int): The number of spins in the system.
    Returns:
        str: A bitstring representing the random state.
    """

    next_state = np.random.randint(0, 2, size=num_spins, dtype=np.int8)
    # s_prime = f"{next_state:0{num_spins}b}"
    s_prime = "".join(str(bit) for bit in next_state)
    return s_prime

get_all_possible_states 

get_all_possible_states(num_spins)

Returns all possible binary strings of length n=num_spins

Paremeters: num_spins: n length of the bitstring Returns: list: A list of all possible binary strings of length num_spins.

Source code in src/qemcmc/utils/helpers.py 
204
205
206
207
208
209
210
211
212
213
214
215
216
def get_all_possible_states(num_spins: int) -> list:
    """
    Returns all possible binary strings of length n=num_spins

    Paremeters:
        num_spins: n length of the bitstring
    Returns:
        list: A list of all possible binary strings of length num_spins.
    """

    num_possible_states = 2 ** (num_spins)
    possible_states = [f"{k:0{num_spins}b}" for k in range(0, num_possible_states)]
    return possible_states