Post Processing Utilities

oqd_heisenberg_ion.common.postprocess

ffs

free_fermion_energy

computes the ground state energy per site associated with a chain of free fermions

N should not be divisible by 4 for this calculation to yield the correct energy

Parameters:

• N (int) –

  number of sites

• phi (float) –

  boundary twist angle

Returns:

• float –

  energy per site

Source code in src/oqd_heisenberg_ion/common/postprocess/ffs.py

def free_fermion_energy(N, phi):
    """
    computes the ground state energy per site associated with a chain of free fermions

    N should not be divisible by 4 for this calculation to yield the correct energy

    Args:
        N (int): number of sites
        phi (float): boundary twist angle

    Returns:
        (float): energy per site
    """
    s_min = -int(N / 2)
    E = 0.0
    for i in range(N):
        s = s_min + i
        eval = -np.cos(2 * np.pi * s / N - phi / N)
        if eval < 0.0:
            E += eval

    # if (s != s_max):
    #    raise Exception("s not equal to s_max")

    return E / N

free_fermion_stiffness

computes the spin stiffness associated with a chain of free fermions

Parameters:

• N (int) –

  number of sites

• phi (float) –

  boundary twist angle

Returns:

• float –

  spin stiffness

Source code in src/oqd_heisenberg_ion/common/postprocess/ffs.py

def free_fermion_stiffness(N, phi):
    """
    computes the spin stiffness associated with a chain of free fermions

    Args:
        N (int): number of sites
        phi (float): boundary twist angle

    Returns:
        (float): spin stiffness
    """

    E_phi = free_fermion_energy(N, phi)
    E_zero = free_fermion_energy(N, 0.0)
    E_minus_phi = free_fermion_energy(N, -phi)

    rho = (N**2) * ((E_phi + E_minus_phi - 2.0 * E_zero) / (phi**2))

    return rho

swt

gamma_k

computes the coefficient gamma corresponding to the kth Fourier mode in linear spin wave theory for the power law decay spin-1/2 XY model For details, see: https://arxiv.org/pdf/2601.20058

Parameters:

• N (int) –

  number of sites

• k (int) –

  Fourier mode

• theta (float) –

  twist angle

• alpha (float) –

  power law decay exponent in interaction couplings

Returns:

• float –

  gamma_k

Source code in src/oqd_heisenberg_ion/common/postprocess/swt.py

def gamma_k(N, k, theta, alpha):
    """
    computes the coefficient gamma corresponding to the kth Fourier mode in linear spin wave theory for the power law decay spin-1/2 XY model
    For details, see: https://arxiv.org/pdf/2601.20058

    Args:
        N (int): number of sites
        k (int): Fourier mode
        theta (float): twist angle
        alpha (float): power law decay exponent in interaction couplings

    Returns:
        (float): gamma_k
    """

    gamma = 0.0
    num_terms = int((N - 1) / 2)
    for r in range(1, num_terms + 1):
        gamma += (1.0 / (r**alpha)) * np.cos(2.0 * np.pi * r * k / N) * np.cos(r * theta)

    return gamma

E_0_LSW

computes the ground state energy in linear spin wave theory for the power law decay spin-1/2 XY model. For details, see: https://arxiv.org/pdf/2601.20058

Parameters:

• N (int) –

  number of sites

• alpha (float) –

  power law decay exponent

• J (float) –

  Energy scale

• theta (float, default: 0.0 ) –

  Boundary twist angle. Defaults to 0.0.

Returns:

• float –

  ground state energy in linear spin wave theory

Source code in src/oqd_heisenberg_ion/common/postprocess/swt.py

def E_0_LSW(N, alpha, J, theta=0.0):
    """
    computes the ground state energy in linear spin wave theory for the power law decay spin-1/2 XY model.
    For details, see: https://arxiv.org/pdf/2601.20058

    Args:
        N (int): number of sites
        alpha (float): power law decay exponent
        J (float): Energy scale
        theta (float, optional): Boundary twist angle. Defaults to 0.0.

    Returns:
        (float): ground state energy in linear spin wave theory
    """

    gamma_0 = gamma_k(N, 0, theta, alpha)
    energy = (-3.0 / 2.0) * N

    for j in range(N):
        gamma_j = gamma_k(N, j, theta, alpha)
        energy += np.sqrt(1.0 - gamma_j / gamma_0) + gamma_j / (2.0 * gamma_0)

    energy *= J * gamma_0 / 2.0

    return energy

