### Pymablock Workflow - Getting Started

Source: https://github.com/quantum-tinkerer/pymablock/blob/main/docs/source/tutorial/getting_started.md

This snippet demonstrates the basic workflow of Pymablock, including importing the library, defining an unperturbed Hamiltonian (H_0), and a random perturbation (H_1). It also shows how to use the `block_diagonalize` function to define the perturbative series.

```APIDOC
## Pymablock Workflow

### Description
This section outlines the initial steps for using Pymablock, including importing necessary libraries, defining Hamiltonian matrices, and applying perturbation theory.

### Method
Python code examples using `pymablock` and `numpy`.

### Endpoint
N/A (Local Python execution)

### Parameters
N/A

### Request Example
```python
import pymablock
import numpy as np

# Diagonal unperturbed Hamiltonian
H_0 = np.diag([-1., -1., 1., 1.])

# Random Hermitian matrix as a perturbation
def random_hermitian(n):
    H = np.random.randn(n, n) + 1j * np.random.randn(n, n)
    H += H.conj().T
    return H

H_1 = 0.2 * random_hermitian(4)

hamiltonian = [H_0, H_1]

H_tilde, U, U_adjoint = pymablock.block_diagonalize(hamiltonian, subspace_indices=[0, 0, 1, 1])
```

### Response
#### Success Response (200)
- **H_tilde** (numpy.ndarray) - The transformed Hamiltonian.
- **U** (numpy.ndarray) - The transformation matrix.
- **U_adjoint** (numpy.ndarray) - The adjoint of the transformation matrix.

#### Response Example
```json
{
  "H_tilde": "[[...]]",
  "U": "[[...]]",
  "U_adjoint": "[[...]]"
}
```

### Notes
- The `subspace_indices` argument is crucial for defining the energy subspaces.
- The `block_diagonalize` function defines the perturbative series but does not perform all computations immediately.
```

--------------------------------

### Install Pymablock

Source: https://github.com/quantum-tinkerer/pymablock/blob/main/docs/source/index.md

Instructions for installing the Pymablock package using common Python package managers like mamba, conda, or pip.

```bash
mamba install pymablock -c conda-forge
```

```bash
pip install pymablock
```

--------------------------------

### Initialize Symbolic Hamiltonian Parameters and Operators

Source: https://github.com/quantum-tinkerer/pymablock/blob/main/docs/source/tutorial/andreev_supercurrent.md

Imports necessary libraries and defines symbolic variables for the quantum dot system, including charging energy, onsite energies, and fermion operators for the dot and superconductors.

```python
from IPython.display import display
import numpy as np
import sympy
from sympy import exp, I, symbols, Symbol, Eq
from sympy.physics.quantum.fermion import FermionOp
from sympy.physics.quantum import Dagger
from pymablock import block_diagonalize

U, N = symbols(r"U N", positive=True)
xis = xi_L, xi_R = symbols(r"\xi_L \xi_R", positive=True)
Gammas = Gamma_L, Gamma_R = symbols(r"\Gamma_L \Gamma_R", positive=True)
ts = t_L_complex, t_R = Symbol("t_Lc", real=False), Symbol("t_R", real=True)

c_up, c_down = FermionOp(r'c_{d, \uparrow}'), FermionOp(r'c_{d, \downarrow}')
d_ups = FermionOp(r'd_{L, \uparrow}'), FermionOp(r'd_{R, \uparrow}')
d_downs = FermionOp(r'd_{L, \downarrow}'), FermionOp(r'd_{R, \downarrow}')

def n(op):
    return Dagger(op) * op
```

--------------------------------

### Initialize SymPy environment for symbolic physics

Source: https://github.com/quantum-tinkerer/pymablock/blob/main/docs/source/tutorial/bilayer_graphene.md

Imports necessary libraries including NumPy and SymPy modules required for symbolic matrix manipulation and vector calculus in physical modeling.

```python
import numpy as np
from sympy import symbols, Matrix, sqrt, Eq, exp, I, pi, Add, MatAdd
from sympy.physics.vector import ReferenceFrame
import sympy
```

--------------------------------

### Defining Quantum Dot and Superconductor Lattice with Kwant

Source: https://github.com/quantum-tinkerer/pymablock/blob/main/docs/source/tutorial/induced_gap.md

Sets up a square lattice for a quantum dot and superconductor system using Kwant. Defines Pauli matrices and initializes a Kwant builder, specifying lattice parameters and onsite potentials for normal and superconducting regions.

```python
sigma_z = ta.array([[1, 0], [0, -1]], float)
sigma_x = ta.array([[0, 1], [1, 0]], float)

syst = kwant.Builder()
lat = kwant.lattice.square(norbs=2)
L, W = 200, 40

def normal_onsite(site, mu_n, t):
    return (-mu_n + 4 * t) * sigma_z

def sc_onsite(site, mu_sc, Delta, t):
    return (-mu_sc + 4 * t) * sigma_z + Delta * sigma_x

syst[lat.shape((lambda pos: abs(pos[1]) < W and abs(pos[0]) < L), (0, 0))] = normal_onsite
syst[lat.shape((lambda pos: abs(pos[1]) < W and abs(pos[0]) < L / 3), (0, 0))] = sc_onsite
```

--------------------------------

### Apply Second Quantization and Spin Operators

Source: https://context7.com/quantum-tinkerer/pymablock/llms.txt

Shows the usage of bosonic operators and Pauli matrices for defining Hamiltonians in the Jaynes-Cummings model context.

