Factorial Design.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "040856e7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from itertools import product,combinations\n",
    "import matplotlib.pyplot as plt\n",
    "from math import log\n",
    "import scipy.stats"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "160554a9",
   "metadata": {},
   "source": [
    "### major functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "a7046422",
   "metadata": {},
   "outputs": [],
   "source": [
    "def design_matrix(n):\n",
    "    D=[]\n",
    "    for i in product([-1,+1],repeat=n):\n",
    "        D.append(list(i))\n",
    "    return D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "fc56033f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def effect_matrix(d):\n",
    "    dm=np.flip(d,1)\n",
    "    (r,c)=np.shape(dm)\n",
    "    x=''\n",
    "    for i in range(c):\n",
    "        x+=str(i)\n",
    "    labels=[]\n",
    "    for i in range(c-1):\n",
    "        for j in combinations(x,i+2):\n",
    "            labels.append(list(map(int,j)))\n",
    "    for i in range(len(labels)):\n",
    "        I=np.full(r,1)\n",
    "        for j in labels[i]:\n",
    "            I=I*dm[:,j]\n",
    "        dm=np.insert(dm,np.shape(dm)[1],I,1)\n",
    "    return (dm,labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "c629b440",
   "metadata": {},
   "outputs": [],
   "source": [
    "def effects(e,R):\n",
    "    avg=np.mean(R,axis=0)\n",
    "    e=e.astype('float64')\n",
    "    for i in range(len(e)):\n",
    "        for j in range(len(e[i])):\n",
    "            e[i][j]=e[i][j]*avg[i]\n",
    "    return e"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "0fa06414",
   "metadata": {},
   "outputs": [],
   "source": [
    "def eff_det(E,labels,R):\n",
    "    AVG=np.mean(R)\n",
    "    e=[(\"AVG\",AVG)]\n",
    "    labels2=[]\n",
    "    for i in labels[-1]:\n",
    "        labels2.append([i])\n",
    "    labels2.extend(labels)\n",
    "    (r,c)=np.shape(E)\n",
    "    for i in range(c):\n",
    "        e.append((l2str(labels2[i]),np.sum(E[:,i])/(2**(len(labels[-1])-1))))\n",
    "    if len(R)==2:\n",
    "        N=2*len(R[0])\n",
    "        s2=np.sum((R[0]-R[1])**2)/(2*len(R[0]))\n",
    "        V_eff=(4/N)*s2\n",
    "        sd_eff=np.sqrt(V_eff)\n",
    "        V_ave=(s2/N)\n",
    "        sd_ave=np.sqrt(V_ave)\n",
    "    return (e,sd_eff,sd_ave)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "62ab3d08",
   "metadata": {},
   "outputs": [],
   "source": [
    "def significants(eff_ints,sd_ave,sd_eff,a):\n",
    "    t_crit=scipy.stats.t.ppf(q=1-a,df=len(eff_ints))\n",
    "    sig=[]\n",
    "    if eff_ints[0][1]/sd_ave>=t_crit:\n",
    "        sig.append(eff_ints[0])\n",
    "    for i in range(1,len(eff_ints)):\n",
    "        if eff_ints[i][1]/sd_eff>=t_crit:\n",
    "            sig.append(eff_ints[i])\n",
    "    return sig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "cb84680e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def seq_eval(sig,seq=''): #,extremes):\n",
    "    if seq:\n",
    "        val=sig[0][1]\n",
    "        for i in range(len(seq)):\n",
    "            if seq[i]==\"+\":\n",
    "                val+=(sig[i+1][1]/2)*1 #extremes[i][1]\n",
    "            elif seq[i]==\"-\":\n",
    "                val+=(sig[i+1][1]/2)*(-1) #extremes[i][0]\n",
    "    return val            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "0f0d253b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def model_eval(sig,seq_list): #,extremes):\n",
    "    L=[]\n",
    "    for seq in seq_list:\n",
    "        L.append((seq,seq_eval(sig,seq)))\n",
    "    return L  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "1ab8bec4",
   "metadata": {},
   "outputs": [],
   "source": [
    "#all main effects considered\n",
    "def model(dm,eff_ints,order=1): #,extremes):\n",
    "    n=log(len(dm),2)\n",
    "    est=[]\n",
    "    for i in range(len(dm)):\n",
    "        val=eff_ints[0][1]\n",
    "        for j in range(len(dm[i])):\n",
    "            val+=((eff_ints[j+1][1])/2)*dm[i][j] #extremes[j][1/0]\n",
    "        est.append(val)\n",
    "    return (est)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "a72f1b31",
   "metadata": {},
   "outputs": [],
   "source": [
    "def resi(A,R):\n",
    "    avg=np.mean(R,axis=0)\n",
    "    res=[]\n",
    "    for i in range(len(avg)):\n",
    "        res.append(avg[i]-A[i])\n",
    "    return res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "44b26dc8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def prob_plot(A,label,title=''):\n",
    "    m=len(A)\n",
    "    Pi=[(100/m)*(i-0.5) for i in range(1,m+1)]\n",
    "    A.sort()\n",
    "    plt.scatter(A,Pi)\n",
    "    plt.axvline(x=0, c=\"black\")\n",
    "    plt.axhline(y=0, c=\"black\")\n",
    "    plt.xlabel(\"Probability\")\n",
    "    plt.ylabel(label)\n",
    "    plt.title(label+\" Vs \"+ \"Probability\"+title)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1dba6faf",
   "metadata": {},
   "source": [
    "### minor functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "16850c50",
   "metadata": {},
   "outputs": [],
   "source": [
    "def l2str(A):\n",
    "    x=''\n",
    "    for i in A:\n",
    "        x+=str(i)\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "6eae0711",
   "metadata": {},
   "outputs": [],
   "source": [
    "def dm2strL(dm):\n",
    "    L=[]\n",
    "    for i in range(len(dm)):\n",
    "        x=''\n",
    "        for j in range(len(dm[i])):\n",
    "            if dm[i][j]==1:\n",
    "                x+='+'\n",
    "            elif dm[i][j]==-1:\n",
    "                x+='-'\n",
    "        L.append(x)\n",
    "    return (L)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c3c86653",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "f6e9a968",
   "metadata": {},
   "source": [
    "### Model formation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "fe00adb8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('AVG', 6.679375),\n",
       " ('0', 5.34125),\n",
       " ('1', -1.0387500000000003),\n",
       " ('2', 5.37875),\n",
       " ('01', -1.1112500000000005),\n",
       " ('02', 6.601249999999999),\n",
       " ('12', 1.4412500000000001),\n",
       " ('012', -1.2112500000000002)]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R=np.array([[7.28,5.03,2.64,2.30,2.54,19.25,5.21,14.41],\n",
    "            [4.54,4.07,4.02,2.04,1.53,13.35,4.31,14.35]])\n",
    "\n",
    "n=3\n",
    "A=design_matrix(n)\n",
    "A=np.array(A)\n",
    "E=effect_matrix(A)\n",
    "E1=effects(E[0],R)\n",
    "(eff_ints,sd_eff,sd_ave)=eff_det(E1,E[1],R)\n",
    "eff_ints"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "97b5b759",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([[-1, -1, -1,  1,  1,  1, -1],\n",
       "        [ 1, -1, -1, -1, -1,  1,  1],\n",
       "        [-1,  1, -1, -1,  1, -1,  1],\n",
       "        [ 1,  1, -1,  1, -1, -1, -1],\n",
       "        [-1, -1,  1,  1, -1, -1,  1],\n",
       "        [ 1, -1,  1, -1,  1, -1, -1],\n",
       "        [-1,  1,  1, -1, -1,  1, -1],\n",
       "        [ 1,  1,  1,  1,  1,  1,  1]]),\n",
       " [[0, 1], [0, 2], [1, 2], [0, 1, 2]])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E #effect matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56ce45e2",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "25a2a3ea",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.8573660609681257"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sd_eff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "d6633c53",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.42868303048406287"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sd_ave"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "50a1fe90",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('AVG', 6.679375), ('0', 5.34125), ('2', 5.37875), ('02', 6.601249999999999)]"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a=0.05\n",
    "sig=significants(eff_ints,sd_ave,sd_eff,a)\n",
    "sig"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd9cdbe7",
   "metadata": {},
   "source": [
    "### Residuals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "3d142ed1",
   "metadata": {},
   "outputs": [],
   "source": [
    "B=design_matrix(len(sig)-1)\n",
    "seq_list=dm2strL(B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "cf13cc9d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('---', -1.9812499999999988),\n",
       " ('--+', 4.620000000000001),\n",
       " ('-+-', 3.3975000000000013),\n",
       " ('-++', 9.998750000000001),\n",
       " ('+--', 3.36),\n",
       " ('+-+', 9.96125),\n",
       " ('++-', 8.73875),\n",
       " ('+++', 15.34)]"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L=model_eval(sig,seq_list)\n",
    "L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "0c96dcb1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'\\nresiduals=resi(A,R) #fill manually\\nresiduals\\n\\n'"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "residuals=resi(A,R) #fill manually\n",
    "residuals\n",
    "\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "073740e8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"\\nprob_plot(residuals,'Residuals')\\n\\n\""
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "prob_plot(residuals,'Residuals')\n",
    "\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "4d678f06",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[4.071249999999999,\n",
       " -2.6675000000000004,\n",
       " 2.5299999999999994,\n",
       " -4.008750000000001,\n",
       " -5.145,\n",
       " 3.741250000000001,\n",
       " -1.3812499999999996,\n",
       " 2.8599999999999994]"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#all main effects considered\n",
    "M=model(A,eff_ints) #,extremes)\n",
    "residuals=resi(M,R)\n",
    "residuals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "284427a1",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjMAAAHFCAYAAAAHcXhbAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAABCFklEQVR4nO3deVxUZf//8feAyO4kJgwkKi6lRN6KJmmZS2rumaWWWmpampaRdZve3oV4m371LvOb5lJ3LqmlraZW3FomWWpYbplLd4XLnSBusSkocH5/+GO+jgMKIzIceT0fj3nkuc6Zcz4zR2feXde5zlgMwzAEAABgUh7uLgAAAOBqEGYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWZwWYsXL5bFYrE/qlSpotDQUD300EP6z3/+c82OO2nSJFkslhJtW7duXQ0ZMuSa1VLaeopz/PhxVa1aVQ899FCx22RkZMjPz0+9evW6qmMVdd5q1aqloUOH6o8//riqfV/KYrHoqaeeKrP9HTx4UBaLRa+88soVty18nQcPHrS3DRkyRHXr1nXY7tK/I0ePHtWkSZO0c+fOq673nnvu0ciRI116rsVi0aRJk+zLGzdulMVi0caNG6+6rrLSrl07tWvXzi3HPnjwoLp3766goCBZLBbFxsZKknbs2KG2bdvKarXKYrFo1qxZZX7sqVOnatWqVS4998UXX1R0dLQKCgrKtigUq4q7C4A5LFq0SI0aNVJOTo6+++47vfzyy/r666+1f/9+Va9evcyPN3z4cHXp0qXM9+tONWvWVK9evbRq1SqdPn26yPdtxYoVOnv2rIYNG1Ymxyw8b2fPntU333yjadOmKTExUT/99JP8/f3L5Bju1L17d23ZskWhoaGX3e6TTz5RtWrV7MtHjx5VfHy86tatq6ZNm7p8/E8//VTfffed3nnnHZf3UdHNnTvXbcd+9tln9f3332vhwoWy2Wz28/zYY48pOztbK1asUPXq1Z3Ca1mYOnWqHnzwQfXu3bvUz33++ec1Z84cLVmyREOHDi3z2uCMMIMSiYqKUosWLSRd+D+1/Px8xcXFadWqVdfkH2utWrVUq1atMt+vuw0bNkwfffSRli9fXmRvxsKFCxUSEqLu3buXyfEuPm/t27dXfn6+/vGPf2jVqlUaOHBgkc85c+aM/Pz8yuT411rNmjVVs2bNK27XrFmza3L8qVOn6v7779dNN910TfZfEURGRrrt2Hv27FHLli2dAsWePXv0+OOPq2vXru4p7AqsVqsGDRqk//mf/9GQIUOuulcXV8YwE1xS+AV57Ngxh/YffvhBvXr1UlBQkHx8fNSsWTO9//77DtucOXNGzz//vCIiIuTj46OgoCC1aNFC7733nn2booZ1zp8/r3Hjxslms8nPz0933XWXkpKSnGorbkioqCGJlStXqnPnzgoNDZWvr68aN26s8ePHKzs7+4rvwYYNG9SuXTvVqFFDvr6+ql27th544AGdOXOm2Ofce++9qlWrlhYtWuS0bt++ffr+++/16KOPqkqVC/+fcfz4cT3xxBMKDw+Xt7e3atasqTvvvFNffvnlFesryh133CFJOnTokKQLQzIBAQH66aef1LlzZwUGBuqee+6RJJ06dUqjRo3STTfdpKpVq6pevXqaOHGicnNzi9z3ggULdPPNN8vb21uRkZFasWKFw/rjx49r1KhRioyMVEBAgIKDg9WhQwdt2rSpyP0VFBTo5ZdfVu3ateXj46MWLVroq6++ctimqHNalIuHmTZu3Kjbb79dkjR06FD7UNykSZO0dOlSWSwWbdmyxWkfkydPlpeXl44ePSrpwlBHUlKSHnnkkat6na4ofN0bNmzQ448/rho1aqhatWp69NFHlZ2drdTUVPXr10833HCDQkND9fzzz+v8+fMO+4iPj1dMTIyCgoJUrVo1RUdH6+2339alvz186TDTxcOAM2fOVEREhAICAtSqVStt3bq1RPWnpqZqxIgRqlWrlqpWraqIiAjFx8crLy9P0v8Nt/3666/64osv7Oeo8HXn5eVp3rx59vaS7rdQbm6uJk+erMaNG8vHx0c1atRQ+/bttXnzZkkXhv+ys7O1ZMkS+zEK34OSfH5J0iOPPKJffvlFX3/9dYneE1wdembgkuTkZEnSzTffbG/7+uuv1aVLF8XExGj+/PmyWq1asWKF+vfvrzNnzti/TMaOHaulS5dqypQpatasmbKzs7Vnzx6dPHnyssd8/PHH9c477+j5559Xp06dtGfPHvXp00eZmZkuv47//Oc/6tatm2JjY+Xv76/9+/dr+vTpSkpK0oYNG4p9XuFYfps2bbRw4ULdcMMN+uOPP5SQkKBz584V27Ph4eGhIUOGaMqUKdq1a5f+8pe/2NcVBpzHHnvM3vbII49o+/btevnll3XzzTfrzz//1Pbt26/4XhXn119/lSSH3oxz586pV69eGjFihMaPH6+8vDzl5OSoffv2+u233xQfH68mTZpo06ZNmjZtmnbu3KnPPvvMYb+rV6/W119/rcmTJ8vf319z587Vww8/rCpVqujBBx+UdCEcSVJcXJxsNpuysrL0ySefqF27dvrqq6+crsuYM2eO6tSpo1mzZqmgoEAzZsxQ165dlZiYqFatWrn0+iUpOjpaixYt0tChQ/X3v//d3gtWq1YtBQcHa9y4cXrjjTccjpGXl6cFCxbo/vvvV1hYmCRp7dq18vT01N133+2w/9K+zqsxfPhw9enTRytWrNCOHTv0t7/9TXl5eTpw4ID69OmjJ554Ql9++aWmT5+usLAwjR071v7cgwcPasSIEapdu7YkaevWrXr66af1xx9/6KWXXrrisd944w01atTIfr3Kiy++qG7duik5OVlWq7XY56Wmpqply5by8PDQSy+9pPr162vLli2aMmWKDh48qEWLFik6OlpbtmzR/fffr/r169uvn4qIiNCWLVvUqlUrPfjgg3ruuedKtV/pwrns2rWrNm3apNjYWHXo0EF5eXnaunWrDh8+rNatW2vLli3q0KGD2rdvrxdffFGS7MOUJf38at68uQICAvTZZ5+pQ4cOV3w/cZUM4DIWLVpkSDK2bt1qnD9/3sjMzDQSEhIMm81m3H333cb58+ft2zZq1Mho1qyZQ5thGEaPHj2M0NBQIz8/3zAMw4iKijJ69+592ePGxcUZF//13LdvnyHJePbZZx22W758uSHJGDx4cLHPvfS1JCcnF3nMgoIC4/z580ZiYqIhydi1a1ex+/zwww8NScbOnTsv+zqK8vvvvxsWi8UYM2aMve38+fOGzWYz7rzzTodtAwICjNjY2FIfo6jztnbtWqNmzZpGYGCgkZqaahiGYQwePNiQZCxcuNDh+fPnzzckGe+//75D+/Tp0w1Jxrp16+xtkgxfX1/7Pg3DMPLy8oxGjRoZDRo0KLbGvLw84/z588Y999xj3H///fb25ORkQ5IRFhZmnD171t6ekZFhBAUFGR07dnR6nRef08GDBxt16tRxOFadOnUc/o5s27bNkGQsWrTIqa64uDijatWqxrFjx+xtK1euNCQZiYmJ9rauXbsajRo1Kvb1Xel1GsaF9y4uLs6+/PXXXxuSjK+//vqy+yx83U8//bRDe+/evQ1JxsyZMx3amzZtakRHRxe7v/z8fOP8+fPG5MmTjRo1ahgFBQX2dW3btjXatm1rXy48P7fddpuRl5dnb09KSjIkGe+9995lax8xYoQREBBgHDp0yKH9lVdeMSQZP//8s72tTp06Rvfu3Z32IckYPXq0S/t95513DEnGW2+9ddk6/f39Hf7OFCrJ51ehO++804iJiSnRtrg6DDOhRO644w55eXkpMDBQXbp0UfXq1fXpp5/ah0N+/fVX7d+/334dRl5env3RrVs3paSk6MCBA5Kkli1b6osvvtD48eO1ceNGnT179orHL+yqvfQ6j379+tlrcMXvv/+uAQMGyGazydPTU15eXmrbtq2kC8M+xWnatKmqVq2qJ554QkuWLNHvv/9e4mNGRESoffv2Wr58uc6dOydJ+uKLL5SamurQKyNdeK8WL16sKVOmaOvWrU5DBVdy8Xnr0aOHbDabvvjiC4WEhDhs98ADDzgsb9iwQf7+/vZelUKFvWuXDvfcc889Dvv09PRU//799euvv+q///2vvX3+/PmKjo6Wj4+PqlSpIi8vL3311VdFvtd9+vSRj4+PfTkwMFA9e/bUN998o/z8/FK9D6Xx5JNPSpLeeuste9ucOXN02223OfTCHD16VMHBwUXuozSv82r06NHDYblx48aS5HTNVePGje1Di4U2bNigjh07ymq12v/uv/TSSzp58qTS0tKueOzu3bvL09PTvtykSRNJcjrOpdauXav27dsrLCzM4XOi8PqXxMTEKx77avb7xRdfyMfHx+nfWkmV5vMrODi4zGcPomiEGZTIO++8o23btmnDhg0aMWKE9u3bp4cffti+vvDameeff15eXl4Oj1GjRkmSTpw4IUl6/fXX9cILL2jVqlVq3769goKC1Lt378tO9S7swrXZbA7tVapUUY0aNVx6TVlZWWrTpo2+//57TZkyRRs3btS2bdv08ccfS9JlP6Tq16+vL7/8UsHBwRo9erTq16+v+vXr63//939LdOxhw4bp5MmTWr16taQLQ0wBAQHq16+fw3YrV67U4MGD9a9//UutWrVSUFCQHn30UaWmppboOIXnbceOHTp69Kh2796tO++802EbPz8/h5k+0oX322azOV17FBwcrCpVqjh1qV96Xi5uK9x25syZevLJJxUTE6OPPvpIW7du1bZt29SlS5ci3+vi9nnu3DllZWWV4NW7JiQkRP3799eCBQuUn5+v3bt3a9OmTU4XbJ89e9YhbBUq7eu8GkFBQQ7LVatWLbY9JyfHvpyUlKTOnTtLuhDavvvuO23btk0TJ060v7YrufTfnbe3d4mee+zYMa1Zs8bpc+LWW2+V9H+fE6VV0v0eP35cYWFh8vBw7euvNJ9fPj4+ZX7OUTSumUGJNG7c2GlWzL/+9S99+OGHevDBB3XjjTdKkiZMmKA+ffoUuY9bbrlFkuTv76/4+HjFx8fr2LFj9v/L6dmzp/bv31/kcws/OFNTUx1mjuTl5Tl9sRZ+weTm5to/YCXnD8kNGzbo6NGj2rhxo703RpL+/PPPK74fktSmTRu1adNG+fn5+uGHHzR79mzFxsYqJCTksveSkS70OlSvXl0LFy5U27ZttXbtWj366KMKCAhw2O7GG2/UrFmzNGvWLB0+fFirV6/W+PHjlZaWpoSEhCvWePF5K05RF0vXqFFD33//vQzDcFiflpamvLw8+/kuVFS4KmwrPHfLli1Tu3btNG/ePIftirvmqbh9Vq1a1el9KmvPPPOMli5dqk8//VQJCQm64YYbnHoFb7zxRvv1MRcr7et0hxUrVsjLy0tr1651CGSu3lelNG688UY1adJEL7/8cpHrC69Julb7rVmzpr799lsVFBS4FGhK8/l16tQpp38ruDbomYFLZsyYoerVq+ull15SQUGBbrnlFjVs2FC7du1SixYtinwEBgY67SckJERDhgzRww8/rAMHDhQ7E6jwosnly5c7tL///vtOMxUK7zmxe/duh/Y1a9Y4LBd+SV8ceKQLs3JKw9PTUzExMXrjjTckSdu3b7/ic3x8fDRgwACtW7dO06dP1/nz56/Y7V27dm099dRT6tSpU4mOcTXuueceZWVlOX25Fd5PpXDGU6GvvvrKYWZbfn6+Vq5cqfr169un2FssFqf3evfu3UXOHJKkjz/+2KE3ITMzU2vWrFGbNm0chjdccaVehObNm6t169aaPn26li9friFDhjjdl6dRo0ZFDi+W9nW6Q+GNFC9+H8+ePaulS5de82P36NFDe/bsUf369Yv8nHA1zJR0v127dlVOTo4WL1582f15e3tfsVflSp9fv//+u1untlcm9MzAJdWrV9eECRM0btw4vfvuuxo0aJAWLFigrl276t5779WQIUN000036dSpU9q3b5+2b9+uDz74QJIUExOjHj16qEmTJqpevbr27dunpUuXqlWrVsXOAmrcuLEGDRqkWbNmycvLSx07dtSePXv0yiuvOA2RdOvWTUFBQRo2bJgmT56sKlWqaPHixTpy5IjDdq1bt1b16tU1cuRIxcXFycvLS8uXL9euXbuu+Prnz5+vDRs2qHv37qpdu7ZycnK0cOFCSVLHjh1L9B4OGzZMb7zxhmbOnKlGjRqpdevWDuvT09PVvn17DRgwQI0aNVJgYKC2bdumhISEYnu/ysqjjz6qN954Q4MHD9bBgwd122236dtvv9XUqVPVrVs3p9d44403qkOHDnrxxRfts5n279/vMD27R48e+sc//qG4uDi1bdtWBw4c0OTJkxUREeEUSKULIbFTp04aO3asCgoKNH36dGVkZCg+Pv6qX1/9+vXl6+ur5cuXq3HjxgoICFBYWJjDF+kzzzyj/v37y2Kx2IdKL9auXTstXLhQv/zyi8OsvtK+Tnfo3r27Zs6cqQEDBuiJJ57QyZMn9corrziFsGth8uTJWr9+vVq3bq0xY8bolltuUU5Ojg4ePKjPP/9c8+fPd+keUyXd78MPP6xFixZp5MiROnDggNq3b6+CggJ9//33aty4sb1X9bbbbtPGjRu1Zs0ahYaGKjAwULfcckuJP79Onjyp//znP3r66afL7L3DZbj7CmRUbIWzJrZt2+a07uzZs0bt2rWNhg0b2mc17Nq1y+jXr58RHBxseHl5GTabzejQoYMxf/58+/PGjx9vtGjRwqhevbrh7e1t1KtXz3j22WeNEydO2LcpakZSbm6u8dxzzxnBwcGGj4+PcccddxhbtmxxmqliGBdmVrRu3drw9/c3brrpJiMuLs7417/+5TTzZfPmzUarVq0MPz8/o2bNmsbw4cON7du3O810ubSeLVu2GPfff79Rp04dw9vb26hRo4bRtm1bY/Xq1aV6f5s1a2ZIMmbMmOG0Licnxxg5cqTRpEkTo1q1aoavr69xyy23GHFxcUZ2dvZl93u583axwYMHG/7+/kWuO3nypDFy5EgjNDTUqFKlilGnTh1jwoQJRk5OjsN2+v8zS+bOnWvUr1/f8PLyMho1amQsX77cYbvc3Fzj+eefN2666SbDx8fHiI6ONlatWuU0+6hwtsz06dON+Ph4o1atWkbVqlWNZs2aGf/+97+LfJ2lnc1kGIbx3nvvGY0aNTK8vLycZhUV1uvt7W106dKlyPcnPT3dCAgIcDp3JX2dhe/d1cxmuvT8Fv49PX78uEN7Ued54cKFxi233GL/Nzht2jTj7bffdno/i5vN9M9//tOprqLex6IcP37cGDNmjBEREWF4eXkZQUFBRvPmzY2JEycaWVlZ9u1KM5upNPs9e/as8dJLLxkNGzY0qlatatSoUcPo0KGDsXnzZvs2O3fuNO68807Dz8/PkGR/D0ry+WUYhvH2228bXl5eDrP8cO1YDOOSOyQBALRmzRr16tVLn332mbp161bkNk8//bS++uor/fzzz9zlFQ7atGmj2rVrOw2N49ogzADARfbu3atDhw7pmWeekb+/v7Zv315sUDl27Jhuvvlmvf32207T2FF5ffPNN+rcubP27t2revXqubucSoELgAHgIqNGjVKvXr1UvXp1vffee5ftcQkJCdHy5cuZfgsHJ0+e1DvvvEOQKUf0zAAAAFOjZwYAAJgaYQYAAJgaYQYAAJjadX/TvIKCAh09elSBgYFMnQQAwCQMw1BmZmaJfkvrug8zR48eVXh4uLvLAAAALjhy5MgV7wp93YeZwt8DOnLkiNNt7wHgSrKzs+0/c3D06FGn32gCcG1kZGQoPDy8yN/1u9R1H2YKh5aqVatGmAFQahf/GGO1atUIM0A5K8klIlwADAAATI0wAwAATI0wAwAATI0wAwAATI0wAwAATI0wAwAATI0wAwAATI0wAwAATI0wAwAATO26vwMwAAAoe/kFhpKSTyktM0fBgT5qGREkTw/3/KAzYQYAAJRKwp4Uxa/Zq5T0HHtbqNVHcT0j1SUqtNzrYZgJAACUWMKeFD25bLtDkJGk1PQcPblsuxL2pJR7TYQZAABQIvkFhuLX7JVRxLrCtvg1e5VfUNQW1w5hBgAAlEhS8imnHpmLGZJS0nOUlHyq/IoSYQYAAJRQWmbxQcaV7coKYQYAAJRIcKBPmW5XVggzAACgRFpGBCnU6qPiJmBbdGFWU8uIoPIsizADAABKxtPDoriekZLkFGgKl+N6Rpb7/WYIMwAAoMS6RIVq3qBo2ayOQ0k2q4/mDYp2y31muGkeAAAolS5RoeoUaeMOwAAAwLw8PSxqVb+Gu8uQxDATAAAwOcIMAAAwNcIMAAAwNcIMAAAwNcIMAAAwNcIMAAAwNcIMAAAwNcIMAAAwNcIMAAAwNcIMAAAwNcIMAAAwNcIMAAAwNcIMAAAwNcIMAAAwNcIMAAAwNcIMAAAwNcIMAAAwNcIMAAAwNcIMAAAwNcIMAAAwNcIMAAAwNcIMAAAwNcIMAAAwNcIMAAAwNcIMAAAwNcIMAAAwNcIMAAAwNcIMAAAwNcIMAAAwNcIMAAAwNcIMAAAwNcIMAAAwNcIMAAAwtSruLgAAAJS//AJDScmnlJaZo+BAH7WMCJKnh8XdZbmEMAMAQCWTsCdF8Wv2KiU9x94WavVRXM9IdYkKdWNlrmGYCQCASiRhT4qeXLbdIchIUmp6jp5ctl0Je1LcVJnrCDMAAFQS+QWG4tfslVHEusK2+DV7lV9Q1BYVF2EGAIBKIin5lFOPzMUMSSnpOUpKPlV+RZUBt4aZvLw8/f3vf1dERIR8fX1Vr149TZ48WQUFBfZtDMPQpEmTFBYWJl9fX7Vr104///yzG6sGAMCc0jKLDzKubFdRuDXMTJ8+XfPnz9ecOXO0b98+zZgxQ//85z81e/Zs+zYzZszQzJkzNWfOHG3btk02m02dOnVSZmamGysHAMB8ggN9ynS7isKtYWbLli2677771L17d9WtW1cPPvigOnfurB9++EHShV6ZWbNmaeLEierTp4+ioqK0ZMkSnTlzRu+++647SwcAwHRaRgQp1Oqj4iZgW3RhVlPLiKDyLOuquTXM3HXXXfrqq6/0yy+/SJJ27dqlb7/9Vt26dZMkJScnKzU1VZ07d7Y/x9vbW23bttXmzZuL3Gdubq4yMjIcHgAAQPL0sCiuZ6QkOQWawuW4npGmu9+MW8PMCy+8oIcffliNGjWSl5eXmjVrptjYWD388MOSpNTUVElSSEiIw/NCQkLs6y41bdo0Wa1W+yM8PPzavggAAEykS1So5g2Kls3qOJRks/po3qBoU95nxq03zVu5cqWWLVumd999V7feeqt27typ2NhYhYWFafDgwfbtLBbHhGgYhlNboQkTJmjs2LH25YyMDAINAAAX6RIVqk6RNu4AXBb++te/avz48XrooYckSbfddpsOHTqkadOmafDgwbLZbJIu9NCEhv5fUkxLS3PqrSnk7e0tb2/va188AAAm5ulhUav6NdxdRplw6zDTmTNn5OHhWIKnp6d9anZERIRsNpvWr19vX3/u3DklJiaqdevW5VorAAComNzaM9OzZ0+9/PLLql27tm699Vbt2LFDM2fO1GOPPSbpwvBSbGyspk6dqoYNG6phw4aaOnWq/Pz8NGDAAHeWDgAAKgi3hpnZs2frxRdf1KhRo5SWlqawsDCNGDFCL730kn2bcePG6ezZsxo1apROnz6tmJgYrVu3ToGBgW6sHAAAVBQWwzDM9QMMpZSRkSGr1ar09HRVq1bN3eUAMJns7GwFBARIkrKysuTv7+/mioDKoTTf3/w2EwAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMDXCDAAAMLUq7i4AAIDrRX6BoaTkU0rLzFFwoI9aRgTJ08Pi7rKue4QZAADKQMKeFMWv2auU9Bx7W6jVR3E9I9UlKtSNlV3/GGYCAOAqJexJ0ZPLtjsEGUlKTc/Rk8u2K2FPipsqqxwIMwAAXIX8AkPxa/bKKGJdYVv8mr3KLyhqC5QFwgwAAFchKfmUU4/MxQxJKek5Sko+VX5FVTKEGQAArkJaZvFBxpXtUHqEGQAArkJwoE+ZbofSI8wAAHAVWkYEKdTqo+ImYFt0YVZTy4ig8iyrUiHMAABwFTw9LIrrGSlJToGmcDmuZyT3m7mGCDMAAFylLlGhmjcoWjar41CSzeqjeYOiuc/MNcZN8wAAKANdokLVKdLGHYDdgDADAEAZ8fSwqFX9Gu4uo9JhmAkAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJgaYQYAAJia28PMH3/8oUGDBqlGjRry8/NT06ZN9eOPP9rXG4ahSZMmKSwsTL6+vmrXrp1+/vlnN1YMAAAqEreGmdOnT+vOO++Ul5eXvvjiC+3du1evvvqqbrjhBvs2M2bM0MyZMzVnzhxt27ZNNptNnTp1UmZmpvsKBwAAFUYVdx58+vTpCg8P16JFi+xtdevWtf/ZMAzNmjVLEydOVJ8+fSRJS5YsUUhIiN59912NGDGivEsGAAAVjFt7ZlavXq0WLVqob9++Cg4OVrNmzfTWW2/Z1ycnJys1NVWdO3e2t3l7e6tt27bavHlzkfvMzc1VRkaGwwMAAFy/3Bpmfv/9d82bN08NGzbUv//9b40cOVJjxozRO++8I0lKTU2VJIWEhDg8LyQkxL7uUtOmTZPVarU/wsPDr+2LAAAAbuXWMFNQUKDo6GhNnTpVzZo104gRI/T4449r3rx5DttZLBaHZcMwnNoKTZgwQenp6fbHkSNHrln9AADA/dwaZkJDQxUZGenQ1rhxYx0+fFiSZLPZJMmpFyYtLc2pt6aQt7e3qlWr5vAAAADXL7eGmTvvvFMHDhxwaPvll19Up04dSVJERIRsNpvWr19vX3/u3DklJiaqdevW5VorAAComNw6m+nZZ59V69atNXXqVPXr109JSUl688039eabb0q6MLwUGxurqVOnqmHDhmrYsKGmTp0qPz8/DRgwwJ2lAwCACsKtYeb222/XJ598ogkTJmjy5MmKiIjQrFmzNHDgQPs248aN09mzZzVq1CidPn1aMTExWrdunQIDA91YOQAAqCgshmEY7i7iWsrIyJDValV6ejrXzwAotezsbAUEBEiSsrKy5O/v7+aKgMqhNN/fbv85AwAAgKvh1mEmAAAqkvwCQ0nJp5SWmaPgQB+1jAiSp0fRtwJBxUGYAQBAUsKeFMWv2auU9Bx7W6jVR3E9I9UlKtSNleFKGGYCAFR6CXtS9OSy7Q5BRpJS03P05LLtStiT4qbKUBKEGQBApZZfYCh+zV4VNRumsC1+zV7lF1zX82VMjTADAKjUkpJPOfXIXMyQlJKeo6TkU+VXFEqFMAMAqNTSMosPMq5sh/JHmAEAVGrBgT5luh3KH2EGAFCptYwIUqjVR8VNwLbowqymlhFB5VkWSoEwAwCo1Dw9LIrrGSlJToGmcDmuZyT3m6nAyiTM5Ofna+fOnTp9+nRZ7A4AgHLVJSpU8wZFy2Z1HEqyWX00b1A095mp4Fy6aV5sbKxuu+02DRs2TPn5+Wrbtq02b94sPz8/rV27Vu3atSvjMgEAuLa6RIWqU6SNOwCbkEth5sMPP9SgQYMkSWvWrFFycrL279+vd955RxMnTtR3331XpkUCAFAePD0salW/hrvLQCm5NMx04sQJ2Ww2SdLnn3+uvn376uabb9awYcP0008/lWmBAAAAl+NSmAkJCdHevXuVn5+vhIQEdezYUZJ05swZeXp6lmmBAAAAl+PSMNPQoUPVr18/hYaGymKxqFOnTpKk77//Xo0aNSrTAgEAAC7HpTAzadIkRUVF6ciRI+rbt6+8vb0lSZ6enho/fnyZFggAAHA5FsMwrutfzsrIyJDValV6erqqVavm7nIAmEx2drYCAgIkSVlZWfL393dzRUDlUJrv7xL3zLz++uslLmDMmDEl3hYAAOBqlDjMvPbaayXazmKxEGYAAEC5KXGYSU5OvpZ1AAAAuITfZgIAAKbm0mwmSfrvf/+r1atX6/Dhwzp37pzDupkzZ151YQAAACXhUpj56quv1KtXL0VEROjAgQOKiorSwYMHZRiGoqOjy7pGAACAYrk0zDRhwgQ999xz2rNnj3x8fPTRRx/pyJEjatu2rfr27VvWNQIAABTLpTCzb98+DR48WJJUpUoVnT17VgEBAZo8ebKmT59epgUCAABcjkthxt/fX7m5uZKksLAw/fbbb/Z1J06cKJvKAAAASsCla2buuOMOfffdd4qMjFT37t313HPP6aefftLHH3+sO+64o6xrBAAAKJZLYWbmzJnKysqSdOF3mrKysrRy5Uo1aNCgxDfXAwAAKAsuhZl69erZ/+zn56e5c+eWWUEAAAClwU3zAACAqbnUM+Ph4SGLxVLs+vz8fJcLAgAAKA2Xwswnn3zisHz+/Hnt2LFDS5YsUXx8fJkUBgAAUBIuhZn77rvPqe3BBx/UrbfeqpUrV2rYsGFXXRgAAEBJlOk1MzExMfryyy/LcpcAAACXVWZh5uzZs5o9e7Zq1apVVrsEAAC4IpeGmapXr+5wAbBhGMrMzJSfn5+WLVtWZsUBAABciUth5rXXXnMIMx4eHqpZs6ZiYmJUvXr1MisOAADgSlwKM0OGDCnjMgAAAFxT4jCze/fuEu+0SZMmLhUDAABQWiUOM02bNpXFYpFhGJLETfMAAECFUOLZTMnJyfr999+VnJysjz/+WBEREZo7d6527NihHTt2aO7cuapfv74++uija1kvAACAgxL3zNSpU8f+5759++r1119Xt27d7G1NmjRReHi4XnzxRfXu3btMiwQAACiOS/eZ+emnnxQREeHUHhERob179151UQAAACXlUphp3LixpkyZopycHHtbbm6upkyZosaNG5dZcQAAAFfi0tTs+fPnq2fPngoPD9df/vIXSdKuXbtksVi0du3aMi0QAADgclwKMy1btlRycrKWLVum/fv3yzAM9e/fXwMGDJC/v39Z1wignOUXGEpKPqW0zBwFB/qoZUSQPD2Kn8EIAO7kUpiRJD8/Pz3xxBNlWQuACiBhT4ri1+xVSvr/DSOHWn0U1zNSXaJC3VgZABStxGFm9erV6tq1q7y8vLR69erLbturV6+rLgxA+UvYk6Inl22XcUl7anqOnly2XfMGRRNoAFQ4FqPwLnhX4OHhodTUVAUHB8vDo/jrhi0WS4W6aV5GRoasVqvS09NVrVo1d5cDVFj5BYbumr7BoUfmYhZJNquPvn2hQ6UacsrOzlZAQIAkKSsri6F0oJyU5vu7xLOZCgoKFBwcbP9zcY+KFGQAlFxS8qlig4wkGZJS0nOUlHyq/IoCgBJwaWp2Uf7888+y2hUAN0jLLD7IuLIdAJQXl8LM9OnTtXLlSvty3759FRQUpJtuukm7du0qs+IAlJ/gQJ8y3Q4AyotLYWbBggUKDw+XJK1fv15ffvmlEhIS1LVrV/31r38t0wIBlI+WEUEKtfqouKthLLowq6llRFB5lgUAV+TS1OyUlBR7mFm7dq369eunzp07q27duoqJiSnTAgGUD08Pi+J6RurJZdtlkRxmNBUGnLiekZXq4l8A5uBSz0z16tV15MgRSVJCQoI6duwoSTIMgwuAARPrEhWqeYOiZbM6DiXZrD5MywZQYbnUM9OnTx8NGDBADRs21MmTJ9W1a1dJ0s6dO9WgQYMyLRBA+eoSFapOkTbuAAzANFwKM6+99prq1q2rI0eOaMaMGfZ7MKSkpGjUqFFlWiCA8ufpYVGr+jXcXQYAlEiJb5pnVtw0D8DV4KZ5gHtck5vmXWrp0qW66667FBYWpkOHDkmSZs2apU8//dTVXQIAAJSaS2Fm3rx5Gjt2rLp27ao///zTftHvDTfcoFmzZrlUyLRp02SxWBQbG2tvMwxDkyZNUlhYmHx9fdWuXTv9/PPPLu0fAABcn1wKM7Nnz9Zbb72liRMnytPT097eokUL/fTTT6Xe37Zt2/Tmm2+qSZMmDu0zZszQzJkzNWfOHG3btk02m02dOnVSZmamK2UDAIDrkEthJjk5Wc2aNXNq9/b2VnZ2dqn2lZWVpYEDB+qtt95S9erV7e2GYWjWrFmaOHGi+vTpo6ioKC1ZskRnzpzRu+++60rZAADgOuRSmImIiNDOnTud2r/44gs1bty4VPsaPXq0unfvbr9XTaHk5GSlpqaqc+fO9jZvb2+1bdtWmzdvdqVsAABwHXJpavZf//pXjR49Wjk5OTIMQ0lJSXrvvfc0depUvf322yXez4oVK/Tjjz/qhx9+cFqXmpoqSQoJCXFoDwkJsV9wXJTc3Fzl5ubalzMyMkpcDwAAMB+XwszQoUOVl5encePG6cyZMxowYIBuuukmzZ49W23atCnRPo4cOaJnnnlG69atk49P8T9cZ7E43qjLMAyntotNmzZN8fHxJXshAADA9Fyemv3444/r0KFDSktLU2pqqpKSkrRjx44S3wH4xx9/VFpampo3b64qVaqoSpUqSkxM1Ouvv64qVarYe2QKe2gKpaWlOfXWXGzChAlKT0+3Pwp/dgEAAFyfShVm/vzzTw0cOFA1a9ZUWFiYXn/9dQUFBemNN95QgwYNtHXrVi1cuLBE+7rnnnv0008/aefOnfZHixYtNHDgQO3cuVP16tWTzWbT+vXr7c85d+6cEhMT1bp162L36+3trWrVqjk8AADA9atUw0x/+9vf9M0332jw4MFKSEjQs88+q4SEBOXk5Ojzzz9X27ZtS7yvwMBARUVFObT5+/urRo0a9vbY2FhNnTpVDRs2VMOGDTV16lT5+flpwIABpSkbAABcx0oVZj777DMtWrRIHTt21KhRo9SgQQPdfPPNLt8o70rGjRuns2fPatSoUTp9+rRiYmK0bt06BQYGXpPjAQAA8ynVbzN5eXnp0KFDCgsLkyT5+fkpKSnJqYelIuG3mQBcDX6bCXCPa/bbTAUFBfLy8rIve3p68g8bAAC4VamGmQzD0JAhQ+Tt7S1JysnJ0ciRI50Czccff1x2FQIAAFxGqcLM4MGDHZYHDRpUpsUAAACUVqnCzKJFi65VHQAAAC5x+aZ5AAAAFQFhBgAAmBphBgAAmBphBgAAmBphBgAAmBphBgAAmBphBgAAmBphBgAAmBphBgAAmBphBgAAmBphBgAAmBphBgAAmBphBgAAmBphBgAAmBphBgAAmBphBgAAmFoVdxcAVFT5BYaSkk8pLTNHwYE+ahkRJE8Pi7vLAgBcgjADFCFhT4ri1+xVSnqOvS3U6qO4npHqEhXqxsoAAJdimAm4RMKeFD25bLtDkJGk1PQcPblsuxL2pLipMgBAUQgzwEXyCwzFr9kro4h1hW3xa/Yqv6CoLQAA7kCYAS6SlHzKqUfmYoaklPQcJSWfKr+iAACXRZgBLpKWWXyQcWU7AMC1R5gBLhIc6FOm2wEArj3CDHCRlhFBCrX6qLgJ2BZdmNXUMiKoPMsCAFwGYQa4iKeHRXE9IyXJKdAULsf1jOR+MwBQgRBmgEt0iQrVvEHRslkdh5JsVh/NGxTNfWYAoILhpnlAEbpEhapTpI07AAOACRBmgGJ4eljUqn4Nd5cBALgChpkAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpuTXMTJs2TbfffrsCAwMVHBys3r1768CBAw7bGIahSZMmKSwsTL6+vmrXrp1+/vlnN1UMAAAqGreGmcTERI0ePVpbt27V+vXrlZeXp86dOys7O9u+zYwZMzRz5kzNmTNH27Ztk81mU6dOnZSZmenGygEAQEVhMQzDcHcRhY4fP67g4GAlJibq7rvvlmEYCgsLU2xsrF544QVJUm5urkJCQjR9+nSNGDHiivvMyMiQ1WpVenq6qlWrdq1fAoDrTHZ2tgICAiRJWVlZ8vf3d3NFQOVQmu/vCnXNTHp6uiQpKChIkpScnKzU1FR17tzZvo23t7fatm2rzZs3F7mP3NxcZWRkODwAAMD1q8KEGcMwNHbsWN11112KioqSJKWmpkqSQkJCHLYNCQmxr7vUtGnTZLVa7Y/w8PBrW/h1KL/A0JbfTurTnX9oy28nlV9QYTrvAABwUsXdBRR66qmntHv3bn377bdO6ywWi8OyYRhObYUmTJigsWPH2pczMjIINKWQsCdF8Wv2KiU9x94WavVRXM9IdYkKdWNlAAAUrUL0zDz99NNavXq1vv76a9WqVcvebrPZJMmpFyYtLc2pt6aQt7e3qlWr5vBAySTsSdGTy7Y7BBlJSk3P0ZPLtithT4qbKgMAoHhuDTOGYeipp57Sxx9/rA0bNigiIsJhfUREhGw2m9avX29vO3funBITE9W6devyLve6ll9gKH7NXhU1oFTYFr9mL0NOAIAKx63DTKNHj9a7776rTz/9VIGBgfYeGKvVKl9fX1ksFsXGxmrq1Klq2LChGjZsqKlTp8rPz08DBgxwZ+nXnaTkU049MhczJKWk5ygp+ZRa1a9RfoUBAHAFbg0z8+bNkyS1a9fOoX3RokUaMmSIJGncuHE6e/asRo0apdOnTysmJkbr1q1TYGBgOVd7fUvLLD7IuLIdAADlxa1hpiS3uLFYLJo0aZImTZp07QuqxIIDfcp0OwAAykuFuAAY7tcyIkihVh8VPUdMsujCrKaWEUHlWRYAAFdEmIEkydPDoriekZLkFGgKl+N6RsrTo7i4AwCAexBmYNclKlTzBkXLZnUcSrJZfTRvUDT3mQEAVEgV5qZ5qBi6RIWqU6RNScmnlJaZo+DAC0NL9MgAACoqwgyceHpYmH4NADANhpkAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpVXF3AWaVX2AoKfmU0jJzFBzoo5YRQfL0sLi7LAAAKh3CjAsS9qQofs1epaTn2NtCrT6K6xmpLlGhbqwMAIDKh2GmUkrYk6Inl213CDKSlJqeoyeXbVfCnhQ3VQYAQOVEmCmF/AJD8Wv2yihiXWFb/Jq9yi8oagsAAHAtEGZKISn5lFOPzMUMSSnpOUpKPlV+RQEAUMkRZkohLbP4IOPKdgAA4OoRZkohONCnTLcDAABXjzBTCi0jghRq9VFxE7AtujCrqWVEUHmWBQBApUaYKQVPD4viekZKklOgKVyO6xnJ/WYAAChHhJlS6hIVqnmDomWzOg4l2aw+mjcomvvMAABQzrhpngu6RIWqU6SNOwADAFABEGZc5OlhUav6NdxdBgAAlR7DTAAAwNQIMwAAwNQIMwAAwNQIMwAAwNQIMwAAwNQIMwAAwNQIMwAAwNQIMwAAwNQIMwAAwNQIMwAAwNQIMwAAwNRMEWbmzp2riIgI+fj4qHnz5tq0aZO7SwIAABVEhQ8zK1euVGxsrCZOnKgdO3aoTZs26tq1qw4fPuzu0gAAQAVgMQzDcHcRlxMTE6Po6GjNmzfP3ta4cWP17t1b06ZNu+LzMzIyZLVadfToUVWrVu1algrgOpSdna2QkBBJ0rFjx+Tv7+/mioDKISMjQ2FhYUpPT7/i93eVcqrJJefOndOPP/6o8ePHO7R37txZmzdvLvI5ubm5ys3NtS9nZGRIksLCwq5doQAqhcJQA6BiqdDDTCdOnFB+fr7TB0hISIhSU1OLfM60adNktVrtj/Dw8PIoFQAAuEmF7pkpZLFYHJYNw3BqKzRhwgSNHTvWvpyRkaHw8HCGmQC4hGEmwD0Kh5lKokKHmRtvvFGenp5OvTBpaWnFdvd6e3vL29vbqd3f358PIQBXhc8RoPzk5+eXeNsKPcxUtWpVNW/eXOvXr3doX79+vVq3bu2mqgAAQEVSoXtmJGns2LF65JFH1KJFC7Vq1UpvvvmmDh8+rJEjR7q7NAAAUAFU+DDTv39/nTx5UpMnT1ZKSoqioqL0+eefq06dOu4uDQAAVAAV/j4zV6vwPjMlmacOAJfKzs5WQECAJCkrK4trZoByUprv7wp9zQwAAMCVEGYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpEWYAAICpVfifM7hahTc4zsjIcHMlAMwoOzvb/ueMjIxS/ZIvANcVfm+X5IcKrvswk5mZKUkKDw93cyUAzC4sLMzdJQCVTmZmpqxW62W3ue5/m6mgoEBHjx5VYGCgLBaLu8sxlYyMDIWHh+vIkSP8rpWbcA4qBs5DxcB5qBjK6zwYhqHMzEyFhYXJw+PyV8Vc9z0zHh4eqlWrlrvLMLVq1arxweFmnIOKgfNQMXAeKobyOA9X6pEpxAXAAADA1AgzAADA1AgzKJa3t7fi4uLk7e3t7lIqLc5BxcB5qBg4DxVDRTwP1/0FwAAA4PpGzwwAADA1wgwAADA1wgwAADA1wgwAADA1wgxKJTc3V02bNpXFYtHOnTvdXU6lcvDgQQ0bNkwRERHy9fVV/fr1FRcXp3Pnzrm7tOve3LlzFRERIR8fHzVv3lybNm1yd0mVyrRp03T77bcrMDBQwcHB6t27tw4cOODusiq1adOmyWKxKDY21t2lSCLMoJTGjRvH79O4yf79+1VQUKAFCxbo559/1muvvab58+frb3/7m7tLu66tXLlSsbGxmjhxonbs2KE2bdqoa9euOnz4sLtLqzQSExM1evRobd26VevXr1deXp46d+7s8COgKD/btm3Tm2++qSZNmri7FDumZqPEvvjiC40dO1YfffSRbr31Vu3YsUNNmzZ1d1mV2j//+U/NmzdPv//+u7tLuW7FxMQoOjpa8+bNs7c1btxYvXv31rRp09xYWeV1/PhxBQcHKzExUXfffbe7y6lUsrKyFB0drblz52rKlClq2rSpZs2a5e6y6JlByRw7dkyPP/64li5dKj8/P3eXg/8vPT1dQUFB7i7junXu3Dn9+OOP6ty5s0N7586dtXnzZjdVhfT0dEni774bjB49Wt27d1fHjh3dXYqD6/6HJnH1DMPQkCFDNHLkSLVo0UIHDx50d0mQ9Ntvv2n27Nl69dVX3V3KdevEiRPKz89XSEiIQ3tISIhSU1PdVFXlZhiGxo4dq7vuuktRUVHuLqdSWbFihX788Uf98MMP7i7FCT0zldikSZNksVgu+/jhhx80e/ZsZWRkaMKECe4u+bpU0vNwsaNHj6pLly7q27evhg8f7qbKKw+LxeKwbBiGUxvKx1NPPaXdu3frvffec3cplcqRI0f0zDPPaPny5fLx8XF3OU64ZqYSO3HihE6cOHHZberWrauHHnpIa9ascfjwzs/Pl6enpwYOHKglS5Zc61KvayU9D4UfIEePHlX79u0VExOjxYsXy8OD/ye5Vs6dOyc/Pz998MEHuv/+++3tzzzzjHbu3KnExEQ3Vlf5PP3001q1apW++eYbRUREuLucSmXVqlW6//775enpaW/Lz8+XxWKRh4eHcnNzHdaVN8IMrujw4cPKyMiwLx89elT33nuvPvzwQ8XExKhWrVpurK5y+eOPP9S+fXs1b95cy5Ytc+uHR2URExOj5s2ba+7cufa2yMhI3XfffVwAXE4Mw9DTTz+tTz75RBs3blTDhg3dXVKlk5mZqUOHDjm0DR06VI0aNdILL7zg9iE/rpnBFdWuXdthOSAgQJJUv359gkw5Onr0qNq1a6fatWvrlVde0fHjx+3rbDabGyu7vo0dO1aPPPKIWrRooVatWunNN9/U4cOHNXLkSHeXVmmMHj1a7777rj799FMFBgbar1eyWq3y9fV1c3WVQ2BgoFNg8ff3V40aNdweZCTCDGAa69at06+//qpff/3VKUTSwXrt9O/fXydPntTkyZOVkpKiqKgoff7556pTp467S6s0CqfFt2vXzqF90aJFGjJkSPkXhAqHYSYAAGBqXDkIAABMjTADAABMjTADAABMjTADAABMjTADAABMjTADAABMjTADAABMjTADwO0mTZqkpk2bXvV+LBaLVq1aVez6gwcPymKxaOfOnZKkjRs3ymKx6M8//5QkLV68WDfccMNV1wGgfBFmAJTKkCFD7L/m7eXlpXr16un5559Xdna2u0u7ovDwcPtdfIvSv39//fLLL/blsgpZAK4tfs4AQKl16dJFixYt0vnz57Vp0yYNHz5c2dnZ9tvOFzp//ry8vLzcVKUzT0/Py/6Ola+vL7/1A5gQPTMASs3b21s2m03h4eEaMGCABg4cqFWrVtl7MhYuXKh69erJ29tbhmHo8OHDuu+++xQQEKBq1aqpX79+OnbsmNN+FyxYoPDwcPn5+alv37724R9J2rZtmzp16qQbb7xRVqtVbdu21fbt2532kZKSoq5du8rX11cRERH64IMP7OsuHWa61MXDTIsXL1Z8fLx27dpl74lavHixHnvsMfXo0cPheXl5ebLZbFq4cGHp30wAV40wA+Cq+fr66vz585KkX3/9Ve+//74++ugje2jo3bu3Tp06pcTERK1fv16//fab+vfv77CPwuetWbNGCQkJ2rlzp0aPHm1fn5mZqcGDB2vTpk3aunWrGjZsqG7duikzM9NhPy+++KIeeOAB7dq1S4MGDdLDDz+sffv2lfo19e/fX88995xuvfVWpaSkKCUlRf3799fw4cOVkJCglJQU+7aff/65srKy1K9fv1IfB8DVY5gJwFVJSkrSu+++q3vuuUeSdO7cOS1dulQ1a9aUJK1fv167d+9WcnKywsPDJUlLly7Vrbfeqm3btun222+XJOXk5GjJkiX2XwSfPXu2unfvrldffVU2m00dOnRwOO6CBQtUvXp1JSYmOvSU9O3bV8OHD5ck/eMf/9D69es1e/ZszZ07t1Svy9fXVwEBAapSpYrD0FTr1q11yy23aOnSpRo3bpykC7/e3LdvXwUEBJTqGADKBj0zAEpt7dq1CggIkI+Pj1q1aqW7775bs2fPliTVqVPHHmQkad++fQoPD7cHGUmKjIzUDTfc4NBjUrt2bXuQkaRWrVqpoKBABw4ckCSlpaVp5MiRuvnmm2W1WmW1WpWVlaXDhw871NaqVSunZVd6Zi5n+PDhWrRokb2uzz77TI899liZHgNAydEzA6DU2rdvr3nz5snLy0thYWEOF/n6+/s7bGsYhiwWi9M+imsvVLiu8L9DhgzR8ePHNWvWLNWpU0fe3t5q1aqVzp07d8V6L3ccVzz66KMaP368tmzZoi1btqhu3bpq06ZNmR4DQMnRMwOg1Pz9/dWgQQPVqVPnirOVIiMjdfjwYR05csTetnfvXqWnp6tx48b2tsOHD+vo0aP25S1btsjDw0M333yzJGnTpk0aM2aMunXrpltvvVXe3t46ceKE0/G2bt3qtNyoUSOXXmfVqlWVn5/v1F6jRg317t1bixYt0qJFizR06FCX9g+gbNAzA+Ca6tixo5o0aaKBAwdq1qxZysvL06hRo9S2bVu1aNHCvp2Pj48GDx6sV155RRkZGRozZoz69etnv16lQYMGWrp0qVq0aKGMjAz99a9/LXIa9QcffKAWLVrorrvu0vLly5WUlKS3337bpdrr1q2r5ORk7dy5U7Vq1VJgYKC8vb0lXRhq6tGjh/Lz8zV48GCX9g+gbNAzA+CaKrwrb/Xq1XX33XerY8eOqlevnlauXOmwXYMGDdSnTx9169ZNnTt3VlRUlMNFuwsXLtTp06fVrFkzPfLIIxozZoyCg4OdjhcfH68VK1aoSZMmWrJkiZYvX67IyEiXan/ggQfUpUsXtW/fXjVr1tR7771nX9exY0eFhobq3nvvVVhYmEv7B1A2LIZhGO4uAgDM5syZMwoLC9PChQvVp08fd5cDVGoMMwFAKRQUFCg1NVWvvvqqrFarevXq5e6SgEqPMAMApXD48GFFRESoVq1aWrx4sapU4WMUcDeGmQAAgKlxATAAADA1wgwAADA1wgwAADA1wgwAADA1wgwAADA1wgwAADA1wgwAADA1wgwAADA1wgwAADC1/wd6vCAT8qiWGwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "prob_plot(residuals,'Residuals',title=\"(all main effects)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d797f54f",
   "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}