Skip to content
Pasqal Documentation

Results are limited to the current section : emulationtools

Getting started with emu-mps

# standard library
from typing import Any
import math


# Pulser and others libraries
import numpy as np
import matplotlib.pyplot as plt
import pulser
from pulser.devices import AnalogDevice


# emu-mps
from emu_mps import (
    MPS,
    MPSConfig,
    MPSBackend,
    BitStrings,
    Fidelity,
    Occupation,
)

Pulser and emu-mps

Section titled “Pulser and emu-mps”

Emu-mps is a Pulser backend, designed to emulate the dynamics of programmable arrays of neutral atoms and here are the steps on how to use it:

  1. using Pulser we will first create the atomic Register
  2. again with Pulser we will then generate the Pulse sequence that produces the AFM state for the created register
  3. we create the configuration object (MPSConfig) and we fill it with the needed observables and time step dt
  4. we instantiate the backend class (MPSBackend) and run the simulation
  5. we inspect the results obtained

Register and Sequence creation

Section titled “Register and Sequence creation”
def afm_sequence_from_register(
    reg: pulser.Register,
    Omega_max: float,
    delta_0: float,
    delta_f: float,
    t_rise: int,
    t_fall: int,
    factor_sweep: int,
    device: Any = pulser.devices.MockDevice,
) -> pulser.Sequence:
    """Sequence that creates AntiFerromagnetic State (AFM) for 1d chain of atoms using pulser.
    This function constructs a sequence of pulses to transition a system of qubits
    distributed in a 1D chain (represented by `reg`) into an AFM state using a specified device.
    The sequence consists of three phases: a rise, a sweep, and a fall.
    For more information, check this [arxiv article](https://arxiv.org/abs/1711.01185)."""


    t_sweep = (delta_f - delta_0) / (2 * np.pi * 10) * 1000 * factor_sweep
    rise = pulser.Pulse.ConstantDetuning(
        pulser.waveforms.RampWaveform(t_rise, 0.0, Omega_max), delta_0, 0.0
    )
    sweep = pulser.Pulse.ConstantAmplitude(
        Omega_max, pulser.waveforms.RampWaveform(int(t_sweep), delta_0, delta_f), 0.0
    )
    fall = pulser.Pulse.ConstantDetuning(
        pulser.waveforms.RampWaveform(t_fall, Omega_max, 0.0), delta_f, 0.0
    )


    seq = pulser.Sequence(reg, device)
    seq.declare_channel("ising_global", "rydberg_global")
    seq.add(rise, "ising_global")
    seq.add(sweep, "ising_global")
    seq.add(fall, "ising_global")


    return seq




def square_perimeter_points(L: int) -> np.ndarray:
    """
    Calculate the coordinates of the points located on the perimeter of a square of size L.
    The square is centered at the origin (0, 0) with sides parallel to the axes.
    The points are ordered starting from the bottom-left corner and moving
    counter-clockwise around the square. The order is important when measuaring the bitstrings


    Args:
        L (int): The length of the side of the square. L should be a positive integer.


    Returns:
        np.ndarray: An array of shape (4*L-4, 2) containing the coordinates of the perimeter points.


    Example:
        >>> square_perimeter_points(3)
        array([[-1, -1],
               [-1,  0],
               [-1,  1],
               [ 0,  1],
               [ 1,  1],
               [ 1,  0],
               [ 1, -1],
               [ 0, -1]])
    """
    pairOrodd = L % 2
    toGrid = int(math.floor(L / 2))
    if pairOrodd == 0:
        axis = list(range(-toGrid, toGrid, 1))
    else:
        axis = list(range(-toGrid, toGrid + 1, 1))
    coord = []
    for i in axis:  # from left, first column of the perimeter
        coord.append([axis[0], i])


    for i in axis[1:-1]:
        coord.append([i, axis[-1]])


    for i in reversed(axis):
        coord.append([axis[-1], i])


    for i in reversed(axis[1:-1]):
        coord.append([i, axis[0]])


    return np.array(coord)
Omega_max = 2 * 2 * np.pi
delta_0 = -6 * Omega_max / 2
delta_f = 1 * Omega_max / 2
t_rise = 500
t_fall = 1500
sweep_factor = 2


square_length = 3
R_interatomic = AnalogDevice.rydberg_blockade_radius(Omega_max / 2)
coords = R_interatomic * square_perimeter_points(square_length)
reg = pulser.Register.from_coordinates(coords,prefix='q')
reg.draw(blockade_radius=R_interatomic, draw_graph=True, draw_half_radius=True)
seq = afm_sequence_from_register(
    reg, Omega_max, delta_0, delta_f, t_rise, t_fall, sweep_factor, AnalogDevice
)
seq.draw("input")

Using emu-mps as backend

Section titled “Using emu-mps as backend”

As mentioned earlier, to run a simulation with the emu-mps backend we need to provide as input a Pulse sequence - which we have just created - and a configuration object.

We still have to create the configuration for the emu-mps backend. This is done via an instantiation of the configuration class MPSConfig which contains all the observables that we wish to measure and the time step chosen for the simulation, along with various algorithm specific parameters that are explained in the documentation.

We start by setting a bigger discretization time than the default one provided ($dt=10$). For simplicity, we will measure the observables only at the final time of the simulation.

We also fix the basis along which the measurements will be done. For the details regarding the conventions used we refer to the Pulser documentation (external).