```python
from sympy import symbols, simplify
from sympy.physics.quantum.boson import BosonOp
from sympy.physics.quantum import Dagger, pauli
from pymablock import block_diagonalize

omega_r, omega_q, g = symbols(r'\omega_r \omega_q g', real=True)
a = BosonOp("a")
H_0 = omega_r * Dagger(a) * a + omega_q * pauli.SigmaZ("s") / 2
H_p = g * (pauli.SigmaPlus("s") * a + pauli.SigmaMinus("s") * Dagger(a))
H_tilde, *_ = block_diagonalize([H_0, H_p], symbols=[g])
H_4 = simplify(H_tilde[0, 0, 4])
```

--------------------------------

### Second Quantization with Bosonic Operators

Source: https://context7.com/quantum-tinkerer/pymablock/llms.txt

Illustrates Pymablock's support for second-quantized Hamiltonians using SymPy's quantum mechanics module, specifically the Jaynes-Cummings model. It demonstrates block diagonalization and extraction of energy corrections.

```APIDOC
## Second Quantization with Bosonic Operators

Pymablock supports second-quantized Hamiltonians using SymPy's quantum mechanics module. This example shows the Jaynes-Cummings model.

```python
from sympy import Matrix, symbols, simplify
from sympy.physics.quantum.boson import BosonOp
from sympy.physics.quantum import Dagger
from pymablock import block_diagonalize

# Define parameters
omega_r, omega_q, g = symbols(r'\omega_r \omega_q g', real=True)

# Bosonic annihilation operator for cavity mode
a = BosonOp("a")

# Jaynes-Cummings Hamiltonian in spin basis
# Diagonal: cavity + qubit energies
H_0 = Matrix([
    [omega_r * Dagger(a) * a + omega_q / 2, 0],
    [0, omega_r * Dagger(a) * a - omega_q / 2]
])

# Off-diagonal: spin-boson coupling
H_p = Matrix([
    [0, g * a],
    [g * Dagger(a), 0]
])

# Block diagonalize
H_tilde, U, U_adjoint = block_diagonalize([H_0, H_p], symbols=[g])

# Get second-order correction (dispersive shift)
H_2 = simplify(H_tilde[0, 0, 2])
print("Second-order correction (contains number operator N_a):")
print(H_2)

# Get fourth-order correction
H_4 = simplify(H_tilde[0, 0, 4])
print("Fourth-order correction:")
print(H_4)

# Convert number operators back to creation/annihilation form
H_2_expanded = H_2[0, 0].doit()
print("Expanded form:", H_2_expanded)
```
```

--------------------------------

### Constructing Hamiltonian Matrices in Python

Source: https://github.com/quantum-tinkerer/pymablock/blob/main/docs/source/tutorial/dispersive_shift.md

This code constructs the Hamiltonian matrices H_0 and H_p in a truncated Hilbert space using Pymablock. It defines creation and annihilation operators for the transmon and resonator, and then combines them to form the total Hamiltonian, including anharmonic terms and coupling.

```python
N = 4  # Number of levels for each boson
a = sympy.zeros(N, N)
for i in range(N-1):
    a[i, i+1] = sympy.sqrt(i+1)
n = sympy.diag(*[i for i in range(N)])

a_t = sympy.KroneckerProduct(a, sympy.eye(N))
a_r = sympy.KroneckerProduct(sympy.eye(N), a)

H_0 = (
    -omega_t * Dagger(a_t) * a_t + omega_r * Dagger(a_r) * a_r
    + alpha * Dagger(a_t)**2 * a_t**2 / 2
)
H_p = (
    -g * (Dagger(a_t) - a_t) * (Dagger(a_r) - a_r)
)
H = H_0 + H_p
```

--------------------------------

### Initialize Hamiltonian and Mask for Selective Diagonalization

Source: https://github.com/quantum-tinkerer/pymablock/blob/main/docs/source/tutorial/getting_started.md

This snippet demonstrates creating a random Hamiltonian and defining a custom boolean mask to select specific off-diagonal terms for preservation during the diagonalization process.

```python
import numpy as np
import matplotlib.pyplot as plt

N = 16
H_0 = np.diag(np.arange(N) + 0.2 * np.random.randn(N))
H_1 = 0.1 * np.random.randn(N, N)
H_1 += H_1.T

smiley_binary = np.array([[1, 1, 0, 0, 0, 1, 1], [1, 1, 0, 0, 0, 1, 1], [0, 0, 0, 1, 0, 0, 0], [0, 0, 0, 1, 0, 0, 0], [1, 0, 0, 0, 0, 0, 1], [0, 1, 1, 1, 1, 1, 0]], dtype=bool)
mask = np.zeros((N, N), dtype=bool)
mask[1:smiley_binary.shape[0] + 1, -smiley_binary.shape[1] - 1:-1] = smiley_binary
mask = ~(mask | mask.T)
np.fill_diagonal(mask, False)
```

--------------------------------

### Kernel Polynomial Method (KPM)

Source: https://github.com/quantum-tinkerer/pymablock/blob/main/docs/source/documentation/pymablock.md

Documentation for KPM functionalities, including greens_function and rescale.

```APIDOC
## Kernel Polynomial Method (KPM)

### Description
Implementation of the Kernel Polynomial Method.

### Module
pymablock.kpm

### Functions
- greens_function
- rescale

### Inheritance
- show-inheritance: True
```

--------------------------------

### Multivariate Perturbation Theory (Python)

Source: https://context7.com/quantum-tinkerer/pymablock/llms.txt