E_0_MF

computes the ground state energy in mean field theory for the power law decay spin-1/2 XY model. For details, see: https://arxiv.org/pdf/2601.20058

Parameters:

• N (int) –

  number of sites

• alpha (float) –

  power law decay exponent

• J (float) –

  Energy scale

• theta (float, default: 0.0 ) –

  Boundary twist angle. Defaults to 0.0.

Returns:

• float –

  ground state energy in mean field theory

Source code in src/oqd_heisenberg_ion/common/postprocess/swt.py

def E_0_MF(N, alpha, J, theta=0.0):
    """
    computes the ground state energy in mean field theory for the power law decay spin-1/2 XY model.
    For details, see: https://arxiv.org/pdf/2601.20058

    Args:
        N (int): number of sites
        alpha (float): power law decay exponent
        J (float): Energy scale
        theta (float, optional): Boundary twist angle. Defaults to 0.0.

    Returns:
        (float): ground state energy in mean field theory
    """

    gamma_0 = gamma_k(N, 0, theta, alpha)

    energy = -J * N * gamma_0 / 4.0

    return energy

E_0_LSW_NN

computes the ground state energy in linear spin wave theory for the nearest neighbor spin-1/2 XY model. For details, see: https://arxiv.org/pdf/2601.20058

Parameters:

• N (int) –

  number of sites

• J (float) –

  Energy scale

Returns:

• float –

  ground state energy in linear spin wave theory

Source code in src/oqd_heisenberg_ion/common/postprocess/swt.py

def E_0_LSW_NN(N, J):
    """
    computes the ground state energy in linear spin wave theory for the nearest neighbor spin-1/2 XY model.
    For details, see: https://arxiv.org/pdf/2601.20058

    Args:
        N (int): number of sites
        J (float): Energy scale

    Returns:
        (float): ground state energy in linear spin wave theory
    """

    energy = (-3.0 / 2.0) * N

    for j in range(N):
        gamma_j = np.cos(2.0 * np.pi * j / N)
        energy += np.sqrt(1.0 - gamma_j) + gamma_j / (2.0)

    energy *= J / 2.0

    return energy

rho_2

computes the spin stiffness for the power law decay spin-1/2 XY model in linear spin wave theory For details, see: https://arxiv.org/pdf/2601.20058

Parameters:

• N (int) –

  number of sites

• alpha (float) –

  power law decay exponent

• theta (float) –

  twist angle

• J_1 (float) –

  energy scale

• J_2 (float, default: 1.0 ) –

  energy scale. Defaults to 1.0. Should typically be set equal to J_1

Returns:

• float –

  spin wave theory spin stiffness

Source code in src/oqd_heisenberg_ion/common/postprocess/swt.py

def rho_2(N, alpha, theta, J_1, J_2=1.0):
    """
    computes the spin stiffness for the power law decay spin-1/2 XY model in linear spin wave theory
    For details, see: https://arxiv.org/pdf/2601.20058

    Args:
        N (int): number of sites
        alpha (float): power law decay exponent
        theta (float): twist angle
        J_1 (float): energy scale
        J_2 (float, optional): energy scale. Defaults to 1.0. Should typically be set equal to J_1

    Returns:
        (float): spin wave theory spin stiffness
    """

    energy_theta = E_0_LSW(N, alpha, J_1, theta)
    energy_0 = E_0_LSW(N, alpha, J_2, 0.0)
    rho_N_alpha = 2.0 * (energy_theta - energy_0) / (N * (theta**2))

    return rho_N_alpha

rho_mf

computes the spin stiffness for the power law decay spin-1/2 XY model in mean field theory For details, see: https://arxiv.org/pdf/2601.20058

Parameters:

• N (int) –

  number of sites

• alpha (float) –

  power law decay exponent

• theta (float) –

  twist angle

• J (float) –

  energy scale

Returns:

