{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "folder_path1 = 'data'\n",
    "folder_path2 = 'figures'\n",
    "folder_path3 = 'models'\n",
    "\n",
    "if not os.path.exists(folder_path1):\n",
    "    os.makedirs(folder_path1)\n",
    "if not os.path.exists(folder_path2):\n",
    "    os.makedirs(folder_path2)\n",
    "if not os.path.exists(folder_path3):\n",
    "    os.makedirs(folder_path3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from observable import jackknife,jackknife_stats,grab,compute_autocorrelation,compute_ess"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "L = 16\n",
    "def wilson_loop(lattice, x, y):\n",
    "    U1 = np.exp(1j * lattice[x, y, 0])\n",
    "    U2 = np.exp(1j * lattice[x, (y + 1) % L, 0])\n",
    "    U3 = np.exp(1j * lattice[(x + 1) % L, y, 1])\n",
    "    U4 = np.exp(1j * lattice[x, y, 1])\n",
    "    return np.real(U1 * U3 * np.conj(U2) * np.conj(U4))\n",
    "\n",
    "def calculate_angle(array1,xp,yp):\n",
    "    right = array1[xp % L, yp % L][0]\n",
    "    up = array1[(xp + 1) % L, yp % L][1]\n",
    "    down = array1[xp % L, yp % L][1]\n",
    "    left = array1[xp % L, (yp + 1) % L][0]\n",
    "    return right + up - left - down\n",
    "\n",
    "def project_angle(angle):\n",
    "    return np.remainder(angle+np.pi, (2 * np.pi)) - np.pi\n",
    "\n",
    "def calculate_topological_charge(array):\n",
    "    top_charge = 0\n",
    "    for i in range(L):\n",
    "        for j in range(L):\n",
    "            angle = calculate_angle(array, i, j)\n",
    "            top_charge += project_angle(angle)\n",
    "\n",
    "    return round(top_charge / (2 * np.pi))\n",
    "\n",
    "\n",
    "\n",
    "def topological_suspectiblity(array):\n",
    "    top_charge = 0\n",
    "    for i in range(L):\n",
    "        for j in range(L):\n",
    "            angle = calculate_angle(array, i, j)\n",
    "            top_charge += project_angle(angle)\n",
    "\n",
    "    return round(top_charge / (2 * np.pi))**2\n",
    "\n",
    "def plaq(lat):\n",
    "    Wilson_loop_values=[]\n",
    "    lattice_size=lat.shape[2]\n",
    "    for conf in range(lat.shape[0]):\n",
    "        wilson_sum = 0\n",
    "        lattice = lat[conf]\n",
    "        for x in range(lattice_size):\n",
    "            for y in range(lattice_size):\n",
    "                wilson_sum += wilson_loop(lattice, x, y)\n",
    "        Wilson_loop_values.append(wilson_sum / (lattice_size ** 2))\n",
    "    return Wilson_loop_values\n",
    "\n",
    "def wilson_loopij(array, x_len, t_len):\n",
    "    Wilson_loop_values=[]\n",
    "    for conf in range(array.shape[0]):\n",
    "        init = np.zeros((L,L))\n",
    "        for x in range(x_len):\n",
    "            init += np.roll(array[conf,:,:,0],-x,axis=0)\n",
    "        for t in range(t_len):\n",
    "            init += np.roll(np.roll(array[conf,:,:,1],-x_len,axis=0),-t,axis=1)\n",
    "        for x in range(x_len):\n",
    "            init -= np.roll(np.roll(array[conf,:,:,0],(-x_len,-t_len),axis=(0,1)),x+1,axis=0)   \n",
    "        for t in range(t_len):\n",
    "            init -= np.roll(np.roll(array[conf,:,:,1],-t_len,axis=1),t+1,axis=1)\n",
    "        wil_l = np.mean(np.cos(init),axis=(0,1))\n",
    "        Wilson_loop_values.append(wil_l)\n",
    "    return Wilson_loop_values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1 HMC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from tqdm import tqdm\n",
    "import math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mc import get_action, get_drift, get_hamiltonian, hmc,gauge_cooling\n",
    "beta = 1\n",
    "chains = 10240*3\n",
    "nk = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "load from file\n"
     ]
    }
   ],
   "source": [
    "import copy\n",
    "from concurrent.futures import ProcessPoolExecutor\n",
    "from tqdm import tqdm \n",
    "\n",
    "filename = 'data/cfgs_L{}_beta{}_{}k.npy'.format(L,beta,nk)\n",
    "\n",
    "def runhmc(chains):\n",
    "    local_cfgs = []\n",
    "    local_acc = []\n",
    "    for chain_idx in tqdm(range(chains), desc=\"Chain Progress\"):\n",
    "        phi = np.zeros((L,L,2))\n",
    "        S = get_action(phi, beta)\n",
    "        for i in range(all_steps):\n",
    "            phi, S, accepted = hmc(phi, S, beta)\n",
    "            local_acc.append(accepted)\n",
    "            phi = gauge_cooling(phi)\n",
    "            if i % eq_step == 0 and i > therm_step:\n",
    "                local_cfgs.append(copy.deepcopy(phi))\n",
    "    return local_cfgs\n",
    "\n",
    "if os.path.exists(filename):\n",
    "    print('load from file') \n",
    "    cfgs = np.load(filename,allow_pickle=True)\n",
    "\n",
    "else:\n",
    "    \n",
    "    print('generating new cfgs')\n",
    "    therm_step = 200\n",
    "    eq_step = 128\n",
    "    all_steps = (therm_step + eq_step)\n",
    "    \n",
    "    # prepare training data-set\n",
    "    chainlist = list(chains for i in range(nk))\n",
    "    max_workers = 0  # Set this variable to control the number of workers\n",
    "\n",
    "    if max_workers > 1:\n",
    "        # Use multiprocessing when more than 1 worker is specified\n",
    "        with ProcessPoolExecutor(max_workers=max_workers) as executor:\n",
    "            results = list(executor.map(runhmc, chainlist))\n",
    "    else:\n",
    "        \n",
    "        # Execute serially without multiprocessing\n",
    "        results = [runhmc(chains) for i in range(nk)]\n",
    "        \n",
    "\n",
    "    # Accumulating and reshaping results\n",
    "    all_cfgs = [cfg for sublist in results for cfg in sublist]\n",
    "    cfgs = np.array(all_cfgs).reshape(-1, L, L, 2)\n",
    "    print(cfgs.shape)\n",
    "\n",
    "    np.save(filename, cfgs)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "cfgs=gauge_cooling(cfgs[:10240*3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gauge_trans(lattice,num):\n",
    "    gauge = (2 * np.random.rand(num, L, L) - 1)*np.pi\n",
    "    transed = np.random.rand(num, L, L,2)\n",
    "    transed[:,:,:,0] = lattice[:num,:,:,0] + gauge - np.roll(gauge,-1,axis=1)\n",
    "    transed[:,:,:,1] = lattice[:num,:,:,1] + gauge - np.roll(gauge,-1,axis=2)\n",
    "    return transed\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "cool_steps = cfgs.shape[0]\n",
    "temp = copy.deepcopy(cfgs)\n",
    "\n",
    "transfered_lattices = []\n",
    "en_times = 4\n",
    "for _ in range(en_times):\n",
    "    temp = gauge_trans(temp, cool_steps)\n",
    "    cooled_temp = gauge_cooling(temp)\n",
    "    transfered_lattices.append(cooled_temp)\n",
    "\n",
    "\n",
    "cfgs = gauge_cooling(cfgs)\n",
    "\n",
    "cfgs = np.concatenate((cfgs, *transfered_lattices), axis=0)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Diffusion Models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1 preparations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using device: cuda\tNVIDIA GeForce RTX 4080 Laptop GPU\n"
     ]
    }
   ],
   "source": [
    "import torch\n",
    "import torch.nn as nn\n",
    "import torch.nn.functional as F\n",
    "import numpy as np\n",
    "\n",
    "if torch.cuda.is_available():\n",
    "\n",
    "  device = 'cuda'\n",
    "  float_dtype = np.float32 # single\n",
    "  torch.set_default_device(device)\n",
    "else:\n",
    "  device = 'cpu'\n",
    "  float_dtype = np.float64 # double\n",
    "  torch.set_default_device(device)\n",
    "print(f\"Using device: {device}\\t\" + (f\"{torch.cuda.get_device_name(0)}\" if torch.cuda.is_available() else \"CPU\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.1.1 normalization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "torch.Size([153600, 2, 16, 16])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from norm import norm, renorm\n",
    "\n",
    "cfgs_nm,cfgs_min,cfgs_max = norm(cfgs)\n",
    "\n",
    "nm = 1\n",
    "\n",
    "configs = torch.from_numpy(cfgs_nm).float().permute(0, 3, 1, 2).to(device) # torch tensor\n",
    "configs.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.1.2 data set "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "from torch.utils.data import Dataset, DataLoader,TensorDataset\n",
    "\n",
    "n_epochs =   250#@param {'type':'integer'}\n",
    "## size of a mini-batch\n",
    "batch_size =  256 #@param {'type':'integer'}\n",
    "## learning rate\n",
    "lr=1e-3 #@param {'type':'number'}\n",
    "\n",
    "num_train = (1+en_times)*cool_steps\n",
    "\n",
    "dataset = TensorDataset(configs[:num_train], configs[:num_train])\n",
    "data_loader = DataLoader(dataset, batch_size = batch_size, shuffle=True, generator=torch.Generator(device))#"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.1.3 stochastic differential equations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "import functools\n",
    "from sde import marginal_prob_std, diffusion_coeff\n",
    "\n",
    "sigma =  25.0#@param {'type':'number'} default = 25\n",
    "marginal_prob_std_fn = functools.partial(marginal_prob_std, sigma=sigma, device=device)\n",
    "diffusion_coeff_fn = functools.partial(diffusion_coeff, sigma=sigma, device=device)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Building ScoreNet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "import tqdm\n",
    "from torch.optim import Adam\n",
    "from loss import loss_fn\n",
    "import torchvision.transforms as transforms\n",
    "from unet import GaussianFourierProjection,Dense,ScoreNet\n",
    "\n",
    "score_model = torch.nn.DataParallel(ScoreNet(marginal_prob_std=marginal_prob_std_fn))\n",
    "score_model = score_model.to(device)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3 Training"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"\\n#@title Training (double click to expand or collapse)\\nn_epochs = 250\\noptimizer = Adam(score_model.parameters(), lr=lr)\\ntqdm_epoch = tqdm.trange(n_epochs)\\nfor epoch in tqdm_epoch:\\n  avg_loss = 0.\\n  num_items = 0\\n  for x, y in data_loader:\\n    x = x.to(device)   \\n    loss = loss_fn(score_model, x, marginal_prob_std_fn)\\n    optimizer.zero_grad()\\n    loss.backward()\\n    optimizer.step()\\n    avg_loss += loss.item() * x.shape[0]\\n    num_items += x.shape[0]\\n  # Print the averaged training loss so far.\\n  tqdm_epoch.set_description('Average Loss: {:5f}'.format(avg_loss / num_items))\\n  # Update the checkpoint after each epoch of training.\\n  if epoch>230:\\n    torch.save(score_model.state_dict(), 'models/noise_ckpt_L{}_beta{}_num{}_nm{}_epoch{}.pth'.format(L,beta,num_train,nm,epoch))\\n\""
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "#@title Training (double click to expand or collapse)\n",
    "n_epochs = 250\n",
    "optimizer = Adam(score_model.parameters(), lr=lr)\n",
    "tqdm_epoch = tqdm.trange(n_epochs)\n",
    "for epoch in tqdm_epoch:\n",
    "  avg_loss = 0.\n",
    "  num_items = 0\n",
    "  for x, y in data_loader:\n",
    "    x = x.to(device)   \n",
    "    loss = loss_fn(score_model, x, marginal_prob_std_fn)\n",
    "    optimizer.zero_grad()\n",
    "    loss.backward()\n",
    "    optimizer.step()\n",
    "    avg_loss += loss.item() * x.shape[0]\n",
    "    num_items += x.shape[0]\n",
    "  # Print the averaged training loss so far.\n",
    "  tqdm_epoch.set_description('Average Loss: {:5f}'.format(avg_loss / num_items))\n",
    "  # Update the checkpoint after each epoch of training.\n",
    "  if epoch>230:\n",
    "    torch.save(score_model.state_dict(), 'models/noise_ckpt_L{}_beta{}_num{}_nm{}_epoch{}.pth'.format(L,beta,num_train,nm,epoch))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.4 Sampling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<All keys matched successfully>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model_beta = 1\n",
    "model_L = 16\n",
    "ckpt_1 = torch.load(\n",
    "    'models/noise_ckpt_L{}_beta{}_num153600_nm1_epoch249.pth'.format(model_L, model_beta),\n",
    "    map_location=device,\n",
    "    weights_only=True \n",
    ")\n",
    "\n",
    "score_model.load_state_dict(ckpt_1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\ProgramData\\anaconda3\\Lib\\site-packages\\torch\\utils\\_device.py:79: UserWarning: To copy construct from a tensor, it is recommended to use sourceTensor.clone().detach() or sourceTensor.clone().detach().requires_grad_(True), rather than torch.tensor(sourceTensor).\n",
      "  return func(*args, **kwargs)\n",
      "100%|██████████| 30/30 [26:06<00:00, 52.22s/it]\n"
     ]
    }
   ],
   "source": [
    "num_sample = 1024 \n",
    "\n",
    "from sampler import MALA_sampler\n",
    "sampler = MALA_sampler\n",
    "alpha=5e-9\n",
    "ratio=0.9\n",
    "t_0=1\n",
    "size=400\n",
    "num_steps = 100\n",
    "target_L = 16\n",
    "target_beta = 7\n",
    "#samples, accept rate, accept rate at each noise scale, action list, topological charge list\n",
    "samples,rate,rate_traj, s_list,Q_list= sampler(score_model, \n",
    "                marginal_prob_std_fn,\n",
    "                diffusion_coeff_fn, \n",
    "                batch_size = num_sample, \n",
    "                num_steps = num_steps, \n",
    "                device=device,\n",
    "                alpha=alpha,ratio=ratio,mh_steps=400,t_mh=0.001,\n",
    "                L =16,start=t_0,size=size,beta=target_beta/model_beta)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.5 Evaluation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_sample=np.array(samples)\n",
    "samples=all_sample.reshape(-1,2,L,L)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1024, 2, 16, 16)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "all_sample.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "samples_re = renorm(samples,cfgs_min,cfgs_max)\n",
    "cfgs_df = samples_re\n",
    "cfgs_test = np.transpose(cfgs_df,(0, 2, 3, 1))\n",
    "cfgs = cfgs[:1024]\n",
    "cfgs_plaq=wilson_loopij(cfgs,1,1)\n",
    "cfgs_test_plaq = wilson_loopij(cfgs_test,1,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 122536 characters]",
      "text/plain": [
       "<Figure size 1600x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(16, 6))\n",
    "\n",
    "plt.subplot(1, 2, 1)\n",
    "\n",
    "plt.plot(np.arange(cfgs_test.shape[0]), cfgs_test_plaq)\n",
    "plt.title(\"Average 1x1 Wilson Loop vs MC Steps\")\n",
    "plt.xlabel(\"MC Steps\")\n",
    "plt.ylabel(\"Average 1x1 Wilson Loop\")\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plaq Hist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 65220 characters]",
      "text/plain": [
       "<Figure size 800x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import scienceplots\n",
    "from matplotlib.ticker import MultipleLocator\n",
    "plt.style.use(['science', 'notebook','grid'])\n",
    "plt.rcParams['figure.dpi'] = 100\n",
    "plt.rcParams['savefig.dpi'] = 300\n",
    "plt.rcParams['font.size'] = 16\n",
    "plt.rcParams['xtick.labelsize'] = 14  # 设置横坐标字体大小\n",
    "plt.rcParams['ytick.labelsize'] = 14  # 设置纵坐标字体大小\n",
    "\n",
    "plt.figure(figsize=(8, 6))\n",
    "\n",
    "\n",
    "\n",
    "Xrange=(0.85,1)\n",
    "plt.hist(cfgs_plaq, bins=40,range=Xrange, histtype='step',linewidth=2,edgecolor='blue',color='blue',alpha=1,density=True,label=r'HMC $\\beta = 7$')\n",
    "plt.hist(cfgs_test_plaq, bins=40,range=Xrange, edgecolor='black',color='orange',alpha=1,density=True, label=r'DM $\\beta = 7$')\n",
    "\n",
    "ax = plt.gca()  \n",
    "ax.yaxis.set_major_locator(MultipleLocator(5)) \n",
    "ax.yaxis.set_minor_locator(MultipleLocator(2.5)) \n",
    "\n",
    "\n",
    "plt.legend()\n",
    "\n",
    "plt.xlabel(r'$\\langle W_{1 \\times 1} \\rangle$', fontsize=16)\n",
    "plt.ylabel(r'Frequency density', fontsize=16)\n",
    "plt.savefig('beta_7_Wilson_loop_step.pdf', format='pdf', bbox_inches='tight', pad_inches=0.1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Wilson loops"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.92519704427842975 0.00620574976761201\n",
      "0.73706997364954174 0.03379398011476174\n",
      "0.49554813620652594 0.07255679981067424\n",
      "0.28336750812313477 0.10592245818325528\n"
     ]
    }
   ],
   "source": [
    "for i in range(1,5):\n",
    "    print(np.mean(wilson_loopij(cfgs_test[:],i,i)),np.std(wilson_loopij(cfgs_test[:],i,i)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9258483702727198 0.00652802527457644\n",
      "0.7384604609967504 0.03066239752436074\n",
      "0.5132801244094619 0.06622602082429403\n",
      "0.3158275714518931 0.09397229272901812\n"
     ]
    }
   ],
   "source": [
    "for i in range(1,5):\n",
    "    print(np.mean(wilson_loopij(cfgs[:],i,i)),np.std(wilson_loopij(cfgs[:],i,i)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Topological Charge Hist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "top_charge_history=[]\n",
    "top_charge_history2=[]\n",
    "for conf in range(cfgs_test.shape[0]):\n",
    "    top_charge_history.append(calculate_topological_charge(cfgs_test[conf]))\n",
    "    top_charge_history2.append(calculate_topological_charge(cfgs[conf]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 51524 characters]",
      "text/plain": [
       "<Figure size 800x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scienceplots\n",
    "\n",
    "\n",
    "plt.style.use(['science', 'notebook','grid'])\n",
    "plt.rcParams['figure.dpi'] = 100\n",
    "plt.rcParams['savefig.dpi'] = 300\n",
    "plt.rcParams['font.size'] = 16\n",
    "plt.rcParams['xtick.labelsize'] = 14\n",
    "plt.rcParams['ytick.labelsize'] = 14  \n",
    "\n",
    "fig_width = 8 \n",
    "fig_height = 6  \n",
    "fig, ax = plt.subplots(figsize=(fig_width, fig_height), dpi=100)\n",
    "\n",
    "bin_range = range(-5, 6)\n",
    "\n",
    "hist2, bins2 = np.histogram(top_charge_history2, bins=bin_range, density=True)\n",
    "hist, bins = np.histogram(top_charge_history, bins=bin_range, density=True)\n",
    "\n",
    "ax.bar(bins2[:-1] - 0.2, hist2, width=0.5, edgecolor='black', color='blue', alpha=1, label=r'HMC $\\beta=7$')\n",
    "ax.bar(bins[:-1] + 0.2, hist, width=0.5, edgecolor='black', color='orange', alpha=1, label=r'DM $\\beta=7$')\n",
    "\n",
    "ax.set_xlabel(r'Topological Charge $Q$', fontsize=16)\n",
    "ax.set_ylabel(r'Frequency density', fontsize=16)\n",
    "ax.set_xticks(range(-5, 6, 1))\n",
    "\n",
    "ax.set_xlim(-5.5, 5.5)\n",
    "#ax.set_ylim(-0.001, 0.3)\n",
    "plt.legend()\n",
    "\n",
    "\n",
    "ax = plt.gca()  # 获取当前的坐标轴\n",
    "\n",
    "\n",
    "\n",
    "plt.legend()\n",
    "\n",
    "ax = plt.gca()  \n",
    "ax.yaxis.set_major_locator(MultipleLocator(0.2)) \n",
    "ax.yaxis.set_minor_locator(MultipleLocator(0.05))  \n",
    "plt.savefig('beta_7_Q_new.pdf', format='pdf', bbox_inches='tight', pad_inches=0.1)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Topological Suspectibility"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.004505157470703125 0.0005369903717622893\n"
     ]
    }
   ],
   "source": [
    "suspect1 = []\n",
    "for conf in range(cfgs_test.shape[0]):\n",
    "    suspect1.append(topological_suspectiblity(cfgs_test[conf]))\n",
    "mean, err = jackknife_stats(suspect1)\n",
    "print(mean/L**2,err/L**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0002365112304687499 3.198772317967985e-05\n"
     ]
    }
   ],
   "source": [
    "suspect2 = []\n",
    "for conf in range(cfgs.shape[0]):\n",
    "    suspect2.append(topological_suspectiblity(cfgs[conf]))\n",
    "mean2, err2 = jackknife_stats(suspect2)\n",
    "print(mean2/L**2,err2/L**2)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