Shows how Pymablock handles multiple perturbation parameters using a dictionary input format. Each key in the dictionary represents the order of each perturbation parameter, allowing for complex, multivariate expansions of the effective Hamiltonian.

```python
import numpy as np
from pymablock import block_diagonalize

# Unperturbed Hamiltonian
H_00 = np.diag([0., 1., 2., 5.])

# Perturbations with different orders in two parameters (lambda_1, lambda_2)
H_10 = 0.1 * np.random.randn(4, 4)  # Linear in lambda_1
H_10 = H_10 + H_10.T  # Make Hermitian

H_01 = 0.1 * np.random.randn(4, 4)  # Linear in lambda_2
H_01 = H_01 + H_01.T

H_20 = 0.05 * np.random.randn(4, 4)  # Quadratic in lambda_1
H_20 = H_20 + H_20.T

# Dictionary format: {(order_lambda1, order_lambda2): perturbation}
hamiltonian = {
    (0, 0): H_00,
    (1, 0): H_10,
    (0, 1): H_01,
    (2, 0): H_20
}

# Block diagonalize with two subspaces
H_tilde, U, U_adjoint = block_diagonalize(
    hamiltonian,
    subspace_indices=[0, 0, 1, 1]
)

# Query specific multivariate order: lambda_1^2 * lambda_2^1
H_AA_21 = H_tilde[0, 0, 2, 1]
print("Correction at order lambda_1^2 * lambda_2:")
print(H_AA_21)

# Get all terms up to certain orders
H_AA_terms = H_tilde[0, 0, :3, :2]  # Up to lambda_1^2, lambda_2^1
print("Shape of result:", H_AA_terms.shape)
```

--------------------------------

### Finalizing and Plotting Kwant System

Source: https://github.com/quantum-tinkerer/pymablock/blob/main/docs/source/tutorial/induced_gap.md

Finalizes the Kwant system builder to create a finalized system object and generates a plot to visualize the quantum dot and superconductor structure. The plot uses color to distinguish between quantum dots (blue) and the superconductor (red).

```python
kwant.plot(
    syst,
    fig_size=(10, 6),
    site_color=(lambda site: abs(site.pos[0]) < L / 3),
    colorbar=False,
    cmap="seismic",
    hop_lw=0,
)

syst = syst.finalized()
f"The system has {len(syst.sites)} sites."
```

--------------------------------

### Importing Libraries for Quantum Tinkering

Source: https://github.com/quantum-tinkerer/pymablock/blob/main/docs/source/tutorial/induced_gap.md

Imports essential Python libraries for numerical computation, quantum mechanics simulations, and plotting. Includes specialized libraries like tinyarray, scipy, numpy, kwant, and pymablock.

```python
import tinyarray as ta
import matplotlib.backends
import scipy.linalg
from scipy.sparse.linalg import eigsh
import numpy as np
import kwant
import matplotlib.pyplot as plt
color_cycle = ["#5790fc", "#f89c20", "#e42536"]

from pymablock import block_diagonalize
```

--------------------------------

### Implicit Mode for Large Sparse Systems

Source: https://context7.com/quantum-tinkerer/pymablock/llms.txt

Demonstrates Pymablock's implicit solver for large sparse Hamiltonians, focusing on a small subspace of interest. This avoids full diagonalization by using iterative methods and direct solvers like MUMPS.

```python
import numpy as np
from scipy import sparse
from pymablock import block_diagonalize
from scipy.sparse.linalg import eigsh

N = 1000
diagonal = np.linspace(-10, 10, N)
H_0 = sparse.diags(diagonal, format='csr')

H_1 = 0.01 * sparse.random(N, N, density=0.01, format='csr')
H_1 = H_1 + H_1.T

_, low_energy_vecs = eigsh(H_0, k=5, which='SM')

H_tilde, U, U_adjoint = block_diagonalize(
    [H_0, H_1],
    subspace_eigenvectors=[low_energy_vecs],
    direct_solver=True
)

H_eff_2 = H_tilde[0, 0, 2]
print(f"Second-order correction shape: {H_eff_2.shape}")
print(f"Effective 5x5 Hamiltonian computed from {N}x{N} system")
```

--------------------------------

### Verify 2x2 Block-Diagonalization with Random Hamiltonian (Python)

Source: https://github.com/quantum-tinkerer/pymablock/blob/main/docs/source/radius.md

This code snippet verifies the correctness of algorithms by comparing results for a 2x2 block-diagonalization case. It generates a random Hamiltonian and sets up the necessary parameters for comparison. Dependencies include NumPy.

```python
import numpy as np

np.random.seed(0)
N = 10
H_0 = np.diag(np.arange(N))
H_1 = np.random.rand(N, N) + 1j * np.random.rand(N, N)
H_1 += H_1.conj().T
```

--------------------------------

### Symbolic Hamiltonians with SymPy

Source: https://context7.com/quantum-tinkerer/pymablock/llms.txt

Demonstrates how Pymablock can work with symbolic matrices using SymPy for analytical perturbation theory calculations. It shows block diagonalization of a symbolic Hamiltonian and extraction of energy corrections.

```APIDOC
## Symbolic Hamiltonians with SymPy

Pymablock works seamlessly with symbolic matrices, enabling analytical perturbation theory calculations.

```python
import sympy
from sympy import symbols, Matrix, sqrt, Rational
from pymablock import block_diagonalize

# Define symbolic parameters
E_a, E_b, V, g = symbols('E_a E_b V g', real=True)

# Symbolic unperturbed Hamiltonian (diagonal)
H_0 = Matrix([
    [E_a, 0],
    [0, E_b]
])