• float –

  mean field theory spin stiffness

Source code in src/oqd_heisenberg_ion/common/postprocess/swt.py

def rho_mf(N, alpha, theta, J):
    """
    computes the spin stiffness for the power law decay spin-1/2 XY model in mean field theory
    For details, see: https://arxiv.org/pdf/2601.20058

    Args:
        N (int): number of sites
        alpha (float): power law decay exponent
        theta (float): twist angle
        J (float): energy scale

    Returns:
        (float): mean field theory spin stiffness
    """

    energy_theta = E_0_MF(N, alpha, J, theta)
    energy_0 = E_0_MF(N, alpha, J, 0.0)
    rho_N_alpha = 2.0 * (energy_theta - energy_0) / (N * (theta**2))

    return rho_N_alpha

utils

combine_different_runs

concatenates two data sets

Parameters:

• data_1 (ndarray[float]) –

  T1 x 1 array containing the first data set

• data_2 (ndarray[float]) –

  T2 x 1 array containing the second data set

• drop_samples_1 (int) –

  number of samples to drop from the top in data1

• drop_samples_2 (int) –

  number of samples to drop from the top in data2

Returns:

• ndarray[float] –

  (T3 x 1) array containing the concatenated data set

Source code in src/oqd_heisenberg_ion/common/postprocess/utils.py

def combine_different_runs(data_1, data_2, drop_samples_1, drop_samples_2):
    """
    concatenates two data sets

    Args:
        data_1 (numpy.ndarray[float]): T1 x 1 array containing the first data set
        data_2 (numpy.ndarray[float]): T2 x 1 array containing the second data set
        drop_samples_1 (int): number of samples to drop from the top in data1
        drop_samples_2 (int): number of samples to drop from the top in data2

    Returns:
        (numpy.ndarray[float]): (T3 x 1) array containing the concatenated data set
    """

    data_3 = []
    data_1 = data_1[drop_samples_1:]
    data_2 = data_2[drop_samples_2:]

    for i in range(len(data_1)):
        data_3.append(data_1[i])

    for i in range(len(data_2)):
        data_3.append(data_2[i])

    return data_3

statistics_binning

Uses binning to post process QMC estimator data

Parameters:

• arr (ndarray[float]) –

  T x 1 array of QMC data where T is the number of simulation steps

• auto_corr_drop (int) –

  number of points to remove for auto-correlation

• eq_drop (int) –

  number of points to remove for equilibration

Returns:

• tuple[float, float] –

  mean and standard error of QMC data

Source code in src/oqd_heisenberg_ion/common/postprocess/utils.py

def statistics_binning(arr, auto_corr_drop, eq_drop):
    """
    Uses binning to post process QMC estimator data

    Args:
        arr (numpy.ndarray[float]): T x 1 array of QMC data where T is the number of simulation steps
        auto_corr_drop (int): number of points to remove for auto-correlation
        eq_drop (int): number of points to remove for equilibration

    Returns:
        (tuple[float,float]): mean and standard error of QMC data
    """

    arr2 = arr[eq_drop:]
    arr3 = arr2[0::auto_corr_drop]
    workingNdim = int(math.log(len(arr3)) / math.log(2))
    mean = np.mean(arr3)
    standardError = max_error_binning(arr3, workingNdim - 6)
    return mean, standardError

error_propagation

computes statistical standard error of array

Parameters:

• data (ndarray[float]) –

  T x 1 array containing the data

Returns:

• float –

  standard error

Source code in src/oqd_heisenberg_ion/common/postprocess/utils.py

def error_propagation(data):
    """
    computes statistical standard error of array

    Args:
        data (numpy.ndarray[float]): T x 1 array containing the data

    Returns:
        (float): standard error
    """

    ndim = len(data)
    error = np.std(data, ddof=0) / np.sqrt(ndim)

    return error

max_error_binning

helper function containing the binning error logic

Parameters:

• data (ndarray[float]) –

  T x 1 array containing the data

• dim (int) –

  number of binning levels

Raises:

• Exception –

  if dim is less than 1 because of insufficient data

Returns:

• float –

  binning error

Source code in src/oqd_heisenberg_ion/common/postprocess/utils.py

