{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fermion Gaussian State (FGS) Simulator\n",
    "\n",
    "This tutorial demonstrates how to use the Fermion Gaussian State (FGS) simulator implemented in `tensorcircuit-ng`. The FGS simulator allows for efficient simulation of non-interacting fermionic systems, which is particularly useful for studying free fermions on lattices.\n",
    "\n",
    "## Introduction\n",
    "\n",
    "The FGS simulator efficiently handles systems governed by quadratic Hamiltonians of the form:\n",
    "$$H = \\sum_{ij} H_{ij} c_i^\\dagger c_j + \\frac{1}{2} \\sum_{ij} (H_{ij}^{(2)} c_i^\\dagger c_j^\\dagger + h.c.)$$\n",
    "\n",
    "Instead of working with the full $2^N$-dimensional Hilbert space, the FGS simulator uses the correlation matrix formalism which scales polynomially with system size.\n",
    "\n",
    "## Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "('complex128', 'float64')"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import tensorcircuit as tc\n",
    "\n",
    "# Set the backend (using numpy for this tutorial)\n",
    "tc.set_backend(\"numpy\")\n",
    "tc.set_dtype(\"complex128\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating an FGS Simulator Instance\n",
    "\n",
    "We can initialize an FGS simulator in several ways:\n",
    "\n",
    "1. By specifying occupied sites in a product state\n",
    "2. From the groudn state of a given Hamiltonian\n",
    "3. Directly with the alpha matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initialized FGS with filled sites [0, 2]\n",
      "Initialized FGS from Hamiltonian ground states\n",
      "Alpha matrix shape: (8, 4)\n"
     ]
    }
   ],
   "source": [
    "# Method 1: Initialize with occupied sites\n",
    "# Create a 4-site system with sites 0 and 2 occupied\n",
    "sim1 = tc.FGSSimulator(L=4, filled=[0, 2])\n",
    "print(\"Initialized FGS with filled sites [0, 2]\")\n",
    "\n",
    "# Method 2: Initialize from a Hamiltonian (ground state)\n",
    "# Create a simple hopping Hamiltonian\n",
    "L = 4\n",
    "hc = np.zeros([2 * L, 2 * L])\n",
    "# Add hopping terms between neighboring sites\n",
    "for i in range(L - 1):\n",
    "    # chi * (c_i^\\dagger c_j + h.c.)\n",
    "    hc[i, i + 1] = 1.0\n",
    "    hc[i + L + 1, i + L] = -1.0\n",
    "\n",
    "sim2 = tc.FGSSimulator(L=4, hc=hc)\n",
    "print(\"Initialized FGS from Hamiltonian ground states\")\n",
    "\n",
    "# Check the alpha matrix of the first simulator\n",
    "alpha = sim1.get_alpha()\n",
    "print(f\"Alpha matrix shape: {alpha.shape}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Single-particle Green's Functions\n",
    "\n",
    "We can compute correlation functions, which are related to single-particle Green's functions:\n",
    "$$C_{ij} = \\langle c_i^\\dagger c_j \\rangle$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Correlation matrix shape: (8, 8)\n",
      "Occupation numbers:\n",
      "Site 0: 1.000\n",
      "Site 1: 0.000\n",
      "Site 2: 1.000\n",
      "Site 3: 0.000\n",
      "\n",
      "Selected off-diagonal correlations:\n",
      "<c_0 c_0^†> = 0.000+0.000j\n",
      "<c_0 c_1^†> = 0.000+0.000j\n",
      "<c_2 c_3^†> = 0.000+0.000j\n"
     ]
    }
   ],
   "source": [
    "# Get the correlation matrix\n",
    "cmatrix = sim1.get_cmatrix()\n",
    "print(f\"Correlation matrix shape: {cmatrix.shape}\")\n",
    "\n",
    "# Check occupation numbers (diagonal elements)\n",
    "print(\"Occupation numbers:\")\n",
    "for i in range(sim1.L):\n",
    "    print(f\"Site {i}: {1-cmatrix[i, i].real:.3f}\")\n",
    "\n",
    "# Compute off-diagonal correlations\n",
    "print(\"\\nSelected off-diagonal correlations:\")\n",
    "print(f\"<c_0 c_0^†> = {sim1.expectation_2body(0, 0):.3f}\")\n",
    "print(f\"<c_0 c_1^†> = {sim1.expectation_2body(0, 1):.3f}\")\n",
    "print(f\"<c_2 c_3^†> = {sim1.expectation_2body(2, 3):.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Time Evolution\n",
    "\n",
    "The FGS simulator supports evolution under quadratic Hamiltonians:\n",
    "1. Hopping terms: $\\chi c_i^\\dagger c_j + h.c.$\n",
    "2. Chemical potential terms: $\\chi c_i^\\dagger c_i$\n",
    "3. Superconducting pairing terms: $\\chi c_i^\\dagger c_j^\\dagger + h.c.$\n",
    "\n",
    "Let's demonstrate hopping:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initial state:\n",
      "Site 0 occupation: 1.000\n",
      "Site 1 occupation: 0.000\n",
      "Site 2 occupation: 0.000\n",
      "Site 3 occupation: 0.000\n",
      "\n",
      "After hopping evolution:\n",
      "Site 0 occupation: 0.000\n",
      "Site 1 occupation: 1.000\n",
      "Site 2 occupation: 0.000\n",
      "Site 3 occupation: 0.000\n"
     ]
    }
   ],
   "source": [
    "# Create a new simulator\n",
    "sim = tc.FGSSimulator(L=4, filled=[0])\n",
    "\n",
    "print(\"Initial state:\")\n",
    "cmatrix_init = sim.get_cmatrix()\n",
    "for i in range(sim.L):\n",
    "    print(f\"Site {i} occupation: {1-cmatrix_init[i, i].real:.3f}\")\n",
    "\n",
    "# Apply hopping between sites 0 and 1 with strength 1.0 for time π/2\n",
    "# This should transfer the fermion from site 0 to site 1\n",
    "sim.evol_hp(0, 1, np.pi)\n",
    "\n",
    "print(\"\\nAfter hopping evolution:\")\n",
    "cmatrix_final = sim.get_cmatrix()\n",
    "for i in range(sim.L):\n",
    "    print(f\"Site {i} occupation: {1-cmatrix_final[i, i].real:.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Entanglement Measures\n",
    "\n",
    "The FGS simulator can efficiently compute entanglement measures like von Neumann entropy and Renyi entropy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Entanglement entropy of sites [0,1]: -0.000000\n",
      "Entanglement entropy of site [0]: 0.416496\n",
      "Renyi-2 entropy of sites [0,1]: -0.000000\n",
      "Renyi-2 entropy of site [0]: 0.287682\n"
     ]
    }
   ],
   "source": [
    "# Create a simple entangled state\n",
    "sim_ent = tc.FGSSimulator(L=4, filled=[0, 2])\n",
    "# Apply a hopping that creates entanglement\n",
    "sim_ent.evol_hp(0, 1, np.pi / 4)\n",
    "\n",
    "# Compute entanglement entropy for different subsystems\n",
    "# Tracing out sites [2, 3] means we look at the entanglement of sites [0, 1]\n",
    "entropy_01 = sim_ent.entropy([2, 3])\n",
    "print(f\"Entanglement entropy of sites [0,1]: {entropy_01.real:.6f}\")\n",
    "\n",
    "# Tracing out sites [1, 2, 3] means we look at the entanglement of site [0]\n",
    "entropy_0 = sim_ent.entropy([1, 2, 3])\n",
    "print(f\"Entanglement entropy of site [0]: {entropy_0.real:.6f}\")\n",
    "\n",
    "# Compute Renyi entropy (n=2) for the same subsystems\n",
    "renyi_01 = sim_ent.renyi_entropy(2, [2, 3])\n",
    "print(f\"Renyi-2 entropy of sites [0,1]: {renyi_01.real:.6f}\")\n",
    "\n",
    "renyi_0 = sim_ent.renyi_entropy(2, [1, 2, 3])\n",
    "print(f\"Renyi-2 entropy of site [0]: {renyi_0.real:.6f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Measurements and Post-selection\n",
    "\n",
    "The FGS simulator supports both projective measurements and post-selection:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Before measurement:\n",
      "Site 0 occupation probability: 0.854\n",
      "Site 1 occupation probability: 0.146\n",
      "\n",
      "Measurement outcome for site 0: 0.0\n",
      "After measurement:\n",
      "Site 0 occupation probability: 0.000\n",
      "Site 1 occupation probability: 1.000\n",
      "\n",
      "Before post-selection:\n",
      "Site 0 occupation: 0.854\n",
      "After post-selecting site 0 as occupied:\n",
      "Site 0 occupation: 1.000\n"
     ]
    }
   ],
   "source": [
    "# Create a superposition state\n",
    "sim_meas = tc.FGSSimulator(L=2, filled=[0])\n",
    "# Put site 0 in an equal superposition of occupied and unoccupied\n",
    "# This is a simplified example - in practice, creating such states requires specific evolutions\n",
    "sim_meas.evol_hp(0, 1, np.pi / 4)\n",
    "\n",
    "print(\"Before measurement:\")\n",
    "cmatrix = sim_meas.get_cmatrix()\n",
    "for i in range(sim_meas.L):\n",
    "    print(f\"Site {i} occupation probability: {1-cmatrix[i, i].real:.3f}\")\n",
    "\n",
    "# Simulate a measurement on site 0 with a random outcome\n",
    "# In practice, you would use a random number generator\n",
    "# Here we manually specify the outcome for reproducibility\n",
    "outcome = sim_meas.cond_measure(0, status=0.1)  # Should likely result in 0 (unoccupied)\n",
    "print(f\"\\nMeasurement outcome for site 0: {outcome}\")\n",
    "\n",
    "print(\"After measurement:\")\n",
    "cmatrix_post = sim_meas.get_cmatrix()\n",
    "for i in range(sim_meas.L):\n",
    "    print(f\"Site {i} occupation probability: {1-cmatrix_post[i, i].real:.3f}\")\n",
    "\n",
    "# Demonstrate post-selection (conditioning on a specific outcome)\n",
    "sim_post = tc.FGSSimulator(L=2, filled=[0])\n",
    "sim_post.evol_hp(0, 1, np.pi / 4)\n",
    "\n",
    "print(\"\\nBefore post-selection:\")\n",
    "print(f\"Site 0 occupation: {1-sim_post.get_cmatrix()[0, 0].real:.3f}\")\n",
    "\n",
    "# Post-select on site 0 being occupied (keep=1)\n",
    "sim_post.post_select(0, keep=1)\n",
    "print(\"After post-selecting site 0 as occupied:\")\n",
    "print(f\"Site 0 occupation: {1-sim_post.get_cmatrix()[0, 0].real:.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Advanced Example: Kitaev Chain\n",
    "\n",
    "Let's simulate a simple Kitaev chain, which includes both hopping and pairing terms:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Kitaev chain ground state:\n",
      "Site occupation numbers:\n",
      "  Site 0: 0.8166\n",
      "  Site 1: 0.8837\n",
      "  Site 2: 0.8746\n",
      "  Site 3: 0.8746\n",
      "  Site 4: 0.8837\n",
      "  Site 5: 0.8166\n",
      "\n",
      "Entanglement entropy of half the chain: 0.349326\n"
     ]
    }
   ],
   "source": [
    "def kitaev_chain(L, mu, Delta, J):\n",
    "    \"\"\"\n",
    "    Create the Hamiltonian matrix for a Kitaev chain\n",
    "\n",
    "    H = -J Σ (c_i^† c_{i+1} + h.c.) - μ Σ c_i^† c_i + Δ Σ (c_i c_{i+1} + h.c.)\n",
    "    \"\"\"\n",
    "    hc = np.zeros([2 * L, 2 * L], dtype=complex)\n",
    "\n",
    "    # Chemical potential term\n",
    "    for i in range(L):\n",
    "        hc[i, i] = -mu\n",
    "        hc[i + L, i + L] = mu\n",
    "\n",
    "    # Hopping terms\n",
    "    for i in range(L - 1):\n",
    "        hc[i, i + 1] = -J\n",
    "        hc[i + L + 1, i + L] = J\n",
    "\n",
    "    # Pairing terms\n",
    "    for i in range(L - 1):\n",
    "        hc[i, i + 1 + L] = Delta\n",
    "        hc[i + 1, i + L] = -Delta\n",
    "        hc[i + L + 1, i] = Delta\n",
    "        hc[i + L, i + 1] = -Delta\n",
    "\n",
    "    return hc\n",
    "\n",
    "\n",
    "# Parameters for the Kitaev chain\n",
    "L = 6\n",
    "mu = 1.0  # Chemical potential\n",
    "Delta = 0.5  # Pairing amplitude\n",
    "J = 1.0  # Hopping amplitude\n",
    "\n",
    "# Create the Hamiltonian\n",
    "kitaev_hc = kitaev_chain(L, mu, Delta, J)\n",
    "\n",
    "# Initialize the ground state of the Kitaev chain\n",
    "sim_kitaev = tc.FGSSimulator(L=L, hc=kitaev_hc)\n",
    "\n",
    "print(\"Kitaev chain ground state:\")\n",
    "cmatrix = sim_kitaev.get_cmatrix()\n",
    "print(\"Site occupation numbers:\")\n",
    "for i in range(L):\n",
    "    print(f\"  Site {i}: {1-cmatrix[i, i].real:.4f}\")\n",
    "\n",
    "# Calculate entanglement entropy for half the chain\n",
    "half_chain_entropy = sim_kitaev.entropy(list(range(L // 2, L)))\n",
    "print(f\"\\nEntanglement entropy of half the chain: {half_chain_entropy.real:.6f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusion\n",
    "\n",
    "This tutorial demonstrated the main features of the FGS simulator:\n",
    "\n",
    "1. Initialization methods\n",
    "2. Computation of correlation functions\n",
    "3. Time evolution under quadratic Hamiltonians\n",
    "4. Entanglement measures\n",
    "5. Measurements and post-selection\n",
    "6. Application to a physical model (Kitaev chain)\n",
    "\n",
    "The FGS simulator provides an efficient way to study non-interacting fermionic systems, avoiding the exponential cost of full Hilbert space simulation while still capturing important physical properties like entanglement and correlation functions."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "2502",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}