# Symbolic perturbation
H_p = Matrix([
    [0, V],
    [V, 0]
])

# Block diagonalize with g as perturbative parameter
H_tilde, U, U_adjoint = block_diagonalize(
    {sympy.S.One: H_0, g: H_p},  # Dictionary with symbolic keys
)

# Get second-order correction
H_2 = H_tilde[0, 0, 2]
print("Second-order energy correction:")
print(sympy.simplify(H_2))
# Output: V**2/(E_a - E_b)

# Get fourth-order correction
H_4 = H_tilde[0, 0, 4]
print("Fourth-order energy correction:")
print(sympy.simplify(H_4))
```
```

--------------------------------

### Second Quantization Module

Source: https://github.com/quantum-tinkerer/pymablock/blob/main/docs/source/documentation/pymablock.md

Details functionalities for second quantization, including Sylvester equation solvers and operator masking.

```APIDOC
## Second Quantization

### Description
Functions for operations in second quantization.

### Module
pymablock.second_quantization

### Functions
- solve_sylvester_2nd_quant
- apply_mask_to_operator

### Inheritance
- show-inheritance: True
```

--------------------------------

### Spin Operators

Source: https://context7.com/quantum-tinkerer/pymablock/llms.txt

Demonstrates the use of Pauli spin operators within Pymablock for more compact Hamiltonian definitions, applied to the Jaynes-Cummings model.

```APIDOC
## Spin Operators

Pymablock also supports Pauli spin operators for more compact Hamiltonian definitions.

```python
from sympy import symbols, simplify
from sympy.physics.quantum.boson import BosonOp
from sympy.physics.quantum import Dagger, pauli
from pymablock import block_diagonalize

# Parameters
omega_r, omega_q, g = symbols(r'\omega_r \omega_q g', real=True)
a = BosonOp("a")

# Jaynes-Cummings using spin operators (scalar form)
H_0 = omega_r * Dagger(a) * a + omega_q * pauli.SigmaZ("s") / 2
H_p = g * (pauli.SigmaPlus("s") * a + pauli.SigmaMinus("s") * Dagger(a))

# Block diagonalize
H_tilde, *_ = block_diagonalize([H_0, H_p], symbols=[g])

# Fourth-order correction depends on spin number operator N_s
H_4 = simplify(H_tilde[0, 0, 4])
print("Fourth-order correction with spin:")
print(H_4)
```
```

--------------------------------

### Import SymPy Quantum Modules

Source: https://github.com/quantum-tinkerer/pymablock/blob/main/docs/source/tutorial/jaynes_cummings.md

Imports necessary functions from SymPy for quantum mechanics, including matrix operations and bosonic operators. These are foundational for defining quantum Hamiltonians.

```python
from sympy import Matrix, Symbol, symbols, Eq, simplify
from sympy.physics.quantum.boson import BosonOp
from sympy.physics.quantum import Dagger
```

--------------------------------

### Block Diagonalize with Eigenvectors (Python)

Source: https://context7.com/quantum-tinkerer/pymablock/llms.txt

Illustrates how to use `block_diagonalize` when the unperturbed Hamiltonian is not diagonal. It involves diagonalizing the unperturbed Hamiltonian first to obtain eigenvectors, which are then used to define the subspaces for the block diagonalization process.

```python
import numpy as np
from pymablock import block_diagonalize

# Non-diagonal unperturbed Hamiltonian
H_00 = np.array([
    [1.0, 0.2, 0.1],
    [0.2, 2.0, 0.3],
    [0.1, 0.3, 5.0]
])

# Perturbation
H_10 = 0.1 * np.array([
    [0.1, 0.2, 0.05],
    [0.2, 0.1, 0.15],
    [0.05, 0.15, 0.1]
])

# Diagonalize H_00 to get eigenvectors
eigenvalues, eigenvectors = np.linalg.eigh(H_00)

# Split eigenvectors into subspaces (2 states in A, 1 state in B)
subspace_eigenvectors = [eigenvectors[:, :2], eigenvectors[:, 2:]]

# Block diagonalize with eigenvectors
H_tilde, U, U_adjoint = block_diagonalize(
    [H_00, H_10],
    subspace_eigenvectors=subspace_eigenvectors
)

# Get effective Hamiltonian for subspace A up to second order
H_eff_A = H_tilde[0, 0, 0] + H_tilde[0, 0, 1] + H_tilde[0, 0, 2]
print("Effective Hamiltonian for subspace A:")
print(H_eff_A)
```

--------------------------------

### Transforming Operators with BlockDiagonalize

Source: https://github.com/quantum-tinkerer/pymablock/blob/main/docs/source/tutorial/getting_started.md

This code demonstrates how to transform an arbitrary operator (X) into the basis of the effective Hamiltonian using Pymablock. It involves creating BlockSeries representations and applying unitary transformations.

```ipython3
H_0 = random_hermitian(4)
H_1 = 0.1 * random_hermitian(4)
X = random_hermitian(4)
_, evecs = np.linalg.eigh(H_0)

H_tilde, U, U_adjoint = block_diagonalize(
    [H_0, H_1], subspace_eigenvectors=[evecs[:, :2], evecs[:, 2:]]
)
```

```ipython3
X_series = pymablock.operator_to_BlockSeries(
    {(0,): X},
    name="X",
    hermitian=True,  # Not important here, but slightly improves performance.
    subspace_eigenvectors=[evecs[:, :2], evecs[:, 2:]]
)
```

