my_vqe.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Performance Evaluation of CVaR-VQE \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ridefinire time to solution\n",
    "in lidar loro trovano l'ottimo per ogni N numero di istanze parallele\n",
    "carleo non-stoquastic QA con troyer\n",
    "\n",
    "Noi geralizziamo facendo una funzione in 2D oppure numero di shots come parametro t in modo fa fare una famiglia di funzioni\n",
    "\n",
    "Cosa vogliamo: resolurce to solution, moltiplico N_shot e N_step\n",
    "ferromagnetico bias\n",
    "mi servono alcuni for\n",
    "- random su theta\n",
    "- scalo su shots \n",
    "- scalo su step (maxiter)\n",
    "\n",
    "alpha suggerito da loro (0.25)\n",
    "\n",
    "valuto con RTS (TTS dove moltiplico per gli shots) non normalizzato per n_max\n",
    "\n",
    "n.b. green quantum computing paper montangero"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Plot for $p=1$ \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import curve_fit\n",
    "from math import ceil\n",
    "from json import JSONEncoder\n",
    "from matplotlib import pyplot as plt\n",
    "from src.utils import collect_results, NumpyArrayEncoder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qubits = np.arange(4, 22, 2, dtype=int)\n",
    "circ_depth = 1\n",
    "\n",
    "psucc, ts, shots, nfevs = collect_results(qubits, circ_depth)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i, (p,t) in enumerate(zip(psucc, ts)):\n",
    "    try:\n",
    "        last_zero = np.argwhere(np.isclose(p, 0.))\n",
    "        print(f\"N={4+i*2}, t={np.asarray(t)[np.asarray(p)==0.].max()}, shots={np.asarray(shots[i])[np.asarray(p)==0.].max()}, nfev={np.asarray(nfevs[i])[np.asarray(p)==0.].max()}\")\n",
    "    except ValueError:\n",
    "        print(f\"No zeros in succes probability for N={4+i*2}\")\n",
    "    p_t = np.vstack((t, p))\n",
    "    p_t = p_t[:, np.argsort(p_t[0, :])]\n",
    "    plt.plot(p_t[0], p_t[1], \".\", label=f\"N={4+i*2}\")\n",
    "\n",
    "    plt.xscale('log')\n",
    "    plt.xlabel(r\"$t$ (shots x iterations)\")\n",
    "    plt.ylabel(r\"$P$\")\n",
    "    plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def running_mean(x, n):\n",
    "    cumsum = np.cumsum(np.pad(x, int(ceil(n/2)), mode='reflect')) \n",
    "    return (cumsum[n:] - cumsum[:-n]) / float(n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i, (p,t) in enumerate(zip(psucc, ts)):\n",
    "    p_t = np.vstack((t, p))\n",
    "    p_t = p_t[:, np.argsort(p_t[0, :])]\n",
    "    len_window = 12\n",
    "    #p_t[1] = np.convolve(p_t[1], np.ones(len_window) / len_window, mode='same')\n",
    "    p_t[1] = running_mean(p_t[1], len_window)\n",
    "    plt.plot(p_t[0][::5], p_t[1][::5], \"o-\", label=f\"N={4+i*2}\")\n",
    "\n",
    "    plt.xscale('log')\n",
    "    plt.xlabel(r\"$t$ (shots x iterations)\")\n",
    "    plt.ylabel(r\"$P$\")\n",
    "    plt.legend()\n",
    "    plt.savefig(\"p_t_n4-18.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for k, (p,s,nfev) in enumerate(zip(psucc, shots, nfevs)):\n",
    "    len_iter = np.max(nfev)//2+1\n",
    "    p = np.reshape(p, (-1,len_iter))\n",
    "    s = np.reshape(s, (-1,len_iter))\n",
    "    nfev = np.reshape(nfev, (-1,len_iter))\n",
    "\n",
    "    fig, (ax1, ax2) = plt.subplots(1,2, figsize=(16,6))\n",
    "    for i in range(p.shape[0]):\n",
    "        ax1.plot(nfev[i,:], p[i,:], '--o', label=s[i,0])\n",
    "    for j in range(p.shape[1]):\n",
    "        ax2.plot(s[:,j], p[:,j], '--s', label=nfev[0,j])\n",
    "\n",
    "    ax1.set_xlabel('Iterations')\n",
    "    ax1.set_ylabel(r'$p$(success)')\n",
    "    #ax1.set_ylim(0,1.1)\n",
    "\n",
    "    ax2.set_xlabel('Shots')\n",
    "    ax2.set_xscale('log', base=2)\n",
    "    ax2.set_ylabel(r'$p$(success)')\n",
    "    #ax2.set_ylim(0,1.1)\n",
    "\n",
    "    ax1.legend(title='Shots')\n",
    "    ax2.legend(title='Iterations')\n",
    "\n",
    "    ax1.set_title(f'N={4+k*2}')\n",
    "    ax2.set_title(f'N={4+k*2}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_r(p, p_d=0.99):\n",
    "    rt = np.log(1-p_d)\n",
    "    rt /= np.log(1-p)\n",
    "    return rt\n",
    "\n",
    "tts_list = []\n",
    "for p,t in zip(psucc, ts):\n",
    "    p = np.asarray(p)\n",
    "    t = np.asarray(t)\n",
    "    #print(p[:20],t.shape)\n",
    "    mask =  p == 0.\n",
    "    mask += p == 1. # mask += p == 1.\n",
    "    #p[mask] = 0.00005\n",
    "    p = np.delete(p, mask)\n",
    "    t = np.delete(t, mask)\n",
    "\n",
    "    rt = compute_r(p)\n",
    "    tts_list.append(t*rt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i, (tts,t) in enumerate(zip(tts_list, ts)):\n",
    "    if 4+i*2 > 18:\n",
    "        pass#continue\n",
    "    t = np.asarray(t)\n",
    "    mask = np.asarray(psucc[i]) == 0.\n",
    "    mask += np.asarray(psucc[i]) == 1. # mask += np.asarray(psucc[i]) == 1.\n",
    "    t = np.delete(t, mask)\n",
    "\n",
    "    tts_t = np.vstack((t, tts))\n",
    "    tts_t = tts_t[:, np.argsort(tts_t[0, :])]\n",
    "    # roll-mean\n",
    "    len_window = 2\n",
    "    tts_t[1] = running_mean(tts_t[1], len_window)\n",
    "    plt.plot(tts_t[0], tts_t[1], \".\", label=f\"N={4+i*2}\")\n",
    "\n",
    "    plt.xscale('log')\n",
    "    plt.yscale('log')\n",
    "    #plt.xlim(1.5,50000)\n",
    "    #plt.ylim(10000, 500000)\n",
    "    plt.xlabel(r\"$t$ (shots x iterations)\")\n",
    "    plt.ylabel(r\"$STS$\")\n",
    "    plt.legend()\n",
    "    plt.savefig(\"STS-t_roll8.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fun_exp(x, a, b):\n",
    "    return a*2**(b*x)\n",
    "\n",
    "def fun_lin(x, a, b):\n",
    "    return a + b*x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pt_all = []\n",
    "for i, (p,t) in enumerate(zip(psucc, ts)):\n",
    "    p_t = np.vstack((t, p))\n",
    "    p_t = p_t[:, np.argsort(p_t[0, :])]\n",
    "    pt_all.append(p_t)\n",
    "\n",
    "colors = ['lightgreen', 'limegreen', 'forestgreen', 'darkgreen','lightseagreen', 'cornflowerblue' , 'blue', 'navy']\n",
    "min_psucc_lst = [0.99, 0.9, 0.75, 0.5, 0.25]\n",
    "min_qubits = 4 \n",
    "max_qubits = 20 \n",
    "for min_psucc, color in zip(min_psucc_lst, colors):\n",
    "    step_size = []\n",
    "    for i, pt in enumerate(pt_all):\n",
    "        mask = pt[1] >= min_psucc\n",
    "        # check if the p is not reached\n",
    "        if mask.sum() == 0:\n",
    "            continue\n",
    "        step_size.append(pt[0][mask][0])\n",
    "\n",
    "    x_axes = 4 + 2*np.arange(len(step_size))\n",
    "\n",
    "    popt, _ = curve_fit(fun_lin, x_axes, np.log2(step_size))\n",
    "    plt.plot(x_axes, step_size, \"o\", label=f\"p={min_psucc}\", color=color)\n",
    "    plt.plot(np.arange(min_qubits, max_qubits+1), fun_exp(np.arange(min_qubits, max_qubits+1), 2**popt[0], popt[1]), \"--\", label=f\"a={2**popt[0]:.2f}, b={popt[1]:.2f}\", color=color)\n",
    "\n",
    "plt.yscale('log', base=2)\n",
    "plt.ylabel(r\"$t$ (shots x iterations)\")\n",
    "plt.xlabel(r\"$N$\")\n",
    "\n",
    "plt.legend(loc=\"best\")\n",
    "plt.title(r\"Fit with $a \\cdot 2^{b \\cdot N}$\")\n",
    "\n",
    "plt.savefig(\"t_succ.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot for $p=2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import curve_fit\n",
    "from math import ceil\n",
    "from json import JSONEncoder\n",
    "from matplotlib import pyplot as plt\n",
    "from src.utils import collect_results, NumpyArrayEncoder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qubits = np.arange(4, 22, 2, dtype=int)\n",
    "circ_depth = 2\n",
    "\n",
    "psucc, ts, shots, nfevs = collect_results(qubits, circ_depth)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i, (p,t) in enumerate(zip(psucc, ts)):\n",
    "    try:\n",
    "        last_zero = np.argwhere(np.isclose(p, 0.))\n",
    "        print(f\"N={4+i*2}, t={np.asarray(t)[np.asarray(p)==0.].max()}, shots={np.asarray(shots[i])[np.asarray(p)==0.].max()}, nfev={np.asarray(nfevs[i])[np.asarray(p)==0.].max()}\")\n",
    "    except ValueError:\n",
    "        print(f\"No zeros in succes probability for N={4+i*2}\")\n",
    "    p_t = np.vstack((t, p))\n",
    "    p_t = p_t[:, np.argsort(p_t[0, :])]\n",
    "    plt.plot(p_t[0], p_t[1], \".\", label=f\"N={4+i*2}\")\n",
    "\n",
    "    plt.xscale('log')\n",
    "    plt.xlabel(r\"$t$ (shots x iterations)\")\n",
    "    plt.ylabel(r\"$P$\")\n",
    "    plt.legend()\n",
    "    #plt.savefig(\"p_t_n4-18.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def running_mean(x, n):\n",
    "    cumsum = np.cumsum(np.pad(x, int(ceil(n/2)), mode='reflect')) \n",
    "    return (cumsum[n:] - cumsum[:-n]) / float(n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i, (p,t) in enumerate(zip(psucc, ts)):\n",
    "    p_t = np.vstack((t, p))\n",
    "    p_t = p_t[:, np.argsort(p_t[0, :])]\n",
    "    len_window = 12\n",
    "    #p_t[1] = np.convolve(p_t[1], np.ones(len_window) / len_window, mode='same')\n",
    "    p_t[1] = running_mean(p_t[1], len_window)\n",
    "    plt.plot(p_t[0][::5], p_t[1][::5], \"o-\", label=f\"N={4+i*2}\")\n",
    "\n",
    "    plt.xscale('log')\n",
    "    plt.xlabel(r\"$t$ (shots x iterations)\")\n",
    "    plt.ylabel(r\"$P$\")\n",
    "    plt.legend()\n",
    "    plt.savefig(\"p_t_n4-18.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for k, (p,s,nfev) in enumerate(zip(psucc, shots, nfevs)):\n",
    "    len_iter = np.max(nfev)//2+1\n",
    "\n",
    "    p = np.reshape(p, (-1,len_iter))\n",
    "    s = np.reshape(s, (-1,len_iter))\n",
    "    nfev = np.reshape(nfev, (-1,len_iter))\n",
    "\n",
    "    fig, (ax1, ax2) = plt.subplots(1,2, figsize=(16,6))\n",
    "    for i in range(p.shape[0]):\n",
    "        ax1.plot(nfev[i,:], p[i,:], '--o', label=s[i,0])\n",
    "    for j in range(p.shape[1]):\n",
    "        ax2.plot(s[:,j], p[:,j], '--s', label=nfev[0,j])\n",
    "\n",
    "    ax1.set_xlabel('Iterations')\n",
    "    ax1.set_ylabel(r'$p$(success)')\n",
    "    #ax1.set_ylim(0,1.1)\n",
    "\n",
    "    ax2.set_xlabel('Shots')\n",
    "    ax2.set_xscale('log', base=2)\n",
    "    ax2.set_ylabel(r'$p$(success)')\n",
    "    #ax2.set_ylim(0,1.1)\n",
    "\n",
    "    ax1.legend(title='Shots')\n",
    "    ax2.legend(title='Iterations')\n",
    "\n",
    "    ax1.set_title(f'N={4+k*2}')\n",
    "    ax2.set_title(f'N={4+k*2}')\n",
    "\n",
    "    #plt.savefig(f'shots_iter_N{4+k*2}.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_r(p, p_d=0.99):\n",
    "    rt = np.log(1-p_d)\n",
    "    rt /= np.log(1-p)\n",
    "    return rt\n",
    "\n",
    "tts_list = []\n",
    "for p,t in zip(psucc, ts):\n",
    "    p = np.asarray(p)\n",
    "    t = np.asarray(t)\n",
    "    #print(p[:20],t.shape)\n",
    "    mask =  p == 0.\n",
    "    mask += p == 1. # mask += p == 1.\n",
    "    #p[mask] = 0.00005\n",
    "    p = np.delete(p, mask)\n",
    "    t = np.delete(t, mask)\n",
    "\n",
    "    rt = compute_r(p)\n",
    "    tts_list.append(t*rt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i, (tts,t) in enumerate(zip(tts_list, ts)):\n",
    "    if 4+i*2 > 18:\n",
    "        pass#continue\n",
    "    t = np.asarray(t)\n",
    "    mask = np.asarray(psucc[i]) == 0.\n",
    "    mask += np.asarray(psucc[i]) == 1. # mask += np.asarray(psucc[i]) == 1.\n",
    "    t = np.delete(t, mask)\n",
    "\n",
    "    tts_t = np.vstack((t, tts))\n",
    "    tts_t = tts_t[:, np.argsort(tts_t[0, :])]\n",
    "    # roll-mean\n",
    "    len_window = 8\n",
    "    tts_t[1] = running_mean(tts_t[1], len_window)\n",
    "    plt.plot(tts_t[0], tts_t[1], \".\", label=f\"N={4+i*2}\")\n",
    "\n",
    "    plt.xscale('log')\n",
    "    plt.yscale('log')\n",
    "    #plt.xlim(1.5,50000)\n",
    "    #plt.ylim(10000, 500000)\n",
    "    plt.xlabel(r\"$t$ (shots x iterations)\")\n",
    "    plt.ylabel(r\"$STS$\")\n",
    "    plt.legend()\n",
    "    plt.savefig(\"STS-t_roll8.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fun_exp(x, a, b):\n",
    "    return a*2**(b*x)\n",
    "\n",
    "def fun_lin(x, a, b):\n",
    "    return a + b*x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fun_exp(x, a, b):\n",
    "    return a*2**(b*x)\n",
    "\n",
    "def fun_lin(x, a, b):\n",
    "    return a + b*x\n",
    "    \n",
    "pt_all = []\n",
    "for i, (p,t) in enumerate(zip(psucc, ts)):\n",
    "    p_t = np.vstack((t, p))\n",
    "    p_t = p_t[:, np.argsort(p_t[0, :])]\n",
    "    pt_all.append(p_t)\n",
    "\n",
    "colors = ['lightgreen', 'limegreen', 'forestgreen', 'darkgreen','lightseagreen', 'cornflowerblue' , 'blue', 'navy']\n",
    "min_psucc_lst = [0.99, 0.9, 0.75, 0.5, 0.25]\n",
    "min_qubits = 4 \n",
    "max_qubits = 20 \n",
    "for min_psucc, color in zip(min_psucc_lst, colors):\n",
    "    step_size = []\n",
    "    for i, pt in enumerate(pt_all):\n",
    "        mask = pt[1] >= min_psucc\n",
    "        # check if the p is not reached\n",
    "        if mask.sum() == 0:\n",
    "            continue\n",
    "        step_size.append(pt[0][mask][0])\n",
    "\n",
    "    x_axes = 4 + 2*np.arange(len(step_size))\n",
    "\n",
    "    popt, _ = curve_fit(fun_lin, x_axes, np.log2(step_size))\n",
    "    plt.plot(x_axes, step_size, \"o\", label=f\"p={min_psucc}\", color=color)\n",
    "    plt.plot(np.arange(min_qubits, max_qubits+1), fun_exp(np.arange(min_qubits, max_qubits+1), 2**popt[0], popt[1]), \"--\", label=f\"a={2**popt[0]:.2f}, b={popt[1]:.2f}\", color=color)\n",
    "\n",
    "plt.yscale('log', base=2)\n",
    "plt.ylabel(r\"$t$ (shots x iterations)\")\n",
    "plt.xlabel(r\"$N$\")\n",
    "\n",
    "plt.legend(loc=\"best\")\n",
    "plt.title(r\"Fit with $a \\cdot 2^{b \\cdot N}$\")\n",
    "\n",
    "plt.savefig(\"t_succ.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Simulation with $p=1,2$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qubits = np.arange(4, 22, 2, dtype=int)\n",
    "\n",
    "psucc_1, ts_1, shots_1, nfevs_1 = collect_results(qubits, 1)\n",
    "psucc_2, ts_2, shots_2, nfevs_2 = collect_results(qubits, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "colors = ['lightgreen','forestgreen','tab:green','darkgreen', 'lightseagreen', 'tab:cyan','cornflowerblue' , 'blue', 'navy']\n",
    "plots1 = []\n",
    "plots2 = []\n",
    "\n",
    "for i, (p,t) in enumerate(zip(psucc_1, ts_1)):\n",
    "    if 4+i*2 > 18:\n",
    "        continue\n",
    "    try:\n",
    "        last_zero = np.argwhere(np.isclose(p, 0.))\n",
    "        print(f\"N={4+i*2}, p=1, t={np.asarray(t)[np.asarray(p)==0.].max()}, shots={np.asarray(shots_1[i])[np.asarray(p)==0.].max()}, nfev={np.asarray(nfevs_1[i])[np.asarray(p)==0.].max()}\")\n",
    "    except ValueError:\n",
    "        print(f\"No zeros in success probability for N={4+i*2}\")\n",
    "    # plot for each dimension\n",
    "    p_t_1 = np.vstack((t, p))\n",
    "    p_t_1 = p_t_1[:, np.argsort(p_t_1[0, :])]\n",
    "    plot1, = plt.plot(p_t_1[0], p_t_1[1], \".\", label=f\"{4+i*2}\", color=colors[i])\n",
    "    plots1.append(plot1)\n",
    "\n",
    "    try:\n",
    "        last_zero = np.argwhere(np.isclose(psucc_2[i], 0.))\n",
    "        print(f\"N={4+i*2}, p=2, t={np.asarray(ts_2[i])[np.asarray(psucc_2[i])==0.].max()}, shots={np.asarray(shots_2[i])[np.asarray(psucc_2[i])==0.].max()}, nfev={np.asarray(nfevs_2[i])[np.asarray(psucc_2[i])==0.].max()}\")\n",
    "    except ValueError:\n",
    "        print(f\"No zeros in succes probability for N={4+i*2}\")\n",
    "    p_t_2 = np.vstack((ts_2[i], psucc_2[i]))\n",
    "    p_t_2 = p_t_2[:, np.argsort(p_t_2[0, :])]\n",
    "    plot2, = plt.plot(p_t_2[0]*2, p_t_2[1], \"o\", label=f\"{4+i*2}\", color=colors[i], fillstyle='none')\n",
    "    plots2.append(plot2)\n",
    "\n",
    "    plt.xscale('log')\n",
    "    plt.xlabel(r\"$t$ (shots x iterations x depth)\")\n",
    "    plt.ylabel(r\"$P$\")\n",
    "\n",
    "# Create a legend for the first line.\n",
    "first_legend = plt.legend(handles=plots1, loc='upper right')\n",
    "# Add the legend manually to the current Axes.\n",
    "ax = plt.gca().add_artist(first_legend)\n",
    "\n",
    "# Create another legend for the second line.\n",
    "plt.legend(handles=plots2, loc='upper left')\n",
    "\n",
    "    #plt.savefig(\"p_t_n4-18.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "colors = ['lightgreen','forestgreen','tab:green','darkgreen', 'lightseagreen', 'tab:cyan','cornflowerblue' , 'blue', 'navy']\n",
    "plots1 = []\n",
    "plots2 = []\n",
    "\n",
    "for i, (p,t) in enumerate(zip(psucc_1, ts_1)):\n",
    "    if 4+i*2 > 18:\n",
    "        continue\n",
    "    p_t1 = np.vstack((t, p))\n",
    "    p_t1 = p_t1[:, np.argsort(p_t1[0, :])]\n",
    "    len_window = 12\n",
    "    #p_t1[1] = np.convolve(p_t1[1], np.ones(len_window) / len_window, mode='same')\n",
    "    p_t1[1] = running_mean(p_t1[1], len_window)\n",
    "    plot1, = plt.plot(p_t1[0][::5], p_t1[1][::5], \"o-\", label=f\"{4+i*2}\", color=colors[i])\n",
    "    plots1.append(plot1)\n",
    "\n",
    "    p_t2 = np.vstack((ts_2[i], psucc_2[i]))\n",
    "    p_t2 = p_t2[:, np.argsort(p_t2[0, :])]\n",
    "    p_t2[1] = running_mean(p_t2[1], len_window)\n",
    "    plot2, = plt.plot(p_t2[0][::5]*2, p_t2[1][::5], \"o-\", label=f\"{4+i*2}\", fillstyle='none', color=colors[i])\n",
    "    plots2.append(plot2)\n",
    "\n",
    "    plt.xscale('log')\n",
    "    plt.xlabel(r\"$t$ (shots x iterations x depth)\")\n",
    "    plt.ylabel(r\"$P$\")\n",
    "\n",
    "# Create a legend for the first line.\n",
    "first_legend = plt.legend(handles=plots1, loc='upper right', title=r'$p=1$')\n",
    "# Add the legend manually to the current Axes.\n",
    "ax = plt.gca().add_artist(first_legend)\n",
    "\n",
    "# Create another legend for the second line.\n",
    "plt.legend(handles=plots2, loc='upper left', title=r'$p=2$')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fun_exp(x, a, b):\n",
    "    return a*2**(b*x)\n",
    "\n",
    "def fun_lin(x, a, b):\n",
    "    return a + b*x\n",
    "    \n",
    "pt_all = []\n",
    "for i, (p,t) in enumerate(zip(psucc_1, ts_1)):\n",
    "    p_t = np.vstack((t, p))\n",
    "    p_t = p_t[:, np.argsort(p_t[0, :])]\n",
    "    pt_all.append(p_t)\n",
    "pt_all_2 = []\n",
    "for i, (p,t) in enumerate(zip(psucc_2, ts_2)):\n",
    "    p_t = np.vstack((t, p))\n",
    "    p_t = p_t[:, np.argsort(p_t[0, :])]\n",
    "    pt_all_2.append(p_t)\n",
    "\n",
    "\n",
    "colors = ['lightgreen', 'limegreen', 'forestgreen', 'darkgreen','lightseagreen', 'cornflowerblue' , 'blue', 'navy']\n",
    "min_psucc_lst = [0.99, 0.9, 0.75, 0.5, 0.25]\n",
    "min_qubits = 4 \n",
    "max_qubits = 20 \n",
    "for min_psucc, color in zip(min_psucc_lst, colors):\n",
    "    step_size = []\n",
    "    for i, pt in enumerate(pt_all):\n",
    "        mask = pt[1] >= min_psucc\n",
    "        # check if the p is not reached\n",
    "        if mask.sum() == 0:\n",
    "            continue\n",
    "        step_size.append(pt[0][mask][0])\n",
    "\n",
    "    x_axes = 4 + 2*np.arange(len(step_size))\n",
    "\n",
    "    popt, _ = curve_fit(fun_lin, x_axes, np.log2(step_size))\n",
    "    plt.plot(x_axes, step_size, \"o\", label=f\"p={min_psucc}\", color=color)\n",
    "    plt.plot(np.arange(min_qubits, max_qubits+1), fun_exp(np.arange(min_qubits, max_qubits+1), 2**popt[0], popt[1]), \"--\", label=f\"a={2**popt[0]:.2f}, b={popt[1]:.2f}\", color=color)\n",
    "\n",
    "    step_size_2 = []\n",
    "    for i, pt_2 in enumerate(pt_all_2):\n",
    "        mask = pt_2[1] >= min_psucc\n",
    "        # check if the p is not reached\n",
    "        if mask.sum() == 0:\n",
    "            continue\n",
    "        step_size_2.append(pt_2[0][mask][0])\n",
    "    x_axes_2 = 4 + 2*np.arange(len(step_size_2))\n",
    "    popt, _ = curve_fit(fun_lin, x_axes_2, np.log2(step_size_2))\n",
    "    plt.plot(x_axes_2, step_size_2, \"o\", label=f\"p={min_psucc}\", color=color, fillstyle='none')\n",
    "    plt.plot(np.arange(min_qubits, max_qubits+1), fun_exp(np.arange(min_qubits, max_qubits+1), 2**popt[0], popt[1]), \"-\", label=f\"a={2**popt[0]:.2f}, b={popt[1]:.2f}\", color=color)\n",
    "\n",
    "plt.yscale('log', base=2)\n",
    "plt.ylabel(r\"$t$ (shots x iterations)\")\n",
    "plt.xlabel(r\"$N$\")\n",
    "\n",
    "#plt.legend(loc=\"best\")\n",
    "plt.title(r\"Fit with $a \\cdot 2^{b \\cdot N}$\")\n",
    "\n",
    "#plt.savefig(\"t_succ.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Simulation with noise $p=2$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import curve_fit\n",
    "from math import ceil\n",
    "from json import JSONEncoder\n",
    "from matplotlib import pyplot as plt\n",
    "from src.utils import collect_results, NumpyArrayEncoder\n",
    "import json\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qubits = [4,6,8,10,12]\n",
    "p = 2\n",
    "# load exact simulations\n",
    "psucc, ts, shots, nfevs = collect_results(qubits, p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ts_noise = []\n",
    "shots_noise = []\n",
    "nfevs_noise = []\n",
    "psucc_noise = []\n",
    "for qubit in qubits:\n",
    "    dir_path = f\"results/N{qubit}/noise/\"\n",
    "    print(f\"directory: {dir_path}\")\n",
    "    p = []\n",
    "    t = []\n",
    "    s = []\n",
    "    it = []\n",
    "    # for each number of qubits\n",
    "    # we have several number of shots and iterations\n",
    "    for filename in sorted(os.listdir(dir_path)):\n",
    "        filename = dir_path + filename\n",
    "        with open(filename, \"r\") as file:\n",
    "            # print(filename)\n",
    "            data = json.load(file)\n",
    "        ever_found = []\n",
    "        # for each shot and iteration param\n",
    "        # we randomized the initial point\n",
    "        # to estimate the right probability\n",
    "        for run in data:\n",
    "            ever_found.append(run[\"ever_found\"])\n",
    "        # compute p('found minimum')\n",
    "        p.append(\n",
    "            np.mean(np.asarray(ever_found, dtype=np.float128))\n",
    "        )  # float128 to avoid to many zeros\n",
    "        # maxiter*shots = actual number of iteration\n",
    "        t.append(run[\"shots\"] * run[\"nfev\"])\n",
    "        s.append(run[\"shots\"])\n",
    "        it.append(run[\"nfev\"])\n",
    "    # update list for each number of qubits\n",
    "    ts_noise.append(t)\n",
    "    psucc_noise.append(p)\n",
    "    shots_noise.append(s)\n",
    "    nfevs_noise.append(it)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "colors = ['lightgreen','forestgreen','tab:green','darkgreen', 'lightseagreen', 'tab:cyan','cornflowerblue' , 'blue', 'navy']\n",
    "plots1 = []\n",
    "plots2 = []\n",
    "\n",
    "for i, (p,t) in enumerate(zip(psucc, ts)):\n",
    "    if 4+i*2 > 18:\n",
    "        continue\n",
    "    try:\n",
    "        last_zero = np.argwhere(np.isclose(p, 0.))\n",
    "        print(f\"N={4+i*2}, p=2, t={np.asarray(t)[np.asarray(p)==0.].max()}, shots={np.asarray(shots[i])[np.asarray(p)==0.].max()}, nfev={np.asarray(nfevs[i])[np.asarray(p)==0.].max()}\")\n",
    "    except ValueError:\n",
    "        print(f\"No zeros in success probability for N={4+i*2}\")\n",
    "    # plot for each dimension\n",
    "    p_t = np.vstack((t, p))\n",
    "    p_t = p_t[:, np.argsort(p_t[0, :])]\n",
    "    plot1, = plt.plot(p_t[0]*2, p_t[1], \".\", label=f\"{4+i*2}\", color=colors[i], alpha=0.8)\n",
    "    plots1.append(plot1)\n",
    "\n",
    "for i, (p,t) in enumerate(zip(psucc_noise, ts_noise)):\n",
    "    try:\n",
    "        last_zero = np.argwhere(np.isclose(psucc_noise[i], 0.))\n",
    "        print(f\"N={4+i*2}, p=2, t={np.asarray(ts_noise[i])[np.asarray(psucc_noise[i])==0.].max()}, shots={np.asarray(shots_noise[i])[np.asarray(psucc_noise[i])==0.].max()}, nfev={np.asarray(nfevs_noise[i])[np.asarray(psucc_noise[i])==0.].max()}\")\n",
    "    except ValueError:\n",
    "        print(f\"No zeros in succes probability for N={4+i*2}\")\n",
    "    p_t_noise = np.vstack((ts_noise[i], psucc_noise[i]))\n",
    "    p_t_noise = p_t_noise[:, np.argsort(p_t_noise[0, :])]\n",
    "    plot2, = plt.plot(p_t_noise[0]*2, p_t_noise[1], \".\", label=f\"{4+i*2}\", color=colors[i], fillstyle='none')\n",
    "    plots2.append(plot2)\n",
    "\n",
    "    plt.xscale('log')\n",
    "    plt.xlabel(r\"$t$ (shots x iterations x depth)\")\n",
    "    plt.ylabel(r\"$P$\")\n",
    "\n",
    "# Create a legend for the first line.\n",
    "first_legend = plt.legend(handles=plots1, loc='lower right', title=\"Exact\")\n",
    "#Add the legend manually to the current Axes.\n",
    "ax = plt.gca().add_artist(first_legend)\n",
    "\n",
    "plt.title(r\"$p=2$\")\n",
    "#Create another legend for the second line.\n",
    "plt.legend(handles=plots2, loc='upper left', title=f\"Noise\\nT1=50\\u03bcs\\nT2=70\\u03bcs\")\n",
    "\n",
    "plt.savefig(\"noise_exact.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "colors = ['lightgreen','tab:green','darkgreen', 'lightseagreen','cornflowerblue' , 'blue', 'navy']\n",
    "plots1 = []\n",
    "plots2 = []\n",
    "\n",
    "def running_mean(x, n):\n",
    "    cumsum = np.cumsum(np.pad(x, int(ceil(n/2)), mode='reflect')) \n",
    "    return (cumsum[n:] - cumsum[:-n]) / float(n)\n",
    "\n",
    "for i, (p,t) in enumerate(zip(psucc, ts)):\n",
    "    if 4+i*2 > 18:\n",
    "        continue\n",
    "    try:\n",
    "        last_zero = np.argwhere(np.isclose(p, 0.))\n",
    "        print(f\"N={4+i*2}, p=2, t={np.asarray(t)[np.asarray(p)==0.].max()}, shots={np.asarray(shots[i])[np.asarray(p)==0.].max()}, nfev={np.asarray(nfevs[i])[np.asarray(p)==0.].max()}\")\n",
    "    except ValueError:\n",
    "        print(f\"No zeros in success probability for N={4+i*2}\")\n",
    "    # plot for each dimension\n",
    "    p_t = np.vstack((t, p))\n",
    "    p_t = p_t[:, np.argsort(p_t[0, :])]\n",
    "    len_window = 12\n",
    "    p_t[1] = running_mean(p_t[1], len_window)\n",
    "    plot1, = plt.plot(p_t[0][::5]*2, p_t[1][::5], \"-o\", label=f\"{4+i*2}\", color=colors[i], alpha=0.8)\n",
    "    plt.vlines(2**(4+i*2), 0, 1, colors=colors[i], linestyles='dashed', alpha=0.6, linewidth=2.5)\n",
    "    plots1.append(plot1)\n",
    "\n",
    "for i, (p,t) in enumerate(zip(psucc_noise, ts_noise)):\n",
    "    try:\n",
    "        last_zero = np.argwhere(np.isclose(psucc_noise[i], 0.))\n",
    "        print(f\"N={4+i*2}, p=2, t={np.asarray(ts_noise[i])[np.asarray(psucc_noise[i])==0.].max()}, shots={np.asarray(shots_noise[i])[np.asarray(psucc_noise[i])==0.].max()}, nfev={np.asarray(nfevs_noise[i])[np.asarray(psucc_noise[i])==0.].max()}\")\n",
    "    except ValueError:\n",
    "        print(f\"No zeros in succes probability for N={4+i*2}\")\n",
    "    p_t_noise = np.vstack((ts_noise[i], psucc_noise[i]))\n",
    "    p_t_noise = p_t_noise[:, np.argsort(p_t_noise[0, :])]\n",
    "    len_window = 12\n",
    "    p_t_noise[1] = running_mean(p_t_noise[1], len_window)\n",
    "    plot2, = plt.plot(p_t_noise[0][::5]*2, p_t_noise[1][::5], \"o-\", label=f\"{4+i*2}\", color=colors[i], fillstyle='none')\n",
    "    plots2.append(plot2)\n",
    "\n",
    "    plt.xscale('log')\n",
    "    plt.xlabel(r\"$t$ (shots x iterations x depth)\")\n",
    "    plt.ylabel(r\"$P$\")\n",
    "\n",
    "# Create a legend for the first line.\n",
    "first_legend = plt.legend(handles=plots1, loc='lower right', title=\"Exact\")\n",
    "#Add the legend manually to the current Axes.\n",
    "ax = plt.gca().add_artist(first_legend)\n",
    "\n",
    "plt.title(r\"$p=2$\")\n",
    "#Create another legend for the second line.\n",
    "plt.legend(handles=plots2, loc='upper left', title=f\"Noise\\nT1=50\\u03bcs\\nT2=70\\u03bcs\")\n",
    "\n",
    "plt.savefig(\"noise_exact_roll.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fun_exp(x, a, b):\n",
    "    return a*2**(b*x)\n",
    "\n",
    "def fun_lin(x, a, b):\n",
    "    return a + b*x\n",
    "    \n",
    "pt_all = []\n",
    "for i, (p,t) in enumerate(zip(psucc, ts)):\n",
    "    p_t = np.vstack((t, p))\n",
    "    p_t = p_t[:, np.argsort(p_t[0, :])]\n",
    "    pt_all.append(p_t)\n",
    "pt_all_2 = []\n",
    "for i, (p,t) in enumerate(zip(psucc_noise, ts_noise)):\n",
    "    p_t = np.vstack((t, p))\n",
    "    p_t = p_t[:, np.argsort(p_t[0, :])]\n",
    "    pt_all_2.append(p_t)\n",
    "\n",
    "\n",
    "colors = ['lightgreen', 'limegreen', 'forestgreen', 'darkgreen','lightseagreen', 'cornflowerblue' , 'blue', 'navy']\n",
    "min_psucc_lst = [0.99, 0.9, 0.5, 0.25]\n",
    "min_qubits = 4 \n",
    "max_qubits = 14\n",
    "plots1 = [] \n",
    "plots2 = []\n",
    "for min_psucc, color in zip(min_psucc_lst, colors[::2]):\n",
    "    step_size = []\n",
    "    for i, pt in enumerate(pt_all):\n",
    "        mask = pt[1] >= min_psucc\n",
    "        # check if the p is not reached\n",
    "        if mask.sum() == 0:\n",
    "            continue\n",
    "        step_size.append(pt[0][mask][0])\n",
    "\n",
    "    x_axes = 4 + 2*np.arange(len(step_size))\n",
    "\n",
    "    popt, _ = curve_fit(fun_lin, x_axes, np.log2(step_size))\n",
    "    plot1, = plt.plot(x_axes, step_size, \"o\", label=f\"p={min_psucc}\", color=color)\n",
    "    plots1.append(plot1)\n",
    "    plot1, = plt.plot(np.arange(min_qubits, max_qubits+1), fun_exp(np.arange(min_qubits, max_qubits+1), 2**popt[0], popt[1]), \"-\", label=f\"a={2**popt[0]:.2f}, b={popt[1]:.2f}\", color=color)\n",
    "    plots1.append(plot1)\n",
    "\n",
    "    step_size_2 = []\n",
    "    for i, pt_2 in enumerate(pt_all_2):\n",
    "        mask = pt_2[1] >= min_psucc\n",
    "        # check if the p is not reached\n",
    "        if mask.sum() == 0:\n",
    "            continue\n",
    "        step_size_2.append(pt_2[0][mask][0])\n",
    "    x_axes_2 = 4 + 2*np.arange(len(step_size_2))\n",
    "    popt, _ = curve_fit(fun_lin, x_axes_2, np.log2(step_size_2))\n",
    "    plot2, = plt.plot(x_axes_2, step_size_2, \"o\", label=f\"p={min_psucc}\", color=color, fillstyle='none')\n",
    "    plots2.append(plot2)\n",
    "    plot2, = plt.plot(np.arange(min_qubits, max_qubits+1), fun_exp(np.arange(min_qubits, max_qubits+1), 2**popt[0], popt[1]), \"--\", label=f\"a={2**popt[0]:.2f}, b={popt[1]:.2f}\", color=color)\n",
    "    plots2.append(plot2)\n",
    "\n",
    "plt.yscale('log', base=2)\n",
    "plt.ylabel(r\"$t$ (shots x iterations)\")\n",
    "plt.xlabel(r\"$N$\")\n",
    "\n",
    "\n",
    "# Create a legend for the first line.\n",
    "first_legend = plt.legend(handles=plots1, loc='lower right', title=\"Exact\", fontsize=\"x-small\")\n",
    "#Add the legend manually to the current Axes.\n",
    "ax = plt.gca().add_artist(first_legend)\n",
    "\n",
    "plt.title(r\"$p=2$\")\n",
    "#Create another legend for the second line.\n",
    "plt.legend(handles=plots2, loc='upper left', title=\"Noise\\n\" + r\"$T_1=50$\" + f\"\\u03bcs\" r\" $T_2=70$\" + f\"\\u03bcs\", fontsize=\"x-small\")\n",
    "\n",
    "plt.title(r\"Fit with $a \\cdot 2^{b \\cdot N}$\")\n",
    "\n",
    "plt.savefig(\"t_scale.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fun_exp(x, a, b):\n",
    "    return a*2**(x**b)\n",
    "\n",
    "def fun_lin(x, a, b):\n",
    "    return a + x**b\n",
    "    \n",
    "pt_all = []\n",
    "for i, (p,t) in enumerate(zip(psucc, ts)):\n",
    "    p_t = np.vstack((t, p))\n",
    "    p_t = p_t[:, np.argsort(p_t[0, :])]\n",
    "    pt_all.append(p_t)\n",
    "pt_all_2 = []\n",
    "for i, (p,t) in enumerate(zip(psucc_noise, ts_noise)):\n",
    "    p_t = np.vstack((t, p))\n",
    "    p_t = p_t[:, np.argsort(p_t[0, :])]\n",
    "    pt_all_2.append(p_t)\n",
    "\n",
    "\n",
    "colors = ['lightgreen', 'limegreen', 'forestgreen', 'darkgreen','lightseagreen', 'cornflowerblue' , 'blue', 'navy']\n",
    "min_psucc_lst = [0.99, 0.9, 0.5, 0.25]\n",
    "min_qubits = 4 \n",
    "max_qubits = 20\n",
    "plots1 = [] \n",
    "plots2 = []\n",
    "for min_psucc, color in zip(min_psucc_lst, colors[::2]):\n",
    "    step_size = []\n",
    "    for i, pt in enumerate(pt_all):\n",
    "        mask = pt[1] >= min_psucc\n",
    "        # check if the p is not reached\n",
    "        if mask.sum() == 0:\n",
    "            continue\n",
    "        step_size.append(pt[0][mask][0])\n",
    "\n",
    "    x_axes = 4 + 2*np.arange(len(step_size))\n",
    "\n",
    "    popt, _ = curve_fit(fun_lin, x_axes, np.log2(step_size))\n",
    "    plot1, = plt.plot(x_axes, step_size, \"o\", label=f\"p={min_psucc}\", color=color)\n",
    "    plots1.append(plot1)\n",
    "    plot1, = plt.plot(np.arange(min_qubits, max_qubits+1), fun_exp(np.arange(min_qubits, max_qubits+1), 2**popt[0], popt[1]), \"-\", label=f\"a={2**popt[0]:.2f}, b={popt[1]:.2f}\", color=color)\n",
    "    plots1.append(plot1)\n",
    "\n",
    "    step_size_2 = []\n",
    "    for i, pt_2 in enumerate(pt_all_2):\n",
    "        mask = pt_2[1] >= min_psucc\n",
    "        # check if the p is not reached\n",
    "        if mask.sum() == 0:\n",
    "            continue\n",
    "        step_size_2.append(pt_2[0][mask][0])\n",
    "    x_axes_2 = 4 + 2*np.arange(len(step_size_2))\n",
    "    popt, _ = curve_fit(fun_lin, x_axes_2, np.log2(step_size_2))\n",
    "    plot2, = plt.plot(x_axes_2, step_size_2, \"o\", label=f\"p={min_psucc}\", color=color, fillstyle='none')\n",
    "    plots2.append(plot2)\n",
    "    plot2, = plt.plot(np.arange(min_qubits, max_qubits+1), fun_exp(np.arange(min_qubits, max_qubits+1), 2**popt[0], popt[1]), \"--\", label=f\"a={2**popt[0]:.2f}, b={popt[1]:.2f}\", color=color)\n",
    "    plots2.append(plot2)\n",
    "\n",
    "plt.yscale('log', base=2)\n",
    "plt.ylabel(r\"$t$ (shots x iterations)\")\n",
    "plt.xlabel(r\"$N$\")\n",
    "\n",
    "\n",
    "# Create a legend for the first line.\n",
    "first_legend = plt.legend(handles=plots1, loc='lower right', title=\"Exact\", fontsize=\"small\")\n",
    "#Add the legend manually to the current Axes.\n",
    "ax = plt.gca().add_artist(first_legend)\n",
    "\n",
    "plt.title(r\"$p=2$\")\n",
    "#Create another legend for the second line.\n",
    "plt.legend(handles=plots2, loc='upper left', title=\"Noise\\n\" + r\"$T_1=50$\" + f\"\\u03bcs\" r\" $T_2=70$\" + f\"\\u03bcs\", fontsize=\"small\")\n",
    "\n",
    "plt.title(r\"Fit with $a \\cdot 2^{N^b}$\")\n",
    "\n",
    "plt.savefig(\"t_scale_stretch.png\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## QAOA+ vs VQE\n",
    "\n",
    "Custom implementation of a QAOA with an increased number of parameters, i.e., a number comparable to the parameters of a standard VQE, in order to make a fair comparison with VQE $p=2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import curve_fit\n",
    "from math import ceil\n",
    "from json import JSONEncoder\n",
    "from matplotlib import pyplot as plt\n",
    "from src.utils import collect_results, NumpyArrayEncoder, create_ansatz\n",
    "import json\n",
    "import os\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qc = create_ansatz(qubits=4, circ_depth=2, ansatz_type='qaoa+')\n",
    "#qc.draw() @uncomment if you want see the circuit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qubits = [4,6,8,10,12,14,16,18]\n",
    "p = 2\n",
    "# load exact simulations\n",
    "psucc_vqe, ts_vqe, shots_vqe, nfevs_vqe = collect_results(qubits, p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ts_qaoap = []\n",
    "shots_qaoap = []\n",
    "nfevs_qaoap = []\n",
    "psucc_qaoap = []\n",
    "for qubit in qubits:\n",
    "    dir_path = f\"results/N{qubit}/qaoa+/\" if qubit < 14 else f\"results/N{qubit}/qaoa+_iter80/\"\n",
    "    print(f\"directory: {dir_path}\")\n",
    "    p = []\n",
    "    t = []\n",
    "    s = []\n",
    "    it = []\n",
    "    # for each number of qubits\n",
    "    # we have several number of shots and iterations\n",
    "    for filename in sorted(os.listdir(dir_path)):\n",
    "        filename = dir_path + filename\n",
    "        with open(filename, \"r\") as file:\n",
    "            # print(filename)\n",
    "            data = json.load(file)\n",
    "        ever_found = []\n",
    "        # for each shot and iteration param\n",
    "        # we randomized the initial point\n",
    "        # to estimate the right probability\n",
    "        for run in data:\n",
    "            ever_found.append(run[\"ever_found\"])\n",
    "        # compute p('found minimum')\n",
    "        p.append(\n",
    "            np.mean(np.asarray(ever_found, dtype=np.float128))\n",
    "        )  # float128 to avoid to many zeros\n",
    "        # maxiter*shots = actual number of iteration\n",
    "        t.append(run[\"shots\"] * run[\"nfev\"])\n",
    "        s.append(run[\"shots\"])\n",
    "        it.append(run[\"nfev\"])\n",
    "    # update list for each number of qubits\n",
    "    ts_qaoap.append(t)\n",
    "    psucc_qaoap.append(p)\n",
    "    shots_qaoap.append(s)\n",
    "    nfevs_qaoap.append(it)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "colors = ['lightgreen','forestgreen','tab:green','darkgreen', 'lightseagreen', 'tab:cyan','cornflowerblue' , 'blue', 'navy']\n",
    "plots1 = []\n",
    "plots2 = []\n",
    "\n",
    "for i, (p,t) in enumerate(zip(psucc_vqe, ts_vqe)):\n",
    "    if 4+i*2 > 18:\n",
    "        continue\n",
    "    try:\n",
    "        last_zero = np.argwhere(np.isclose(p, 0.))\n",
    "        print(f\"N={4+i*2}, p=2, t={np.asarray(t)[np.asarray(p)==0.].max()}, shots={np.asarray(shots_vqe[i])[np.asarray(p)==0.].max()}, nfev={np.asarray(nfevs_vqe[i])[np.asarray(p)==0.].max()}\")\n",
    "    except ValueError:\n",
    "        print(f\"No zeros in success probability for N={4+i*2}\")\n",
    "    # plot for each dimension\n",
    "    p_t = np.vstack((t, p))\n",
    "    p_t = p_t[:, np.argsort(p_t[0, :])]\n",
    "    plot1, = plt.plot(p_t[0]*2, p_t[1], \".\", label=f\"{4+i*2}\", color=colors[i], alpha=0.8)\n",
    "    plots1.append(plot1)\n",
    "\n",
    "for i, (p,t) in enumerate(zip(psucc_qaoap, ts_qaoap)):\n",
    "    try:\n",
    "        last_zero = np.argwhere(np.isclose(psucc_qaoap[i], 0.))\n",
    "        print(f\"N={4+i*2}, p=2, t={np.asarray(ts_qaoap[i])[np.asarray(psucc_qaoap[i])==0.].max()}, shots={np.asarray(shots_qaoap[i])[np.asarray(psucc_qaoap[i])==0.].max()}, nfev={np.asarray(nfevs_qaoap[i])[np.asarray(psucc_qaoap[i])==0.].max()}\")\n",
    "    except ValueError:\n",
    "        print(f\"No zeros in succes probability for N={4+i*2}\")\n",
    "    p_t_qaoap = np.vstack((ts_qaoap[i], psucc_qaoap[i]))\n",
    "    p_t_qaoap = p_t_qaoap[:, np.argsort(p_t_qaoap[0, :])]\n",
    "    plot2, = plt.plot(p_t_qaoap[0]*2, p_t_qaoap[1], \".\", label=f\"{4+i*2}\", color=colors[i], fillstyle='none')\n",
    "    plots2.append(plot2)\n",
    "\n",
    "    plt.xscale('log')\n",
    "    plt.xlabel(r\"$t$ (shots x iterations x depth)\")\n",
    "    plt.ylabel(r\"$P$\")\n",
    "\n",
    "# Create a legend for the first line.\n",
    "first_legend = plt.legend(handles=plots1, loc='lower right', title=\"VQE\")\n",
    "#Add the legend manually to the current Axes.\n",
    "ax = plt.gca().add_artist(first_legend)\n",
    "\n",
    "plt.title(r\"$p=2$\")\n",
    "#Create another legend for the second line.\n",
    "plt.legend(handles=plots2, loc='upper left', title=f\"QAOA+\")\n",
    "\n",
    "#plt.savefig(\"qaoa_vqe.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "colors = ['lightgreen','tab:green','darkgreen', 'lightseagreen','cornflowerblue' , 'blue', 'navy']\n",
    "plots1 = []\n",
    "plots2 = []\n",
    "\n",
    "def running_mean(x, n):\n",
    "    cumsum = np.cumsum(np.pad(x, int(ceil(n/2)), mode='reflect')) \n",
    "    return (cumsum[n:] - cumsum[:-n]) / float(n)\n",
    "\n",
    "for i, (p,t) in enumerate(zip(psucc_vqe, ts_vqe)):\n",
    "    if 4+i*2 > 18:\n",
    "        continue\n",
    "    try:\n",
    "        last_zero = np.argwhere(np.isclose(p, 0.))\n",
    "        print(f\"N={4+i*2}, p=2, t={np.asarray(t)[np.asarray(p)==0.].max()}, shots={np.asarray(shots_vqe[i])[np.asarray(p)==0.].max()}, nfev={np.asarray(nfevs_vqe[i])[np.asarray(p)==0.].max()}\")\n",
    "    except ValueError:\n",
    "        print(f\"No zeros in success probability for N={4+i*2}\")\n",
    "    # plot for each dimension\n",
    "    p_t = np.vstack((t, p))\n",
    "    p_t = p_t[:, np.argsort(p_t[0, :])]\n",
    "    len_window = 12\n",
    "    p_t[1] = running_mean(p_t[1], len_window)\n",
    "    plot1, = plt.plot(p_t[0][::5]*2, p_t[1][::5], \"-o\", label=f\"{4+i*2}\", color=colors[i], alpha=0.8)\n",
    "    plt.vlines(2**(4+i*2), 0, 1, colors=colors[i], linestyles='dashed', alpha=0.6, linewidth=2.5)\n",
    "    plots1.append(plot1)\n",
    "\n",
    "for i, (p,t) in enumerate(zip(psucc_qaoap, ts_qaoap)):\n",
    "    try:\n",
    "        last_zero = np.argwhere(np.isclose(psucc_qaoap[i], 0.))\n",
    "        print(f\"N={4+i*2}, p=2, t={np.asarray(ts_qaoap[i])[np.asarray(psucc_qaoap[i])==0.].max()}, shots={np.asarray(shots_qaoap[i])[np.asarray(psucc_qaoap[i])==0.].max()}, nfev={np.asarray(nfevs_qaoap[i])[np.asarray(psucc_qaoap[i])==0.].max()}\")\n",
    "    except ValueError:\n",
    "        print(f\"No zeros in succes probability for N={4+i*2}\")\n",
    "    p_t_qaoap = np.vstack((ts_qaoap[i], psucc_qaoap[i]))\n",
    "    p_t_qaoap = p_t_qaoap[:, np.argsort(p_t_qaoap[0, :])]\n",
    "    len_window = 12\n",
    "    p_t_qaoap[1] = running_mean(p_t_qaoap[1], len_window)\n",
    "    plot2, = plt.plot(p_t_qaoap[0][::5]*2, p_t_qaoap[1][::5], \"o-\", label=f\"{4+i*2}\", color=colors[i], fillstyle='none')\n",
    "    plots2.append(plot2)\n",
    "\n",
    "    plt.xscale('log')\n",
    "    plt.xlabel(r\"$t$ (shots x iterations x depth)\")\n",
    "    plt.ylabel(r\"$P$\")\n",
    "\n",
    "# Create a legend for the first line.\n",
    "first_legend = plt.legend(handles=plots1, loc='lower right', title=\"Exact\")\n",
    "#Add the legend manually to the current Axes.\n",
    "ax = plt.gca().add_artist(first_legend)\n",
    "\n",
    "plt.title(r\"$p=2$\")\n",
    "#Create another legend for the second line.\n",
    "plt.legend(handles=plots2, loc='upper left', title=f\"QAOA+\")\n",
    "\n",
    "plt.savefig(\"noise_exact_roll.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fun_exp(x, a, b):\n",
    "    return a*2**(b*x)\n",
    "\n",
    "def fun_lin(x, a, b):\n",
    "    return a + b*x\n",
    "    \n",
    "pt_all = []\n",
    "for i, (p,t) in enumerate(zip(psucc_vqe, ts_vqe)):\n",
    "    p_t = np.vstack((t, p))\n",
    "    p_t = p_t[:, np.argsort(p_t[0, :])]\n",
    "    pt_all.append(p_t)\n",
    "pt_all_2 = []\n",
    "for i, (p,t) in enumerate(zip(psucc_qaoap, ts_qaoap)):\n",
    "    p_t = np.vstack((t, p))\n",
    "    p_t = p_t[:, np.argsort(p_t[0, :])]\n",
    "    pt_all_2.append(p_t)\n",
    "\n",
    "\n",
    "colors = ['lightgreen', 'limegreen', 'forestgreen', 'darkgreen','lightseagreen', 'cornflowerblue' , 'blue', 'navy']\n",
    "min_psucc_lst = [0.99, 0.9, 0.5, 0.25]\n",
    "min_qubits = 4 \n",
    "max_qubits = 16\n",
    "plots1 = [] \n",
    "plots2 = []\n",
    "\n",
    "plt.plot(np.arange(min_qubits, max_qubits+1), fun_exp(np.arange(min_qubits, max_qubits+1), 1, 1), color='red', alpha=0.7)\n",
    "\n",
    "for min_psucc, color in zip(min_psucc_lst, colors[::2]):\n",
    "    step_size = []\n",
    "    for i, pt in enumerate(pt_all):\n",
    "        mask = pt[1] >= min_psucc\n",
    "        # check if the p is not reached\n",
    "        if mask.sum() == 0:\n",
    "            continue\n",
    "        step_size.append(pt[0][mask][0]*2)\n",
    "\n",
    "    x_axes = 4 + 2*np.arange(len(step_size))\n",
    "\n",
    "    popt, _ = curve_fit(fun_lin, x_axes, np.log2(step_size))\n",
    "    plot1, = plt.plot(x_axes, step_size, \"o\", label=f\"p={min_psucc}\", color=color)\n",
    "    plots1.append(plot1)\n",
    "    plot1, = plt.plot(np.arange(min_qubits, max_qubits+1), fun_exp(np.arange(min_qubits, max_qubits+1), 2**popt[0], popt[1]), \"-\", label=f\"a={2**popt[0]:.2f}, b={popt[1]:.2f}\", color=color)\n",
    "    plots1.append(plot1)\n",
    "\n",
    "    step_size_2 = []\n",
    "    for i, pt_2 in enumerate(pt_all_2):\n",
    "        mask = pt_2[1] >= min_psucc\n",
    "        # check if the p is not reached\n",
    "        if mask.sum() == 0:\n",
    "            continue\n",
    "        step_size_2.append(pt_2[0][mask][0]*2)\n",
    "    x_axes_2 = 4 + 2*np.arange(len(step_size_2))\n",
    "    popt, _ = curve_fit(fun_lin, x_axes_2, np.log2(step_size_2))\n",
    "    plot2, = plt.plot(x_axes_2, step_size_2, \"o\", label=f\"p={min_psucc}\", color=color, fillstyle='none')\n",
    "    plots2.append(plot2)\n",
    "    plot2, = plt.plot(np.arange(min_qubits, max_qubits+1), fun_exp(np.arange(min_qubits, max_qubits+1), 2**popt[0], popt[1]), \"--\", label=f\"a={2**popt[0]:.2f}, b={popt[1]:.2f}\", color=color)\n",
    "    plots2.append(plot2)\n",
    "\n",
    "plt.yscale('log', base=2)\n",
    "plt.ylabel(r\"$t$ (shots x iterations)\")\n",
    "plt.xlabel(r\"$N$\")\n",
    "\n",
    "\n",
    "# Create a legend for the first line.\n",
    "first_legend = plt.legend(handles=plots1, loc='upper left', title=\"VQE\", fontsize=\"x-small\")\n",
    "#Add the legend manually to the current Axes.\n",
    "ax = plt.gca().add_artist(first_legend)\n",
    "\n",
    "plt.title(r\"$p=2$\")\n",
    "#Create another legend for the second line.\n",
    "plt.legend(handles=plots2, loc='lower right', title=\"QAOA+\", fontsize=\"x-small\")\n",
    "\n",
    "plt.title(r\"Fit with $a \\cdot 2^{b \\cdot N}$\")\n",
    "\n",
    "#plt.savefig(\"t_scale.png\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### Comparison QAOA+ and VQE for iteration\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for k, (p_qaoa, s_qaoa, nfev_qaoa, p_vqe, s_vqe, nfev_vqe) in enumerate(zip(psucc_qaoap, shots_qaoap, nfevs_qaoap, psucc_vqe, shots_vqe, nfevs_vqe)):\n",
    "    len_iter_qaoa = np.max(nfev_qaoa)//2+1\n",
    "    p_qaoa = np.reshape(p_qaoa, (-1,len_iter_qaoa))\n",
    "    s_qaoa = np.reshape(s_qaoa, (-1,len_iter_qaoa))\n",
    "    nfev_qaoa = np.reshape(nfev_qaoa, (-1,len_iter_qaoa))\n",
    "\n",
    "    len_iter_vqe = np.max(nfev_vqe)//2+1\n",
    "    p_vqe = np.reshape(p_vqe, (-1,len_iter_vqe))\n",
    "    s_vqe = np.reshape(s_vqe, (-1,len_iter_vqe))\n",
    "    nfev_vqe = np.reshape(nfev_vqe, (-1,len_iter_vqe))\n",
    "\n",
    "    fig, (ax1, ax2) = plt.subplots(1,2, figsize=(16,6))\n",
    "    for i in range(p_qaoa.shape[0]):\n",
    "        ax1.plot(nfev_qaoa[i,:], p_qaoa[i,:], '--o', label=s_qaoa[i,0])\n",
    "    for i in range(p_vqe.shape[0]):\n",
    "        ax2.plot(nfev_vqe[i,:], p_vqe[i,:], '--o', label=s_vqe[i,0])\n",
    "\n",
    "    ax1.set_xlabel('Iterations')\n",
    "    ax1.set_ylabel(r'$F(success \\leq t)$')\n",
    "    #ax1.set_ylim(0,1.1)\n",
    "\n",
    "    ax2.set_xlabel('Iterations')\n",
    "    ax2.set_ylabel(r'$F(success \\leq t)$')\n",
    "    #ax2.set_ylim(0,1.1)\n",
    "\n",
    "    ax1.legend(title='Shots')\n",
    "    ax2.legend(title='Shots')\n",
    "\n",
    "    ax1.set_title(f'QAOA+ N={4+k*2}')\n",
    "    ax2.set_title(f'VQE N={4+k*2}')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### Comparison QAOA+ and VQE for shots\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for k, (p_qaoa, s_qaoa, nfev_qaoa, t_qaoa, p_vqe, s_vqe, nfev_vqe, t_vqe) in enumerate(zip(psucc_qaoap, shots_qaoap, nfevs_qaoap, ts_qaoap, psucc_vqe, shots_vqe, nfevs_vqe, ts_vqe)):\n",
    "    len_iter_qaoa = np.max(nfev_qaoa)//2+1\n",
    "    p_qaoa = np.reshape(p_qaoa, (-1,len_iter_qaoa))\n",
    "    s_qaoa = np.reshape(s_qaoa, (-1,len_iter_qaoa))\n",
    "    nfev_qaoa = np.reshape(nfev_qaoa, (-1,len_iter_qaoa))\n",
    "    t_qaoa = np.reshape(t_qaoa, (-1,len_iter_qaoa))\n",
    "\n",
    "    len_iter_vqe = np.max(nfev_vqe)//2+1\n",
    "    p_vqe = np.reshape(p_vqe, (-1,len_iter_vqe))\n",
    "    s_vqe = np.reshape(s_vqe, (-1,len_iter_vqe))\n",
    "    nfev_vqe = np.reshape(nfev_vqe, (-1,len_iter_vqe))\n",
    "    t_vqe = np.reshape(t_vqe, (-1,len_iter_vqe))\n",
    "\n",
    "    fig, (ax1, ax2) = plt.subplots(1,2, figsize=(16,6))\n",
    "    for j in range(p_qaoa.shape[1]):\n",
    "        ax1.plot(s_qaoa[:,j], p_qaoa[:,j], '--s', label=nfev_qaoa[0,j])\n",
    "    for j in range(p_vqe.shape[1]):\n",
    "        ax2.plot(s_vqe[:,j], p_vqe[:,j], '--s', label=nfev_vqe[0,j])\n",
    "\n",
    "    ax1.set_xlabel('Shots')\n",
    "    ax1.set_xscale('log', base=2)\n",
    "    ax1.set_ylabel(r'$F(success \\leq t)$')\n",
    "    #ax1.set_ylim(0,1.1)\n",
    "\n",
    "    ax2.set_xlabel('Shots')\n",
    "    ax2.set_xscale('log', base=2)\n",
    "    ax2.set_ylabel(r'$F(success \\leq t)$')\n",
    "    #ax2.set_ylim(0,1.1)\n",
    "\n",
    "    ax1.legend(title='Iterations')\n",
    "    ax2.legend(title='Iterations')\n",
    "\n",
    "    ax1.set_title(f'QAOA+ N={4+k*2}')\n",
    "    ax2.set_title(f'VQE N={4+k*2}')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### Comparison QAOA+ and VQE for $t$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for k, (p_qaoa, s_qaoa, nfev_qaoa, t_qaoa, p_vqe, s_vqe, nfev_vqe, t_vqe) in enumerate(zip(psucc_qaoap, shots_qaoap, nfevs_qaoap, ts_qaoap, psucc_vqe, shots_vqe, nfevs_vqe, ts_vqe)):\n",
    "    len_iter_qaoa = np.max(nfev_qaoa)//2+1 if 4+k*2 < 14 else np.max(nfev_qaoa)//5+1\n",
    "    p_qaoa = np.reshape(p_qaoa, (-1,len_iter_qaoa))\n",
    "    s_qaoa = np.reshape(s_qaoa, (-1,len_iter_qaoa))\n",
    "    nfev_qaoa = np.reshape(nfev_qaoa, (-1,len_iter_qaoa))\n",
    "    t_qaoa = np.reshape(t_qaoa, (-1,len_iter_qaoa))*2\n",
    "\n",
    "    len_iter_vqe = np.max(nfev_vqe)//2+1\n",
    "    p_vqe = np.reshape(p_vqe, (-1,len_iter_vqe))\n",
    "    s_vqe = np.reshape(s_vqe, (-1,len_iter_vqe))\n",
    "    nfev_vqe = np.reshape(nfev_vqe, (-1,len_iter_vqe))\n",
    "    t_vqe = np.reshape(t_vqe, (-1,len_iter_vqe))*2\n",
    "\n",
    "    fig, (ax1, ax2) = plt.subplots(1,2, figsize=(16,6))\n",
    "    for i in range(p_qaoa.shape[0]):\n",
    "        ax1.plot(t_qaoa[i,:], p_qaoa[i,:], '--o', label=s_qaoa[i,0])\n",
    "    ax1.plot(np.arange(0,2**(7+k*2),2), 1 - (1-1/(2**(4+k*2)))**np.arange(0,2**(7+k*2),2), '--', color=\"black\", label=\"Random sampling\")\n",
    "    for i in range(p_vqe.shape[0]):\n",
    "        ax2.plot(t_vqe[i,:], p_vqe[i,:], '--o', label=s_vqe[i,0])\n",
    "    ax2.plot(np.arange(0,2**(7+k*2),2), 1 - (1-1/(2**(4+k*2)))**np.arange(0,2**(7+k*2),2), '--', color=\"black\", label=\"Random sampling\")\n",
    "    #ax1.plot(np.arange(0, 2**(2*(4+k*2))), 1 - (1-1/(2**(4+k*2)))**np.arange(0, 2**(2*(4+k*2))), '--', color=\"black\", label=\"Random sampling\")\n",
    "\n",
    "    ax1.set_xlabel('t (shots x iterations)')\n",
    "    ax1.set_ylabel(r'$F(success \\leq t)$')\n",
    "    ax1.set_xscale('log', base=2)\n",
    "    #ax1.set_yscale('log', base=2)\n",
    "    #ax1.set_ylim(0,1.1)\n",
    "\n",
    "    ax2.set_xlabel('t (shots x iterations)')\n",
    "    ax2.set_ylabel(r'$F(success \\leq t)$')\n",
    "    ax2.set_xscale('log', base=2)\n",
    "    #ax2.set_yscale('log', base=2)\n",
    "    #ax2.set_ylim(0,1.1)\n",
    "\n",
    "    ax1.legend(title='Shots')\n",
    "    ax2.legend(title='Shots')\n",
    "\n",
    "    ax1.set_title(f'QAOA+ N={4+k*2}')\n",
    "    ax2.set_title(f'VQE N={4+k*2}')\n",
    "\n",
    "    # if 4+k*2 == 14:\n",
    "    #     plt.savefig(\"qaoa_vqe.png\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### QAOA+ $p=1$ vs VQE $p=1$ and QAOA+ $p=2$\n",
    "\n",
    "Comparison only for size $N=12$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dir_path = f\"results/N12/qaoa+_p1/\"\n",
    "print(f\"directory: {dir_path}\")\n",
    "p_qaoa = []\n",
    "t_qaoa = []\n",
    "s_qaoa = []\n",
    "nfev_qaoa = []\n",
    "# for each number of qubits\n",
    "# we have several number of shots and iterations\n",
    "for filename in sorted(os.listdir(dir_path)):\n",
    "    filename = dir_path + filename\n",
    "    with open(filename, \"r\") as file:\n",
    "        # print(filename)\n",
    "        data = json.load(file)\n",
    "    ever_found = []\n",
    "    # for each shot and iteration param\n",
    "    # we randomized the initial point\n",
    "    # to estimate the right probability\n",
    "    for run in data:\n",
    "        ever_found.append(run[\"ever_found\"])\n",
    "    # compute p('found minimum')\n",
    "    p_qaoa.append(\n",
    "        np.mean(np.asarray(ever_found, dtype=np.float128))\n",
    "    )  # float128 to avoid to many zeros\n",
    "    # maxiter*shots = actual number of iteration\n",
    "    t_qaoa.append(run[\"shots\"] * run[\"nfev\"])\n",
    "    s_qaoa.append(run[\"shots\"])\n",
    "    nfev_qaoa.append(run[\"nfev\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 4\n",
    "p_vqe, s_vqe, nfev_vqe, t_vqe = psucc_1[k], shots_1[k], nfevs_1[k], ts_1[k]\n",
    "\n",
    "\n",
    "len_iter_qaoa = np.max(nfev_qaoa)//2+1 if 4+k*2 < 14 else np.max(nfev_qaoa)//5+1\n",
    "p_qaoa = np.reshape(p_qaoa, (-1,len_iter_qaoa))\n",
    "s_qaoa = np.reshape(s_qaoa, (-1,len_iter_qaoa))\n",
    "nfev_qaoa = np.reshape(nfev_qaoa, (-1,len_iter_qaoa))\n",
    "t_qaoa = np.reshape(t_qaoa, (-1,len_iter_qaoa))\n",
    "\n",
    "len_iter_vqe = np.max(nfev_vqe)//2+1\n",
    "p_vqe = np.reshape(p_vqe, (-1,len_iter_vqe))\n",
    "s_vqe = np.reshape(s_vqe, (-1,len_iter_vqe))\n",
    "nfev_vqe = np.reshape(nfev_vqe, (-1,len_iter_vqe))\n",
    "t_vqe = np.reshape(t_vqe, (-1,len_iter_vqe))\n",
    "\n",
    "fig, (ax1, ax2) = plt.subplots(1,2, figsize=(16,6))\n",
    "for i in range(p_qaoa.shape[0]):\n",
    "    ax1.plot(t_qaoa[i,:], p_qaoa[i,:], '--o', label=s_qaoa[i,0])\n",
    "ax1.plot(np.arange(0,2**(7+k*2),2), 1 - (1-1/(2**(4+k*2)))**np.arange(0,2**(7+k*2),2), '--', color=\"black\", label=\"Random sampling\")\n",
    "for i in range(p_vqe.shape[0]):\n",
    "    ax2.plot(t_vqe[i,:], p_vqe[i,:], '--o', label=s_vqe[i,0])\n",
    "ax2.plot(np.arange(0,2**(7+k*2),2), 1 - (1-1/(2**(4+k*2)))**np.arange(0,2**(7+k*2),2), '--', color=\"black\", label=\"Random sampling\")\n",
    "#ax1.plot(np.arange(0, 2**(2*(4+k*2))), 1 - (1-1/(2**(4+k*2)))**np.arange(0, 2**(2*(4+k*2))), '--', color=\"black\", label=\"Random sampling\")\n",
    "\n",
    "ax1.set_xlabel('t (shots x iterations)')\n",
    "ax1.set_ylabel(r'$F(success \\leq t)$')\n",
    "ax1.set_xscale('log', base=2)\n",
    "#ax1.set_yscale('log', base=2)\n",
    "#ax1.set_ylim(0,1.1)\n",
    "\n",
    "ax2.set_xlabel('t (shots x iterations)')\n",
    "ax2.set_ylabel(r'$F(success \\leq t)$')\n",
    "ax2.set_xscale('log', base=2)\n",
    "#ax2.set_yscale('log', base=2)\n",
    "#ax2.set_ylim(0,1.1)\n",
    "\n",
    "ax1.legend(title='Shots')\n",
    "ax2.legend(title='Shots')\n",
    "\n",
    "ax1.set_title(f'QAOA+ N={4+k*2}')\n",
    "ax2.set_title(f'VQE N={4+k*2}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k=4\n",
    "p_qaoa_p2, s_qaoa_p2, nfev_qaoa_p2, t_qaoa_p2 = psucc_qaoap[k], shots_qaoap[k], nfevs_qaoap[k], ts_qaoap[k]\n",
    "\n",
    "\n",
    "len_iter_qaoa = np.max(nfev_qaoa)//2+1 if 4+k*2 < 14 else np.max(nfev_qaoa)//5+1\n",
    "p_qaoa = np.reshape(p_qaoa, (-1,len_iter_qaoa))\n",
    "s_qaoa = np.reshape(s_qaoa, (-1,len_iter_qaoa))\n",
    "nfev_qaoa = np.reshape(nfev_qaoa, (-1,len_iter_qaoa))\n",
    "t_qaoa = np.reshape(t_qaoa, (-1,len_iter_qaoa))\n",
    "\n",
    "len_iter_vqe = np.max(nfev_qaoa_p2)//2+1\n",
    "p_qaoa_p2 = np.reshape(p_qaoa_p2, (-1,len_iter_vqe))\n",
    "s_qaoa_p2 = np.reshape(s_qaoa_p2, (-1,len_iter_vqe))\n",
    "nfev_qaoa_p2 = np.reshape(nfev_qaoa_p2, (-1,len_iter_vqe))\n",
    "t_qaoa_p2 = np.reshape(t_qaoa_p2, (-1,len_iter_vqe))\n",
    "\n",
    "fig, (ax1, ax2) = plt.subplots(1,2, figsize=(16,6))\n",
    "for i in range(p_qaoa.shape[0]):\n",
    "    ax1.plot(t_qaoa[i,:], p_qaoa[i,:], '--o', label=s_qaoa[i,0])\n",
    "ax1.plot(np.arange(0,2**(7+k*2),2), 1 - (1-1/(2**(4+k*2)))**np.arange(0,2**(7+k*2),2), '--', color=\"black\", label=\"Random sampling\")\n",
    "for i in range(p_qaoa_p2.shape[0]):\n",
    "    ax2.plot(t_qaoa_p2[i,:]*2, p_qaoa_p2[i,:], '--o', label=s_qaoa_p2[i,0])\n",
    "ax2.plot(np.arange(0,2**(7+k*2),2), 1 - (1-1/(2**(4+k*2)))**np.arange(0,2**(7+k*2),2), '--', color=\"black\", label=\"Random sampling\")\n",
    "#ax1.plot(np.arange(0, 2**(2*(4+k*2))), 1 - (1-1/(2**(4+k*2)))**np.arange(0, 2**(2*(4+k*2))), '--', color=\"black\", label=\"Random sampling\")\n",
    "\n",
    "ax1.set_xlabel('t (shots x iterations)')\n",
    "ax1.set_ylabel(r'$F(success \\leq t)$')\n",
    "ax1.set_xscale('log', base=2)\n",
    "#ax1.set_yscale('log', base=2)\n",
    "#ax1.set_ylim(0,1.1)\n",
    "\n",
    "ax2.set_xlabel('t (shots x iterations)')\n",
    "ax2.set_ylabel(r'$F(success \\leq t)$')\n",
    "ax2.set_xscale('log', base=2)\n",
    "#ax2.set_yscale('log', base=2)\n",
    "#ax2.set_ylim(0,1.1)\n",
    "\n",
    "ax1.legend(title='Shots')\n",
    "ax2.legend(title='Shots')\n",
    "\n",
    "ax1.set_title(f'QAOA+ p=1 N={4+k*2}')\n",
    "ax2.set_title(f'QAOA+ p=2 N={4+k*2}')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## QAOA $p=2$ vs VQE\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import curve_fit\n",
    "from math import ceil\n",
    "from json import JSONEncoder\n",
    "from matplotlib import pyplot as plt\n",
    "from src.utils import collect_results, NumpyArrayEncoder, create_ansatz\n",
    "import json\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qubits = [8,10,12,14,16,18]\n",
    "p = 2\n",
    "# load exact simulations\n",
    "psucc_vqe, ts_vqe, shots_vqe, nfevs_vqe = collect_results(qubits, p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ts_qaoa = []\n",
    "shots_qaoa = []\n",
    "nfevs_qaoa = []\n",
    "psucc_qaoa = []\n",
    "for qubit in qubits:\n",
    "    dir_path = f\"results/N{qubit}/qaoa_p2/\" \n",
    "    print(f\"directory: {dir_path}\")\n",
    "    p = []\n",
    "    t = []\n",
    "    s = []\n",
    "    it = []\n",
    "    # for each number of qubits\n",
    "    # we have several number of shots and iterations\n",
    "    for filename in sorted(os.listdir(dir_path)):\n",
    "        filename = dir_path + filename\n",
    "        with open(filename, \"r\") as file:\n",
    "            # print(filename)\n",
    "            data = json.load(file)\n",
    "        ever_found = []\n",
    "        # for each shot and iteration param\n",
    "        # we randomized the initial point\n",
    "        # to estimate the right probability\n",
    "        for run in data:\n",
    "            ever_found.append(run[\"ever_found\"])\n",
    "        # compute p('found minimum')\n",
    "        p.append(\n",
    "            np.mean(np.asarray(ever_found, dtype=np.float128))\n",
    "        )  # float128 to avoid to many zeros\n",
    "        # maxiter*shots = actual number of iteration\n",
    "        t.append(run[\"shots\"] * run[\"nfev\"])\n",
    "        s.append(run[\"shots\"])\n",
    "        it.append(run[\"nfev\"])\n",
    "    # update list for each number of qubits\n",
    "    ts_qaoa.append(t)\n",
    "    psucc_qaoa.append(p)\n",
    "    shots_qaoa.append(s)\n",
    "    nfevs_qaoa.append(it)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "colors = ['lightgreen','forestgreen','tab:green','darkgreen', 'lightseagreen', 'tab:cyan','cornflowerblue' , 'blue', 'navy']\n",
    "plots1 = []\n",
    "plots2 = []\n",
    "\n",
    "for i, (p,t) in enumerate(zip(psucc_vqe, ts_vqe)):\n",
    "    if 4+i*2 > 18:\n",
    "        continue\n",
    "    try:\n",
    "        last_zero = np.argwhere(np.isclose(p, 0.))\n",
    "        print(f\"N={4+i*2}, p=2, t={np.asarray(t)[np.asarray(p)==0.].max()}, shots={np.asarray(shots_vqe[i])[np.asarray(p)==0.].max()}, nfev={np.asarray(nfevs_vqe[i])[np.asarray(p)==0.].max()}\")\n",
    "    except ValueError:\n",
    "        print(f\"No zeros in success probability for N={4+i*2}\")\n",
    "    # plot for each dimension\n",
    "    p_t = np.vstack((t, p))\n",
    "    p_t = p_t[:, np.argsort(p_t[0, :])]\n",
    "    plot1, = plt.plot(p_t[0]*2, p_t[1], \".\", label=f\"{4+i*2}\", color=colors[i], alpha=0.8)\n",
    "    plots1.append(plot1)\n",
    "\n",
    "for i, (p,t) in enumerate(zip(psucc_qaoa, ts_qaoa)):\n",
    "    try:\n",
    "        last_zero = np.argwhere(np.isclose(psucc_qaoa[i], 0.))\n",
    "        print(f\"N={4+i*2}, p=2, t={np.asarray(ts_qaoa[i])[np.asarray(psucc_qaoa[i])==0.].max()}, shots={np.asarray(shots_qaoa[i])[np.asarray(psucc_qaoa[i])==0.].max()}, nfev={np.asarray(nfevs_qaoa[i])[np.asarray(psucc_qaoa[i])==0.].max()}\")\n",
    "    except ValueError:\n",
    "        print(f\"No zeros in succes probability for N={4+i*2}\")\n",
    "    p_t_qaoa = np.vstack((ts_qaoa[i], psucc_qaoa[i]))\n",
    "    p_t_qaoa = p_t_qaoa[:, np.argsort(p_t_qaoa[0, :])]\n",
    "    plot2, = plt.plot(p_t_qaoa[0]*2, p_t_qaoa[1], \".\", label=f\"{4+i*2}\", color=colors[i], fillstyle='none')\n",
    "    plots2.append(plot2)\n",
    "\n",
    "    plt.xscale('log')\n",
    "    plt.xlabel(r\"$t$ (shots x iterations x depth)\")\n",
    "    plt.ylabel(r\"$P$\")\n",
    "\n",
    "# Create a legend for the first line.\n",
    "first_legend = plt.legend(handles=plots1, loc='lower right', title=\"VQE\")\n",
    "#Add the legend manually to the current Axes.\n",
    "ax = plt.gca().add_artist(first_legend)\n",
    "\n",
    "plt.title(r\"$p=2$\")\n",
    "#Create another legend for the second line.\n",
    "plt.legend(handles=plots2, loc='upper left', title=f\"QAOA\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fun_exp(x, a, b):\n",
    "    return a*2**(b*x)\n",
    "\n",
    "def fun_lin(x, a, b):\n",
    "    return a + b*x\n",
    "    \n",
    "pt_all = []\n",
    "for i, (p,t) in enumerate(zip(psucc_vqe, ts_vqe)):\n",
    "    p_t = np.vstack((t, p))\n",
    "    p_t[0] *= 2 \n",
    "    p_t = p_t[:, np.argsort(p_t[0, :])]\n",
    "    pt_all.append(p_t)\n",
    "pt_all_2 = []\n",
    "for i, (p,t) in enumerate(zip(psucc_qaoa, ts_qaoa)):\n",
    "    p_t = np.vstack((t, p))\n",
    "    p_t[0] *= 2 \n",
    "    p_t = p_t[:, np.argsort(p_t[0, :])]\n",
    "    pt_all_2.append(p_t)\n",
    "\n",
    "\n",
    "colors = ['lightgreen', 'limegreen', 'forestgreen', 'darkgreen','lightseagreen', 'cornflowerblue' , 'blue', 'navy']\n",
    "min_psucc_lst = [0.99, 0.9, 0.5, 0.25]\n",
    "min_qubits = 8 \n",
    "max_qubits = 18\n",
    "plots1 = [] \n",
    "plots2 = []\n",
    "\n",
    "plt.plot(np.arange(min_qubits, max_qubits+1), fun_exp(np.arange(min_qubits, max_qubits+1), 1, 1), color='red', linewidth=3.5, linestyle=\"dashed\")\n",
    "\n",
    "for min_psucc, color in zip(min_psucc_lst, colors[::2]):\n",
    "    step_size = []\n",
    "    for i, pt in enumerate(pt_all):\n",
    "        mask = pt[1] >= min_psucc\n",
    "        # check if the p is not reached\n",
    "        if mask.sum() == 0:\n",
    "            continue\n",
    "        step_size.append(pt[0][mask][0])\n",
    "\n",
    "    x_axes = 8 + 2*np.arange(len(step_size))\n",
    "    popt, _ = curve_fit(fun_lin, x_axes, np.log2(step_size))\n",
    "    plot1, = plt.plot(x_axes, step_size, \"o\", label=f\"p={min_psucc}\", color=color)\n",
    "    plots1.append(plot1)\n",
    "    plot1, = plt.plot(np.arange(min_qubits, max_qubits+1), fun_exp(np.arange(min_qubits, max_qubits+1), 2**popt[0], popt[1]), \"-\", label=f\"a={2**popt[0]:.2f}, b={popt[1]:.2f}\", color=color)\n",
    "    plots1.append(plot1)\n",
    "\n",
    "    step_size_2 = []\n",
    "    for i, pt_2 in enumerate(pt_all_2):\n",
    "        mask = pt_2[1] >= min_psucc\n",
    "        # check if the p is not reached\n",
    "        if mask.sum() == 0:\n",
    "            continue\n",
    "        step_size_2.append(pt_2[0][mask][0])\n",
    "\n",
    "    x_axes_2 = 8 + 2*np.arange(len(step_size_2))\n",
    "    popt, _ = curve_fit(fun_lin, x_axes_2, np.log2(step_size_2))\n",
    "    plot2, = plt.plot(x_axes_2, step_size_2, \"o\", label=f\"p={min_psucc}\", color=color, fillstyle='none')\n",
    "    plots2.append(plot2)\n",
    "    plot2, = plt.plot(np.arange(min_qubits, max_qubits+1), fun_exp(np.arange(min_qubits, max_qubits+1), 2**popt[0], popt[1]), \"--\", label=f\"a={2**popt[0]:.2f}, b={popt[1]:.2f}\", color=color)\n",
    "    plots2.append(plot2)\n",
    "\n",
    "plt.yscale('log', base=2)\n",
    "plt.ylabel(r\"$t$ (shots x iterations)\")\n",
    "plt.xlabel(r\"$N$\")\n",
    "\n",
    "\n",
    "# Create a legend for the first line.\n",
    "first_legend = plt.legend(handles=plots1, loc='upper left', title=\"VQE\", fontsize=\"x-small\")\n",
    "#Add the legend manually to the current Axes.\n",
    "ax = plt.gca().add_artist(first_legend)\n",
    "\n",
    "plt.title(r\"$p=2$\")\n",
    "#Create another legend for the second line.\n",
    "plt.legend(handles=plots2, loc='upper right', title=\"QAOA\", fontsize=\"x-small\")\n",
    "\n",
    "plt.title(r\"Fit with $a \\cdot 2^{b \\cdot N}$\")\n",
    "\n",
    "plt.savefig(\"scale_qaoap2.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_108 = fun_exp(108, 166, 0.39) / 2\n",
    "print(f\"{t_108:.2e} circuit repetions (shots x iterations)\")\n",
    "print(f\"{t_108*10e-9:.2e}s\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for k, (p_qaoa, s_qaoa, nfev_qaoa, t_qaoa, p_vqe, s_vqe, nfev_vqe, t_vqe) in enumerate(zip(psucc_qaoa, shots_qaoa, nfevs_qaoa, ts_qaoa, psucc_vqe, shots_vqe, nfevs_vqe, ts_vqe)):\n",
    "    len_iter_qaoa = np.max(nfev_qaoa)//2+1\n",
    "    #print(len_iter_qaoa, np.max(nfev_qaoa), len(s_qaoa), 4+k*2)\n",
    "    p_qaoa = np.reshape(p_qaoa, (-1,len_iter_qaoa))\n",
    "    s_qaoa = np.reshape(s_qaoa, (-1,len_iter_qaoa))\n",
    "    nfev_qaoa = np.reshape(nfev_qaoa, (-1,len_iter_qaoa))\n",
    "    t_qaoa = np.reshape(t_qaoa, (-1,len_iter_qaoa))*2\n",
    "\n",
    "    len_iter_vqe = np.max(nfev_vqe)//2+1\n",
    "    p_vqe = np.reshape(p_vqe, (-1,len_iter_vqe))\n",
    "    s_vqe = np.reshape(s_vqe, (-1,len_iter_vqe))\n",
    "    nfev_vqe = np.reshape(nfev_vqe, (-1,len_iter_vqe))\n",
    "    t_vqe = np.reshape(t_vqe, (-1,len_iter_vqe))*2\n",
    "\n",
    "    fig, (ax1, ax2) = plt.subplots(1,2, figsize=(16,6))\n",
    "    for i in range(p_qaoa.shape[0]):\n",
    "        ax1.plot(t_qaoa[i,:], p_qaoa[i,:], '--o', label=s_qaoa[i,0])\n",
    "    ax1.plot(np.arange(0, 2**(8+k*2)), 1 - (1-1/(2**(4+k*2)))**np.arange(0, 2**(8+k*2)), '--', color=\"black\", label=\"Random sampling\")\n",
    "    for i in range(p_vqe.shape[0]):\n",
    "        ax2.plot(t_vqe[i,:], p_vqe[i,:], '--o', label=s_vqe[i,0])\n",
    "    ax2.plot(np.arange(0, 2**(8+k*2)), 1 - (1-1/(2**(4+k*2)))**np.arange(0, 2**(8+k*2)), '--', color=\"black\", label=\"Random sampling\")\n",
    "\n",
    "    ax1.set_xlabel('t (shots x iterations)')\n",
    "    ax1.set_ylabel(r'$F(success \\leq t)$')\n",
    "    ax1.set_xscale('log', base=2)\n",
    "    #ax1.set_yscale('log', base=2)\n",
    "    #ax1.set_ylim(0,1.1)\n",
    "\n",
    "    ax2.set_xlabel('t')\n",
    "    ax2.set_ylabel(r'$F(success \\leq t)$')\n",
    "    ax2.set_xscale('log', base=2)\n",
    "    #ax2.set_yscale('log', base=2)\n",
    "    #ax2.set_ylim(0,1.1)\n",
    "\n",
    "    ax1.legend(title='Shots')\n",
    "    ax2.legend(title='Shots')\n",
    "\n",
    "    ax1.set_title(f'QAOA N={4+k*2}')\n",
    "    ax2.set_title(f'VQE N={4+k*2}')\n",
    "\n",
    "    #if 4+k*2 == 14:\n",
    "    #    plt.savefig(\"qaoa_vqe.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## QAOA $p=3$ vs QAOA $p=2$\n",
    "\n",
    "Comparison between QAOA (classic version) with different depth."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import curve_fit\n",
    "from math import ceil\n",
    "from json import JSONEncoder\n",
    "from matplotlib import pyplot as plt\n",
    "from src.utils import collect_results, NumpyArrayEncoder, create_ansatz\n",
    "import json\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qubits = [4,6,8,10,12,14,16,18]\n",
    "p = 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ts_qaoa2 = []\n",
    "shots_qaoa2 = []\n",
    "nfevs_qaoa2 = []\n",
    "psucc_qaoa2 = []\n",
    "for qubit in qubits:\n",
    "    dir_path = f\"results/N{qubit}/qaoa_p2/\"\n",
    "    print(f\"directory: {dir_path}\")\n",
    "    p = []\n",
    "    t = []\n",
    "    s = []\n",
    "    it = []\n",
    "    # for each number of qubits\n",
    "    # we have several number of shots and iterations\n",
    "    for filename in sorted(os.listdir(dir_path)):\n",
    "        filename = dir_path + filename\n",
    "        with open(filename, \"r\") as file:\n",
    "            # print(filename)\n",
    "            data = json.load(file)\n",
    "        ever_found = []\n",
    "        # for each shot and iteration param\n",
    "        # we randomized the initial point\n",
    "        # to estimate the right probability\n",
    "        for run in data:\n",
    "            ever_found.append(run[\"ever_found\"])\n",
    "        # compute p('found minimum')\n",
    "        p.append(\n",
    "            np.mean(np.asarray(ever_found, dtype=np.float128))\n",
    "        )  # float128 to avoid to many zeros\n",
    "        # maxiter*shots = actual number of iteration\n",
    "        t.append(run[\"shots\"] * run[\"nfev\"])\n",
    "        s.append(run[\"shots\"])\n",
    "        it.append(run[\"nfev\"])\n",
    "    # update list for each number of qubits\n",
    "    ts_qaoa2.append(t)\n",
    "    psucc_qaoa2.append(p)\n",
    "    shots_qaoa2.append(s)\n",
    "    nfevs_qaoa2.append(it)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ts_qaoa3 = []\n",
    "shots_qaoa3 = []\n",
    "nfevs_qaoa3 = []\n",
    "psucc_qaoa3 = []\n",
    "for qubit in qubits:\n",
    "    dir_path = f\"results/N{qubit}/qaoa_p3/\"\n",
    "    print(f\"directory: {dir_path}\")\n",
    "    p = []\n",
    "    t = []\n",
    "    s = []\n",
    "    it = []\n",
    "    # for each number of qubits\n",
    "    # we have several number of shots and iterations\n",
    "    for filename in sorted(os.listdir(dir_path)):\n",
    "        filename = dir_path + filename\n",
    "        with open(filename, \"r\") as file:\n",
    "            # print(filename)\n",
    "            data = json.load(file)\n",
    "        ever_found = []\n",
    "        # for each shot and iteration param\n",
    "        # we randomized the initial point\n",
    "        # to estimate the right probability\n",
    "        for run in data:\n",
    "            ever_found.append(run[\"ever_found\"])\n",
    "        # compute p('found minimum')\n",
    "        p.append(\n",
    "            np.mean(np.asarray(ever_found, dtype=np.float128))\n",
    "        )  # float128 to avoid to many zeros\n",
    "        # maxiter*shots = actual number of iteration\n",
    "        t.append(run[\"shots\"] * run[\"nfev\"])\n",
    "        s.append(run[\"shots\"])\n",
    "        it.append(run[\"nfev\"])\n",
    "    # update list for each number of qubits\n",
    "    ts_qaoa3.append(t)\n",
    "    psucc_qaoa3.append(p)\n",
    "    shots_qaoa3.append(s)\n",
    "    nfevs_qaoa3.append(it)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "colors = ['lightgreen','forestgreen','tab:green','darkgreen', 'lightseagreen', 'tab:cyan','cornflowerblue' , 'blue', 'navy']\n",
    "plots1 = []\n",
    "plots2 = []\n",
    "\n",
    "for i, (p,t) in enumerate(zip(psucc_qaoa2, ts_qaoa2)):\n",
    "    if 4+i*2 > 18:\n",
    "        continue\n",
    "    try:\n",
    "        last_zero = np.argwhere(np.isclose(p, 0.))\n",
    "        print(f\"N={4+i*2}, p=2, t={np.asarray(t)[np.asarray(p)==0.].max()}, shots={np.asarray(shots_qaoa2[i])[np.asarray(p)==0.].max()}, nfev={np.asarray(nfevs_qaoa2[i])[np.asarray(p)==0.].max()}\")\n",
    "    except ValueError:\n",
    "        print(f\"No zeros in success probability for N={4+i*2}\")\n",
    "    # plot for each dimension\n",
    "    p_t = np.vstack((t, p))\n",
    "    p_t = p_t[:, np.argsort(p_t[0, :])]\n",
    "    plot1, = plt.plot(p_t[0]*2, p_t[1], \".\", label=f\"{4+i*2}\", color=colors[i], alpha=0.8)\n",
    "    plots1.append(plot1)\n",
    "\n",
    "for i, (p,t) in enumerate(zip(psucc_qaoa3, ts_qaoa3)):\n",
    "    try:\n",
    "        last_zero = np.argwhere(np.isclose(psucc_qaoa3[i], 0.))\n",
    "        print(f\"N={4+i*2}, p=2, t={np.asarray(ts_qaoa3[i])[np.asarray(psucc_qaoa3[i])==0.].max()}, shots={np.asarray(shots_qaoa3[i])[np.asarray(psucc_qaoa3[i])==0.].max()}, nfev={np.asarray(nfevs_qaoa3[i])[np.asarray(psucc_qaoa3[i])==0.].max()}\")\n",
    "    except ValueError:\n",
    "        print(f\"No zeros in succes probability for N={4+i*2}\")\n",
    "    p_t_qaoa3 = np.vstack((ts_qaoa3[i], psucc_qaoa3[i]))\n",
    "    p_t_qaoa3 = p_t_qaoa3[:, np.argsort(p_t_qaoa3[0, :])]\n",
    "    plot2, = plt.plot(p_t_qaoa3[0]*2, p_t_qaoa3[1], \".\", label=f\"{4+i*2}\", color=colors[i], fillstyle='none')\n",
    "    plots2.append(plot2)\n",
    "\n",
    "    plt.xscale('log')\n",
    "    plt.xlabel(r\"$t$ (shots x iterations x depth)\")\n",
    "    plt.ylabel(r\"$P$\")\n",
    "\n",
    "# Create a legend for the first line.\n",
    "first_legend = plt.legend(handles=plots1, loc='lower right', title=r\"QAOA $p=2$\")\n",
    "#Add the legend manually to the current Axes.\n",
    "ax = plt.gca().add_artist(first_legend)\n",
    "\n",
    "plt.title(r\"$p=2$\")\n",
    "#Create another legend for the second line.\n",
    "plt.legend(handles=plots2, loc='upper left', title=r\"QAOA $p=3$\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### Clean results\n",
    "\n",
    "Clean repeated $t$ and consider only the best one, i.e., the one with the high success porbability."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_psucc_qaoa2 = []\n",
    "new_ts_qaoa2 = []\n",
    "for i, (ps,ts) in enumerate(zip(psucc_qaoa2, ts_qaoa2)):\n",
    "    # plot for each dimension\n",
    "    p_t = np.vstack((ts, ps))\n",
    "    p_t = p_t[:, np.argsort(p_t[0, :])]\n",
    "    new_p = [p_t[1,0]]\n",
    "    new_t = [p_t[0,0]]\n",
    "    new_len = 0\n",
    "    for j, (t,p) in enumerate(p_t.T):\n",
    "        if new_t[new_len] == t:\n",
    "            if new_p[new_len] > p:\n",
    "                continue\n",
    "            new_p[new_len] = p\n",
    "        else:\n",
    "            new_t.append(t)\n",
    "            new_p.append(p)\n",
    "            new_len += 1\n",
    "    new_psucc_qaoa2.append(new_p)\n",
    "    new_ts_qaoa2.append(new_t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fun_exp(x, a, b):\n",
    "    return a*2**(b*x)\n",
    "\n",
    "def fun_lin(x, a, b):\n",
    "    return a + b*x\n",
    "    \n",
    "pt_all = []\n",
    "for i, (p,t) in enumerate(zip(psucc_qaoa2, ts_qaoa2)):\n",
    "    p_t = np.vstack((t, p))\n",
    "    p_t[0] *= 2 \n",
    "    p_t = p_t[:, np.argsort(p_t[0, :])]\n",
    "    pt_all.append(p_t)\n",
    "pt_all_2 = []\n",
    "for i, (p,t) in enumerate(zip(psucc_qaoa3, ts_qaoa3)):\n",
    "    p_t = np.vstack((t, p))\n",
    "    p_t[0] *= 3 \n",
    "    p_t = p_t[:, np.argsort(p_t[0, :])]\n",
    "    pt_all_2.append(p_t)\n",
    "\n",
    "\n",
    "colors = ['lightgreen', 'limegreen', 'forestgreen', 'darkgreen','lightseagreen', 'cornflowerblue' , 'blue', 'navy']\n",
    "min_psucc_lst = [0.99, 0.9, 0.5, 0.25]\n",
    "min_qubits = 4 \n",
    "max_qubits = 18\n",
    "plots1 = [] \n",
    "plots2 = []\n",
    "\n",
    "plt.plot(np.arange(min_qubits, max_qubits+1), fun_exp(np.arange(min_qubits, max_qubits+1), 1, 1), color='red', linewidth=3.5, linestyle=\"dashed\")\n",
    "\n",
    "for min_psucc, color in zip(min_psucc_lst, colors[::2]):\n",
    "    step_size = []\n",
    "    for i, pt in enumerate(pt_all):\n",
    "        mask = pt[1] >= min_psucc\n",
    "        # check if the p is not reached\n",
    "        if mask.sum() == 0:\n",
    "            continue\n",
    "        step_size.append(pt[0][mask][0])\n",
    "\n",
    "    x_axes = 4 + 2*np.arange(len(step_size))\n",
    "\n",
    "    popt, _ = curve_fit(fun_lin, x_axes, np.log2(step_size))\n",
    "    plot1, = plt.plot(x_axes, step_size, \"o\", label=f\"p={min_psucc}\", color=color)\n",
    "    plots1.append(plot1)\n",
    "    plot1, = plt.plot(np.arange(min_qubits, max_qubits+1), fun_exp(np.arange(min_qubits, max_qubits+1), 2**popt[0], popt[1]), \"-\", label=f\"a={2**popt[0]:.2f}, b={popt[1]:.2f}\", color=color)\n",
    "    plots1.append(plot1)\n",
    "\n",
    "    step_size_2 = []\n",
    "    for i, pt_2 in enumerate(pt_all_2):\n",
    "        mask = pt_2[1] >= min_psucc\n",
    "        # check if the p is not reached\n",
    "        if mask.sum() == 0:\n",
    "            continue\n",
    "        step_size_2.append(pt_2[0][mask][0])\n",
    "    x_axes_2 = 4 + 2*np.arange(len(step_size_2))\n",
    "    popt, _ = curve_fit(fun_lin, x_axes_2, np.log2(step_size_2))\n",
    "    plot2, = plt.plot(x_axes_2, step_size_2, \"o\", label=f\"p={min_psucc}\", color=color, fillstyle='none')\n",
    "    plots2.append(plot2)\n",
    "    plot2, = plt.plot(np.arange(min_qubits, max_qubits+1), fun_exp(np.arange(min_qubits, max_qubits+1), 2**popt[0], popt[1]), \"--\", label=f\"a={2**popt[0]:.2f}, b={popt[1]:.2f}\", color=color)\n",
    "    plots2.append(plot2)\n",
    "\n",
    "plt.yscale('log', base=2)\n",
    "plt.ylabel(r\"$t$ (shots x iterations)\")\n",
    "plt.xlabel(r\"$N$\")\n",
    "\n",
    "\n",
    "# Create a legend for the first line.\n",
    "first_legend = plt.legend(handles=plots1, loc='upper left', title=r\"QAOA $p=2$\", fontsize=\"x-small\")\n",
    "#Add the legend manually to the current Axes.\n",
    "ax = plt.gca().add_artist(first_legend)\n",
    "\n",
    "plt.title(r\"$p=2$\")\n",
    "#Create another legend for the second line.\n",
    "plt.legend(handles=plots2, loc='lower right', title=r\"QAOA $p=3$\", fontsize=\"x-small\")\n",
    "\n",
    "plt.title(r\"Fit with $a \\cdot 2^{b \\cdot N}$\")\n",
    "\n",
    "#plt.savefig(\"scale_qaoap2.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for k, (p_qaoa2, s_qaoa2, nfev_qaoa2, t_qaoa2, p_qaoa3, s_qaoa3, nfev_qaoa3, t_qaoa3) in enumerate(zip(psucc_qaoa2, shots_qaoa2, nfevs_qaoa2, ts_qaoa2, psucc_qaoa3, shots_qaoa3, nfevs_qaoa3, ts_qaoa3)):\n",
    "    len_iter_qaoa2 = np.max(nfev_qaoa2)//2+1\n",
    "    p_qaoa2 = np.reshape(p_qaoa2, (-1,len_iter_qaoa2))\n",
    "    s_qaoa2 = np.reshape(s_qaoa2, (-1,len_iter_qaoa2))\n",
    "    nfev_qaoa2 = np.reshape(nfev_qaoa2, (-1,len_iter_qaoa2))\n",
    "    t_qaoa2 = np.reshape(t_qaoa2, (-1,len_iter_qaoa2)) * 2\n",
    "\n",
    "    len_iter_qaoa3 = np.max(nfev_qaoa3)//2+1\n",
    "    p_qaoa3 = np.reshape(p_qaoa3, (-1,len_iter_qaoa3))\n",
    "    s_qaoa3 = np.reshape(s_qaoa3, (-1,len_iter_qaoa3))\n",
    "    nfev_qaoa3 = np.reshape(nfev_qaoa3, (-1,len_iter_qaoa3))\n",
    "    t_qaoa3 = np.reshape(t_qaoa3, (-1,len_iter_qaoa3)) * 3\n",
    "\n",
    "    fig, (ax1, ax2) = plt.subplots(1,2, figsize=(16,6))\n",
    "    for i in range(p_qaoa2.shape[0]):\n",
    "        ax1.plot(t_qaoa2[i,:], p_qaoa2[i,:], '--o', label=s_qaoa2[i,0])\n",
    "    ax1.plot(np.arange(0, 2**(8+k*2)), 1 - (1-1/(2**(4+k*2)))**np.arange(0, 2**(8+k*2)), '--', color=\"black\", label=\"Random sampling\")\n",
    "    for i in range(p_qaoa3.shape[0]):\n",
    "        ax2.plot(t_qaoa3[i,:], p_qaoa3[i,:], '--o', label=s_qaoa3[i,0])\n",
    "    ax2.plot(np.arange(0, 2**(8+k*2)), 1 - (1-1/(2**(4+k*2)))**np.arange(0, 2**(8+k*2)), '--', color=\"black\", label=\"Random sampling\")\n",
    "\n",
    "    ax1.set_xlabel('t (shots x iterations)')\n",
    "    ax1.set_ylabel(r'$F(success \\leq t)$')\n",
    "    ax1.set_xscale('log', base=2)\n",
    "    #ax1.set_yscale('log', base=2)\n",
    "    #ax1.set_ylim(0,1.1)\n",
    "\n",
    "    ax2.set_xlabel('t')\n",
    "    ax2.set_ylabel(r'$F(success \\leq t)$')\n",
    "    ax2.set_xscale('log', base=2)\n",
    "    #ax2.set_yscale('log', base=2)\n",
    "    #ax2.set_ylim(0,1.1)\n",
    "\n",
    "    ax1.legend(title='Shots')\n",
    "    ax2.legend(title='Shots')\n",
    "\n",
    "    ax1.set_title(f'QAOA p=2 N={4+k*2}')\n",
    "    ax2.set_title(f'QAOA p=3 N={4+k*2}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Test circuits with noise\n",
    "\n",
    "Implement RealAmplitude ansatz and look for reliable noise setting. I load the noise from a qiskit class, from the actual real divece and I compare them with the exact simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from math import pi\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from qiskit.circuit.library import  RealAmplitudes\n",
    "from qiskit.providers.fake_provider import FakeMumbaiV2, FakeVigoV2, FakeMumbai \n",
    "from qiskit.providers.aer import AerProvider, AerSimulator\n",
    "from qiskit_aer.noise import NoiseModel\n",
    "from qiskit import IBMQ, transpile\n",
    "from qiskit.visualization import plot_histogram\n",
    "from qiskit.providers.aer.noise import thermal_relaxation_error, depolarizing_error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### Fake Device\n",
    "\n",
    "Loads from qiskit class a shot of a real device."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mumbai has 27 qubots and quantum volume 128\n",
    "device_backend = FakeMumbai()\n",
    "noise_model = NoiseModel.from_backend(device_backend)\n",
    "sim_fakenoise = AerSimulator(\n",
    "            method='statevector',\n",
    "            max_parallel_threads=1,\n",
    "            noise_model = noise_model, \n",
    "        )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### Custom noise\n",
    "\n",
    "Implement custom noise in two different ways, setting damping time $T_1$ and dephasing time $T_2$ or depolarizing error $p_1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# T1 and T2 values for all qubits\n",
    "# T1s = np.random.normal(50e3, 10e3) # Sampled from normal distribution mean 50 microsec\n",
    "# T2s = np.random.normal(70e3, 10e3)  # Sampled from normal distribution mean 50 microsec\n",
    "\n",
    "t1 = 50e3\n",
    "t2 = 70e3\n",
    "\n",
    "# Truncate random T2s <= T1s\n",
    "#T2s = np.array([min(T2s[j], 2 * T1s[j]) for j in range(20)])\n",
    "\n",
    "# Instruction times (in nanoseconds)\n",
    "time_u1 = 0   # virtual gate\n",
    "time_u2 = 50  # (single X90 pulse)\n",
    "time_u3 = 100 # (two X90 pulses)\n",
    "time_cx = 300\n",
    "time_reset = 1000  # 1 microsecond\n",
    "time_measure = 1000 # 1 microsecond\n",
    "\n",
    "\n",
    "# QuantumError objects\n",
    "errors_reset = thermal_relaxation_error(t1, t2, time_reset)\n",
    "errors_measure = thermal_relaxation_error(t1, t2, time_measure)\n",
    "errors_u1  = thermal_relaxation_error(t1, t2, time_u1)\n",
    "errors_u2  = thermal_relaxation_error(t1, t2, time_u2)\n",
    "errors_u3  = thermal_relaxation_error(t1, t2, time_u3)\n",
    "errors_cx = thermal_relaxation_error(t1, t2, time_cx).expand(thermal_relaxation_error(t1, t2, time_cx))\n",
    "# Add errors to noise model\n",
    "noise_thermal = NoiseModel()\n",
    "noise_thermal.add_all_qubit_quantum_error(errors_reset, \"reset\")\n",
    "noise_thermal.add_all_qubit_quantum_error(errors_measure, \"measure\")\n",
    "noise_thermal.add_all_qubit_quantum_error(errors_u1, \"u1\")\n",
    "noise_thermal.add_all_qubit_quantum_error(errors_u2, \"u2\")\n",
    "noise_thermal.add_all_qubit_quantum_error(errors_u3, \"u3\")\n",
    "\n",
    "noise_thermal.add_all_qubit_quantum_error(errors_cx, \"cx\")\n",
    "\n",
    "\n",
    "sim_customnoise = AerSimulator(\n",
    "            method='statevector',\n",
    "            max_parallel_threads=1,\n",
    "            noise_model=noise_thermal, \n",
    "        )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# depolarizing error\n",
    "\n",
    "prob_1 = 1e-2 / 10. #1e-3 -> 1e-4 (best scenario)\n",
    "prob_2 = prob_1 * 10. ##1e-2 -> 1e-3 (best) \n",
    "error_1 = depolarizing_error(prob_1, 1)\n",
    "error_2 = depolarizing_error(prob_2, 2)\n",
    "\n",
    "# Add errors to noise model\n",
    "noise_model = NoiseModel()\n",
    "noise_model.add_all_qubit_quantum_error(error_1, ['u1', 'u2', 'u3', 'h'])\n",
    "noise_model.add_all_qubit_quantum_error(error_2, ['cx'])\n",
    "\n",
    "sim_customnoise2 = AerSimulator(\n",
    "            method='statevector',\n",
    "            max_parallel_threads=1,\n",
    "            noise_model=noise_model, \n",
    "        )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### Real device noise \n",
    "\n",
    "Load the noise on the fly from current online devices, the maximum qubits number is thus limited to 7. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get a real backend from a real provider\n",
    "provider = IBMQ.load_account()\n",
    "# define noise from online ibm computer\n",
    "backend = provider.get_backend('ibmq_manila')\n",
    "noise_model = NoiseModel.from_backend(backend)\n",
    "# generate a simulator that mimics the real quantum system with the latest calibration results\n",
    "sim_realnoise = AerSimulator(\n",
    "            method='statevector',\n",
    "            max_parallel_threads=1,\n",
    "            noise_model = noise_model, \n",
    "        )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### Exact simulation\n",
    "\n",
    "We compare the noisy simulator with the exact one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "exact_sim = AerSimulator(\n",
    "            method='statevector',\n",
    "            max_parallel_threads=1,\n",
    "        )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qubits = 4\n",
    "circ_depth = 2\n",
    "seed = 42\n",
    "\n",
    "\n",
    "# create the circuit\n",
    "# standard VQE ansatz\n",
    "qc = RealAmplitudes(\n",
    "    qubits,\n",
    "    reps=circ_depth,\n",
    "    insert_barriers=True,\n",
    "    entanglement=\"circular\",\n",
    ")\n",
    "# measure all qubits at the end of the circuit\n",
    "qc.measure_all()\n",
    "\n",
    "# define generator for initial point\n",
    "rng = np.random.default_rng(seed=seed)\n",
    "params =  np.ones(qc.num_parameters) * pi \n",
    "#params = np.random.uniform(-2 * pi, 2 * pi, qc.num_parameters)\n",
    "\n",
    "qc = qc.assign_parameters(params)\n",
    "\n",
    "print(qc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Execute the circuit with fixed params\n",
    "job = exact_sim.run(transpile(qc, exact_sim))\n",
    "# Grab results from the job\n",
    "result = job.result().get_counts()\n",
    "# return counts\n",
    "plot_histogram(result, title='Exact') #filename='exact_sim.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Execute the circuit with fixed params\n",
    "job = sim_realnoise.run(transpile(qc, sim_realnoise))\n",
    "# Grab results from the job\n",
    "result = job.result().get_counts()\n",
    "# return counts\n",
    "plot_histogram(result, title='Real device Manila')#filename='realnoise.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Execute the circuit with fixed params\n",
    "job = sim_fakenoise.run(transpile(qc, sim_fakenoise))\n",
    "# Grab results from the job\n",
    "result = job.result().get_counts()\n",
    "# return counts\n",
    "plot_histogram(result, title='Fake device Mumbai') #filename='fakenoise.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Execute the circuit with fixed params\n",
    "job = sim_customnoise2.run(transpile(qc, sim_customnoise2))\n",
    "# Grab results from the job\n",
    "result = job.result().get_counts()\n",
    "# return counts\n",
    "plot_histogram(result, title=f'Depolarizing noise {prob_1*1e2}%') #filename='customnoise2.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Execute the circuit with fixed params\n",
    "job = sim_customnoise.run(transpile(qc, sim_customnoise))\n",
    "# Grab results from the job\n",
    "result = job.result().get_counts()\n",
    "# return counts\n",
    "plot_histogram(result, title=f'Damping {t1//1e3}\\u03bcs  Dephasing {t2//1e3}\\u03bcs',)# filename='customnoise.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Use Pandas\n",
    "\n",
    "It's time to use Pandas."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import itertools\n",
    "import os\n",
    "import math\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import json\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sizes = np.arange(4,22,2)\n",
    "circ_depths = np.asarray([1,2,3,4])\n",
    "algs = np.asarray([\"qaoa\", \"qaoa+\", \"vqe\"])\n",
    "\n",
    "model_info = [(algo, p) for p in circ_depths for algo in algs] "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "run_info = {\"t\": [], \"shots\": [], \"iter\": [], \"p_succ\": [], \"qubits\": [], \"circ_depth\": [], \"algo\": [], \"succ_true\": [], \"succ_false\": [], 'p_succ_final': []}\n",
    "for n in sizes:\n",
    "    for (algo, p) in model_info:\n",
    "        dir_path = f\"results/N{n}/{algo}_p{p}/\"\n",
    "        if not os.path.isdir(dir_path):\n",
    "            print(f\"{dir_path} does not exist\")\n",
    "            continue\n",
    "        for filename in sorted(os.listdir(dir_path)):\n",
    "            filename = dir_path + filename\n",
    "            with open(filename, \"r\") as file:\n",
    "                # data are a list with #initial_points elements\n",
    "                runs = json.load(file)\n",
    "            ever_founds = []\n",
    "            success_final = []\n",
    "            iterations = [] \n",
    "            # for each shot and iteration param\n",
    "            # we randomized the initial point\n",
    "            # to estimate the right probability\n",
    "            for r in runs:\n",
    "                ever_founds.append(r[\"ever_found\"])\n",
    "                iterations.append(r[\"nfev\"])\n",
    "                success_final.append(r[\"success\"])\n",
    "            # compute p('found minimum')\n",
    "            #run_info[\"p_succ\"].append([ever_founds])\n",
    "            run_info[\"p_succ\"].append(\n",
    "                np.mean(np.asarray(ever_founds, dtype=np.float128))\n",
    "            )  # float128 to avoid to many zeros\n",
    "            run_info[\"p_succ_final\"].append(\n",
    "                np.mean(np.asarray(success_final, dtype=np.float128))\n",
    "            )  # float128 to avoid to many zeros\n",
    "            run_info[\"succ_true\"].append(np.count_nonzero(ever_founds))\n",
    "            run_info[\"succ_false\"].append(len(ever_founds) - np.count_nonzero(ever_founds))\n",
    "\n",
    "            run_info[\"iter\"].append(int(math.ceil(np.mean(np.asarray(iterations)))))\n",
    "            # maxiter*shots = actual number of iteration\n",
    "            run_info[\"t\"].append(int(r[\"shots\"] * run_info[\"iter\"][-1]))\n",
    "            run_info[\"shots\"].append(int(r[\"shots\"]))\n",
    "            run_info[\"qubits\"].append(int(n))\n",
    "            run_info[\"circ_depth\"].append(int(p))\n",
    "            run_info[\"algo\"].append(str(algo))\n",
    "        \n",
    "jobdf = pd.DataFrame.from_dict(run_info)\n",
    "jobdf.to_csv(f\"data_0pi_iter1.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "jobdf.head(10)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### VQE with noise\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import itertools\n",
    "import os\n",
    "import math\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import json\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sizes = np.arange(4,14,2)\n",
    "circ_depths = [2]\n",
    "algs = [\"vqe\"]\n",
    "\n",
    "model_info = [(n, algo, p) for p in circ_depths for algo in algs for n in sizes] "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "run_info = {\"t\": [], \"shots\": [], \"iter\": [], \"p_succ\": [], \"qubits\": [], \"circ_depth\": [], \"algo\": [], \"succ_true\": [], \"succ_false\": [], \"noise\": []}\n",
    "for (n, algo, p) in model_info:\n",
    "    dir_path = f\"results/N{n}/noise/\"\n",
    "    if not os.path.isdir(dir_path):\n",
    "        print(f\"{dir_path} does not exist\")\n",
    "        continue\n",
    "    for filename in sorted(os.listdir(dir_path)):\n",
    "        filename = dir_path + filename\n",
    "        with open(filename, \"r\") as file:\n",
    "            # data are a list with #initial_points elements\n",
    "            runs = json.load(file)\n",
    "        ever_founds = []\n",
    "        iterations = [] \n",
    "        # for each shot and iteration param\n",
    "        # we randomized the initial point\n",
    "        # to estimate the right probability\n",
    "        for r in runs:\n",
    "            ever_founds.append(r[\"ever_found\"])\n",
    "            iterations.append(r[\"nfev\"])\n",
    "        # compute p('found minimum')\n",
    "        #run_info[\"p_succ\"].append([ever_founds])\n",
    "        run_info[\"p_succ\"].append(\n",
    "            np.mean(np.asarray(ever_founds, dtype=np.float128))\n",
    "          )  # float128 to avoid to many zeros\n",
    "        run_info[\"succ_true\"].append(np.count_nonzero(ever_founds))\n",
    "        run_info[\"succ_false\"].append(len(ever_founds) - np.count_nonzero(ever_founds))\n",
    "\n",
    "        run_info[\"iter\"].append(int(math.ceil(np.mean(np.asarray(iterations)))))\n",
    "        # maxiter*shots = actual number of iteration\n",
    "        run_info[\"t\"].append(int(r[\"shots\"] * run_info[\"iter\"][-1]))\n",
    "        run_info[\"shots\"].append(int(r[\"shots\"]))\n",
    "        run_info[\"qubits\"].append(int(n))\n",
    "        run_info[\"circ_depth\"].append(int(p))\n",
    "        run_info[\"algo\"].append(str(algo))\n",
    "        run_info[\"noise\"].append(bool(True))\n",
    "        \n",
    "jobdf = pd.DataFrame.from_dict(run_info)\n",
    "jobdf.to_csv(f\"datanoise.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "jobdf.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for shot, data in jobdf[(jobdf['circ_depth'] == 2) & (jobdf['algo'] == 'qaoa')].groupby(\"shots\"):\n",
    "    print(f\"{shot} {len(data['p_succ'])}\", [p=='True' for p in data[\"p_succ\"]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(np.arange(0,2**(5+jobdf['qubits'][0]),2), 1 - (1-1/(2**jobdf['qubits'][0]))**np.arange(0,2**(5+jobdf['qubits'][0]),2), '--', color=\"black\", label=\"Random sampling\")\n",
    "for shot, data in jobdf[(jobdf['circ_depth'] == 2) & (jobdf['algo'] == 'qaoa')].groupby(\"shots\"):\n",
    "    if shot > 64:\n",
    "        continue\n",
    "    #plt.plot(data[\"t\"], data[\"p_succ\"].apply(np.mean, 1, dtype=np.float128), \"-o\", label=f\"{shot}\")\n",
    "    plt.plot(data[\"t\"]*jobdf['circ_depth'].mean(), data[\"p_succ\"], \"-o\", label=f\"{shot}\")\n",
    "    plt.xscale(\"log\", base=2)\n",
    "plt.title(f\"$N=8$ QAOA $p=2$\")\n",
    "plt.legend(title='Shots')\n",
    "\n",
    "plt.savefig(\"qaoa_ferro.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "jobdf.to_csv(\"data_n6.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Comparison for Spin Glasses\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import itertools\n",
    "import os\n",
    "import math\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import json\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sizes = [16] \n",
    "circ_depths = [1,2,3,4]\n",
    "algs = [\"qaoa\", \"qaoa+\", \"vqe\"]\n",
    "seeds = np.arange(1, 31)\n",
    "\n",
    "model_info = [(n, algo, p, seed) for p in circ_depths for algo in algs for seed in seeds for n in sizes]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "run_info = {\"t\": [], \"shots\": [], \"iter\": [], \"p_succ\": [], \"qubits\": [], \"circ_depth\": [], \"algo\": [], \"seed\": []}\n",
    "for (n, algo, p, seed) in model_info:\n",
    "    dir_path = f\"results/N{n}/spinglass/{algo}p{p}_{seed}/\"\n",
    "    if not os.path.isdir(dir_path):\n",
    "        #print(f\"{dir_path} does not exist\")\n",
    "        continue\n",
    "    for filename in sorted(os.listdir(dir_path)):\n",
    "        filename = dir_path + filename\n",
    "        with open(filename, \"r\") as file:\n",
    "            # data are a list with #initial_points elements\n",
    "            runs = json.load(file)\n",
    "        ever_founds = []\n",
    "        iterations = [] \n",
    "        # for each shot and iteration param\n",
    "        # we randomized the initial point\n",
    "        # to estimate the right probability\n",
    "        for r in runs:\n",
    "            ever_founds.append(bool(r[\"ever_found\"]))\n",
    "            iterations.append(r[\"nfev\"])\n",
    "        # compute p('found minimum')\n",
    "        #run_info[\"p_succ\"].append(ever_founds)\n",
    "        run_info[\"p_succ\"].append(\n",
    "            np.mean(ever_founds, dtype=np.float128 )\n",
    "        )  # float128 to avoid to many zeros\n",
    "        run_info[\"iter\"].append(int(math.ceil(np.mean(iterations))))\n",
    "        # maxiter*shots = actual number of iteration\n",
    "        run_info[\"t\"].append(r[\"shots\"] * run_info[\"iter\"][-1])\n",
    "        run_info[\"shots\"].append(r[\"shots\"])\n",
    "        run_info[\"qubits\"].append(n)\n",
    "        run_info[\"circ_depth\"].append(p)\n",
    "        run_info[\"algo\"].append(algo)\n",
    "        run_info[\"seed\"].append(seed)\n",
    "        \n",
    "jobs_df = pd.DataFrame(run_info)\n",
    "#jobs_df.to_csv(f\"data_n{n}_spin_glass.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "jobs_df.describe()\n",
    "jobs_df.to_csv(f\"dataspinglass-n16.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "jobs_df[(jobs_df['circ_depth'] == 2) & (jobs_df['qubits'] == 16)].head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_df = pd.read_csv(f\"data_spinglass-ALL.csv\", index_col=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "complete_df = pd.concat([all_df, jobs_df], ignore_index=True, verify_integrity=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "complete_df.drop('Unnamed: 0', axis=1, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "complete_df.to_csv(f\"dataspinglass-NEW.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#fig, (axa, axb, axc) = plt.subplots(1, 3, figsize=(20.4, 5), constrained_layout=True)\n",
    "for shot, data in jobs_df[jobs_df['algo'] == 'qaoa'].groupby('shots'):\n",
    "    print(data['p_succ'].groupby('seed'))\n",
    "    #print(data[\"p_succ\"])\n",
    "    break\n",
    "    # plt.xscale(\"log\")\n",
    "    # plt.title(f\"$N=6$ QAOA $p=2$\")\n",
    "    # plt.legend(title='Shots')\n",
    "#     axa.plot(data[\"t\"], data[\"p_succ\"].apply(np.mean, 1, dtype=np.float128), \"-o\", label=f\"{shot}\")\n",
    "#     axa.plot(data[\"t\"], data[\"p_succ\"].apply(np.mean, 1, dtype=np.float128), \"-o\", label=f\"{shot}\")\n",
    "#     axa.plot(data[\"t\"], data[\"p_succ\"].apply(np.mean, 1, dtype=np.float128), \"-o\", label=f\"{shot}\")\n",
    "    \n",
    "#     plt.xscale(\"log\")\n",
    "\n",
    "# axa.title(f\"$N={n}$ {algo} $p={data.loc[0, ['circ_depth']]}$\")\n",
    "# axb.title(f\"$N={n}$ {algo} $p={data.loc[0, ['circ_depth']]}$\")\n",
    "# axc.title(f\"$N={n}$ {algo} $p={data.loc[0, ['circ_depth']]}$\")\n",
    "\n",
    "# plt.legend(title='Shots')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for _, df in jobs_df.groupby([\"shots\", \"algo\"]):\n",
    "    print(df['p_succ'].mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Comparison $N=14$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import json\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 14 \n",
    "circ_depths = np.asarray([1,2,3])\n",
    "algs = np.asarray([\"qaoa\", \"qaoa+\", \"vqe\"])\n",
    "\n",
    "model_info = [(n, algo, p) for p in circ_depths for algo in algs]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "run_info = {\"t\": [], \"shots\": [], \"iter\": [], \"p_succ\": [], \"qubits\": [], \"circ_depth\": [], \"algo\": []}\n",
    "#for (n, algo, p) in model_info:\n",
    "\n",
    "for (n, algo, p) in model_info:\n",
    "    dir_path = f\"results/N{n}/{algo}_p{p}/\"\n",
    "    if not os.path.isdir(dir_path):\n",
    "        print(f\"{dir_path} does not exist\")\n",
    "        continue\n",
    "    for filename in sorted(os.listdir(dir_path)):\n",
    "        filename = dir_path + filename\n",
    "        with open(filename, \"r\") as file:\n",
    "            # data are a list with #initial_points elements\n",
    "            runs = json.load(file)\n",
    "        ever_founds = []\n",
    "        iterations = [] \n",
    "        # for each shot and iteration param\n",
    "        # we randomized the initial point\n",
    "        # to estimate the right probability\n",
    "        for r in runs:\n",
    "            ever_founds.append(r[\"ever_found\"])\n",
    "            iterations.append(r[\"nfev\"])\n",
    "        # compute p('found minimum')\n",
    "        #run_info[\"p_succ\"].append(ever_founds)\n",
    "        run_info[\"p_succ\"].append(\n",
    "            np.mean(ever_founds, dtype=np.float128 )\n",
    "        )  # float128 to avoid to many zeros\n",
    "        run_info[\"iter\"].append(np.mean(iterations))\n",
    "        # maxiter*shots = actual number of iteration\n",
    "        run_info[\"t\"].append(r[\"shots\"] * run_info[\"iter\"][-1])\n",
    "        run_info[\"shots\"].append(r[\"shots\"])\n",
    "        run_info[\"qubits\"].append(n)\n",
    "        run_info[\"circ_depth\"].append(p)\n",
    "        run_info[\"algo\"].append(algo)\n",
    "        \n",
    "jobs_df = pd.DataFrame(run_info)\n",
    "#jobs_df.to_csv(f\"data_n{n}_spin_glass.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "jobs_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for shot, data in jobs_df[(jobs_df['circ_depth'] == 1) & (jobs_df['algo'] == 'vqe')].groupby(\"shots\"):\n",
    "    plt.plot(data[\"t\"], data[\"p_succ\"], \"-o\", label=f\"{shot}\")\n",
    "    plt.plot()\n",
    "    plt.xscale(\"log\", base=2)\n",
    "\n",
    "x_axes = np.arange(0, 2**(3+n))\n",
    "plt.plot(x_axes, 1 - (1-1/(2**n))**x_axes, '--', color=\"black\", label=\"Random sampling\")\n",
    "plt.title(f\"$N=14$ VQE $p=1$\")\n",
    "plt.legend(title='Shots')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for shot, data in jobs_df[(jobs_df['circ_depth'] == 2) & (jobs_df['algo'] == 'vqe')].groupby(\"shots\"):\n",
    "    plt.plot(data[\"t\"]*2, data[\"p_succ\"], \"-o\", label=f\"{shot}\")\n",
    "    plt.plot()\n",
    "    plt.xscale(\"log\", base=2)\n",
    "\n",
    "x_axes = np.arange(0, 2**(3+n))\n",
    "plt.plot(x_axes, 1 - (1-1/(2**n))**x_axes, '--', color=\"black\", label=\"Random sampling\")\n",
    "plt.title(f\"$N=14$ VQE $p=2$\")\n",
    "plt.legend(title='Shots')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for shot, data in jobs_df[(jobs_df['circ_depth'] == 2) & (jobs_df['algo'] == 'qaoa')].groupby(\"shots\"):\n",
    "    plt.plot(data[\"t\"]*2, data[\"p_succ\"], \"-o\", label=f\"{shot}\")\n",
    "    plt.plot()\n",
    "    plt.xscale(\"log\", base=2)\n",
    "\n",
    "x_axes = np.arange(0, 2**(3+n))\n",
    "plt.plot(x_axes, 1 - (1-1/(2**n))**x_axes, '--', color=\"black\", label=\"Random sampling\")\n",
    "plt.title(f\"$N=14$ QAOA $p=2$\")\n",
    "plt.legend(title='Shots')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for shot, data in jobs_df[(jobs_df['circ_depth'] == 2) & (jobs_df['algo'] == 'qaoa+')].groupby(\"shots\"):\n",
    "    plt.plot(data[\"t\"]*2, data[\"p_succ\"], \"-o\", label=f\"{shot}\")\n",
    "    plt.plot()\n",
    "    plt.xscale(\"log\", base=2)\n",
    "\n",
    "x_axes = np.arange(0, 2**(3+n))\n",
    "plt.plot(x_axes, 1 - (1-1/(2**n))**x_axes, '--', color=\"black\", label=\"Random sampling\")\n",
    "plt.title(f\"$N=14$ QAOA+ $p=2$\")\n",
    "plt.legend(title='Shots')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "jobs_df[jobs_df[\"algo\"] == \"qaoa+\"].head(20)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# TODO\n",
    "\n",
    "- [x] nel plot dello scaling esponenziale della prob di successo inserire anche altre famiglie di p, i. e., p=[0.5, 0.6, 0.7, 0.8, 0.9]\n",
    "- [x] studiare i rumori mediani degli ultimi device di IBM e capire quali sono i valori, a quel punto inserire il rumore nel codice o con una classe prefissata capendo come estrarre i valori del rumore oppure settandoli a mano\n",
    "- [x] verificare la differenza tra i risultati del circuito rumoroso e quello del simulatore esatto per una taglia ragionevole e stampare i risultati per Guglielmo\n",
    "- [x] $p=2$ servono 10k per alcuni circuiti, da 14 in poi. VALUTARE BENE!! (non serve lo scaling si vede)\n",
    "- [x] d'ora in poi i tempi della cpu vanno normalizzati per numero di cpu e mediati per le iterazioni, l'obiettivo è avere un t/CPU medio\n",
    "- [x] provare modello frustrato $\\pm1$ (scritto ma ancora non ho fatto girare)\n",
    "- [x] inizializzare con hadamart H gate su QAOA-inspired\n",
    "- [ ] plot in cui faccio vedere lo shape dell'energia e, sopra, anche alcuni passi dell'ottimizzatore come istogramma degli output, così da visualizzare che se non becca il ground state all'inizio ha poche speranze di beccarlo alla fine.  \n",
    "- [x] QAOA $p=2$ pochi shots più iterazioni (N=14,16,18)\n",
    "- [ ] Spin glass (QUANTE REALIZZAZIONI di $J_{ij}$ 100?) solo su size $N=14, 16$ e VQE, QAOA $p=2$, e QAOA+. OBC connect.\n",
    "- [ ] Spin glass con campo locale random $h_i$ connectivity [?]. \n",
    "- [ ] Figure: fig1 shape E vs x e diversi circuiti; fig2 vqe vs qaoa $p=2$ con $F(success)$ due panel; fig3 scaling VQE (con e senza rumore) e QAOA $p=2$ due panel; fig4 spin glass\n",
    "- è\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Test QAOA exact\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import functools\n",
    "\n",
    "import numpy as np\n",
    "import qiskit\n",
    "import qiskit_aer\n",
    "from scipy import sparse\n",
    "\n",
    "from src import utils"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qubits = 4\n",
    "dim = 1\n",
    "type_ising = 'ferro'\n",
    "h_field = -0.05\n",
    "seed = 42\n",
    "verbose = 1\n",
    "circ_depth = 1\n",
    "type_ansatz = 'qaoa'\n",
    "initial_points = [1]\n",
    "backend = \"statevector\"\n",
    "measure = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# hamiltonian is defined with +\n",
    "# following http://spinglass.uni-bonn.de/ notation\n",
    "ising, global_min = utils.create_ising1d(qubits, dim, type_ising, h_field, seed)\n",
    "print(ising)\n",
    "print(f\"J: {ising.adj_dict}\\nh: {ising.h_field}\\n\")\n",
    "# create the quantum ansatz\n",
    "ansatz = utils.create_ansatz(\n",
    "    qubits,\n",
    "    circ_depth,\n",
    "    type_ansatz,\n",
    "    ising,\n",
    "    measure,\n",
    "    verbose=verbose,\n",
    ")\n",
    "# define generator for initial points\n",
    "rng = np.random.default_rng(seed=seed)\n",
    "# define different starting points\n",
    "if len(initial_points) == 1:\n",
    "    initial_points.insert(0, 0)\n",
    "thetas0 = utils.get_init_points(initial_points, ansatz.num_parameters, rng)[0]\n",
    "\n",
    "def kron(gate_lst):\n",
    "    return functools.reduce(np.kron, gate_lst)\n",
    "\n",
    "def pauli_z(i, n):\n",
    "    if i < 0 or i >= n or n < 1:\n",
    "        raise ValueError(\"Bad value of i and/or n.\")\n",
    "    pauli_z_lst = [\n",
    "        np.array([[1, 0], [0, -1]]) if j == i else np.eye(2) for j in range(n)\n",
    "    ]\n",
    "    return kron(pauli_z_lst)\n",
    "\n",
    "q_hamiltonian = sparse.coo_array(\n",
    "    sum(\n",
    "        [\n",
    "            ising.adj_matrix[i, j]\n",
    "            * pauli_z(i, ising.spins)\n",
    "            @ pauli_z(j, ising.spins)\n",
    "            for i in range(ising.spins)\n",
    "            for j in range(ising.spins)\n",
    "        ]\n",
    "    )\n",
    "    + sum([ising.h_field[i] * pauli_z(i, ising.spins) for i in range(ising.spins)])\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "simulator = qiskit_aer.AerSimulator(method=backend)\n",
    "\n",
    "def update_ansatz(ansatz, parameters: np.ndarray) -> qiskit.QuantumCircuit:\n",
    "    # assign current thetas to the parametric circuit\n",
    "    # in_place=False to not overwrite the parametric circuit\n",
    "    circuit = ansatz.assign_parameters(parameters)\n",
    "    return circuit\n",
    "\n",
    "circuit = update_ansatz(ansatz, thetas0)\n",
    "circuit.save_statevector()\n",
    "\n",
    "circuit = qiskit.transpile(circuit, simulator)\n",
    "\n",
    "job = simulator.run(circuit)\n",
    "# Grab results from the job\n",
    "result = job.result()\n",
    "result_statevector = result.get_statevector(circuit).data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = q_hamiltonian @ result_statevector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q_hamiltonian.dot(result_statevector).conjugate().dot(result_statevector)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(np.real(result_statevector)**2 + np.imag(result_statevector)**2).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.vdot(result_statevector,result_statevector)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exact simulation\n",
    "No shot noise\n",
    "\n",
    "\n",
    "### Test qiskit implementation\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import qiskit\n",
    "import qiskit_aer\n",
    "\n",
    "from src import utils, vqe\n",
    "from src import ising"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_qaoa_ansatz(\n",
    "    num_qubits: int,\n",
    "    circ_depth: int,\n",
    "    hamiltonian: ising.Ising,\n",
    "    increase_params: bool = False,\n",
    ") -> qiskit.QuantumCircuit:\n",
    "    \"\"\"Create a quantum ansatz inspired by the problem Hamiltonian.\n",
    "    Setting increase_params=False the ansatz becomes the typical one of QAOA, with 2*circ_depth parameters.\n",
    "    Otherwise, increase_params=True increases the parameters to (2*num_qubits - 1)*circ_depth.\n",
    "    See also https://qiskit.org/textbook/ch-applications/qaoa.html.\n",
    "\n",
    "    Args:\n",
    "        num_qubits (int): Number of qubits.\n",
    "        circ_depth (int): Depth of the circuit, minimum 1.\n",
    "        hamiltonian (ising.Ising): Problem Hamiltonian.\n",
    "        increase_params (bool, optional): Allows betas and gammas to change in each layer. Defaults to False.\n",
    "\n",
    "    Returns:\n",
    "        qiskit.QuantumCircuit: QAOA or QAOA+ ansatz.\n",
    "    \"\"\"\n",
    "    # define circuit\n",
    "    name = f\"QAOA+ p={circ_depth}\" if increase_params else f\"QAOA p={circ_depth}\"\n",
    "    qc = qiskit.QuantumCircuit(num_qubits, name=name)\n",
    "    # set initial entalgled state\n",
    "    for i in range(num_qubits):\n",
    "        qc.h(i)\n",
    "    if increase_params:\n",
    "        # create the parameters for the circuit\n",
    "        thetas = qiskit.circuit.ParameterVector(\n",
    "            \"$\\\\theta$\", (2 * num_qubits) * circ_depth\n",
    "        )\n",
    "    else:\n",
    "        # calssical QAOA ansatz\n",
    "        betas = qiskit.circuit.ParameterVector(\"$\\\\beta$\", circ_depth)\n",
    "        gammas = qiskit.circuit.ParameterVector(\"$\\\\gamma$\", circ_depth)\n",
    "    # add circ_depth layers\n",
    "    for i in range(circ_depth):\n",
    "        #add R_zz parametric gates\n",
    "        for j, j_coupling in enumerate(hamiltonian.adj_dict.values()):\n",
    "            if increase_params:\n",
    "                qc.rzz(\n",
    "                    j_coupling * 2*thetas[(2 * num_qubits) * i + j],\n",
    "                    j,\n",
    "                    (j + 1) % hamiltonian.spins,\n",
    "                )\n",
    "            else:\n",
    "                qc.rzz(j_coupling * 2 * gammas[i], j, (j + 1) % hamiltonian.spins)\n",
    "        qc.barrier()\n",
    "        #add Rz parametric gates\n",
    "        for j, h_j in enumerate(hamiltonian.h_field):\n",
    "            if increase_params:\n",
    "                qc.rz(h_j * 2 * thetas[(2 * num_qubits) * i + num_qubits - 1 + j], j)\n",
    "            else:\n",
    "                qc.rz(h_j * 2 * gammas[i], j)\n",
    "        # add Rx parametric gates\n",
    "        for j in range(num_qubits):\n",
    "            if increase_params:\n",
    "                qc.rx(-2 * thetas[(2 * num_qubits) * (i + 1) - 1], j)\n",
    "            else:\n",
    "                qc.rx(-2 * betas[i], j)\n",
    "        # do not put barrier in the last iteration\n",
    "        if i == circ_depth - 1:\n",
    "            continue\n",
    "    return qc\n",
    "\n",
    "\n",
    "def create_ansatz(\n",
    "    qubits: int,\n",
    "    circ_depth: int,\n",
    "    ansatz_type: str,\n",
    "    hamiltonian: ising.Ising,\n",
    "    measure: bool,\n",
    "    verbose: int = 0,\n",
    "):\n",
    "    if ansatz_type == \"vqe\":\n",
    "        # standard VQE ansatz\n",
    "        qc = qiskit.circuit.library.RealAmplitudes(\n",
    "            qubits,\n",
    "            reps=circ_depth,\n",
    "            insert_barriers=True,\n",
    "            entanglement=\"circular\",\n",
    "            name=f\"VQE p={circ_depth}\",\n",
    "        )\n",
    "    elif ansatz_type == \"qaoa\" or ansatz_type == \"qaoa+\":\n",
    "        increase_params = True if ansatz_type == \"qaoa+\" else False\n",
    "        qc = create_qaoa_ansatz(qubits, circ_depth, hamiltonian, increase_params)\n",
    "    else:\n",
    "        raise NotImplementedError(f\"Ansatz type {ansatz_type} not found\")\n",
    "    if measure:\n",
    "        qc.barrier()\n",
    "        qc.measure_all()\n",
    "    if verbose > 0 and qubits < 8:\n",
    "        print(qc)\n",
    "    return qc\n",
    "\n",
    "\n",
    "def update_ansatz(ansatz, parameters: np.ndarray) -> qiskit.QuantumCircuit:\n",
    "    # assign current thetas to the parametric circuit\n",
    "    # in_place=False to not overwrite the parametric circuit\n",
    "    circuit = ansatz.assign_parameters(parameters)\n",
    "    return circuit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_ising1d(\n",
    "    spins: int,\n",
    "    dim: int,\n",
    "    type_ising: str,\n",
    "    h_field: float,\n",
    "    seed: int,\n",
    "):\n",
    "    # hamiltonian is defined with +\n",
    "    # following http://spinglass.uni-bonn.de/ notation\n",
    "    # and D-Wave https://docs.dwavesys.com/docs/latest/c_gs_2.html#underlying-quantum-physics\n",
    "    J, h = utils.get_ising_params(spins, h_field, type_ising, seed)\n",
    "    adj_dict = {}\n",
    "    for i in range(spins):\n",
    "        if i == spins - 1:\n",
    "            continue\n",
    "        adj_dict[(i, (i + 1) % spins)] = J[i]\n",
    "    # class devoted to set the couplings and get the energy\n",
    "    model = ising.Ising(spins, dim=dim, adj_dict=adj_dict, h_field=h)\n",
    "    # exact enumeration would be too expensive for large sizes\n",
    "    if spins < 22 and (type_ising == \"binary\" or type_ising == \"random\"):\n",
    "        min_eng = utils.compute_ising_min(model)\n",
    "    elif type_ising == \"ferro\":\n",
    "        min_eng = model.energy(-np.ones(model.spins))\n",
    "    else:\n",
    "        min_eng = -np.inf\n",
    "    return model, min_eng"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qubits = 3\n",
    "dim = 1\n",
    "type_ising = 'ferro'\n",
    "h_field = 0.05\n",
    "seed = 42\n",
    "\n",
    "circ_depth = 1\n",
    "type_ansatz = 'qaoa'\n",
    "verbose = 1\n",
    "initial_points = [1]\n",
    "shots = [None]\n",
    "backend = 'automatic'\n",
    "measure = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# hamiltonian is defined with +\n",
    "# following http://spinglass.uni-bonn.de/ notation\n",
    "hamiltonian, global_min = create_ising1d(qubits, dim, type_ising, h_field, seed)\n",
    "print(hamiltonian)\n",
    "print(f\"J: {hamiltonian.adj_dict}\\nh: {hamiltonian.h_field}\\n\")\n",
    "# create the quantum ansatz\n",
    "ansatz = create_ansatz(\n",
    "    qubits,\n",
    "    circ_depth,\n",
    "    type_ansatz,\n",
    "    hamiltonian,\n",
    "    measure=measure,\n",
    "    verbose=verbose,\n",
    ")\n",
    "# define generator for initial points\n",
    "rng = np.random.default_rng(seed=seed)\n",
    "# define different starting points\n",
    "if len(initial_points) == 1:\n",
    "    initial_points.insert(0, 0)\n",
    "#thetas0 = utils.get_init_points(initial_points, ansatz.num_parameters, rng)[0]\n",
    "thetas0 = [1, 1]\n",
    "\n",
    "simulator = qiskit_aer.AerSimulator(method=backend)\n",
    "\n",
    "circuit = update_ansatz(ansatz, thetas0)\n",
    "circuit.save_statevector()\n",
    "\n",
    "circuit = qiskit.transpile(circuit, simulator)\n",
    "\n",
    "job = simulator.run(circuit)\n",
    "# Grab results from the job\n",
    "result = job.result()\n",
    "statevector = result.get_statevector(circuit)\n",
    "\n",
    "print(statevector)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### Simulation with Matrix\n",
    "\n",
    "Using qutip and Mazzola's code\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import qutip\n",
    "from qutip import qobj"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def my_kron(arr_list):\n",
    "    \"\"\"Takes a list of numpy arrays [A1, A2, ...] and computes their tensor product A1 (x) A2 (x) ...\"\"\"\n",
    "    if len(arr_list) == 1:\n",
    "        return arr_list[0]\n",
    "    else:\n",
    "        return np.kron(arr_list[0], my_kron(arr_list[1:]))\n",
    "def X(i, n):\n",
    "    if i<0 or i>=n or n<1:\n",
    "        raise ValueError('Bad value of i and/or n.')\n",
    "    X_list = [np.array([[0,1],[1,0]]) if j==i else np.eye(2) for j in range(n)]\n",
    "    return my_kron(X_list)\n",
    "def Y(i, n):\n",
    "    if i<0 or i>=n or n<1:\n",
    "        raise ValueError('Bad value of i and/or n.')\n",
    "    X_list = [np.array([[0,-1j],[1j,0]]) if j==i else np.eye(2) for j in range(n)]\n",
    "    return my_kron(X_list)\n",
    "def Z(i, n):\n",
    "    if i<0 or i>=n or n<1:\n",
    "        raise ValueError('Bad value of i and/or n.')\n",
    "    X_list = [np.array([[1,0],[0,-1]]) if j==i else np.eye(2) for j in range(n)]\n",
    "    return my_kron(X_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ansatz_matrixform_expectation(thetas, n_qubits, h_prob, h_mix, psi_in, n_blocks):\n",
    "    # initialize the circuit operations (in matrix form) - make the ansatz\n",
    "    thetas_counter = 0\n",
    "    psi_circuit = psi_in\n",
    "    for block in range(n_blocks):\n",
    "        temp = -1.j*(-1)*thetas[thetas_counter]*h_prob\n",
    "        psi_circuit = temp.expm()*psi_circuit\n",
    "        thetas_counter += 1\n",
    "        temp = -1.j*(-1)*thetas[thetas_counter]*h_mix\n",
    "        psi_circuit = temp.expm()*psi_circuit\n",
    "        thetas_counter += 1\n",
    "    # measure the expectation value\n",
    "    val = qutip.expect(-1*h_prob, psi_circuit)\n",
    "    \n",
    "    #print(\"C\", val)\n",
    "    \n",
    "    return val , psi_circuit\n",
    "\n",
    "def initstate_Had(L):\n",
    "  psi_initial = qutip.basis(2,0)+qutip.basis(2,1)\n",
    "  for i in range(L-1):\n",
    "      psi_initial = qutip.tensor(psi_initial, qutip.basis(2,0)+qutip.basis(2,1))\n",
    "  return psi_initial/psi_initial.norm()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Llist = [3]\n",
    "\n",
    "for L in Llist:\n",
    "    Mshots=16\n",
    "    Gamma=0\n",
    "    Jcoupl=1\n",
    "    h_Q = np.zeros(L)-0.05*np.ones(L)\n",
    "    J_Q = np.diag(np.ones(L-1), k=1)\n",
    "    #J_Q[0,L-1]=1 #PBC\n",
    "    LambdaReg=L+1\n",
    "    H_z = sum([J_Q[i,j]*Z(i,L) @ Z(j,L) for i in range(L) for j in range(L)]) + sum([h_Q[i]*Z(i,L) for i in range(L)])\n",
    "    H_x = sum([X(i,L) for i in range(L)])\n",
    "\n",
    "    #print(H_z)\n",
    "    H = -Gamma*H_x-Jcoupl*H_z\n",
    "\n",
    "    dims = [[2] * L, [2] * L]\n",
    "    mazzola_hamiltonian = qobj.Qobj(H, dims=dims)\n",
    "    h_prob = qobj.Qobj(H_z, dims=dims)\n",
    "    #print(h_prob.expm())\n",
    "    h_mix = qobj.Qobj(H_x, dims=dims)\n",
    "\n",
    "    thetas=[1,1] #[6.15816242, 2.378081]\n",
    "    n_blocks=1 \n",
    "\n",
    "    psi_in = initstate_Had(L)\n",
    "\n",
    "    value, psi = ansatz_matrixform_expectation(thetas, L, h_prob, h_mix, psi_in,n_blocks)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Guglielmo's loss {value}\")\n",
    "print(f\"Guglielmo's psi \\n {psi.data}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "loss = np.real(\n",
    "    hamiltonian.quantum_hamiltonian.dot(statevector)\n",
    "    .conjugate()\n",
    "    .dot(statevector)\n",
    ")\n",
    "print(f\"my loss {loss}\")\n",
    "print(f\"my psi \\n{statevector.data[::-1]}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"hamiltonian: {thetas} circuit:{thetas0}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "cvar",
   "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.10.8"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "b3c625637230a5b7ad254c94d5c84587d6b3e71b3826eabd556557ff2223426f"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}