dt = 5.0
eval_times = [0.0, 1.0]


basis = ("r","g")

BitString counts

Section titled “BitString counts”

It samples the evolved state at desired time steps, returning the bitStrings in a counter.

sampling_times = 1000
bitstrings = BitStrings(evaluation_times=eval_times, num_shots=sampling_times)

Fidelity

Section titled “Fidelity”

The fidelity is computed as $\langle \psi_{evolved} | \phi_{given} \rangle^2$ where

  • $\psi_{evolved}$ is the system state at desired steps
  • $\phi_{given}$ is the state onto which we want to project the system state.

In this tutorial we will compute the fidelity against one of the two dominant antiferromagnetic states $\phi_{given} = |rgrgrgrg>$

nqubits = len(seq.register.qubit_ids)


afm_string_pure = {"rgrgrgrg": 1.0} # |10101010> in ground-rydberg basis


afm_mps_state = MPS.from_state_amplitudes(
    eigenstates=basis, amplitudes=afm_string_pure
)
fidelity_mps_pure = Fidelity(evaluation_times=eval_times, state=afm_mps_state,tag_suffix="1")

Occupation of the excited states

Section titled “Occupation of the excited states”

It is computed as $\langle \psi_{evolved} |\frac{(1+Z_i)}{2}|\psi_{evolved}\rangle$ and often informally referred to as the magnetization at each atom site.

duration = seq.get_duration()
step = int(dt)
occ_eval_times = [t / duration for t in range(0, duration+step, step)]
density = Occupation(evaluation_times=occ_eval_times)

Configuration object for emu-mps

Section titled “Configuration object for emu-mps”
mpsconfig = MPSConfig(
    dt=dt,
    precision=1.0e-9,
    observables=[
        bitstrings,
        fidelity_mps_pure,
        density,
    ],
    log_level=10,
)

Finally, we run

Section titled “Finally, we run”

We instantiate the backend and run the simulation. The MPSBackend class takes the sequence and config objects as arguments.

sim = MPSBackend(seq, config=mpsconfig)
results = sim.run()

Inspecting the result object

Section titled “Inspecting the result object”

In the following lines, we are going to give brief code examples of how you can get the information from the results object

results.get_result_tags()

Computation statistics

Section titled “Computation statistics”

During the simulation the sequence is discretized into steps according to dt. For each step results.statistics contains:

  • $\chi$: maximum bond dimension used.
  • $|\Psi|$: MPS memory footprint.
  • RSS: peak memory allocation.
  • Δt: wall-clock time for that step.

Example:

# print statistics of step 1
statistics = results.statistics
print(f"Number of discretized steps:\n {len(statistics)}")
print("Step 1 statistics\n",statistics[1])

Bitstrings analysis

Section titled “Bitstrings analysis”

Below, we retrieve the bitstrings computed. We observe that they were indeed computed only at a single time $ns = 3400$, and we find those that were sampled the highest number of times with their relative counting.

results.get_result_times(bitstrings)
bitstrings_final = results.get_result(bitstrings, 1.0)


max_val = max(bitstrings_final.values())  # max number of counts in the bitstring
max_string = [key for key, value in bitstrings_final.items() if value == max_val]
print(
    "The most frequent bitstring is {} which was sampled {} times".format(
        max_string, max_val
    )
)


filtered_counts = [count for count in bitstrings_final.values() if count > 20]
filtered_bitstrings = [
    bitstring for bitstring, count in bitstrings_final.items() if count > 20
]
x_labels = range(len(filtered_bitstrings))
with plt.style.context("seaborn-v0_8-darkgrid"):
    fig, ax = plt.subplots()
    ax.bar(x_labels, filtered_counts, color="teal", alpha=0.8)
    ax.set_xlabel("Bitstrings")
    ax.set_ylabel("Counts")
    ax.set_title("Histogram of Bitstring Counts (Counts > 20)")
    ax.set_xticks(x_labels)
    ax.set_xticklabels(filtered_bitstrings, rotation="vertical")
    ax.spines["top"].set_visible(False)
    ax.spines["right"].set_visible(False)
    plt.tight_layout()
    plt.show()

Fidelity analysis

Section titled “Fidelity analysis”

Here we compute the fidelity of the system against the two different AFM state realizations defined above.

fidelity_pure = results.get_result(fidelity_mps_pure, 1.0)


print(
    "The probability of the system being in the sate |rgrgrgr> is equal to {} ".format(
        fidelity_pure,
    )
)

Evolution of the state in time

Section titled “Evolution of the state in time”

Here, we plot the time evolution of the magnetization of the system sites, and we observe how the system slowly reaches the AFM state.

magnetization_values = np.array(list(results.occupation))
magnetization_times = results.get_result_times(density)
# rescaling the time
real_times = []
for time in magnetization_times:
    real_times.append(time * seq.get_duration())
fig, ax = plt.subplots(figsize=(8, 4), layout="constrained")


num_time_points, positions = magnetization_values.shape
x, y = np.meshgrid(np.arange(num_time_points), np.arange(1, positions + 1))
im = plt.pcolormesh(real_times, y, magnetization_values.T, shading="auto")
ax.set_xlabel("Time [ns]")
ax.set_ylabel("Qubit")
ax.set_title("State Density")
ax.set_yticks(np.arange(1, positions + 1))
cb = fig.colorbar(im, ax=ax)