Z2gauge_frame.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "6795344a",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: redefinition of constant X. This may fail, cause incorrect answers, or produce other errors.\n",
      "WARNING: redefinition of constant Y. This may fail, cause incorrect answers, or produce other errors.\n",
      "WARNING: redefinition of constant Z. This may fail, cause incorrect answers, or produce other errors.\n",
      "WARNING: redefinition of constant II. This may fail, cause incorrect answers, or produce other errors.\n"
     ]
    }
   ],
   "source": [
    "using LinearAlgebra\n",
    "using Combinatorics\n",
    "using Statistics\n",
    "using PyPlot\n",
    "using BenchmarkTools\n",
    "using Profile\n",
    "\n",
    "const X = [0, 1, 1, 0]\n",
    "const Y = [0, 1im, -1im, 0]\n",
    "const Z = [1, 0, 0, -1]\n",
    "const II = [1, 0, 0, 1];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "7ba90395",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "z2gauge_2nd_frame_op (generic function with 1 method)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function n_choose_k(n, k)\n",
    "    if k==0\n",
    "        return [Bool[0 for _ in 1:n]]\n",
    "    elseif k==n\n",
    "        return [Bool[1 for _ in 1:n]]\n",
    "    else\n",
    "        a = n_choose_k(n-1, k)\n",
    "        b = n_choose_k(n-1, k-1)\n",
    "        return [\n",
    "            [Bool[0, i...] for i in a]..., \n",
    "            [Bool[1, i...] for i in b]...\n",
    "        ]\n",
    "    end\n",
    "end\n",
    "\n",
    "\n",
    "function charge_q_indices(q)\n",
    "    indices = Int[]\n",
    "    \n",
    "    tmp = n_choose_k(4, q)\n",
    "    for a in tmp\n",
    "        idx = sum([2^(i-1) * a[i] for i in 1:n]) + 1\n",
    "        push!(indices, idx)\n",
    "    end\n",
    "    \n",
    "    tmp = n_choose_k(4, 4-q)\n",
    "    for a in tmp\n",
    "        idx = sum([2^(i-1) * a[i] for i in 1:n]) + 1\n",
    "        push!(indices, idx)\n",
    "    end\n",
    "    \n",
    "    sort!(indices)\n",
    "    return indices\n",
    "end\n",
    "\n",
    "\n",
    "# \"\"\"\n",
    "# Hilbert space:\n",
    "# 1 \\otimes 1bar \\otimes 2 \\otimes 2bar\n",
    "# \"\"\"\n",
    "function charged_paired_state(n, perm, q1, q2)\n",
    "    d = 2^n\n",
    "    result = zeros(Int, d^4)\n",
    "    indices1 = charge_q_indices(q1)\n",
    "    indices2 = charge_q_indices(q2)\n",
    "    for a1 in indices1, a2 in indices2\n",
    "        if perm==false\n",
    "            idx = (a1-1) * d^0 + (a1-1) * d^1 + (a2-1) * d^2 + (a2-1) * d^3\n",
    "        else\n",
    "            idx = (a1-1) * d^0 + (a2-1) * d^1 + (a2-1) * d^2 + (a1-1) * d^3\n",
    "        end\n",
    "        result[idx + 1] += 1\n",
    "    end\n",
    "    return result\n",
    "end\n",
    "\n",
    "\n",
    "function z2gauge_2nd_frame_op()\n",
    "    d = 2^4\n",
    "    part1 = zeros(Float32, d^4 * d^4)\n",
    "    part2 = zeros(Float32, d^4 * d^4)\n",
    "    \n",
    "    for q1 in 0:3, q2 in 0:3\n",
    "        if q1==q2\n",
    "            continue\n",
    "        end\n",
    "        tmp1 = charged_paired_state(n, false, q1, q2)\n",
    "        tmp2 = charged_paired_state(n, true, q1, q2)\n",
    "        part1 += (kron(tmp1, tmp1) + kron(tmp2, tmp2)) / (binomial(n, q1) * binomial(n, q2))\n",
    "    end\n",
    "    \n",
    "    for q in 0:n\n",
    "        dq = binomial(n, q)\n",
    "        tmp1 = charged_paired_state(n, false, q, q)\n",
    "        tmp2 = charged_paired_state(n, true, q, q)\n",
    "        if dq==1\n",
    "            part2 += kron(tmp1, tmp1)\n",
    "        else\n",
    "            part2 += (kron(tmp1, tmp1) + kron(tmp2, tmp2)) / (dq^2 - 1)\n",
    "            part2 -= (kron(tmp1, tmp2) + kron(tmp2, tmp1)) / ((dq^2 - 1) * dq)\n",
    "        end\n",
    "    end\n",
    "    result = part1 + part2\n",
    "    result = reshape(result, d^4, d^4)\n",
    "    return result\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "9d420df6",
   "metadata": {},
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "OutOfMemoryError()",
     "output_type": "error",
     "traceback": [
      "OutOfMemoryError()",
      "",
      "Stacktrace:",
      " [1] Array at .\\boot.jl:406 [inlined]",
      " [2] Array at .\\boot.jl:415 [inlined]",
      " [3] zeros at .\\array.jl:525 [inlined]",
      " [4] zeros at .\\array.jl:521 [inlined]",
      " [5] z2gauge_2nd_frame_op() at .\\In[7]:60",
      " [6] top-level scope at In[8]:1",
      " [7] include_string(::Function, ::Module, ::String, ::String) at .\\loading.jl:1091"
     ]
    }
   ],
   "source": [
    "q = z2gauge_2nd_frame_op()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8e923b0",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.5.3",
   "language": "julia",
   "name": "julia-1.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.5.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}