Fourier_Ptychography_for_Phase_Retrieval.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ptychography using gradient descent\n",
    "### Chengyu Wang, Duke University\n",
    "### David J. Brady, University of Arizona\n",
    "### Timothy J. Schulz, Michigan Technological University\n",
    "\n",
    "This script implements the 2D Fourier ptychography and a phase retrieval algorithm using gradient descent as described in \"Photon-limited bounds for phase retrieval.\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import tensorflow as tf\n",
    "import numpy as np\n",
    "from numpy import matlib\n",
    "import random\n",
    "tf.compat.v1.disable_eager_execution()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Error metrics\n",
    "\n",
    "If $\\underline{F}$ is an $N$-element optical field, and an estimate of this field is $\\widehat{\\underline{F}}$.\n",
    "One way to address this is to optimize the MSE over a global phase:\n",
    "\\begin{equation} \n",
    "MSE = \\underset{\\psi}{argmin} \\frac{1}{N} \\left\\| e^{j \\psi} \\underline{F} - \\widehat{\\underline{F}} \\right\\|^2. \n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def minAngMSE(xt,xest):\n",
    "    mse=(np.abs(np.vdot(xt,xt))+np.abs(np.vdot(xest,xest))-2*np.abs(np.vdot(xt,xest)));\n",
    "    return mse/xt.shape[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Forward model\n",
    "For the forward model, $N_x$ is the number of pixels in each dimension of the random gaussian signal, $L_x$ is the size of the subaperture/selected window, and $M_x$ is the size of the each sampld frame. The subaperture moves circularly with step size $\\Delta$. The total number of frames is $(\\frac{N_x}{\\Delta})^2$. We consider $N_x$ and $\\Delta$ to be powers of 2.\n",
    "\n",
    "A phase shift can be added to the window using the argument *is_mask*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Dimensions\n",
    "N = 256\n",
    "L = 64\n",
    "pad_size = 32\n",
    "M = L + pad_size*2\n",
    "pitch = 16 # Delta\n",
    "total_frame = int(N/pitch)**2\n",
    "\n",
    "is_mask = False ## add a phase shift\n",
    "is_noisy = True  ## add Poisson noise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Phase mask\n",
    "Define the phase mask. 3 different examples are shown here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "## random phase mask\n",
    "mask = np.exp(1j * np.random.random([L,L]) * 2 * np.pi)\n",
    "\n",
    "# ## periodic phase mask\n",
    "# K = np.pi/L\n",
    "# mesh = np.arange(L)\n",
    "# xlocations,ylocations = np.meshgrid(mesh,mesh);\n",
    "# mask = np.exp(1j * np.cos(K*xlocations) * np.cos(K*ylocations)  * 2*np.pi)\n",
    "\n",
    "# ## MURA mask\n",
    "# def qres(p):\n",
    "#     qr = np.zeros([p-1])\n",
    "#     for i in range(p-1):\n",
    "#         qr[np.mod((i+1)**2,p)-1]=1\n",
    "#     return qr\n",
    "\n",
    "# def gen_mura(p):\n",
    "#     ## p is a prime number\n",
    "#     mask = np.zeros([p,p])\n",
    "#     mask[1:,:] = 1\n",
    "#     qr = qres(p)\n",
    "#     for i in range(1,p):\n",
    "#         for j in range(1,p):\n",
    "#             if qr[i-1] == 1 and qr[j-1] == 1:\n",
    "#                 mask[i,j] = 1\n",
    "#             elif qr[i-1] == 0 and qr[j-1] == 0:\n",
    "#                 mask[i,j] = 1\n",
    "#             else:\n",
    "#                 mask[i,j] = 0\n",
    "#     return mask\n",
    "# mask = gen_mura(L)*2-1 + 0 * 1j"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "## selected windows\n",
    "## xlocations/ylocations save the cordinates of the top left pixels of subaperture/selected windows\n",
    "mesh = np.arange(0,N,pitch)\n",
    "xlocations,ylocations = np.meshgrid(mesh,mesh);\n",
    "xlocations = xlocations.astype(int).flatten()\n",
    "ylocations = ylocations.astype(int).flatten()\n",
    "\n",
    "## Energy constraint\n",
    "## Count how many times each pixel (in Fourier domain) is measured\n",
    "## The energy will be eqaully divided to all measurements\n",
    "counts = np.zeros([2*N,2*N])\n",
    "for i in range(total_frame):\n",
    "     counts[xlocations[i]:xlocations[i]+L,ylocations[i]:ylocations[i]+L] += 1\n",
    "counts = counts[0:N,0:N] + counts[0:N,N:2*N] + counts[N:2*N,0:N] + counts[N:2*N,N:2*N]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define the tensorflow graph to compute gradient"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "## define the tensorflow model\n",
    "X_FT_EST = tf.compat.v1.placeholder(tf.complex128, shape=(None,L,L))\n",
    "paddings = tf.constant([[0,0],[pad_size, pad_size], [pad_size, pad_size]])\n",
    "MEASUREMENT = tf.compat.v1.placeholder(tf.float64, shape=(None,M,M))\n",
    "if is_mask:\n",
    "    MASK = tf.constant(mask)\n",
    "    X_FT_EST_PAD = tf.pad(tf.multiply(X_FT_EST,MASK),paddings)\n",
    "else:\n",
    "    X_FT_EST_PAD = tf.pad(X_FT_EST,paddings)\n",
    "FORWARD = tf.abs(tf.signal.ifft2d(tf.signal.ifftshift(X_FT_EST_PAD,[1,2]))) * M\n",
    "loss = tf.reshape(tf.cast(tf.math.reduce_sum(tf.math.square(abs(MEASUREMENT)-abs(FORWARD))),tf.complex128),[1,1])\n",
    "# loss = tf.reshape(tf.cast(tf.compat.v1.losses.mean_squared_error(abs(MEASUREMENT),abs(FORWARD)),tf.complex128),[1,1])\n",
    "weight = [[1 + 1j*0]]\n",
    "gradient = tf.gradients(loss, X_FT_EST, grad_ys = weight)\n",
    "\n",
    "## session\n",
    "sess = tf.compat.v1.Session()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulation\n",
    "\n",
    "In phase retirevel, we initialize the gradient descent with a simple projection method. The projection method itself is sensitive to noise, resulting in poor MSE, but it approximates the groundtruth faster than gradient descent method. The phase retrieval implemented in this script has two stages:\n",
    "\n",
    "- First a projection method approximates the grountruth\n",
    "- Second the gradient descent algorithm improves the accuracy.\n",
    "\n",
    "FYI: Because of the circular connections, the if ... elif ... elif ... else part is used to consider different locations of the winsow."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_simulation(XT,num_projection_loop = 50, num_gd_iteration = 50,lr = 0.5,return_error = False,verbose = True):\n",
    "    \n",
    "    ## XT: signal\n",
    "    ## num_projection_loop: number of loops for GS projection\n",
    "    ## num_gd_iteration: number of iterations for gradient descent\n",
    "    ## lr: learning rate for the gradient descent\n",
    "    ## return_error: return the reconstruction error\n",
    "    ## verbose: print the system configurations and the reconstruction error\n",
    "    \n",
    "    num_trial = XT.shape[0]\n",
    "    MSE = np.zeros([num_trial])\n",
    "    for ct in range(num_trial):\n",
    "        \n",
    "        ## signal\n",
    "        xt = XT[ct,:,:]\n",
    "        objectFT = np.fft.fftshift(np.fft.fft2(xt))/N/np.sqrt(counts)\n",
    "        objectFT = matlib.repmat(objectFT, 2, 2)\n",
    "\n",
    "        ## measurement\n",
    "        imSeqLowRes = np.zeros([total_frame, M, M]);\n",
    "        for i in range(total_frame):\n",
    "            imSeqLowFT = objectFT[xlocations[i]:xlocations[i] + L,ylocations[i]:ylocations[i] + L]\n",
    "            if is_mask:\n",
    "                imSeqLowFT = np.pad(imSeqLowFT * mask,pad_size)\n",
    "            else:\n",
    "                imSeqLowFT = np.pad(imSeqLowFT,pad_size)\n",
    "            imSeqLowRes[i,:,:] = np.abs(np.fft.ifft2(np.fft.ifftshift(imSeqLowFT))) * M\n",
    "        if is_noisy:\n",
    "            imSeqLowRes = np.sqrt(np.random.poisson(np.power(np.abs(imSeqLowRes),2)))\n",
    "\n",
    "        ## Projection\n",
    "        imageRecover = np.random.randn(N,N);\n",
    "        imageRecoverFT = np.fft.fftshift(np.fft.fft2(imageRecover))\n",
    "        seq = list(range(total_frame))\n",
    "        for loop in range(num_projection_loop):\n",
    "            random.shuffle(seq)\n",
    "            for i in seq:\n",
    "                if not xlocations[i] > N - L and not ylocations[i] > N-L:\n",
    "                    if is_mask:\n",
    "                        imLowRes = np.fft.ifft2(np.fft.ifftshift(np.pad(imageRecoverFT[xlocations[i]:xlocations[i]+L,ylocations[i]:ylocations[i]+L] * mask,pad_size)))\n",
    "                    else:\n",
    "                        imLowRes = np.fft.ifft2(np.fft.ifftshift(np.pad(imageRecoverFT[xlocations[i]:xlocations[i]+L,ylocations[i]:ylocations[i]+L],pad_size)))\n",
    "\n",
    "                    imLowRes = imSeqLowRes[i,:,:] * np.exp(1j * np.angle(imLowRes))\n",
    "                    inverse = np.fft.fftshift(np.fft.fft2(imLowRes)) / M\n",
    "                    if is_mask:\n",
    "                        imageRecoverFT[xlocations[i]:xlocations[i]+L,ylocations[i]:ylocations[i]+L] = inverse[pad_size:pad_size+L,pad_size:pad_size+L]/mask\n",
    "                    else:\n",
    "                        imageRecoverFT[xlocations[i]:xlocations[i]+L,ylocations[i]:ylocations[i]+L] = inverse[pad_size:pad_size+L,pad_size:pad_size+L]\n",
    "                else:\n",
    "                    imageRecoverFTRep = matlib.repmat(imageRecoverFT, 2, 2)\n",
    "                    if is_mask:\n",
    "                        imLowRes = np.fft.ifft2(np.fft.ifftshift(np.pad(imageRecoverFTRep[xlocations[i]:xlocations[i]+L,ylocations[i]:ylocations[i]+L] * mask,pad_size)))\n",
    "                    else:\n",
    "                        imLowRes = np.fft.ifft2(np.fft.ifftshift(np.pad(imageRecoverFTRep[xlocations[i]:xlocations[i]+L,ylocations[i]:ylocations[i]+L],pad_size)))\n",
    "                    imLowRes = imSeqLowRes[i,:,:] * np.exp(1j * np.angle(imLowRes))\n",
    "                    inverse = np.fft.fftshift(np.fft.fft2(imLowRes)) / M\n",
    "                    if is_mask:\n",
    "                        inverse = inverse[pad_size:pad_size+L,pad_size:pad_size+L]/mask\n",
    "                    else:\n",
    "                        inverse = inverse[pad_size:pad_size+L,pad_size:pad_size+L]\n",
    "                    if xlocations[i] > N - L and not ylocations[i] > N-L:\n",
    "                        imageRecoverFT[xlocations[i]:N,ylocations[i]:ylocations[i]+L] = inverse[0:N-xlocations[i],:]\n",
    "                        imageRecoverFT[0:L+xlocations[i]-N,ylocations[i]:ylocations[i]+L] = inverse[N-xlocations[i]:,:]\n",
    "                    elif not xlocations[i] > N - L and ylocations[i] > N-L:\n",
    "                        imageRecoverFT[xlocations[i]:xlocations[i]+L,ylocations[i]:N] = inverse[:,0:N-ylocations[i]]\n",
    "                        imageRecoverFT[xlocations[i]:xlocations[i]+L,0:L+ylocations[i]-N] = inverse[:,N-ylocations[i]:]\n",
    "                    else:\n",
    "                        imageRecoverFT[xlocations[i]:N,ylocations[i]:N] = inverse[0:N-xlocations[i],0:N-ylocations[i]]\n",
    "                        imageRecoverFT[0:L+xlocations[i]-N,0:L+ylocations[i]-N] = inverse[N-xlocations[i]:N,N-ylocations[i]:]\n",
    "                        imageRecoverFT[xlocations[i]:N,0:L+ylocations[i]-N] = inverse[0:N-xlocations[i],N-ylocations[i]:]\n",
    "                        imageRecoverFT[0:L+xlocations[i]-N,ylocations[i]:N] = inverse[N-xlocations[i]:,0:N-ylocations[i]]\n",
    "#             if verbose:\n",
    "#                 imageRecover=np.fft.ifft2(np.fft.ifftshift(imageRecoverFT*np.sqrt(counts))) * N\n",
    "#                 print(loop,minAngMSE(xt.flatten(),imageRecover.flatten()),end='\\r',flush=True)\n",
    "\n",
    "        ## Gradient descent\n",
    "        for iter_ in range(num_gd_iteration):\n",
    "            imageRecoverFTRep = matlib.repmat(imageRecoverFT, 2, 2)\n",
    "            batch_tensor = np.zeros([total_frame ,L, L]) * 1j\n",
    "            for i in range(total_frame):\n",
    "                batch_tensor[i,:,:] = imageRecoverFTRep[xlocations[i]:xlocations[i]+L,ylocations[i]:ylocations[i]+L]\n",
    "            gradient_sum = np.zeros([N,N]) * 1j\n",
    "            gradient_tensor = sess.run(gradient,feed_dict={X_FT_EST:batch_tensor,MEASUREMENT:imSeqLowRes})[0]\n",
    "            for i in range(total_frame):\n",
    "                if not xlocations[i] > N - L and not ylocations[i] > N-L:\n",
    "                    gradient_sum[xlocations[i]:xlocations[i]+L,ylocations[i]:ylocations[i]+L] += gradient_tensor[i,:,:]\n",
    "                elif xlocations[i] > N - L and not ylocations[i] > N-L:\n",
    "                    gradient_sum[xlocations[i]:N,ylocations[i]:ylocations[i]+L] += gradient_tensor[i,0:N-xlocations[i],:]\n",
    "                    gradient_sum[0:L+xlocations[i]-N,ylocations[i]:ylocations[i]+L] += gradient_tensor[i,N-xlocations[i]:,:]\n",
    "                elif not xlocations[i] > N - L and ylocations[i] > N-L:\n",
    "                    gradient_sum[xlocations[i]:xlocations[i]+L,ylocations[i]:N] += gradient_tensor[i,:,0:N-ylocations[i]]\n",
    "                    gradient_sum[xlocations[i]:xlocations[i]+L,0:L+ylocations[i]-N] += gradient_tensor[i,:,N-ylocations[i]:]\n",
    "                else:\n",
    "                    gradient_sum[xlocations[i]:N,ylocations[i]:N] += gradient_tensor[i,0:N-xlocations[i],0:N-ylocations[i]]\n",
    "                    gradient_sum[0:L+xlocations[i]-N,0:L+ylocations[i]-N] += gradient_tensor[i,N-xlocations[i]:N,N-ylocations[i]:]\n",
    "                    gradient_sum[xlocations[i]:N,0:L+ylocations[i]-N] += gradient_tensor[i,0:N-xlocations[i],N-ylocations[i]:]\n",
    "                    gradient_sum[0:L+xlocations[i]-N,ylocations[i]:N] += gradient_tensor[i,N-xlocations[i]:,0:N-ylocations[i]]      \n",
    "            imageRecoverFT = imageRecoverFT - lr * gradient_sum/counts\n",
    "#             if verbose:\n",
    "#                 imageRecover=np.fft.ifft2(np.fft.ifftshift(imageRecoverFT*np.sqrt(counts))) * N\n",
    "#                 print(iter_,minAngMSE(xt.flatten(),imageRecover.flatten()),end='\\r',flush=True)\n",
    "        imageRecover=np.fft.ifft2(np.fft.ifftshift(imageRecoverFT*np.sqrt(counts))) * N\n",
    "        MSE[ct] = minAngMSE(xt.flatten(),imageRecover.flatten())\n",
    "    if verbose:\n",
    "        print('photons = %d, pitch = %d,L = %d, M = %d, MSE = %f.'%(photons,pitch,L,M,np.mean(MSE)))\n",
    "    if return_error:\n",
    "        return np.mean(MSE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "photons = 1000000, pitch = 16,L = 64, M = 128, MSE = 1.109946.\n"
     ]
    }
   ],
   "source": [
    "num_trial = 10 # number of signals in simulation\n",
    "photons = 1e6 # the expected number of photons per element\n",
    "XT = np.random.randn(num_trial,N,N)+1j*np.random.randn(num_trial,N,N);\n",
    "XT = np.sqrt(photons)*XT/np.sqrt(2.)\n",
    "run_simulation(XT,num_projection_loop=100,num_gd_iteration=50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MSE vs L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "photons = 1000000, pitch = 4,L = 16, M = 32, MSE = 1.036618.\n",
      "photons = 1000000, pitch = 8,L = 32, M = 64, MSE = 1.030497.\n",
      "photons = 1000000, pitch = 16,L = 64, M = 128, MSE = 1.023170.\n",
      "photons = 1000000, pitch = 32,L = 128, M = 256, MSE = 1.021987.\n"
     ]
    }
   ],
   "source": [
    "num_trial = 1\n",
    "photons = 1e6\n",
    "XT = np.random.randn(num_trial,N,N)+1j*np.random.randn(num_trial,N,N);\n",
    "XT = np.sqrt(photons)*XT/np.sqrt(2.)\n",
    "\n",
    "MSE = np.zeros([4])\n",
    "L_list = [16,32,64,128]\n",
    "for idx in range(len(L_list)):\n",
    "    ## dimensions\n",
    "    L = L_list[idx]\n",
    "    N = 256\n",
    "    pad_size = int(L/2)\n",
    "    M = 2 * L\n",
    "    pitch = int(L/4)\n",
    "    total_frame = int(N/pitch)**2\n",
    "\n",
    "    is_mask = True\n",
    "    is_noisy = True \n",
    "    \n",
    "    ## random phase mask\n",
    "    mask = np.exp(1j * np.random.random([L,L]) * 2 * np.pi)\n",
    "    \n",
    "    ## selected windows\n",
    "    mesh = np.arange(0,N,pitch)\n",
    "    xlocations,ylocations = np.meshgrid(mesh,mesh);\n",
    "    xlocations = xlocations.astype(int).flatten()\n",
    "    ylocations = ylocations.astype(int).flatten()\n",
    "\n",
    "    ## Energy constraint\n",
    "    counts = np.zeros([2*N,2*N])\n",
    "    for i in range(total_frame):\n",
    "         counts[xlocations[i]:xlocations[i]+L,ylocations[i]:ylocations[i]+L] += 1\n",
    "    counts = counts[0:N,0:N] + counts[0:N,N:2*N] + counts[N:2*N,0:N] + counts[N:2*N,N:2*N]\n",
    "    \n",
    "    ## define the tensorflow model\n",
    "    tf.compat.v1.reset_default_graph()\n",
    "    X_FT_EST = tf.compat.v1.placeholder(tf.complex128, shape=(None,L,L))\n",
    "    paddings = tf.constant([[0,0],[pad_size, pad_size], [pad_size, pad_size]])\n",
    "    MEASUREMENT = tf.compat.v1.placeholder(tf.float64, shape=(None,M,M))\n",
    "    if is_mask:\n",
    "        MASK = tf.constant(mask)\n",
    "        X_FT_EST_PAD = tf.pad(tf.multiply(X_FT_EST,MASK),paddings)\n",
    "    else:\n",
    "        X_FT_EST_PAD = tf.pad(X_FT_EST,paddings)\n",
    "    FORWARD = tf.abs(tf.signal.ifft2d(tf.signal.ifftshift(X_FT_EST_PAD,[1,2]))) * M\n",
    "    loss = tf.reshape(tf.cast(tf.math.reduce_sum(tf.math.square(abs(MEASUREMENT)-abs(FORWARD))),tf.complex128),[1,1])\n",
    "    weight = [[1 + 1j*0]]\n",
    "    gradient = tf.gradients(loss, X_FT_EST, grad_ys = weight)\n",
    "    sess = tf.compat.v1.Session()\n",
    "    \n",
    "    ## simulations\n",
    "    ## num_projection_loop is set based on experience\n",
    "    if L == 16:\n",
    "        MSE[idx] = run_simulation(XT,num_projection_loop = 200,num_gd_iteration = 100,return_error = True)\n",
    "    else:\n",
    "        MSE[idx] = run_simulation(XT,num_projection_loop = int(1000/L),num_gd_iteration = 20,return_error = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEWCAYAAACnlKo3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deZwdZZ3v8c+3t3Q6W4ekA9kgIBBUBhAXcBRscQNGxdkHFxbHiV6XO77U8TqLgtvV6ziLXke4cWRihInOcF1QQcEZW9QxbCIRZDFKQjoLITudvTu/+eOpTqpPTp/TXemT053+vl+v8+pT9VTVeU7X6fOtep6nqhURmJmZFdFQ7wqYmdnY5RAxM7PCHCJmZlaYQ8TMzApziJiZWWEOETMzK8whchRJul7SB+tdj+GQdK2kG+tdj6GQtErSy+tdDxudjuTvT1JIOnWk63QscIhw8Mtnt6Se3GPOSL9ORLwtIj46UtvLPtg7s/qulfT3khqzMn+hjgBJC7Lf860l82+UdO0Iv9aVku6TtENSt6RPSWrKlXdJ2pP7jD5asn6HpH+VtF3SVkk3DfF1+9/j/SXzZ0raJ2nViLzBQ9s9XdI3JT0laYuk70lamCuvyYFL/u9PUqek7pF+jSMhaYmkj1Uof6Gk/yqZ95Js3w26Xq05RA55TURMzj3WjeTG+7/cC67bVKH47IiYDLwMeD3wZ0Vfp56qvMfR4DxJv13j12gD3g3MBM4j7dP3lSzzztxndGFJ2deADcCJwCzg08N9fUln5qZfDzw+zG0MRTtwC7AQOB64G/hmDV7nWPM7wMGDGUnNwGeAu+pWIxwiFUmaIOkfJa3LHv8oaUJWdpWkH5csf/CUNzuquE7SrZJ2Ai8tPdKQ9GpJP5e0TdJ/STorV7ZK0v+StALYWe1LNiIeAX4EnCnpy6Qvkm9lR6zvl/QdSe8qqe8KSb+bPX+2pDuyI8MnJf1VbtEWSUslPS3pIUnPy23jmdkR8ras7LW5shmSvpUdWd8j6WP531n2+3qHpF8Bv8rmfUbSmmyd+yRdkFv+Wkk3S/pqVpefSTq75FdxTva+tmfLtWbrPijpNbltNUvaJOk5lX6vOZ8CPj7EZQuJiOsi4kcRsS8i1gI3AS8ayrqSXgnMB/4iIrZHxP6IuL/aeiW+DFyZm74CWDrMbVQVEXdHxBcjYktE7Af+AVgoacaRbLf/7ELSX2X7dpWkN+TKl2SfwUnAbcCc3FndHEmN2bq/zj5f90man3uJl0v6VfZZ/ydJOpL6FnApuRAB3gvcDjxylOsxgEOksr8GzgfOAc4GXgD8zTDWfz3pi2cKUBo4zwFuAN4KzAD+H3BLf0hlLicdfbRHRG+lF5L0LOAC4P6IeBPwBIfOrj4FfAl4Y275s4G5wHckTQG+D3wXmAOcCvxHbvOvBb7CoSPIz2XbaAa+RfogzwLeBdyUa5r4J2AncALpyyn/BdXvdaSj7mdl0/eQft/HAf8K/Ht/EGQuA/49V/6NrB79/gi4GDgZOAu4Kpu/NP/+SX+Q6yPifknflvSBMnXL+zxwuobQRCjpxdkXzWCPF1fbRuZC4KGSeZ/IviB/IqkzN/984FHgS5I2Z6H9kiG+Tr8bgT/JvkyfBUymylFuFtiDvc/PD/F1LwQ2RMTmYda3nBNIZ3JzSZ+3xbnPIwARsRO4BFhX0vLwHtLf3KXAVODNwK7cqq8Gnk/6XP0R8KpyFZD0+ir7/8ThvilJs0lnbfdn0ydl9fvIcLc14iJi3D+AVUAPsC17fCOb/2vg0txyrwJWZc+vAn5csp0ATs2eLwGWlpQvAT6WPb8O+GhJ+aPAS3J1enOVegewA9ia1fVjQENu/Zfnlm3Nljstm/408Pns+eWk8Cn3GtcC389NPwvYnT2/gNR80pArX5at0wjsBxbmyj6W/51l9b+oynvcSmqy66/L8lxZA7AeuCD3nt+YK/8UcH32fA7wNDA1m74ZeP8QPhsLsno2AW/vf33SF+61NfxMvhnoBmbm5p1HOiCZQPqCfBp4Rla2OKvnnwLNwJ+QPsszh/Ba+ff4fdLn/JOkg6iXk33ma/Q+5wFrgctLPnM3FthWJ9ALTMrN+zfgg2X+/jqB7jJ/f5dV+Ft7ccl2PzDCv4uD9StT9qfAF3PT3wT+uNp6R+PhM5FDXhcR7dnjddm8OcDq3DKrs3lDtaZC2UnAe/NHKKTmiPz2K63f79yImB4Rz4iIv4mIA+UWiog9wFeBN0pqIAXHl7Pi+aQQGsyG3PNdQKtS89ocYE3Ja64mHQV2kL6U8u+h3PsZME/S+yQ9nDVHbQOmkY4sD1s+e91uBv7OSus6OVt2HfAT4PcltZOORIfU8Zzzz8Dx+WaxWpD0OuATwCURsal/fkTcFRFPR8TeiPgS6f1cmhXvJn3ZfzFSU9ZXSL+rITWH5SwlHSDlPx81IamDdBb7+YhYNkKb3RrpTKPfcP5mh/t3MHmYdTsSB5uyss/flIj46lF8/UE5RCpbR/qy73diNg9SM01bf4GkE8qsX+kWyWuAj+eCqz0i2kr+mI7kFsvl1v0S8AZSh+2uiPhpri6nFHiNdcD8LJT6nUg6snyKdFQ4L1eWb18+rJ5Z/8f7SU0F0yOiHdgO5Nue5+eWb8i2P9RBEP1Nen8I/DRSv8OQRcQ+4MPAR0vqNICkCzRwpF/p44IK614MfIHUFPmLalXK1WMFh+/zIp+f/09qQv1NRDxRbWGlfrDB3uf1FdabTgqQWyJiJPuapmd9Hv3yf7N55X43a4BnHGkFJL2hyv4fVnNW1lz7EuCObNbLgOdJ2iBpA/DHwLsl1WVwgkOksmXA3ygNnZwJfIjUjAHwAPBsSedkbfbXDnPbXwDeJuk8JZMk/U7WPzESnqQkGLLQOAD8HQOPMr8NzJb0bqXBBFMknTeE17iLdET2fqWO6k7gNcBXIqKPNFroWkltks4gddRWMoUUPE8BTZI+RGqbznuupN/LzoTeDewFlg+hrgDfAM4F/pziHcZfJjUNXjzYApE6xydXePyo3HqSLiKdHf1+RNxdUtYu6VWSWiU1KXUYX0jqxwL4OukL9MqsT+MPSAH7k2z9ayV1VXtz2VH8RcBbqi2bLf/sCu/zbYO8z6nA94CfRMRgfVEN2Xvtf/QPaFkiaUmVan1YUksW1q8m9aGVehKYIWlabt4/Ax+VdFr2N3mWCnT2R8RNVfZ/pXBuLHnfLcCLgRURsSNb5oPA6aS+w3NI/ZRfAK4ebl1HgkOkso8B95KO8n4B/CybR0Q8RurU+j5pZNGPB9lGWRFxL2k47udI7f4rOdQJPBI+QQrAbZLyw0SXAr/FoTAkIp4GXkEKgA2k9/PSai+QHZm/htQ0tInU+XxFpJFiAO8kNUdtIH35LiN96Q/me6QvxcdIzRB7OLwJ7JukI6+twJuA34s0wqeqiNhNOtI+mRRwAEi6TQNHo1XaRh/pYOK4oSw/TB8k/b5uzR213paVNZM+e0+RftfvIjXBPpbVawtpAMT7SGdvHyC17/c3h80nC5RqIuLeiKjUrHOkfpfUQX11hSP0y0lNdP2P/vpUex8bSJ+NdaRAflvu83hQNm8Z8Jvsb2QO8Pekvo7bSX2NXwQmFn+bhXyAge/7PykZ2ps1aW7of2TL7cw+A0edso4ZGyckXQEsioihjhAaydf+P8AJEVFulNZQ1r+WNHDhjdWWrbCNDwGnH8k2xiJJPwdeFiMzAqousqPyB4Czyh04ZGfCN0bEvNKysUzSL4E/iIhf1rsu5Yz2C7xsBElqI40wGurQyyN9vTOAFtJZ3PNJI0yG1ExSo/ocl9XhTfWqQ71ExDn1rsORys58n1nvehxNWXAuHa0BAjVszpJ0g6SNkh4cpPwMST+VtLekuQVJF0t6VNJKVR+/b0Mg6VWkppAnSddXHA1TSM1GO0kjw/6OOl2ZLOnPSE1jt0XEnfWog9lwRbrw9JP1rkclNWvOknQh6dqLpRFxZpnyWaSRT68jDcv7dDa/kdQm/grS8M17SGPIR20Sm5mNVzU7E8mO9gbt6ImIjRFxD+mCtLwXACsj4jfZ6etXSFcpm5nZKDMa+0TmMnBETjfpSt3DSFoELAJobW197oknDvtuAjZKHDhwgIYGDxYcq7z/xq7HHntsU0R0FF1/NIbIkEXEYtLtHli4cGE8+uijVdaw0aqrq4vOzs56V8MK8v4buyStrr7U4EbjocNaBl7Z3H9vHTMzG2VGY4jcA5wm6eRseNufkK7INDOzUaZmzVmSlpHulDlT6T+IXUO66paIuF7pXlP3km5rcUDSu4FnRcQOSe8kXb3cCNwQEaW3wzYzs1GgZiESEZdXKd/AwJvz5ctuZeA/XzEzs1FoNDZnmZnZGOEQMTOzwhwiZmZWmEPEzMwKc4iYmVlhDhEzMyvMIWJmZoU5RMzMrDCHiJmZFeYQMTOzwhwiZmZWmEPEzMwKc4iYmVlhDhEzMyvMIWJmZoU5RMzMrDCHiJmZFeYQMTOzwhwiZmZWmEPEzMwKc4iYmVlhDhEzMyvMIWJmZoU5RMzMrDCHiJmZFeYQMTOzwhwiZmZWmEPEzMwKc4iYmVlhDhEzMyvMIWJmZoU5RMzMrDCHiJmZFeYQMTOzwhwiZmZWWM1CRNINkjZKenCQckn6rKSVklZIOjdX9ilJD0l6OFtGtaqnmZkVV8szkSXAxRXKLwFOyx6LgOsAJP028CLgLOBM4PnAS2pYTzMzK6hmIRIRdwJbKixyGbA0kuVAu6TZQACtQAswAWgGnqxVPc3MrLimOr72XGBNbrobmBsRP5X0A2A9IOBzEfFwuQ1IWkQ6i6Gjo4Ourq7a1thqpqenx/tvDPP+G7/qGSJlSToVeCYwL5t1h6QLIuJHpctGxGJgMcDChQujs7PzqNXTRlZXVxfef2OX99/4Vc/RWWuB+bnpedm83wWWR0RPRPQAtwEvrEP9zMysinqGyC3AFdkorfOB7RGxHngCeImkJknNpE71ss1ZZmZWXzVrzpK0DOgEZkrqBq4hdZITEdcDtwKXAiuBXcDV2ao3AxcBvyB1sn83Ir5Vq3qamVlxNQuRiLi8SnkA7ygzvw94a63qZWZmI8dXrJuZWWEOETMzK8whYmZmhTlEzMysMIeImZkV5hAxM7PCHCJmZlaYQ8TMzApziJiZWWEOETMzK8whYmZmhTlEzMysMIeImZkV5hAxM7PCHCJmZlaYQ8TMzApziJiZWWEOETMzK8whYmZmhTlEzMysMIeImZkV5hAxM7PCHCJmZlaYQ8TMzApziJiZWWEOETMzK8whYmZmhTlEzMysMIeImZkV5hAxM7PCHCJmZlaYQ8TMzApziJiZWWEOETMzK6xmISLpBkkbJT04SLkkfVbSSkkrJJ2bKztR0u2SHpb0S0kLalVPMzMrrpZnIkuAiyuUXwKclj0WAdflypYCfxsRzwReAGysUR3NzOwINNVqwxFxZ5UziMuApRERwHJJ7ZJmA9OBpoi4I9tOT63qaGZmR6ZmITIEc4E1uenubN48YJukrwEnA98HPhARfaUbkLSIdBZDR0cHXV1dta6z1UhPT4/33xjm/Td+1TNEBtMEXAA8B3gC+CpwFfDF0gUjYjGwGGDhwoXR2dl51CppI6urqwvvv7HL+2/8quforLXA/Nz0vGxeN/DziPhNRPQC3wDOLbO+mZnVWT1D5BbgimyU1vnA9ohYD9wDtEvqyJa7CPhlvSppZmaDq1lzlqRlQCcwU1I3cA3QDBAR1wO3ApcCK4FdwNVZWZ+k9wH/IUnAfcAXalVPMzMrrpajsy6vUh7AOwYpuwM4qxb1MjOzkeMr1s3MrDCHiJmZFeYQMTOzwhwiZmZWmEPEzMwKc4iYmVlhDhEzMyvMIWJmZoVVDBFJb8w9f1FJ2TtrVSkzMxsbqp2JvCf3/P+WlL15hOtiZmZjTLUQ0SDPy02bmdk4Uy1EYpDn5abNzGycqXYDxjMkrSCddTwje042fUpNa2ZmZqNetRB55lGphZmZjUkVQyQiVuenJc0ALgSeiIj7alkxMzMb/aoN8f22pDOz57OBB0mjsr4s6d1HoX5mZjaKVetYPzkiHsyeXw3cERGvAc7DQ3zNzMa9aiGyP/f8ZaR/aUtEPA0cqFWlzMxsbKjWsb5G0ruAbuBc4LsAkiaS/b90MzMbv6qdifwp8GzgKuCPI2JbNv984F9qWC8zMxsDqo3O2gi8rcz8HwA/qFWlzMxsbKgYIpJuqVQeEa8d2eqYmdlYUq1P5IXAGmAZcBe+X5aZmeVUC5ETgFcAlwOvB74DLIuIh2pdMTMzG/0qdqxHRF9EfDciriR1pq8Euvy/RMzMDKqfiSBpAvA7pLORBcBnga/XtlpmZjYWVOtYXwqcSbrI8MO5q9fNzMyqnom8EdgJ/DnwP6WD/eoCIiKm1rBuZmY2ylW7TqTaxYhmZjaOOSTMzKwwh4iZmRXmEDEzs8IcImZmVphDxMzMCnOImJlZYTULEUk3SNooqewFiko+K2mlpBWSzi0pnyqpW9LnalVHMzM7MrU8E1kCXFyh/BLgtOyxCLiupPyjwJ01qZmZmY2ImoVIRNwJbKmwyGXA0kiWA+2SZgNIei5wPHB7repnZmZHruoNGGtoLul/lfTrBuZKehL4O9ItV15eaQOSFpHOYujo6KCrq6s2NbWa6+np8f4bw7z/xq96hshg3g7cGhHduXt1lRURi4HFAAsXLozOzs7a185qoqurC++/scv7b/yqZ4isBebnpudl814IXCDp7cBkoEVST0R8oA51NDOzCuoZIrcA75T0FeA8YHtErAfe0L+ApKuA5zlAzMxGp5qFiKRlQCcwU1I3cA3QDBAR15P+R8mlpP+WuAu4ulZ1MTOz2qhZiETE5VXKA3hHlWWWkIYKm5nZKDQaO9YL2bInWHznr5nb3sac9lbmTp/IzEkTaGio3DlvZmbFHTMh0rMv+N+3PjJgXktjA3PaW5nTPpG57RPTz+np+dz2icxub2VCU2OdamxmNvYdMyFy4tQG7rn2lazbtpu1W3ezbttuurftZt22Pazduosf/WoTTz69h4iB63VMmcCc9onMa5+YzmByYTOvvY2pE5uoNtTYzGy8OmZCBGBqazNTT2jmjBPK/+v3fb0HeHLHHrqzkFm77dDPh9fv4PsPP8ne3gMD1pnU0sjc6RMHnM3My03PmjKBpkbfx9LMxqdjKkSqaWlqYP5xbcw/rq1seUSweee+g2cza0uCZkX3drbs3DdgncYGccLUdAaTwqb1YL9Mf9i0tYyrX7OZjSP+dsuRxMzJE5g5eQJnzWsvu8yufb2piaw/XHJhc/fjW9iwYw99Bwa2mU1vax70TGbu9InMmNTiJjMzG5McIsPU1tLEqbMmc+qsyWXL+w4ET+7Yc/DsZW2uj2bV5p38ZOUmdu7rG7DOhKaGQ30xuT6ZOe2tzGtv44RprbQ0ucnMzEYfh8gIa2wQc7IgeF6Z8ohgx+7eXMDsYt32PQfPaP7z0Y089fTeAetIMCsbADA3dwYzZ9rEg/010yY2H503aGaW4xA5yiQxra2ZaW3NPGtO+QEAe3v7WL9tT26EWXY2s303D67dzu0PPcm+voEDAKZMaDpsAEAazpz6aDqmTKDR18yY2QhziIxCE5oaWTBzEgtmTipbfuBAsGnn3mz48m7WbtvFum2HRp3dt3or23fvH7BOU4OY3d568OxlbmnT2bSJTGzxNTNmNjwOkTGooUHMmtLKrCmtnDO//ACAnr29h/plSoY0L//1Zjbs2ENJ/z8zJrUMeiYzp72V4zwAwMxKOESOUZMnNHH68VM4/fgpZct7+w6wYceebKTZruyMJo06W/lUDz987Cl27x84AKC1ueFgyMwr6ZOZ2z6RE6a10uxrZszGFYfIONXU2MC86W3Mm94GHHdYeUSwbdf+w0aY9Z/N3LF+B5t6Bl4z0yA4fmq528wcOpuZ0uoBAGbHEoeIlSWJ6ZNamD6phTPnTiu7zJ79fazrv7XMgLOZXfx8zTZue3A9+/sGtplNbW067FqZOe0TeWpbH8/ZtZ9pbQ4Zs7HEIWKFtTY3ckrHZE7pGPyamU09ew+/zczW3XRv3c1dj2/h6T29B5f/yPLbaW9rZsGMSSyY0ZYGF8xIAwxOnjHJAWM2CjlErGYaG8TxU1s5fmorzz1petllduzZz7ptu7n1h3czZfYpPL55J6s37+SeVVv55gPrBtww0wFjNvo4RKyu+m+aueH4JjovPGVA2Z79fazZsovHN+1k9eZdhQJmwYw22ttajvK7Mhs/HCI2arU2N3La8VM4rcwIs/6AWbV5F6s27awaMCfNmMTJDhizEecQsTFpOAGzanN6VAuYk2ZM4uSZDhiz4XCI2DGnlgFz0ow2Tp45yQFjlnGI2LhSLWC6t+7i8U27WL15J49vGjxgpk1szjr1HTA2vjlEzDKtzY2cOmsKp84aWsCs3ryrYsAsmNHGAgeMHeMcImZDUClg9vb2jyIbGDD3rtrKLVUCZsHMQ0HjgLGxyCFidoQmNFUPmFWbdh3sf1m1aXgBs2DGJKZPcsDY6OQQMauhIgFz32oHjI0dDhGzOqkeMLsHjCCrGDADroFxwNjR4xAxG4VSwEzm1FmH35esXMCs3lw9YNIoMgeMjSyHiNkYUzRgvvXAugH/iKxcwKTrYhwwNnQOEbNjyFACJj+CbNXmncMOmPa2Zv+HSzvIIWI2ThQJmJ89cXjATG1tyq59SbeIOem4NlY92Us8spHmxgaaG0VzUwMtjQ00NzbQ1KiDz0vLGhscRmOdQ8TMqgZM99bURJYPmPvXbOXbK3IBc/89w37dBkFzYwqVpkZlQdNAS1MWOP3TjQ00Nx2azpel8my6aeB0U74sK28pXbdJNDUcej5wu4det6lBPgMrwyFiZhVNaGrkGR2TeUaZfz62t7ePddv2cOd/Lefs5zyX3r4D7Os7wP6+YH/vAfbnp/uy6d403ds/nStL5SXT2bL7eg/Qs7fv4Hb3Z9vd1/+899B0rRwWdtmZ1WFhVnoGloViSy6gmps0yNlatmxTyXQWiv0Beqj88Ndpbmig4Sid5TlEzKywCU2NnDxzEqunNnLO/PZ6VweAiKDvQAwMmL4D7O8tmT4s3A4Pr9LAOiy8esuEWTa9Z2/fwdc9FKaHArZ/+kBUf09FNDboUODkgmvAdOORB41DxMyOKZJoahRNjTCRxnpXp6oUeAODal9vyXQWWL0H8uFVEk69A8Os97DyMmE2AmdtNQsRSTcArwY2RsSZZcoFfAa4FNgFXBURP5N0DnAdMBXoAz4eEV+tVT3NzOqpsUE0NjTS2lyfwLvxLUe2fsPIVKOsJcDFFcovAU7LHotIwQEpUK6IiGdn6/+jpNFxnmxmZgPU7EwkIu6UtKDCIpcBSyMigOWS2iXNjojHcttYJ2kj0AFsq1VdzcysmHr2icwF1uSmu7N56/tnSHoB0AL8utwGJC0incXQ0dFBV1dXrepqNdbT0+P9N4Z5/41fo7ZjXdJs4MvAlRFRtvcnIhYDiwEWLlwYnZ2dR6+CNqK6urrw/hu7vP/Gr1r2iVSzFpifm56XzUPSVOA7wF9HxPI61M3MzIagniFyC3CFkvOB7RGxXlIL8HVSf8nNdayfmZlVUcshvsuATmCmpG7gGqAZICKuB24lDe9dSRqRdXW26h8BFwIzJF2VzbsqIn5eq7qamVkxtRyddXmV8gDeUWb+jcCNtaqXmZmNnHo2Z5mZ2RjnEDEzs8IcImZmVphDxMzMCnOImJlZYQ4RMzMrzCFiZmaFOUTMzKwwh4iZmRXmEDEzs8IcImZmVphDxMzMCnOImJlZYQ4RMzMrzCFiZmaFOUTMzKwwh4iZmRXmEDEzs8IcImZmVphDxMzMCnOImJlZYQ4RMzMrzCFiZmaFOUTMzKwwh4iZmRXmEDEzs8IcImZmVphDxMzMCnOImJlZYQ4RMzMrzCFiZmaFOUTMzKwwh4iZmRXmEDEzs8IcImZmVljNQkTSDZI2SnpwkHJJ+qyklZJWSDo3V3alpF9ljytrVUczMzsytTwTWQJcXKH8EuC07LEIuA5A0nHANcB5wAuAayRNr2E9zcysoJqFSETcCWypsMhlwNJIlgPtkmYDrwLuiIgtEbEVuIPKYWRmZnXSVMfXngusyU13Z/MGm38YSYtIZzEAewdrOrMxYSawqd6VsMK8/8auhUeycj1D5IhFxGJgMYCkeyPieXWukhXk/Te2ef+NXZLuPZL16zk6ay0wPzc9L5s32HwzMxtl6hkitwBXZKO0zge2R8R64HvAKyVNzzrUX5nNMzOzUaZmzVmSlgGdwExJ3aQRV80AEXE9cCtwKbAS2AVcnZVtkfRR4J5sUx+JiEod9P0Wj+gbsKPN+29s8/4bu45o3ykiRqoiZmY2zviKdTMzK8whYmZmhY3JEBnsliqS3iXpEUkPSfpUvepng5PUKuluSQ9k++nD2fybJD0q6cFs/zbXu65WnqR2STdnf2sPS3phruy9kkLSzHrW0Q4p930p6W+z/bdC0tcltWfzmyV9SdIvsn37l9W2PyZDhDK3VJH0UtJV8GdHxLOBT9ehXlbdXuCiiDgbOAe4OBuddxNwBvBbwETgLfWrolXxGeC7EXEGcDbwMICk+aTRlE/UsW52uCUcftePO4AzI+Is4DGgPyz+EJgQEb8FPBd4q6QFlTY+JkNkkFuq/A/gkxGxN1tm41GvmFWV3eamJ5tszh4REbdmZQHcTbo+yEYZSdOAC4EvAkTEvojYlhX/A/B+wKN1RpFy35cRcXtE9GaTyzn09xbAJElNpIO5fcCOStsfkyEyiNOBCyTdJemHkp5f7wpZeZIaJf0c2Ei6T9pdubJm4E3Ad+tVP6voZOAp4F8k3S/pnyVNknQZsDYiHqhz/Wz43gzclj2/GdgJrCedUX662iUWx1KINAHHAecDfwH8myTVt0pWTkT0RcQ5pKOfF0g6M1f8eeDOiPhRfWpnVTQB5wLXRcRzSF841wJ/BXyojvWyAiT9NdBLak6GdOf0PmAO6YDhvZJOqbSNYylEuoGvZS0idwMHSDeFs1Eqawb5AVl7raRrgA7gPfWslxc9ZZ0AAAGGSURBVFXUDXTnzh5vJoXKycADklaRDg5+JumE+lTRhkLSVcCrgTfEoQsGX0/q79qfdQn8BKh4T7RjKUS+AbwUQNLpQAu+q+ioI6kjNxJkIvAK4BFJbyH9G4DLI+JAPetog4uIDcAaSf13fn0Z8LOImBURCyJiASlozs2WtVFI0sWk/qvXRsSuXNETwEXZMpNILTuPVNrWmLyL7yC3VLkBuCEbxrYPuDKXrjZ6zAa+JKmRdBDzbxHxbUm9wGrgp1kr5Nci4iN1rKcN7l3ATZJagN+Q3bLIRqdBvi//EpgA3JH9vS2PiLcB/0Tq73oIEPAvEbGi4vb9PWtmZkUdS81ZZmZ2lDlEzMysMIeImZkV5hAxM7PCHCJmZlaYQ8SshiT1VF/KbOxyiJiZWWEOETMzK8whYmZmhTlEzMysMIeImZkV5hAxM7PCHCJmtdUmqTv38P9KsWOK7+JrZmaF+UzEzMwKc4iYmVlhDhEzMyvMIWJmZoU5RMzMrDCHiJmZFeYQMTOzwv4bIn3HLOcCoEUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(MSE)\n",
    "plt.xticks([0, 1, 2, 3], ['16', '32', '64','128']) \n",
    "plt.grid()\n",
    "plt.xlim([0,3])\n",
    "plt.ylim([1.0,1.1])\n",
    "plt.title('Fourier Ptychography: N = 256, M = 2L, pitch = L/4')\n",
    "plt.xlabel('L')\n",
    "plt.ylabel('MSE')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}