```ipython3
X_tilde = pymablock.series.cauchy_dot_product(U_adjoint, X_series, U)
```

```ipython3
assert np.allclose(X_tilde[0, 0, 0], X_series[0, 0, 0])
X_tilde[0, 0, 2]
```

--------------------------------

### Construct and Display Full Hamiltonian

Source: https://github.com/quantum-tinkerer/pymablock/blob/main/docs/source/tutorial/jaynes_cummings.md

Constructs the unperturbed (H_0) and perturbative (H_p) parts of the Hamiltonian using SymPy matrices and operators. It then defines a helper function to display the combined Hamiltonian symbolically.

```python
H_0 = Matrix([[wr * Dagger(a) * a + wq / 2, 0], [0, wr * Dagger(a) * a - wq / 2]])
H_p = Matrix([[0,  g * a], [g * Dagger(a), 0]])

def display_eq(title, expr):
    return Eq(Symbol(title), expr, evaluate=False)

display_eq('H', H_0 + H_p)
```

--------------------------------

### General Hamiltonians Definition

Source: https://github.com/quantum-tinkerer/pymablock/blob/main/docs/source/tutorial/getting_started.md

Defines a general Hamiltonian with two perturbative parameters and demonstrates how to represent it using a dictionary. It also includes a hint for efficient perturbation strength evaluation.

```APIDOC
## General Hamiltonians

### Description
Defines a general Hamiltonian with two perturbative parameters $\lambda_1$ and $\lambda_2$ in the form $H = H_{00} + \lambda_1 H_{10} + \lambda_1^2 H_{20} + \lambda_2 H_{01}$. The Hamiltonian is represented as a dictionary where keys are tuples indicating the order of perturbative parameters.

### Parameters
#### Request Body
- **H_00** (numpy.ndarray) - Unperturbed Hamiltonian.
- **H_10** (numpy.ndarray) - Linear term in the first perturbative parameter.
- **H_20** (numpy.ndarray) - Quadratic term in the first perturbative parameter.
- **H_01** (numpy.ndarray) - Linear term in the second perturbative parameter.

### Request Example
```python
H_00 = random_hermitian(7)
H_10 = 0.1 * random_hermitian(7)
H_20 = 0.1 * random_hermitian(7)
H_01 = 0.1 * random_hermitian(7)
hamiltonian = {(0, 0): H_00, (1, 0): H_10, (2, 0): H_20, (0, 1): H_01}
```

### Efficiency Hint
To evaluate $\tilde{H}$, provide the values of the perturbative parameters at the last step. This is more efficient than recomputing the perturbative series if only the perturbation strength varies.
```

--------------------------------

### Block Diagonalization with Three or More Subspaces

Source: https://context7.com/quantum-tinkerer/pymablock/llms.txt

Shows how to perform block diagonalization into three or more distinct subspaces. This allows for the analysis of effective Hamiltonians within each subspace for systems with multiple energy scales.

```python
import numpy as np
from pymablock import block_diagonalize

H_0 = np.diag([0., 0.1, 1., 1.1, 5., 5.1])

H_1 = 0.05 * np.random.randn(6, 6)
H_1 = H_1 + H_1.T

_, evecs = np.linalg.eigh(H_0)
subspace_eigenvectors = [
    evecs[:, :2],
    evecs[:, 2:4],
    evecs[:, 4:]
]

H_tilde, U, U_adjoint = block_diagonalize(
    [H_0, H_1],
    subspace_eigenvectors=subspace_eigenvectors
)

print(f"Block structure: {H_tilde.shape}")

H_AA = H_tilde[0, 0, 0] + H_tilde[0, 0, 2]
H_BB = H_tilde[1, 1, 0] + H_tilde[1, 1, 2]
H_CC = H_tilde[2, 2, 0] + H_tilde[2, 2, 2]

print("Effective H_AA:\n", H_AA)
print("Effective H_BB:\n", H_BB)
print("Effective H_CC:\n", H_CC)
```

--------------------------------

### Subsetting and Querying BlockSeries Entries

Source: https://github.com/quantum-tinkerer/pymablock/blob/main/docs/source/tutorial/getting_started.md

This code illustrates how to create a sub-block of a BlockSeries by slicing and how to query specific terms of the perturbative series. It shows that querying a term triggers computation and caching of intermediate results.

```ipython3
H_tilde_subblocks = H_tilde[:2, :2]
print(type(H_tilde_subblocks))
print(H_tilde_subblocks.shape)
print(H_tilde_subblocks.n_infinite)
```

```ipython3
%time H_tilde[0, 0, 2, 3]
```

```ipython3
f"{len(H_tilde._data.keys())=}, {len(U._data.keys())=}"
```

```ipython3
%time H_tilde[0, 0, 2, 3]
```

```ipython3
%time H_tilde_subblocks[0, 0, 2, 3]
```

--------------------------------

### Define Hamiltonian and Perturbation

Source: https://github.com/quantum-tinkerer/pymablock/blob/main/docs/source/tutorial/dispersive_shift.md

Initializes the transmon-resonator system using SymPy Boson operators and defines the unperturbed Hamiltonian H0 and the coupling perturbation Hp.

```python
from itertools import product
import numpy as np
import sympy
from sympy.physics.quantum import Dagger
from sympy.physics.quantum.boson import BosonOp

symbols = sympy.symbols(r"\omega_{t} \omega_{r} \alpha g", real=True, positive=True)
omega_t, omega_r, alpha, g = symbols
a_t, a_r = BosonOp("a_t"), BosonOp("a_r")

H_0 = (-omega_t * Dagger(a_t) * a_t + omega_r * Dagger(a_r) * a_r + alpha * Dagger(a_t)**2 * a_t**2 / 2)
H_p = (-g * (Dagger(a_t) - a_t) * (Dagger(a_r) - a_r))
```