def max_error_binning(data, dim):
    """
    helper function containing the binning error logic

    Args:
        data (numpy.ndarray[float]): T x 1 array containing the data
        dim (int): number of binning levels

    Raises:
        Exception: if dim is less than 1 because of insufficient data

    Returns:
        (float): binning error
    """

    if dim <= 1:
        raise Exception("Not enough points MC steps were used for the binning method\n")

    error = np.zeros(dim)
    error[0] = error_propagation(data)

    for i in range(1, dim):
        bin_dim = int(len(data) / 2)
        data_bin = np.zeros(bin_dim)

        for j in range(bin_dim):
            data_bin[j] = 0.5 * (data[2 * j] + data[2 * j + 1])
        data = data_bin
        error[i] = error_propagation(data)

    return np.max(error)

ed_energy

computes the equilibrium energy from exact diagonalization calculations at given temperature

Parameters:

• evals (ndarray[float]) –

  T x 1 array containing the energies for each state

• T (float) –

  temperature

Returns:

• float –

  equilibrium energy

Source code in src/oqd_heisenberg_ion/common/postprocess/utils.py

def ed_energy(evals, T):
    """
    computes the equilibrium energy from exact diagonalization calculations at given temperature

    Args:
        evals (numpy.ndarray[float]): T x 1 array containing the energies for each state
        T (float): temperature

    Returns:
        (float): equilibrium energy
    """

    beta = 1.0 / T

    energy = 0.0
    partition = 0.0

    for i in range(len(evals)):
        energy += np.exp(-beta * evals[i]) * evals[i]
        partition += np.exp(-beta * evals[i])

    energy /= partition

    return energy

compute_histogram_from_shot_data

computes the histogram associated with shot data

Parameters:

• N (int) –

  number of sites in the lattice

• num_shots (int) –

  number of data points

• shot_data (ndarray[int]) –

  T x N array of shot data. shot_data[t,k] represents the spin configuration at time t and site k

Returns:

• ndarray[float] –

  (T x 1) array containing the histogram frequences where T is the number of possible bitstrings given N

Source code in src/oqd_heisenberg_ion/common/postprocess/utils.py

def compute_histogram_from_shot_data(N, num_shots, shot_data):
    """
    computes the histogram associated with shot data

    Args:
        N (int): number of sites in the lattice
        num_shots (int): number of data points
        shot_data (numpy.ndarray[int]): T x N array of shot data. shot_data[t,k] represents the spin configuration at time t and site k

    Returns:
        (numpy.ndarray[float]): (T x 1) array containing the histogram frequences where T is the number of possible bitstrings given N
    """

    num_bits = 2**N
    freqs = np.zeros(num_bits)
    for i in range(num_shots):
        bit_str = 0
        for j in range(N):
            signed_spin_config = shot_data[i, j]
            bit_i_j = int(0.5 * (signed_spin_config + np.abs(signed_spin_config)))
            bit_str += (2 ** (j)) * (bit_i_j)
        freqs[bit_str] += 1.0 / num_shots

    return freqs

kl_divergence

computes the KL divergence between two histograms

Parameters:

• hist1 (ndarray[float]) –

  T x 1 array containing the histogram frequences of the first distribution

• hist2 (ndarray[float]) –

  T x 1 array containing the histogram frequences of the second distribution

• num_bits (int) –

  Number of unique measurements in histograms

Returns:

• float –

  KL divergence between data sets

Source code in src/oqd_heisenberg_ion/common/postprocess/utils.py

def kl_divergence(hist1, hist2, num_bits):
    """
    computes the KL divergence between two histograms

    Args:
        hist1 (numpy.ndarray[float]): T x 1 array containing the histogram frequences of the first distribution
        hist2 (numpy.ndarray[float]): T x 1 array containing the histogram frequences of the second distribution
        num_bits (int): Number of unique measurements in histograms

    Returns:
        (float): KL divergence between data sets
    """

    kl_divergence = 0.0
    for j in range(num_bits):
        if hist1[j] != 0 and hist2[j] != 0:
            kl_divergence += hist1[j] * (np.log(hist1[j]) - np.log(hist2[j]))

    return kl_divergence