ex1_XOR_3bits.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "8a4fb02f",
   "metadata": {},
   "source": [
    "### First double multiplexer, then single phase tuner\n",
    "Including even ports template and odd ports template, use odd port nets for test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "28ace8b7",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:absl:No GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)\n"
     ]
    }
   ],
   "source": [
    "import jax.numpy as jnp\n",
    "from jax import grad, jit, vmap, value_and_grad\n",
    "from jax import random\n",
    "from jax import custom_vjp\n",
    "\n",
    "from functools import partial\n",
    "import numpy as np\n",
    "\n",
    "key = random.PRNGKey(0)\n",
    "CTYPE = jnp.complex64\n",
    "RTYPE = jnp.float32"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e65fad7",
   "metadata": {},
   "source": [
    "### Helper functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "c7ffbce8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# this block won't be changed.\n",
    "# this functions is a helper function to get the location; won't be modified\n",
    "def get_ABCD_index(num_ports, parity):\n",
    "    if parity==\"even\":\n",
    "        ind = jnp.arange(0, num_ports-1, 2)\n",
    "    elif parity==\"odd\":\n",
    "        ind = jnp.arange(1, num_ports-1, 2)\n",
    "    else:\n",
    "        raise ValueError(\"parity should be `even` or `odd`, but get:\", parity)\n",
    "    ind_A = (ind, ind)\n",
    "    ind_D = (ind+1, ind+1)\n",
    "    ind_B = (ind, ind+1)\n",
    "    ind_C = (ind+1, ind)\n",
    "    return ind_A, ind_B, ind_C, ind_D\n",
    "\n",
    "# these two functions can be isolated; won't be modified\n",
    "def double_transfer(theta, r1, r2, a, dx):\n",
    "    theta = -(theta + dx)\n",
    "    denom = 1+a**2*jnp.exp(2*1j*theta)*r1**2-2*a*jnp.exp(1j*theta)*r1*r2\n",
    "    numA = r1+a**2*jnp.exp(2*1j*theta)*r1-a*jnp.exp(1j*theta)*r2-a*jnp.exp(1j*theta)*r1**2*r2\n",
    "    numB = -1j*a*jnp.exp(1j*theta)*(-1+r1**2)*jnp.sqrt(1-r2**2)\n",
    "    numC = -1j*a*jnp.exp(1j*theta)*(-1+r1**2)*jnp.sqrt(1-r2**2)\n",
    "    numD = r1+a**2*jnp.exp(2*1j*theta)*r1-a*jnp.exp(1j*theta)*r2-a*jnp.exp(1j*theta)*r1**2*r2\n",
    "    return numA/denom, numB/denom, numC/denom, numD/denom\n",
    "\n",
    "def single_transfer(phi, r, a):\n",
    "    return jnp.exp(1j*jnp.pi+1j*phi)*(a-r*jnp.exp(-1j*phi))/(1.-r*a*jnp.exp(1j*phi))\n",
    "\n",
    "# activation can be separated as well\n",
    "def ring_activation(g, phi_b, alpha, a_act, r_act, xout):\n",
    "    # ring\n",
    "    dphi = 0.5*g*jnp.abs(xout)**2 + 0.5*phi_b\n",
    "    t_ring = jnp.exp(1j*(np.pi+dphi))*(a_act-r_act*jnp.exp(1j*dphi))/(1.0-r_act*a_act*jnp.exp(1j*dphi))\n",
    "    xout = (1j * jnp.sqrt(1.0-alpha) * t_ring) * xout\n",
    "    return xout\n",
    "\n",
    "@partial(jit, static_argnums=(0,))\n",
    "def get_transfer_even(num_ports, dA, dB, dC, dD):\n",
    "    \"\"\"\n",
    "    Get the even depths transfer matrix.\n",
    "    \"\"\"\n",
    "    ind_ABCD = get_ABCD_index(num_ports, \"even\")\n",
    "    T = jnp.eye(num_ports, dtype=CTYPE)\n",
    "    T = T.at[ind_ABCD[0]].set(dA)\n",
    "    T = T.at[ind_ABCD[1]].set(dB)\n",
    "    T = T.at[ind_ABCD[2]].set(dC)\n",
    "    T = T.at[ind_ABCD[3]].set(dD)\n",
    "    return T\n",
    "\n",
    "@partial(jit, static_argnums=(0,))\n",
    "def get_transfer_odd(num_ports, dA, dB, dC, dD):\n",
    "    \"\"\"\n",
    "    Get the odd depths transfer matrix.\n",
    "    \"\"\"\n",
    "    ind_ABCD = get_ABCD_index(num_ports, \"odd\")\n",
    "    T = jnp.eye(num_ports, dtype=CTYPE)\n",
    "    T = T.at[ind_ABCD[0]].set(dA)\n",
    "    T = T.at[ind_ABCD[1]].set(dB)\n",
    "    T = T.at[ind_ABCD[2]].set(dC)\n",
    "    T = T.at[ind_ABCD[3]].set(dD)\n",
    "    return T\n",
    "\n",
    "double_transfer_even_jit = jit(double_transfer)\n",
    "single_transfer_even_jit = jit(single_transfer)\n",
    "double_transfer_odd_jit = jit(double_transfer)\n",
    "single_transfer_odd_jit = jit(single_transfer)\n",
    "ring_activation_jit = jit(ring_activation)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "befb4cd7",
   "metadata": {},
   "source": [
    "### Define network parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "f435adf2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# double ring parameters\n",
    "r1 = 0.99580889\n",
    "r2 = 0.99999208\n",
    "dx = 0.00133703\n",
    "a = 1\n",
    "\n",
    "# single ring parameters\n",
    "rsingle = r1\n",
    "asingle = a\n",
    "\n",
    "# activation function parameters\n",
    "g = np.pi\n",
    "phi_b = 0\n",
    "alpha = 0.1\n",
    "a_act = 0.9\n",
    "r_act = 0.9\n",
    "\n",
    "# initiate around the phase value\n",
    "around = 1e-2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "fbacdc99",
   "metadata": {},
   "outputs": [],
   "source": [
    "# depth = ports\n",
    "num_ports = 3\n",
    "num_depth = num_ports\n",
    "num_layers = 2\n",
    "\n",
    "# for even and odd depths (assume the first depth is even.)\n",
    "num_depth_even = int(jnp.ceil(num_depth/2))\n",
    "num_depth_odd = int(jnp.floor(num_depth/2))\n",
    "num_ports_even = int(jnp.floor(num_ports/2))\n",
    "num_ports_odd = int(jnp.ceil(num_ports/2) - 1)\n",
    "num_total_even = num_layers*num_depth_even*num_ports_even \n",
    "num_total_odd = num_layers*num_depth_odd*num_ports_odd\n",
    "\n",
    "# get the optimization parameters\n",
    "# once the structure is setup, the shape of transfer matrix is setup.\n",
    "key, subkey1, subkey2, subkey3, subkey4 = random.split(key, 5)\n",
    "\n",
    "# our parameters\n",
    "theta_even = random.uniform(subkey1, \n",
    "                            shape=(num_layers, num_depth_even, \n",
    "                                   num_ports_even),\n",
    "                           )*around\n",
    "theta_odd = random.uniform(subkey2, \n",
    "                           shape=(num_layers, num_depth_odd, \n",
    "                                  num_ports_odd),\n",
    "                          )*around\n",
    "\n",
    "phi_even = random.uniform(subkey3, \n",
    "                          shape=(num_layers, num_depth_even, \n",
    "                                 num_ports_even),\n",
    "                         )*around\n",
    "phi_odd = random.uniform(subkey4, \n",
    "                         shape=(num_layers, num_depth_odd, \n",
    "                                num_ports_odd),\n",
    "                        )*around\n",
    "\n",
    "params = {}\n",
    "params['theta_even'] = theta_even\n",
    "params['theta_odd'] = theta_odd\n",
    "params['phi_even'] = phi_even\n",
    "params['phi_odd'] = phi_odd\n",
    "init_params = params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "2e0af520",
   "metadata": {},
   "outputs": [],
   "source": [
    "jit_matmul = jit(jnp.matmul)\n",
    "\n",
    "# let's focus on the double ring only first\n",
    "@jit\n",
    "def clements_evenports(r1, r2, a, dx, rsingle, asingle, \n",
    "                       theta_even, theta_odd, phi_even, phi_odd, xin):\n",
    "    \"\"\"\n",
    "    1-layer clements setup.\n",
    "    theta_even: (num_depth, num_even_devices)\n",
    "    theta_odd: (num_depth, num_odd_devices)\n",
    "    phi_even: (num_depth, num_even_devices)\n",
    "    phi_odd: (num_depth, num_even_devices)\n",
    "    xin: (x_vec_shape, batch_size)\n",
    "    \"\"\"\n",
    "    dAeven, dBeven, dCeven, dDeven = double_transfer_even_jit(theta_even, r1, r2, a, dx)\n",
    "    dAodd, dBodd, dCodd, dDodd = double_transfer_odd_jit(theta_odd, r1, r2, a, dx)\n",
    "    phase_even = single_transfer_even_jit(phi_even, rsingle, asingle)\n",
    "    phase_odd = single_transfer_odd_jit(phi_odd, rsingle, asingle)\n",
    "    # if the ports number is even\n",
    "    for d in range(int(num_depth/2)):\n",
    "        T_even = get_transfer_even(num_ports, \n",
    "                                   dAeven[d]*phase_even[d], dBeven[d]*phase_even[d], \n",
    "                                   dCeven[d], dDeven[d])\n",
    "        T_odd = get_transfer_odd(num_ports, \n",
    "                                 dAodd[d]*phase_odd[d], dBodd[d]*phase_odd[d], \n",
    "                                 dCodd[d], dDodd[d])\n",
    "        xin = jit_matmul(T_odd, jit_matmul(T_even, xin))\n",
    "    return xin\n",
    "\n",
    "@jit\n",
    "def clements_oddports(r1, r2, a, dx, rsingle, asingle, \n",
    "                      theta_even, theta_odd, phi_even, phi_odd, xin):\n",
    "    \"\"\"\n",
    "    1-layer clements setup.\n",
    "    theta_even: (num_depth, num_even_devices)\n",
    "    theta_odd: (num_depth, num_odd_devices)\n",
    "    phi_even: (num_depth, num_even_devices)\n",
    "    phi_odd: (num_depth, num_even_devices)\n",
    "    xin: (x_vec_shape, batch_size)\n",
    "    \"\"\"\n",
    "    dAeven, dBeven, dCeven, dDeven = double_transfer_even_jit(theta_even, r1, r2, a, dx)\n",
    "    dAodd, dBodd, dCodd, dDodd = double_transfer_odd_jit(theta_odd, r1, r2, a, dx)\n",
    "    phase_even = single_transfer_even_jit(phi_even, rsingle, asingle)\n",
    "    phase_odd = single_transfer_odd_jit(phi_odd, rsingle, asingle)\n",
    "    \n",
    "    # if the ports number is odd\n",
    "    for d in range(int(num_depth/2)):\n",
    "        T_even = get_transfer_even(num_ports, \n",
    "                                   dAeven[d]*phase_even[d], dBeven[d]*phase_even[d], \n",
    "                                   dCeven[d], dDeven[d])\n",
    "        T_odd = get_transfer_odd(num_ports, \n",
    "                                 dAodd[d]*phase_odd[d], dBodd[d]*phase_odd[d], \n",
    "                                 dCodd[d], dDodd[d])\n",
    "        xin = jit_matmul(T_odd, jit_matmul(T_even, xin))\n",
    "    T_even = get_transfer_even(num_ports, \n",
    "                               dAeven[d+1]*phase_even[d+1], dBeven[d+1]*phase_even[d+1], \n",
    "                               dCeven[d+1], dDeven[d+1])\n",
    "    xin = jit_matmul(T_even, xin)\n",
    "    return xin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "38c8b997",
   "metadata": {},
   "outputs": [],
   "source": [
    "@jit\n",
    "def network_even(theta_even, theta_odd, phi_even, phi_odd, xin):\n",
    "    xin = clements_evenports(r1, r2, a, dx, rsingle, asingle, \n",
    "                              theta_even[0], theta_odd[0], \n",
    "                              phi_even[0], phi_odd[0], xin)\n",
    "    xin = ring_activation(g, phi_b, alpha, a_act, r_act, xin)\n",
    "    for l in range(1, num_layers):\n",
    "        xin = clements_evenports(r1, r2, a, dx, rsingle, asingle, \n",
    "                                  theta_even[l], theta_odd[l], \n",
    "                                  phi_even[l], phi_odd[l], xin)\n",
    "        xin = ring_activation(g, phi_b, alpha, a_act, r_act, xin)\n",
    "    \n",
    "    xin = jnp.abs(xin)**2\n",
    "    return xin\n",
    "\n",
    "@jit\n",
    "def network_odd(theta_even, theta_odd, phi_even, phi_odd, xin):\n",
    "    xin = clements_oddports(r1, r2, a, dx, rsingle, asingle, \n",
    "                             theta_even[0], theta_odd[0], \n",
    "                             phi_even[0], phi_odd[0], xin)\n",
    "    xin = ring_activation(g, phi_b, alpha, a_act, r_act, xin)\n",
    "    for l in range(1, num_layers):\n",
    "        xin = clements_oddports(r1, r2, a, dx, rsingle, asingle, \n",
    "                                 theta_even[l], theta_odd[l], \n",
    "                                 phi_even[l], phi_odd[l], xin)\n",
    "        xin = ring_activation(g, phi_b, alpha, a_act, r_act, xin)\n",
    "    \n",
    "    xin = jnp.abs(xin)**2\n",
    "    return xin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "b808fc5d",
   "metadata": {},
   "outputs": [],
   "source": [
    "@jit\n",
    "def loss(params, xin, Y):\n",
    "    xout = network_odd(params['theta_even'], params['theta_odd'], \n",
    "                       params['phi_even'], params['phi_odd'], xin)\n",
    "    l = jnp.abs(xout[0,:]-Y[0,:])**2\n",
    "    return jnp.sum(l)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "60a4ee9a",
   "metadata": {},
   "source": [
    "### XOR dataset generation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "b51b8c20",
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_data_cascade_XOR(N, y_scale=1.0, normalize=False):\n",
    "    if N > 8:\n",
    "        print(\"we limit the total xor bits smaller than 8, but got N = :\", N)\n",
    "        raise\n",
    "    Y = np.zeros((1, 2**N))\n",
    "    number_array = np.arange(2**N, dtype=np.uint8)[:,None]\n",
    "    bit_array    = np.unpackbits(number_array, axis=1)\n",
    "    inds = np.array(np.count_nonzero(bit_array, axis=1) % 2,\n",
    "                    dtype=bool)\n",
    "    X = bit_array.astype(np.complex)[:, (8-N):].T\n",
    "    Y[0,inds] = y_scale\n",
    "    X = X + 1e-9\n",
    "    if normalize:\n",
    "        X = X/np.linalg.norm(X, axis=0)\n",
    "        X[:,0] = 1e-9\n",
    "    return jnp.array(X), jnp.array(Y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "99774861",
   "metadata": {},
   "outputs": [],
   "source": [
    "num_bits = 3\n",
    "y_max = 0.2\n",
    "X, Y = generate_data_cascade_XOR(num_bits, y_scale=y_max, normalize=True)\n",
    "X_padded = jnp.pad(X, ((0,num_ports-num_bits),(0,0)), mode='constant', constant_values=1)\n",
    "# X_padded shape: (num_ports, num_samples)\n",
    "# l_test = loss(params, X_padded, Y)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1a3dd55",
   "metadata": {},
   "source": [
    "### Optimizer + train"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "bec76c8f",
   "metadata": {},
   "outputs": [],
   "source": [
    "from jax.experimental import optimizers\n",
    "import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "56a84744",
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_xor(x, y, params, num_epochs, lr=1e-3):\n",
    "    jit_value_grad  = jit(value_and_grad(loss, (0, )))\n",
    "    \n",
    "    @jit\n",
    "    def update(params, x, y, opt_state):\n",
    "        value, grads = jit_value_grad(params, x, y)\n",
    "        opt_state = opt_update(0, grads[0], opt_state)\n",
    "        return get_params(opt_state), opt_state, value\n",
    "\n",
    "    # Defining an optimizer in Jax\n",
    "    opt_init, opt_update, get_params = optimizers.adam(lr)\n",
    "    opt_state = opt_init(params)\n",
    "    \n",
    "    loss_list = []\n",
    "    params, opt_state, l = update(params, x, y, opt_state)\n",
    "    loss_list.append(l)\n",
    "    \n",
    "    for epoch in range(num_epochs):\n",
    "        params, opt_state, l = update(params, x, y, opt_state)\n",
    "        loss_list.append(l)\n",
    "    return loss_list, params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "0d28a2f9",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "training time:  4.255724191665649\n"
     ]
    }
   ],
   "source": [
    "num_epochs = 2000\n",
    "lr = 1e-4\n",
    "start_time = time.time()\n",
    "loss_list, trained_params = run_xor(X_padded, Y, params, num_epochs, lr)\n",
    "print(\"training time: \", time.time()-start_time)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ed4fe09",
   "metadata": {},
   "source": [
    "### Analyze results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "0007dff2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "b7d9d332",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f820dfa32b0>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUUAAADQCAYAAACQqHMZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAgCElEQVR4nO3dd3yV5f3/8dcnJ4sRgoGA7AAJSwGBiBNFcKCCe2CXfrXwxTrab7VtrFbF2jrqqi2/KlaLq1AHVhBkiCgiMmXPBAgSRsJMQva4fn+cO+EkJOGcnHWfcz7PxyMPDlfu+74+9znJO/e8bjHGoJRSyikq2AUopZSdaCgqpZQLDUWllHKhoaiUUi40FJVSyoWGolJKuYgOdgENEZFxwLiEhIQJffr0CXY5Sqkws2bNmsPGmOSGvid2vk4xPT3drF69OthlKKXCjIisMcakN/Q93X1WSikXGopKKeVCQ1EppVxoKCqllIuwCcWD+aX85sP1rNt7PNilKKVCWNiEYlF5JR+uyWHPkaJgl6KUCmG2DEURGSciU/Pz892eJy7auSplFdX+KkspFQFsGYrGmNnGmImJiYluzxNbE4pVGopKqeazZSg2R1y0A4DySg1FpVTzhVEoWluKlVVBrkQpFcrCJhRjHc5V0S1FpZQ3wiYUo6KEGIdQpqGolPJC2IQiOLcWdUtRKeWNsArFuBiHhqJSyithFYqxjig90aKU8kpYhWJcjO4+K6W8E1ah6NxS1FBUSjVfwEJRRHqJyJsi8pG/+tAtRaWUt9wKRRF5S0TyRGRTvfYxIrJdRLJEJKOpZRhjdhlj7vGm2NOJdURRrrf5KaW84O6Dq6YBfwfeqWkQEQcwBbgCyAFWicgswAE8U2/+u40xeV5Xexqx0VE6IIRSyituhaIxZomIpNRrHg5kGWN2AYjIDOB6Y8wzwNjmFiQiE4GJAN27d/do3rhoB8dLKprbtVJKeXVMsQuw1+X/OVZbg0SknYi8BgwRkUcam84YM9UYk26MSU9ObvAJhI0qLq9k/d7jvLV0t0fzKaVUDW9CURpoa/R5qcaYI8aYScaY3tbWpM9t2lcAwAsLtvtj8UqpCOBNKOYA3Vz+3xXY7105Ts0ZZBbg5yN6AhAd1VBeK6XU6XkTiquANBHpKSKxwHhgli+Kas4gswD/d3kf7h3Zm4LSSt79LtsXpSilIoy7l+RMB74D+opIjojcY4ypBO4H5gNbgQ+MMZt9UVRztxSjooTeya0B+MOnPilFKRVh3D37fEcj7XOBuT6tyLnc2cDs9PT0CZ7O27djQu3r/OIKElvG+LI0pVSYC6vb/AAGdk3krgtTANi837MtTaWUsmUoNnf3ucaDo9MA+MuC7RjT6AlxpZQ6hS1DsbknWmoktYrlpqFdWPvDcTbvL/BxdUqpcGbLUPSFx8cOICEumjvfWsns9fsprdBxFpVSpyd23L0UkXHAuNTU1AmZmZnNXs7XOw4xedZmdh0uAqBzYjy9O7Smd3Jreie34ooBZ3JmYryPqlZKhQoRWWOMSW/we3YMxRrp6elm9erVXi2jrLKKhVty2XWoiN2Hi9h56AQ7805QVF6FI0r4841nc/u5nt1jrZQKbU2Foruj5ISsuGgHYwd1rtNmjGH34SJ+89EGnp+3nZuGdiXGEbZHEpRSHojIJBAReiW35t5Le3OkqJxvMg8FuySllE3YMhS9vSTHXZf0SaZFjIMlOw77tR+lVOiwZSh6e0mOu2Kjoxjaoy0rdx/1az9KqdBhy1AMpOEp7dh6sID8Yh2cVimloch5vZIwBmas+iHYpSilbMCWoRioY4oAQ7q3BWB9znG/96WUsj9bhmKgjimC85KdW4Z15Zsdhykp17telIp0tgzFQLt5aFcKyyqZuTYn2KUopYJMQxE4v1cSAI9+suk0Uyqlwp2GIs6LuUf36wBAbkFpkKtRSgWThqLlF5f1BuBNfTyqUhHNlqEYyLPPNYZ2PwOAqUt2UVFVHbB+lVL2YstQDOTZ5xoiwgu3DgZg2c4jAetXKWUvtgzFYLl2YCeSWsUy5cusYJeilAoSDUUXLWId3Htpb1ZmH+WlhTv0+S5KRSANxXpuTe8KwKuLMpm+cm+Qq1FKBZqGYj1tW8byzW8vA+D3n2xkxkq9J1qpSKKh2IBuSS356uGRAGTM3Mj/vuvdIxGUUqFDQ7ERKe1bsfWpMbSOi2b+5lx+/rYGo1KRwJahGIzrFBvSItbBskdGAfDF1lzmbToQ1HqUUv5ny1AMxnWKjWkTH1N7jPHxTzdzvLg8yBUppfzJlqFoN92SWvLvCeeRV1jG7a8vD3Y5Sik/0lB004W925MQH8323EKW79I7XpQKVxqKHlj6u1F0aduC8VOX62g6SoUpDUUPJLaI4c83DQRgxHOLg1yNUsofNBQ9dGmfZFLataS8qpqP1uhI3UqFGw3FZpj3q0tI69Ca38/cyP7jJcEuRynlQxqKzRAf4+Cl286horqaJ2ZtprpaB45QKlxoKDbTwK6JPDAqjYVbcpm5dl+wy1FK+UhAQ1FEbhCRN0TkUxG5MpB9+8OvRqcxtHtbnpy1mUOFZcEuRynlA26Hooi8JSJ5IrKpXvsYEdkuIlkiktHUMowx/zXGTADuAm5vVsU2EhUlPDHuLE6UVfLSwh3BLkcp5QOebClOA8a4NoiIA5gCXA0MAO4QkQEiMlBEPqv31cFl1ses+ULe4G5tGX9uN6av/IHdh4uCXY5Syktuh6IxZglwtF7zcCDLGLPLGFMOzACuN8ZsNMaMrfeVJ07PAZ8bY7733WoE16+v6EN0lPC3RZnBLkUp5SVvjyl2AVyHp86x2hrzAHA5cIuITGpoAhGZKCKrRWT1oUOHvCwvMDq0ieeuC1OYuXafDkqrVIjzNhSlgbZGr08xxrxqjBlmjJlkjHmtkWmmGmPSjTHpycnJXpYXOA9d2ReA5+Zto7C0IsjVKKWay9tQzAG6ufy/K7Dfy2XaZjxFT7SIdTDr/os4VlzBwCcX6EOvlApR3obiKiBNRHqKSCwwHpjlbVF2Gk/RE4O6tuXsLm0AuO/fYXPIVKmI4sklOdOB74C+IpIjIvcYYyqB+4H5wFbgA2PMZv+UGho++cVFJMRHM3fjQRZsPhjscpRSHhI77uaJyDhgXGpq6oTMzNA7o7vtYAFjXvkGgDkPXsxZnUNri1epcCcia4wx6Q19z5a3+YXq7nONfme2YepPhwFw7atL+WCVPj9aqVBhy1AMxRMt9V151pk8d7Nz7MXffryBlIw5VFZVB7kqpdTp2DIUQ31Lscbt53bn5dsH1/5/xPOLKausCmJFSqnTsWUohpMbh3Tly4cuBeBAfil9H5tHSsYcfc6LUjZly1AMh91nV72SW7Ptj2NIToirbRs/dTmvLsokv1gv9FbKTmx59rlGenq6Wb16dbDL8Km8wlI+XrOP5+Ztq207N+UMnrr+bPp3ahPEypSKHE2dfdZQDJLi8kr+PHcr7y0/ea/0iLT2XDe4M9EO4dI+HUhqFRvECpUKXxqKNlZUVsmX2/L417e72X6wkKLykydihqckMaR7W64Z2Il+nRKIi3YEsVKlwkfIhWKoX7zdXNXVhre/y+aTtfvIKyjjYL1nS5/fK4mjReUM65FEt6QWjO7Xkb5nJgSpWqVCV8iFYo1I2FJsSlbeCT7bsJ/KKsPWAwV8k3WY8sq61zp2aduC9glxRAncNzKV0f07UG3AEdXQAEZKKdBQDBvlldV8tT2Pti1jyTlWzBdbc1mVfazR58NcO6gTvdq34kB+KbcM60r/Tm1IbBET4KqVsh8NxTC3ZX8By3YeZv/xUjbkHGf1nmNNTh8XHcUvRqYybnAnUtq1Ikq3KlWE0VCMMMYYDhWW0a51HPuPl/DVjkMs2XGIrLwTDT5HJtYRxS8vT+P7Pce4oHc77rowhRcX7qBz2xb89PweQVgDpfwr5EIxUk+0BEplVTWLtx9iwjun/4MTY10e9ODoVPp0TCA+Rs+Aq9AXcqFYQ7cUAyO/pIKNOfmcKKtkzZ6jHMgv5bMNB5qcxxEljB3Uiapqw+TrzqJd67gmp1fKTjQUVbPlFZaSlXuCpz7bwraDhQ1Ok5wQx+NjBzCke1u6ntEywBUq5TkNReUT5ZXVVFRV8+7yPXybdZilWYdp7Mfn3JQzeGBUGpf0CZ2Hj6nIoaGo/CqvsJQnZ21m7sZTH78QHSUsyxhFhzbxQahMqYZpKKqAOVpUznvL9/D+ij3kFpy8frJ7Ukt6Jbfid2P66cAXKuhCLhT17HP4+HJbLndPq/uH7eLU9rz38/OCVJFSIRiKNXRLMXwUllbw7vI9PD9ve21bjEP48qGRdEvSkzMqsELuwVUq/CTEx/CLkalsfWoMdwzvDkBFlWHE84v526JMSiv0MQ3KHnRLUQWFMYbXl+zi2c9PDrZ745AuvHz7OcErSkUM3VJUtiMiTLq0N8sfGV172c4na/eRkjGHf3y1k6KyyiBXqCKVbikqW8jKO8HlL31dp+2yvsn8YewAeiW3DlJVKlzplqKyvdQOrcl+9loeHJ1W27Z4+yFGvfg1H6zaG8TKVKTRUFS28usr+rD7mWu477LetW2//XgD//efdcErSkUUW4ZiuD3iVHlGRPjNVf3I/NPVjBvcGTh5vHHTPv2ZUP6lxxSV7eUXV3DBs4soth7q9e8J53Fh7/ZBrkqFMj2mqEJaYssYNk++ioev7APAj95YQUrGnCBXpcKVhqIKCSLC/aPS+PjeC2rbXlq4I4gVqXCloahCyrAeSSzLGAXAq4sy+fuXem+88i0NRRVyOrdtwRe/vhSAFxbs0JMvyqc0FFVISu3QmjfvdB4nH/u3pSzLOhzkilS40FBUIWt0/45c3r8jAD/65wpKynVQCeU9DUUV0v55ZzoO67nV/R+fF+RqVDjQUFQhb/Pkq2pfL83U3WjlnYCFooj0F5HXROQjEbk3UP2q8Bcf42DW/RcB8JM3V/CyXqqjvOBWKIrIWyKSJyKb6rWPEZHtIpIlIhlNLcMYs9UYMwm4DWjwSnKlmmtQ17b85Hzn4LV/XaSX6ajmc3dLcRowxrVBRBzAFOBqYABwh4gMEJGBIvJZva8O1jzXAUuBRT5bA6UsT98wsPb1IzM3BLESFcrcCkVjzBLgaL3m4UCWMWaXMaYcmAFcb4zZaIwZW+8rz1rOLGPMhcCPfbkSStVY+ehoAKav3Mun6/YFuRoVirw5ptgFcB3oLsdqa5CIjBSRV0XkdWBuE9NNFJHVIrL60KFDXpSnIlGHhHgeGJUKwC9nrMPOA54oe/ImFKWBtkZ/Ao0xXxljHjTG/K8xZkoT0001xqQbY9KTk5O9KE9Fqoeu7MvQ7m0BPb6oPOdNKOYA3Vz+3xXY7105TjqeovLW9InnIwKvfJHJniNFwS5HhRBvQnEVkCYiPUUkFhgPzPJFUcaY2caYiYmJib5YnIpAcdEO/jPROaLOmFe+CXI1KpS4e0nOdOA7oK+I5IjIPcaYSuB+YD6wFfjAGLPZF0XplqLyheE9kwAoqagiJWMOhwrLglyRCgU68rYKa65PCRzW4ww+mnQBIg0dDleRREfeVhErtUNrXr59MABr9hzj719mBbkiZXe2DEXdfVa+dOOQrtw70vl0wBcX7uD6Kd9SVW3fPSQVXLYMRT3Ronztd2P68fOLewKwfu9x9h8vIa+wlIqq6iBXpuzGlqGolD88NnZA7esRzy9m+J8W8cjMjUGsSNmRLUNRd5+Vv7g++ArgozU5QapE2ZUtQ1F3n5W/DOuRxIQRPYNdhrIxW4aiUv70+2v68+Gkk1uMgycv4ERZZRArUnaioagijohwbkoSr9x+DgD5JRWc/cR89h0vISvvRHCLU0Fny1DUY4oqEG4YUndQp4ue/bL2Qm8VuWwZinpMUQXKit+PPqUtv7iCqmpD9uEifeZLBNLb/FTEyy0o5bw/1x0MvnVcdO1xxlfvGMJ1gzsHozTlJ3qbn1JN6Ngmnh1PX12nzfXEy4PT1zY675b9BRSX60macGLLUNRjiirQYqOjWPf4FY1+P7+k4pS24vJKrnn1myZD01VRWSVFepbb9mwZinpMUQVD25axfHzvhQ1+b/DkBcxY+QNzNhyobduR6zxTvXrPMbeWf9YT8xnyx4XeF6r8KjrYBShlJ8N6nMHTN5zNY//ddMr3MqxbAnu0u5jEFjHcMOVbAKI8GIqsvFLvtbY7DUWl6vnJ+T04XlxO93atGtw1/vE/V9TZnXZnxB29ODx02HL3Walgu39UGtcN7kz/Tm1O+V7944sNHW+sb/Z6nzy+SAWAhqJSTZh9/0W8ffdwt6Y1xrA08zDV1Yb84grW7DnKlMXOQW19PRpPVbUhv6SCA/klTPt2t0+W+dTsLbX1RjJb7j6LyDhgXGpqarBLUREu2hHFpX2SWfzwSC574atGp1udfZRfzljHvuMlpHZoXed2wfsuq/tzvPaHYwzpfoZXdf3xsy1MW5ZNn46t2ZF7givPOpPObVu4Ne+8TQcY2bcD8TGOOu1vWeFav95IY8stRT37rOymZ/tWZD97LQO7NPwzectr37HveAnAKfdP178r5sb/t4yxf/sGYwxHi8o9rmXW+v1MW5YNnDwDXlnl3k0Y6/YeZ9J73zN5tk+eMee2j9fk8LZVs93ZMhSVsqsPJ13AyL7JHs3zs7dWnNK2aV8BPR+Zy9A/LmRjTj7GGKYszuLYaUJyz5GiBk/+VLl5Z1rNdZLZh4vdmr6+5buONOtRDg99uJ4nZnkXxCkZc0jJmEOln0dL11BUygPxMQ6m/c9wtj41hp9d0INnbxp42nlOlyEz1+awcvdR/jJ/+2mvY/zhaMNhVu1mKNZcPtScYFu+6wjjpy4P+sO/iiuq/Lp8DUWlmqFFrIOnrj+b8cO7s/7xK71a1r++zaag9OQlOykZc6huJLTmbjzYYLu7Ibff2sVfmX2UjTkn7xj7NqvhgS/yCkt5ctZmKqqqyS0oBWDnIecuuzEmYM+4cX0//D1cg4aiUl5KbBnDv+46l57tW/H3Hw0hOSHO42VMeKfuwCcljWwNTV/5Q4Pt7obTQx+ur329cMvJgP1808k7dZ79fFvt6ydnbWbasmwWb8urDaO5G53TTp69hbRHPyclYw4vL9zhVv+//Wg9hwrL3JrWVUX1yfVr7A+Gr2goKuUDl/XrwOKHRzJ2UGdWPDKa3c9cw+s/Hdbs5X2wem/tMbTsw0Wnnf5Y0anXShpj+OsXmbUngOp79cus2ovKYxwno+C1r3fy2H83kpIxhy37CwBYvusoBmcYVVqhNM3lxMlfF2U2Wptr/R+szql9Ls6ynYdPG+bF5ZVUVlXzwaq9tW2VGopKhZaoKEFEuOqsM/n6NyNpE3/yyrcrB3R0axmTZ2+pfT3yha9IyZjDpHfXNDr9T95cwTNzt9ZpW5+Tz8tf7OC+978nv7iC9XuPnzLf3dNWAc5deFfvLXdukWYfcR7DfOvb3azKbvoe7/LKavJLKurcyrh4Wx4j613K9Ny8baz94Rg/emMFLyzY3ujyKqqqGfD4fFIf/Zw/fHryJM25f/qCIyc839p0ly3HU3S5TnFCZmbjf4GUChWLtuayMvsoGWP6sedI8SlB4SvrHr+CVnHOEB7x3GIOWscBm3Jpn2S+3nHIo362Pz2Gvo/Nq9PWsU0cuQXOsNr9zDUYA71+P7fB+UekteebzMNc3r8Dd1/ckwt6tUPq3UN+oqySs5+Y3+D8VwzoyBs/a3A4RLc0NZ6iLUOxhg4yq8JVfkkFgycvqNN2Xs8kVuw+GqSKPHO6IP3XXedyYWq7U4KzMc/fMojoKKG8sprb0rsRFSUNvkc1rhjQkbM7J/LD0WIevqoPnRLdu3C9hoaiUjZUWVXNiwt3cPPQrqS0a0m0I4rtBwu56pUlwS4t4JJaxda5kL1zYjyzHriY9Ke/cGv+7Gev9ag/DUWlQkhBaQXvLd/DJWnJVFUbOraJ5/xnTj4uYfPkqzirkd3KSOXLULTlvc9KRbI28TH8YmTd+4+zn72WzNxCRIRWcdG8e89wFm7JpVf7VsTHOMgrLGPqkl1uDVH24q2D61yao+rSUFQqRKR1TKh9PSItmRFpdW83fHB0GtsOFjDmlW/od2YCv7mqLw99uJ6WMQ7255884XLzsK58u/MwczYcoKzeoLe/HJ3W5OU1zfXzi3syun9H7nhjuc+X7Wu6+6xUhMgvrmB/fkmdMSJX7DpC+4Q4qqoNfVxCd9HWXACenrOVGRPPJ7FFDA9MX0tCXDQz1+5rcPlZf7qa1Ec/P6X9+VsGcVt6NwDW7DnKzf/4zuPauye15O6LUnjS5VIlV3pMUSkVVPM3H6RTYjxntoknPtZBrCOK+BgH+4+XsP1gIUuzDtMpMZ5R/TrQK7l1nXmNMRwrrmDoHxcy5qwzmbf5IC/cOphN+/K588IU8gpKqaw2GOO8jnLSyN78+oo+AGzIOc7x4gqyjxTxuHXt4vM3D+K2c7t5VL+GolLKdsoqq4iJiiIqyv1n3LiqrjbNnldPtCilbCcu2nH6iZrQ3EA87XL9slSllApRAQ1FEWklImtEZGwg+1VKKXe5FYoi8paI5InIpnrtY0Rku4hkiUiGG4v6HfBBcwpVSqlAcPeY4jTg78A7NQ0i4gCmAFcAOcAqEZkFOIBn6s1/NzAI2ALEe1eyUkr5j1uhaIxZIiIp9ZqHA1nGmF0AIjIDuN4Y8wxwyu6xiFwGtAIGACUiMtcYc8pgaiIyEZgI0L17dw9WRSmlvOfN2ecuwF6X/+cA5zU2sTHmUQARuQs43FAgWtNNBaZa0x4SkT0e1tUeaHhsdf8LZt+R3n8kr3uk99+cvns09g1vQrGh8+GnvejRGDPN3Q6MMZ49Ng0QkdWNXX/kb8HsO9L7j+R1j/T+fd23N2efcwDXy8i7Avu9K0cppYLLm1BcBaSJSE8RiQXGA7N8U5ZSSgWHu5fkTAe+A/qKSI6I3GOMqQTuB+YDW4EPjDHePe3aN6ZGaN+R3n8kr3uk9+/Tvm1977NSSgWa3uanlFIuwiYUm3F3jafL7yYii0Vkq4hsFpFfWu1Pisg+EVlnfV3jMs8jVj3bReQqH9SQLSIbrX5WW21JIrJQRDKtf8/wR/8i0tdlHdeJSIGI/Mqf69/QnVTNWV8RGWa9b1ki8qrUf2yc+33/RUS2icgGEflERNpa7SkiUuLyHrzmTd9N9O/xe+3j/v/j0ne2iKzzx/o38bsWkM8eY0zIf+G8i2Yn0AuIBdYDA3zcRydgqPU6AdiB80L0J4GHG5h+gFVHHNDTqs/hZQ3ZQPt6bc8DGdbrDOA5f/Vf7/0+iPNaL7+tP3AJMBTY5M36AiuBC3BeRvY5cHUz+74SiLZeP+fSd4rrdPWW43HfTfTv8Xvty/7rff9F4HF/rD+N/64F5LMPly3F2rtrjDHlwAzgel92YIw5YIz53npdiPPkUpcmZrkemGGMKTPG7AayrDp97Xrgbev128ANAeh/NLDTGNPUhfVe92+MWQLUf+anR+srIp2ANsaY74zzt+Qdl3k86tsYs8A4TzACLMd5GVqjmtt3Y/03wafrfrr+ra2t24DpTS3Di/e+sd+1gHz24RKKDd1d01RgeUWctzwOAVZYTfdbu1RvuWzS+6MmAywQ50hDE622jsaYA+D8YQI6+LH/GuOp+wsRqPUHz9e3i/Xa13XcjXPLo0ZPEVkrIl+LyAiXmnzdtyfvtb/WfQSQa4xxfZiLX9a/3u9aQD77cAnFZt1d06yORFoDHwO/MsYUAP8AegPnAAdw7lb4q6aLjDFDgauB+0TkkqZK9UP/iPOa1OuAD62mQK5/k6U10p/P6xCRR4FK4H2r6QDQ3RgzBPg18G8RaeOHvj19r/31GdxB3T+Kfln/Bn7XGp20kX6a1X+4hGJA7q4RkRicH9L7xpiZAMaYXGNMlXHey/0GJ3cRfV6TMWa/9W8e8InVV661m1Czu5Lnr/4tVwPfG2NyrVoCtv4WT9c3h7q7uV7VISJ34hzw5MfWLhnWbtsR6/UanMe0+vi672a81z7tH0BEooGbgP+41OXz9W/od40AffbhEop+v7vGOo7yJrDVGPOSS3snl8luBGrO1s0CxotInIj0BNJwHvRtbv+tRCSh5jXOg/6brH7utCa7E/jUH/27qLOVEKj1d+HR+lq7WYUicr71Gf7MZR6PiMgYnGOCXmeMKXZpTxbnUHqISC+r712+7Ntatkfvta/7t1wObDPG1O6W+nr9G/tdI1CfvTtnokLhC7gG51mqncCjflj+xTg3vTcA66yva4B3gY1W+yygk8s8j1r1bMfNs35N9N8L5xm29cDmmnUE2gGLgEzr3yR/9G8tryVwBEh0afPb+uMM3wNABc6/+vc0Z32BdJwBshPnuKDSzL6zcB67qvn8X7Omvdn6TNYD3wPjvOm7if49fq992b/VPg2YVG9an64/jf+uBeSz1ztalFLKRbjsPiullE9oKCqllAsNRaWUcqGhqJRSLjQUlVLKhYaiihgiMlJEPgt2HcreNBSVUsqFhqKyHRH5iYistMbme11EHCJyQkReFJHvRWSRiCRb054jIsvl5BiHZ1jtqSLyhYist+bpbS2+tYh8JM5xEd93a3w9FVE0FJWtiEh/4Hacg1+cA1QBPwZa4bzneijwNfCENcs7wO+MMYNw3u1R0/4+MMUYMxi4EOfdGeAcceVXOMfg6wVc5OdVUiHGm+c+K+UPo4FhwCprI64Fzhv/qzk5CMF7wEwRSQTaGmO+ttrfBj607hHvYoz5BMAYUwpgLW+lse7bFefI0SnAUr+vlQoZGorKbgR42xjzSJ1GkT/Um66p+1Ob2iUuc3ldhf4OqHp091nZzSLgFhHpALXP5eiB82f1FmuaHwFLjTH5wDGXQU1/CnxtnGPv5YjIDdYy4kSkZSBXQoUu/SupbMUYs0VEHsM5wngUzlFa7gOKgLNEZA2Qj/O4IziHkHrNCr1dwP9Y7T8FXheRp6xl3BrA1VAhTEfJUSFBRE4YY1oHuw4V/nT3WSmlXOiWolJKudAtRaWUcqGhqJRSLjQUlVLKhYaiUkq50FBUSikXGopKKeXi/wPjfxeb5lqutAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUIAAADCCAYAAAA8asvwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAWHElEQVR4nO3df3BUd/3v8ec7kRjFtlhIv19nCSRSbgMtBJOlwpe0irW/BEFsOwWvc1G/QECZzvWOWNreezVW51tvGUc737bArr1tHVvm4ohS229Bb+swlfZbwkgrBGgJSU1Sv/yqbYkYAuR9/9glNyQLOZvsZjec12NmJ3vOfj7nvPdzlhfnnN09a+6OiEiYFeS6ABGRXFMQikjoKQhFJPQUhCISegpCEQk9BaGIhN4Hcl1AKmPGjPGysrJclyEiF5mdO3cedfeS3vPzMgjLysqor6/PdRkicpExs7dSzdehsYiEnoJQREJPQSgioacgvMh1trTQOGcue6++hsY5c+lsacl1SSJ5J1AQmtktZrbfzA6Y2eoUj/9nM3s9edtuZpVB+14s6urq+NmaNX1CJx6PU1dXl7OaXl+4iM6mJjhzhs6mJlqWr8hpTfmkrq6OeDzeZ77GJ39la5v1G4RmVgg8DNwKTAYWmdnkXs2agE+5+1TgfmB9Gn0vCpFIhFEPP0LHwYPdofP6wkXU1tYSiURyVtMHjx6Frq7EjK4uTh48mNOa8kkkEuH+b6ykflZN939eP1uzRuOTx7K1zay/y3CZ2Uzgu+5+c3L6HgB3/5fztP8osNvdI+n2PSsajfpw/PhMw6TJWI/xPO3Ov6/6FkuWLMlZTfWzaig+epRCM7rMONjRwaH77s1pTflE4zP8DGabmdlOd4/2nh/k0DgC9Dyx1Jqcdz7/DPxbun3NbJmZ1ZtZ/ZEjRwKUlX8++PGP02UGwBl3To4Zk/N/UFM3PE3HmDGcdudgRwfvfuPrOa8pn4x8910Kk9uswJ0JxcUanzyXjW0WJAgtxbyUu5FmNptEEN6dbl93X+/uUXePlpT0+eD3sFC69lH+fvnlnHanqbOTvQu+kOuSKCotZddXv8LUN/Yzr7mJk6NG5bqkvFJUVnbOf14nLr88xxVJf7KxzYIEYStQ2mN6LPB270ZmNhWIA/Pd/Vg6fS8WT27ZwidfeZl/X/UtDt13L1+5++6UJ3aHUjwep7a2llgsRiwWo7a2Nuc15ZPtsz/NwY4O3IyOMWO4rb5e45PnsrLN3P2CNxJfwzsIlANFwGvA1b3ajAMOAP+Ubt9Ut+rqah9uYrGYFxQUeCwWu+C8sNeUTzQ+w89gtxlQ7ykyp9/vGrv7aTNbCWwBCoHH3H2PmS1PPr4W+J/AaOARS+yynvbEYW7KvoOL7vzU1tbGunXrzjlXcfZ+W1ubaspDGp/hJ1vbrN93jXNhuL5rLCL5bTDvGouIXNQUhCISegpCEQk9BaGIhJ6CUERCT0EoIqGnIBSR0FMQikjoKQhFJPQUhCISegpCEQk9BaGIhJ6CUERCT0EoIqGnIBSR0FMQikjoKQhFJPQUhCISegpCEQk9BaGIhJ6CUERCT0EoIqGnIBSR0FMQikjoKQhFJPQUhCISegpCEQm9QEFoZreY2X4zO2Bmq1M8XmFmL5vZSTP7Vq/Hms3sT2a2y8zqM1W4iEimfKC/BmZWCDwM3Ai0AjvMbLO7N/Ro9g5wF/CF8yxmtrsfHWStIiJZEWSP8FrggLsfdPdOYAMwv2cDdz/s7juAU1moUUQkq4IEYQRo6THdmpwXlANbzWynmS1LpzgRkaHQ76ExYCnmeRrrmOXub5vZFcBvzWyfu2/rs5JESC4DGDduXBqLFxEZnCB7hK1AaY/pscDbQVfg7m8n/x4GNpE41E7Vbr27R909WlJSEnTxIiKDFiQIdwATzazczIqAhcDmIAs3s5FmdsnZ+8BNwO6BFisikg39Hhq7+2kzWwlsAQqBx9x9j5ktTz6+1sz+EagHLgW6zOy/ApOBMcAmMzu7rqfc/fmsPBMRkQEKco4Qd38OeK7XvLU97v8HiUPm3t4HKgdToIhItumbJSISegpCEQk9BaGIhJ6CUERCT0EoIqGnIBSR0FMQikjoKQhFJPQUhCISegpCEQk9BaGIhJ6CUERCT0EoIqGnIBSR0At0GS4Ryb5Tp07R2tpKR0dHrksZ9oqLixk7diwjRowI1F5BKJInWltbueSSSygrKyN5MWMZAHfn2LFjtLa2Ul5eHqiPDo1F8kRHRwejR49WCA6SmTF69Oi09qwVhCJ5RCGYGemOo4JQRAB49913eeSRR7K+nl/96lc0NDRkfT3pUBCKDEN1dXXE4/E+8+PxOHV1dQNaZrpB6O50dXWlvR4FoYhkRCQSoba29pwwjMfj1NbWEolEBrTM1atX09jYyLRp0/jmN7/JDTfcQFVVFVOmTOHXv/41AM3NzUyaNImvf/3rVFVV0dLSwv33309FRQU33ngjixYtYs2aNQA0NjZyyy23UF1dzXXXXce+ffvYvn07mzdvZtWqVUybNo3GxsbBD0YmuHve3aqrq10kbBoaGtJqH4vFvKCgwGOx2Dn3B6qpqcmvvvpqd3c/deqUv/fee+7ufuTIEZ8wYYJ3dXV5U1OTm5m//PLL7u6+Y8cOr6ys9BMnTvj777/vV155pT/44IPu7v6Zz3zG33jjDXd3f+WVV3z27Nnu7r548WLfuHHjgOsMKtV4AvWeInP08RmRYWrJkiUALF26FIBYLNY9b7DcnXvvvZdt27ZRUFBAW1sbhw4dAmD8+PHMmDEDgJdeeon58+fzoQ99CIDPf/7zALS3t7N9+3buuOOO7mWePHkyI7Vlg4JQRPr4+c9/zpEjR9i5cycjRoygrKys++MoI0eO7G6X2Mnqq6uri1GjRrFr166hKHfQdI5QZJg6e04wFosRi8X6nDNM1yWXXMLx48cBeO+997jiiisYMWIEL774Im+99VbKPjU1NTzzzDN0dHTQ3t7Os88+C8Cll15KeXk5GzduBBKB+dprr/VZT77QHqHIMHQ2BNetW3fO4XBtbS3AgA6RR48ezaxZs7jmmmuYPn06+/btIxqNMm3aNCoqKlL2mT59OvPmzaOyspLx48cTjUa57LLLgMRe5YoVK/j+97/PqVOnWLhwIZWVlSxcuJClS5fy0EMP8Ytf/IIJEyYMYAQyy863a5tL0WjU6+vrc12GyJDau3cvkyZNCtS2rq6OSCTSJ/Di8ThtbW185zvfyUaJKbW3t/ORj3yEEydOcP3117N+/XqqqqqGbP3nk2o8zWynu0d7t9UeocgwdL6gy9SbJelYtmwZDQ0NdHR0sHjx4rwIwXQFCkIzuwX4CVAIxN39gV6PVwD/G6gC7nP3NUH7isjw9tRTT+W6hEHr980SMysEHgZuBSYDi8xscq9m7wB3AWsG0FdEJKeCvGt8LXDA3Q+6eyewAZjfs4G7H3b3HcCpdPuKiORakCCMAC09pluT84II3NfMlplZvZnVHzlyJODiRUQGL0gQprqeTdC3mgP3dff17h5192hJSUnAxYuIDF6QIGwFSntMjwXeDrj8wfQVkWHs97//PXPnzgVg8+bNPPDA+d8nHeglwL773e92X+RhMIIE4Q5gopmVm1kRsBDYHHD5g+krInnozJkzafeZN28eq1evPu/jQ3UtxPPpNwjd/TSwEtgC7AX+j7vvMbPlZrYcwMz+0cxagf8G/HczazWzS8/XN1tPRiRMOltaaJwzl71XX0PjnLl0trT036kfzc3NVFRUsHjxYqZOncrtt9/OiRMnKCsr43vf+x41NTVs3LiRrVu3MnPmTKqqqrjjjjtob28H4Pnnn6eiooKamhp++ctfdi/38ccfZ+XKlQAcOnSIBQsWUFlZSWVlJdu3bz/nEmCrVq0C4MEHH2T69OlMnTr1nM9N/uAHP+Cqq67is5/9LPv37x/0c4aAnyN09+eA53rNW9vj/n+QOOwN1FdEBq9l+Qo6m5qgq4vOpiZalq9gwrO/GfRy9+/fz09/+lNmzZrF1772te49teLiYl566SWOHj3KF7/4RX73u98xcuRIfvjDH/KjH/2Ib3/72yxdupQXXniBK6+8kjvvvDPl8u+66y4+9alPsWnTJs6cOUN7ezsPPPAAu3fv7r5Iw9atW3nzzTd59dVXcXfmzZvHtm3bGDlyJBs2bOCPf/wjp0+fpqqqiurq6kE/Z32zRGSY6mxuhrNXiO7qSkxnQGlpKbNmzQLgy1/+Mg899BBAd7C98sorNDQ0dLfp7Oxk5syZ7Nu3j/LyciZOnNjdd/369X2W/8ILL/Dkk08CUFhYyGWXXcZf//rXc9ps3bqVrVu38olPfAJIfI3vzTff5Pjx4yxYsIAPf/jDQOKQOxMUhCLDVFFZWfceIQUFFJWVZWS5vX/46Oz02ctvuTs33ngjTz/99Dntdu3albEfn3J37rnnnu6LSJz14x//OCs/cKXLcIkMU6VrH6WovBwKCykqL6d07aMZWe6f//xnXn75ZQCefvppampqznl8xowZ/OEPf+DAgQMAnDhxgjfeeIOKigqampq6L7/fOyjPuuGGG3j00UStZ86c4f333+9zaa6bb76Zxx57rPvcY1tbG4cPH+b6669n06ZN/P3vf+f48eM888wzGXnOCkKRYaqotJQJz/6GSXt2M+HZ31BUWtp/pwAmTZrEE088wdSpU3nnnXdYsWLFOY+XlJTw+OOPs2jRIqZOncqMGTPYt28fxcXFrF+/njlz5lBTU8P48eNTLv8nP/kJL774IlOmTKG6upo9e/accwmwVatWcdNNN/GlL32JmTNnMmXKFG6//XaOHz9OVVUVd955J9OmTeO2227juuuuy8hz1mW4RPJEOpfhypbm5mbmzp3L7t27c1pHJqRzGS7tEYpI6CkIRaRbWVnZRbE3mC4FoYiEnoJQJI/k4zn74SjdcVQQiuSJ4uJijh07pjAcJHfn2LFjFBcXB+6jD1SL5ImxY8fS2tqKrsc5eMXFxYwdm/JbvykpCEXyxIgRIygvL891GaGkQ2MRCT0FoYiEnoJQREJPQSgioacgFJHQUxCKSOgpCEUk9BSEIhJ6CkIRCT0FoYiEnoJQREJPQSgioacgFJHQUxCKSOgpCEUk9AIFoZndYmb7zeyAma1O8biZ2UPJx183s6oejzWb2Z/MbJeZ6Tc6RSTv9HthVjMrBB4GbgRagR1mttndG3o0uxWYmLx9Eng0+fes2e5+NGNVi4hkUJA9wmuBA+5+0N07gQ3A/F5t5gNPesIrwCgz+1iGaxURyYogQRgBWnpMtybnBW3jwFYz22lmywZaqIhItgT5zRJLMa/3z2xdqM0sd3/bzK4Afmtm+9x9W5+VJEJyGcC4ceMClCUikhlB9ghbgdIe02OBt4O2cfezfw8Dm0gcavfh7uvdPeru0ZKSkmDVi4hkQJAg3AFMNLNyMysCFgKbe7XZDPyX5LvHM4D33P0vZjbSzC4BMLORwE3A7gzWLyIyaP0eGrv7aTNbCWwBCoHH3H2PmS1PPr4WeA74HHAAOAF8Ndn9H4BNZnZ2XU+5+/MZfxYiIoNg7r1P9+VeNBr1+np95FBEMsvMdrp7tPd8fbNEREJPQSgioacgFJHQUxCKSOgpCEUk9BSEIhJ6CkIRCT0FoYiEnoJQREJPQSgioacgFJHQUxCKSOgpCEUk9BSEIhJ6CkIRCT0FoYiEnoJQREJPQSgioacgFJHQUxCKSOgpCEUk9BSEIhJ6CkIRCT0FoYiEnoJQhlRnSwuNc+ay9+praJwzl86WllyXJP0IwzZTEMqQalm+gs6mJjhzhs6mJlqWr8h1SXklH0MnDNtMQShDqrO5Gbq6EhNdXYlp6ZaPoROGbaYglCFVVFYGBcmXXUFBYlq65WPohGGbBQpCM7vFzPab2QEzW53icTOzh5KPv25mVUH7DlRdXR3xeLzP/Hg8Tl1dXaZWIxlUV1fH9tmfpqi8HAoLKSovp3Tto9pmSXV1dfxt1Kg+oZPL8QnLNus3CM2sEHgYuBWYDCwys8m9mt0KTEzelgGPptF3QCKRCPd/YyX1s2q6z6f8bM0aamtriUQimViFZFgkEuErd9/Niwu+wKQ9u5nw7G94cssWbbOkSCTCbfX1/O2jH+0One2zP53T8QnNNnP3C96AmcCWHtP3APf0arMOWNRjej/wsSB9U92qq6s9iB3/NMv/9J+u8oarKnx3xSTfXFbusVgsUF/JjVgs5gUFBR6Lxc65Lwn5OD75WNNAAfWeInM+ECArI0DPt65agU8GaBMJ2BcAM1tGYm+ScePGBSgLRr77LpgBUODOhOJiPr9kSaC+khtLkttn6dKlAMRise55kp/jk481ZVqQc4SWYp4HbBOkb2Km+3p3j7p7tKSkJEBZiZO4XckgPOPOicsvD9RPRKSnIHuErUBpj+mxwNsB2xQF6Dtg22d/mlENDUwoLqZj9Ghuq6/nf8TjF93/VheTeDxObW0tsVgMgNraWgBts6R8HJ98rCnjUh0v+7nn/z4AHATKSQTba8DVvdrMAf6NxB7gDODVoH1T3YKcI0x1rmK4n7+42GmbXVg+jk8+1jQYDPQcobufNrOVwBagEHjM3feY2fLk42uB54DPAQeAE8BXL9R3sOEN0NbWxrp16875X+ns/ba2tkysQjJM2+zC8nF88rGmbLBESOaXaDTq9fX1uS5DRC4yZrbT3aO95+ubJSISegpCEQm9vDw0NrMjwFtpdBkDHM1SOQORb/VA/tWkei4s3+qB/KtpIPWMd/c+n8/LyyBMl5nVpzruz5V8qwfyrybVc2H5Vg/kX02ZrEeHxiISegpCEQm9iyUI1+e6gF7yrR7Iv5pUz4XlWz2QfzVlrJ6L4hyhiMhgXCx7hCIiAzYsgjDVVa7N7HIz+62ZvZn8+9Ee7e9Jtt1vZjfnuiYzG21mL5pZu5n96xDWc4eZ7TGzLjOL9mqf1TFKp54hGp/HzOywme3uMS9nr6F06snh+OTy9RO4noyMT6ovIOfTjcR3lBuBj/P/L9wwGfhfwOpkm9XAD5P3JyfbfJDExR4agcIc1zQSqAGWA/86hGM0CbgK+D0Q7dE+q2M0gHqyOj7JdVwPVAG7e8zL5WsonXpyNT45ef0MoJ5Bj89w2CO8Fjjg7gfdvRPYAMxP3p5ItnkC+ELy/nxgg7ufdPcmEheCuDaXNbn739z9JaAjw3VcsB533+vu+1O0z/YYpVXPEIwP7r4NeKfX7Jy9htKpJ1fjk8PXT1r1ZGJ8hkMQnu/q1//g7n8BSP69op/2uawp29J9ztkeo6HYBpmQy9dQOvXkm+GyfQMbDkEY+CrXA2w/EEOxjnTk2xjl2/ika7jXn20X3fgMhyA839WvD5nZxwCSfw/30z6XNWVbus8522M0FNsgE3L5GkqnnnwzXLZvYMMhCHcAE82s3MyKgIXA5uRtcbLNYuDXyfubgYVm9kEzKyfxE6Ov5rimbDtfPeeT7TFKt55cyeVrKJ168k2uxid7svEOVBbe0foc8AaJd6fuS84bDfxf4M3k38t7tL8v2XY/cGue1NRM4uRvO4n/UScPQT0Lkus6CRzi3J9WzeoYDaCebI/P08BfgFPJ5f9zLl9DA6gnF+OTy9dPuvUManz0zRIRCb3hcGgsIpJVCkIRCT0FoYiEnoJQREJPQSgioacgFJHQUxCKSOgpCEUk9P4f3W8HQRK7GiMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(5,3))\n",
    "\n",
    "plt.plot(loss_list)\n",
    "plt.yscale('log')\n",
    "plt.xlabel('epoch')\n",
    "\n",
    "test = network_odd(trained_params['theta_even'], trained_params['theta_odd'], \n",
    "                   trained_params['phi_even'], trained_params['phi_odd'], X_padded)\n",
    "\n",
    "plt.figure(figsize=(5,3))\n",
    "plt.plot(jnp.arange(X.shape[1])+1, jnp.abs(Y[0,:]), \"x\", mec=\"k\", mew=1.25, ms=7, label='target')\n",
    "plt.plot(jnp.arange(X.shape[1])+1, jnp.abs(test[0,:]), \"o\", c=\"#d62728\", mew=1, ms=4, label='predicted')\n",
    "xticks = ['000','001','010','011','100','101','110','111']\n",
    "plt.xticks(jnp.arange(int(2**num_bits))+1,xticks)\n",
    "plt.legend()\n",
    "\n",
    "# wl = 1.515e-6\n",
    "# r0 = 3.225e-6\n",
    "# k0 = anp.pi*2/wl\n",
    "# phi_n_factor = 2*anp.pi*r0*k0\n",
    "# n_list = phi_list/2/anp.pi/r0/k0\n",
    "# print('n value', n_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4e3158ac",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a26715b",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}