--------------------------------

### Block Diagonalize Hamiltonian (Python)

Source: https://context7.com/quantum-tinkerer/pymablock/llms.txt

Demonstrates the primary use of `block_diagonalize` with numerical numpy arrays for a simple two-subspace Hamiltonian. It shows how to define the unperturbed and perturbed Hamiltonians, perform the diagonalization, and extract specific orders of the effective Hamiltonian for a subspace.

```python
import numpy as np
from pymablock import block_diagonalize

# Define diagonal unperturbed Hamiltonian with two subspaces
H_0 = np.diag([-1., -1., 1., 1.])

# Define a Hermitian perturbation
H_1 = np.array([
    [0.0, 0.1, 0.05, 0.02],
    [0.1, 0.0, 0.03, 0.04],
    [0.05, 0.03, 0.0, 0.1],
    [0.02, 0.04, 0.1, 0.0]
])

# Perform block diagonalization
# subspace_indices: 0 for subspace A, 1 for subspace B
H_tilde, U, U_adjoint = block_diagonalize(
    [H_0, H_1],
    subspace_indices=[0, 0, 1, 1]
)

# Get second-order correction to the effective Hamiltonian of subspace A
H_AA_2 = H_tilde[0, 0, 2]
print("Second-order correction to H_AA:")
print(H_AA_2)
# Output: 2x2 numpy array with perturbative corrections

# Verify off-diagonal blocks are zero (block-diagonalized)
H_AB_3 = H_tilde[0, 1, 3]
print(f"Off-diagonal block (should be zero): {H_AB_3}")
# Output: zero

# Get multiple orders at once using numpy-style indexing
H_AA_up_to_3 = H_tilde[0, 0, :4]  # Orders 0, 1, 2, 3
print("Orders 0-3 of H_AA:", H_AA_up_to_3)
```

--------------------------------

### Define Floquet Hamiltonian with Pymablock

Source: https://github.com/quantum-tinkerer/pymablock/blob/main/docs/source/tutorial/spin_rwa_floquet.md

This snippet initializes the symbolic representation of a periodically driven spin-1/2 system. It defines the ladder operators, Pauli matrices, and constructs the Floquet Hamiltonian using Pymablock's LadderOp and NumberOperator classes.

```python
from IPython.display import display
import numpy as np
from sympy import Matrix, Symbol, symbols, Eq, simplify, cancel, I
from sympy.physics.quantum import Dagger
from sympy.physics.quantum import pauli
from pymablock.block_diagonalization import block_diagonalize
from pymablock.number_ordered_form import LadderOp, NumberOperator

def display_eq(title, value):
    display(Eq(Symbol(title), value, evaluate=False))

omega_0, Omega, g = symbols('omega_0 Omega g', real=True)

a = LadderOp("a")
N_a = NumberOperator(a)

sigma_x = pauli.SigmaX("s")
sigma_z = pauli.SigmaZ("s")

H_0 = omega_0/2 * sigma_z + Omega * N_a
H_p = g/2 * sigma_x * (a + Dagger(a))

display_eq('H_{Floquet}', H_0 + H_p)
```

--------------------------------

### Perform Symbolic Perturbation Theory with SymPy

Source: https://context7.com/quantum-tinkerer/pymablock/llms.txt

Demonstrates how to define symbolic Hamiltonians and compute higher-order energy corrections using block_diagonalize with SymPy symbols.

```python
import sympy
from sympy import symbols, Matrix, sqrt, Rational
from pymablock import block_diagonalize

E_a, E_b, V, g = symbols('E_a E_b V g', real=True)
H_0 = Matrix([[E_a, 0], [0, E_b]])
H_p = Matrix([[0, V], [V, 0]])
H_tilde, U, U_adjoint = block_diagonalize({sympy.S.One: H_0, g: H_p})
H_2 = H_tilde[0, 0, 2]
print(sympy.simplify(H_2))
```

--------------------------------

### Import Pymablock and NumPy

Source: https://github.com/quantum-tinkerer/pymablock/blob/main/docs/source/tutorial/getting_started.md

Imports the Pymablock library and the NumPy library for numerical operations. These are essential for setting up and manipulating matrices in Pymablock.

```python
import pymablock
import numpy as np
```

--------------------------------

### Defining Perturbation Matrices (barrier, delta_mu)

Source: https://github.com/quantum-tinkerer/pymablock/blob/main/docs/source/tutorial/induced_gap.md

Calculates the sparse matrices representing the barrier perturbation and the asymmetry in dot potentials. These matrices are derived using Kwant's `hamiltonian_submatrix` and `operator.Density` for use in perturbative calculations.

```python
barrier = syst.hamiltonian_submatrix(
    params={**{p: 0 for p in params.keys()}, "t_barrier": 1}, sparse=True
).real
delta_mu = (
    kwant.operator.Density(syst, (lambda site: sigma_z * site.pos[0] / L)).tocoo().real
)
```

--------------------------------

### Perform Full Perturbative Diagonalization

Source: https://github.com/quantum-tinkerer/pymablock/blob/main/docs/source/tutorial/spin_rwa_floquet.md

Demonstrates how to use block_diagonalize to eliminate off-diagonal terms in a Hamiltonian. It retrieves effective Hamiltonian orders and simplifies the resulting symbolic expressions.

```python
H_full, U_full, U_adjoint_full = block_diagonalize(
    [H_0, H_p],
    symbols=[g],
)

display_eq('H_{eff}^{(0)}', H_full[0, 0, 0])
display_eq('H_{eff}^{(1)}', H_full[0, 0, 1])
display_eq('H_{eff}^{(2)}', H_full[0, 0, 2])
display_eq('H_{eff}^{(2)}', simplify(H_full[0, 0, 2]).doit().expand())
```

--------------------------------

### Defining Hoppings and Barrier in Kwant System

Source: https://github.com/quantum-tinkerer/pymablock/blob/main/docs/source/tutorial/induced_gap.md

Defines the hopping parameters between lattice sites and specifically for the tunnel barrier between the quantum dot and superconductor. This involves setting up lambda functions to determine hopping based on site positions and parameters.

```python
syst[lat.neighbors()] = lambda site1, site2, t: -t * sigma_z

def barrier(site1, site2):
    return (abs(site1.pos[0]) - L / 3) * (abs(site2.pos[0]) - L / 3) < 0

syst[(hop for hop in syst.hoppings() if barrier(*hop))] = (
    lambda site1, site2, t_barrier: -t_barrier * sigma_z
)
```

--------------------------------

### Series Module

Source: https://github.com/quantum-tinkerer/pymablock/blob/main/docs/source/documentation/pymablock.md

Documentation for the BlockSeries class and related functions/data in the series module.

```APIDOC
## Series

### Description
Handles series expansions and related operations.

### Classes
- BlockSeries
  - members: (all members of the class)
  - show-inheritance: True
  - class-doc-from: class

### Module
pymablock.series

### Functions
- cauchy_dot_product

### Data
- zero
- one

### Inheritance
- show-inheritance: True
```

--------------------------------

### Calculating Unperturbed Hamiltonian (h_0) with Kwant

Source: https://github.com/quantum-tinkerer/pymablock/blob/main/docs/source/tutorial/induced_gap.md

Computes the sparse Hamiltonian matrix (h_0) for the unperturbed system using Kwant's `hamiltonian_submatrix` function. This is crucial for subsequent analysis involving sparse matrices and low-energy eigenvectors.

```python
params = dict(
    mu_n=0.05,
    mu_sc=0.3,
    Delta=0.05,
    t=1.,
    t_barrier=0.,
)

h_0 = syst.hamiltonian_submatrix(params=params, sparse=True).real
```

--------------------------------

### Define Elimination Rules for Rotating Wave Approximation

Source: https://github.com/quantum-tinkerer/pymablock/blob/main/docs/source/tutorial/spin_rwa_floquet.md

Illustrates how to define custom elimination rules to target specific terms for removal during the perturbative expansion. This is useful for implementing the Rotating Wave Approximation by filtering out counter-rotating terms.

```python
k = Symbol('k', integer=True, positive=True)
sigma_plus = pauli.SigmaPlus("s")
to_eliminate = sigma_plus * a + Dagger(sigma_plus) * Dagger(a) + sigma_z * (a**k + Dagger(a)**k)
```

--------------------------------

### Convert Number-Ordered Forms to SymPy Expressions

Source: https://github.com/quantum-tinkerer/pymablock/blob/main/docs/source/second_quantization.md

Demonstrates how to convert a BlockSeries of number-ordered forms into standard SymPy expressions and how to evaluate them back into their original operator representation.

```python
# Convert number-ordered forms to regular sympy expressions
expr = result.applyfunc(lambda x: x.as_expr())

# Replace number operators with original creation/annihilation operators
original_operators = result.doit()
```

--------------------------------

### Benchmark and Compare Algorithm Convergence

Source: https://github.com/quantum-tinkerer/pymablock/blob/main/docs/source/radius.md

Defines a function to compare the error norms of eigenvalues between Pymablock and Schrieffer-Wolff methods as a function of perturbation order and alpha. It generates plots for error comparison and performance timing.

```python
def compare_schrieffer_wolff_pymablock(H_0, H_1, mask, n_max=100, alpha_max=0.5, title=""):
    H_tilde, *_ = pymablock.block_diagonalize([H_0, H_1], fully_diagonalize={0: mask})
    H_tilde_sw, S = schrieffer_wolff(H_0, H_1, mask)

    alpha = np.linspace(0, alpha_max, 100).reshape(1, 1, 1, -1)
    powers = np.arange(n_max).reshape(-1, 1, 1, 1)

    start_time = time.time()
    H_orders = H_tilde[0, 0, :n_max]
    H_orders[0] = H_0
    time_H_orders = time.time() - start_time

    start_time = time.time()
    H_orders_sw = H_tilde_sw[:n_max]
    time_H_orders_sw = time.time() - start_time

    eigvals_la = np.linalg.eigvalsh(np.cumsum(np.array(list(H_orders))[..., None] * alpha**powers, axis=0).transpose(0, 3, 1, 2))
    eigvals_sw = np.linalg.eigvalsh(np.cumsum(np.array(list(H_orders_sw))[..., None] * alpha**powers, axis=0).transpose(0, 3, 1, 2))
    eigvals_exact = np.linalg.eigvalsh((H_0[..., None] + H_1[..., None] * alpha[0, 0, 0]).transpose(2, 0, 1))[None, ...]

    fig = plt.figure(layout="constrained")
    ax1, ax2 = fig.subplots(2, 1, sharex=True)

    for order in range(n_max):
        la_err = np.linalg.norm(eigvals_la[order] - eigvals_exact, axis=(0, -1))
        sw_err = np.linalg.norm(eigvals_sw[order] - eigvals_exact, axis=(0, -1))
        ratio = np.divide(sw_err, la_err, out=np.zeros_like(sw_err), where=la_err != 0)
        ax1.plot(alpha.flatten()[1:], la_err[1:], c=plt.cm.inferno(order / n_max))
        ax2.plot(alpha.flatten()[1:], ratio[1:], c=plt.cm.inferno(order / n_max))

    ax1.set_ylabel("Error Pymablock")
    ax1.set_title(title)
    ax2.set_ylabel("Error SW / Error Pymablock")
    ax2.set_title(f"Time Pymablock {time_H_orders:.2f}s, Time SW: {time_H_orders_sw:.2f}s")
    ax2.set_xlabel(r"$\alpha$")
    for ax in (ax1, ax2): ax.semilogy()
```

--------------------------------

### Perform Block-Diagonalization with Pymablock

Source: https://github.com/quantum-tinkerer/pymablock/blob/main/docs/source/radius.md

Initializes a boolean mask for a 2x2 block-diagonalization and executes the block_diagonalization process using pymablock and a Schrieffer-Wolff implementation.

```python
mask = np.zeros((N, N), dtype=bool)
mask[: N // 2, N // 2 :] = True
mask[N // 2 :, : N // 2] = True

H_tilde, *_ = pymablock.block_diagonalize([H_0, H_1], fully_diagonalize={0: mask})
H_tilde_sw, S = schrieffer_wolff(H_0, H_1, mask)

n = 5
np.testing.assert_almost_equal(H_tilde_sw[n], H_tilde[0, 0, n])
```

--------------------------------

### General Momentum Grouping with np.mgrid (Python)

Source: https://github.com/quantum-tinkerer/pymablock/blob/main/docs/source/tutorial/bilayer_graphene.md

Provides a more general method to group terms by momentum power using `np.mgrid`. It generates power grids and selects terms based on the sum of powers, which is more scalable than manual array definition.

```python
k_powers = np.mgrid[:4, :4]
k_square = k_powers[..., np.sum(k_powers, axis=0) == 2]
k_cube = k_powers[..., np.sum(k_powers, axis=0) == 3]
```

--------------------------------

### Apply Elimination Rules in Block Diagonalization

Source: https://github.com/quantum-tinkerer/pymablock/blob/main/docs/source/second_quantization.md

Shows how to pass the defined elimination rules to the block_diagonalize function to control term truncation in the effective Hamiltonian.

```python
H_tilde, *_ = block_diagonalize(
    [H_0, H_p],
    fully_diagonalize=elimination_rules,
    symbols=[g],
)
```

--------------------------------

### Specifying Subspaces for Block Diagonalization

Source: https://github.com/quantum-tinkerer/pymablock/blob/main/docs/source/tutorial/getting_started.md

Demonstrates how to specify subspaces by splitting the eigenvectors of the unperturbed Hamiltonian ($H_{00}$) and then using these subspaces to block-diagonalize the full Hamiltonian.

```APIDOC
## Specifying the Subspaces

### Description
This section explains how to define perturbative series by splitting the eigenvectors of $H_{00}$ into groups that define different subspaces. The `block_diagonalize` function is then used to transform the Hamiltonian into these subspaces.

### Method
1. Compute the eigenvectors of $H_{00}$.
2. Split the eigenvectors into distinct subspace sets.
3. Use `pymablock.block_diagonalize` with the Hamiltonian and subspace eigenvectors.

### Parameters
#### Request Body
- **hamiltonian** (dict) - The Hamiltonian dictionary as defined previously.
- **subspace_eigenvectors** (list of numpy.ndarray) - A list where each element is a numpy array containing the eigenvectors for a specific subspace.

### Request Example
```python
# Assuming H_00 and hamiltonian are already defined
_, evecs = np.linalg.eigh(H_00)
subspace_eigenvectors = [evecs[:, :3], evecs[:, 3:6], evecs[:, 6:]]  # Example: three subspaces

H_tilde, U, U_adjoint = block_diagonalize(
  hamiltonian=hamiltonian,
  subspace_eigenvectors=subspace_eigenvectors
)
```

### Important Note
The `block_diagonalize` function transforms all Hamiltonians to the basis of `subspace_vectors`. Consequently, the unperturbed Hamiltonian ($H_{00}$) becomes diagonal in this new basis. The returned `U` is the unitary transformation that achieves this block-diagonalization.

### Response
#### Success Response (200)
- **H_tilde** (dict) - The block-diagonalized Hamiltonian.
- **U** (dict) - The unitary transformation matrix.
- **U_adjoint** (dict) - The adjoint of the unitary transformation matrix.
```

--------------------------------

### Substitute Hamiltonian Parameters for Perturbation

Source: https://github.com/quantum-tinkerer/pymablock/blob/main/docs/source/tutorial/andreev_supercurrent.md

Substitutes symbolic variables in the Hamiltonian with physical parameters and introduces a perturbative parameter delta-phi to facilitate energy derivative calculations.

```python
t_L, phi, dphi = symbols(r"t_L \phi \delta\phi", positive=True)
H = H.subs({t_L_complex: t_L * exp(I * (phi + dphi))})

E_0, E_1, E_2 = symbols(r"E_0 E_1 E_2", positive=True)
E_0_value, E_1_value, E_2_value = [(U * (N - i)**2 / 2).expand() for i in range(3)]
H = H.subs({E_2_value: E_2}).subs({E_1_value: E_1}).subs({E_0_value: E_0})
```