Balance Optimization Subset Selection.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Input initial parameter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 163,
   "metadata": {},
   "outputs": [],
   "source": [
    "samples=40\n",
    "no_bins=4 #bins = (1,2,,,n_bins-1)\n",
    "n_cov=4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create random matrix "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 164,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "mat=np.random.random((samples,n_cov))*10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 165,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_bins_array(initial_mat,no_bins,n_cov):\n",
    "    bmat=np.zeros((len(initial_mat),len(initial_mat[0])))\n",
    "    bins=np.linspace(0,10,no_bins)\n",
    "    for b in range(n_cov):\n",
    "        data=initial_mat[:,b]\n",
    "        digitized=np.digitize(data,bins)\n",
    "        bmat[:,b]=digitized\n",
    "    return bmat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create bin matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 166,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[3. 3. 2. 2.]\n",
      " [3. 1. 2. 1.]\n",
      " [2. 3. 1. 1.]\n",
      " [1. 3. 3. 3.]\n",
      " [1. 1. 3. 2.]\n",
      " [2. 2. 1. 2.]\n",
      " [1. 1. 1. 2.]\n",
      " [3. 2. 2. 3.]\n",
      " [3. 1. 2. 2.]\n",
      " [1. 2. 1. 3.]\n",
      " [3. 1. 3. 1.]\n",
      " [1. 3. 2. 3.]\n",
      " [3. 3. 1. 1.]\n",
      " [2. 2. 1. 1.]\n",
      " [2. 3. 3. 3.]\n",
      " [3. 3. 2. 1.]\n",
      " [1. 3. 3. 3.]\n",
      " [1. 3. 2. 3.]\n",
      " [1. 2. 2. 2.]\n",
      " [1. 3. 1. 1.]\n",
      " [2. 3. 3. 2.]\n",
      " [2. 2. 1. 2.]\n",
      " [2. 2. 1. 1.]\n",
      " [2. 3. 2. 3.]\n",
      " [1. 2. 2. 1.]\n",
      " [2. 2. 1. 2.]\n",
      " [3. 1. 2. 2.]\n",
      " [1. 2. 3. 2.]\n",
      " [1. 3. 1. 3.]\n",
      " [2. 3. 1. 1.]\n",
      " [2. 3. 2. 2.]\n",
      " [2. 2. 2. 3.]\n",
      " [3. 1. 1. 1.]\n",
      " [1. 1. 3. 1.]\n",
      " [2. 2. 3. 1.]\n",
      " [1. 2. 1. 1.]\n",
      " [3. 3. 2. 1.]\n",
      " [2. 1. 3. 2.]\n",
      " [2. 2. 2. 1.]\n",
      " [1. 3. 1. 1.]]\n"
     ]
    }
   ],
   "source": [
    "bins=create_bins_array(mat,no_bins,n_cov)\n",
    "print(bins)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model Formulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 167,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyqubo import Array,Constraint\n",
    "s=Array.create('s',shape=samples,vartype='SPIN')\n",
    "\n",
    "def objective_function(s,bins,no_bins):\n",
    "    H=0\n",
    "    for col in range(len(mat[0])):\n",
    "        for b in range(1,no_bins+1):\n",
    "            M=0\n",
    "            for val,ss in zip(bins[:,col],s):\n",
    "                if val==b:\n",
    "                    M=M+ss\n",
    "            H=H+(M-0)**2\n",
    "    return H-0\n",
    "\n",
    "H=(objective_function(s,bins,no_bins))\n",
    "#M=4\n",
    "#con=sum(ss for ss in s)**2\n",
    "#H=H-M*Constraint(con,label='q-1')\n",
    "model=H.compile()\n",
    "qubo,offset=model.to_qubo()\n",
    "bqm=model.to_dimod_bqm()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### QUBO matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 168,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "       0      1      2      3      4      5      6      7      8      9   \\\n",
      "0  -204.0   16.0    8.0    8.0    8.0    8.0    8.0   16.0   24.0    0.0   \n",
      "1     0.0 -188.0    8.0    0.0    8.0    0.0    8.0   16.0   24.0    0.0   \n",
      "2     0.0    0.0 -240.0    8.0    0.0   16.0    8.0    0.0    0.0    8.0   \n",
      "3     0.0    0.0    0.0 -192.0   16.0    0.0    8.0    8.0    0.0   16.0   \n",
      "4     0.0    0.0    0.0    0.0 -172.0    8.0   24.0    0.0   16.0    8.0   \n",
      "5     0.0    0.0    0.0    0.0    0.0 -212.0   16.0    8.0    8.0   16.0   \n",
      "6     0.0    0.0    0.0    0.0    0.0    0.0 -192.0    0.0   16.0   16.0   \n",
      "7     0.0    0.0    0.0    0.0    0.0    0.0    0.0 -180.0   16.0   16.0   \n",
      "8     0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0 -172.0    0.0   \n",
      "9     0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0 -200.0   \n",
      "10    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0   \n",
      "11    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0   \n",
      "12    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0   \n",
      "13    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0   \n",
      "14    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0   \n",
      "15    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0   \n",
      "16    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0   \n",
      "17    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0   \n",
      "18    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0   \n",
      "19    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0   \n",
      "20    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0   \n",
      "21    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0   \n",
      "22    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0   \n",
      "23    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0   \n",
      "24    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0   \n",
      "25    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0   \n",
      "26    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0   \n",
      "27    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0   \n",
      "28    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0   \n",
      "29    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0   \n",
      "30    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0   \n",
      "31    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0   \n",
      "32    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0   \n",
      "33    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0   \n",
      "34    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0   \n",
      "35    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0   \n",
      "36    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0   \n",
      "37    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0   \n",
      "38    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0   \n",
      "39    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0   \n",
      "\n",
      "    ...       30     31     32     33     34     35     36     37     38  \\\n",
      "0   ...     24.0    8.0    8.0    0.0    0.0    0.0   24.0    8.0    8.0   \n",
      "1   ...      8.0    8.0   24.0   16.0    8.0    8.0   24.0    8.0   16.0   \n",
      "2   ...     16.0    8.0   16.0    8.0   16.0   16.0   16.0    8.0   16.0   \n",
      "3   ...      8.0    8.0    0.0   16.0    8.0    8.0    8.0    8.0    0.0   \n",
      "4   ...      8.0    0.0    8.0   24.0    8.0    8.0    0.0   24.0    0.0   \n",
      "5   ...     16.0   16.0    8.0    0.0   16.0   16.0    0.0   16.0   16.0   \n",
      "6   ...      8.0    0.0   16.0   16.0    0.0   16.0    0.0   16.0    0.0   \n",
      "7   ...      8.0   24.0    8.0    0.0    8.0    8.0   16.0    0.0   16.0   \n",
      "8   ...     16.0    8.0   16.0    8.0    0.0    0.0   16.0   16.0    8.0   \n",
      "9   ...      0.0   16.0    8.0    8.0    8.0   24.0    0.0    0.0    8.0   \n",
      "10  ...      0.0    0.0   24.0   24.0   16.0    8.0   16.0   16.0    8.0   \n",
      "11  ...     16.0   16.0    0.0    8.0    0.0    8.0   16.0    0.0    8.0   \n",
      "12  ...      8.0    0.0   24.0    8.0    8.0   16.0   24.0    0.0    8.0   \n",
      "13  ...      8.0   16.0   16.0    8.0   24.0   24.0    8.0    8.0   24.0   \n",
      "14  ...     16.0   16.0    0.0    8.0   16.0    0.0    8.0   16.0    8.0   \n",
      "15  ...     16.0    8.0   16.0    8.0    8.0    8.0   32.0    0.0   16.0   \n",
      "16  ...      8.0    8.0    0.0   16.0    8.0    8.0    8.0    8.0    0.0   \n",
      "17  ...     16.0   16.0    0.0    8.0    0.0    8.0   16.0    0.0    8.0   \n",
      "18  ...     16.0   16.0    0.0    8.0    8.0   16.0    8.0    8.0   16.0   \n",
      "19  ...      8.0    0.0   16.0   16.0    8.0   24.0   16.0    0.0    8.0   \n",
      "20  ...     24.0    8.0    0.0    8.0   16.0    0.0    8.0   24.0    8.0   \n",
      "21  ...     16.0   16.0    8.0    0.0   16.0   16.0    0.0   16.0   16.0   \n",
      "22  ...      8.0   16.0   16.0    8.0   24.0   24.0    8.0    8.0   24.0   \n",
      "23  ...     24.0   24.0    0.0    0.0    8.0    0.0   16.0    8.0   16.0   \n",
      "24  ...      8.0   16.0    8.0   16.0   16.0   24.0   16.0    0.0   24.0   \n",
      "25  ...     16.0   16.0    8.0    0.0   16.0   16.0    0.0   16.0   16.0   \n",
      "26  ...     16.0    8.0   16.0    8.0    0.0    0.0   16.0   16.0    8.0   \n",
      "27  ...      8.0    8.0    0.0   16.0   16.0   16.0    0.0   16.0    8.0   \n",
      "28  ...      8.0    8.0    8.0    8.0    0.0   16.0    8.0    0.0    0.0   \n",
      "29  ...     16.0    8.0   16.0    8.0   16.0   16.0   16.0    8.0   16.0   \n",
      "30  ...   -224.0   16.0    0.0    0.0    8.0    0.0   16.0   16.0   16.0   \n",
      "31  ...      0.0 -200.0    0.0    0.0   16.0    8.0    8.0    8.0   24.0   \n",
      "32  ...      0.0    0.0 -188.0   16.0    8.0   16.0   16.0    8.0    8.0   \n",
      "33  ...      0.0    0.0    0.0 -188.0   16.0   16.0    8.0   16.0    8.0   \n",
      "34  ...      0.0    0.0    0.0    0.0 -208.0   16.0    8.0   16.0   24.0   \n",
      "35  ...      0.0    0.0    0.0    0.0    0.0 -228.0    8.0    0.0   16.0   \n",
      "36  ...      0.0    0.0    0.0    0.0    0.0    0.0 -220.0    0.0   16.0   \n",
      "37  ...      0.0    0.0    0.0    0.0    0.0    0.0    0.0 -172.0    8.0   \n",
      "38  ...      0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0 -228.0   \n",
      "39  ...      0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0   \n",
      "\n",
      "       39  \n",
      "0     8.0  \n",
      "1     8.0  \n",
      "2    24.0  \n",
      "3    16.0  \n",
      "4     8.0  \n",
      "5     8.0  \n",
      "6    16.0  \n",
      "7     0.0  \n",
      "8     0.0  \n",
      "9    16.0  \n",
      "10    8.0  \n",
      "11   16.0  \n",
      "12   24.0  \n",
      "13   16.0  \n",
      "14    8.0  \n",
      "15   16.0  \n",
      "16   16.0  \n",
      "17   16.0  \n",
      "18    8.0  \n",
      "19   32.0  \n",
      "20    8.0  \n",
      "21    8.0  \n",
      "22   16.0  \n",
      "23    8.0  \n",
      "24   16.0  \n",
      "25    8.0  \n",
      "26    0.0  \n",
      "27    8.0  \n",
      "28   24.0  \n",
      "29   24.0  \n",
      "30    8.0  \n",
      "31    0.0  \n",
      "32   16.0  \n",
      "33   16.0  \n",
      "34    8.0  \n",
      "35   24.0  \n",
      "36   16.0  \n",
      "37    0.0  \n",
      "38    8.0  \n",
      "39 -240.0  \n",
      "\n",
      "[40 rows x 40 columns]\n"
     ]
    }
   ],
   "source": [
    "import re\n",
    "qmat=np.zeros((samples,samples))\n",
    "for key,val in qubo.items():\n",
    "    x=int(int(re.search(r'\\d+', key[0]).group()))\n",
    "    y=int(int(re.search(r'\\d+', key[1]).group()))\n",
    "    qmat[min(x,y)][max(x,y)]=val\n",
    "\n",
    "import pandas as pd\n",
    "q=pd.DataFrame(qmat)\n",
    "print(q)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Optimization "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 169,
   "metadata": {},
   "outputs": [],
   "source": [
    "from dwave.system.samplers import DWaveSampler\n",
    "from dwave.system.composites import EmbeddingComposite\n",
    "import neal\n",
    "sampler=EmbeddingComposite(DWaveSampler())\n",
    "#sampler=neal.SimulatedAnnealingSampler()\n",
    "response=sampler.sample(bqm,num_reads=100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Select one of the result that gives minimum energy as the answer to this problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 170,
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import defaultdict    \n",
    "result=defaultdict(list)\n",
    "for r in response.data(['sample', 'energy']):\n",
    "    result[r.energy].append(r.sample)\n",
    "\n",
    "res=pd.DataFrame()\n",
    "min_energy=min(result.keys())\n",
    "for i in list(result[min_energy]):\n",
    "    temp=pd.DataFrame.from_dict(i,'index').transpose()\n",
    "    res=res.append(temp)\n",
    "res=res.drop_duplicates()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 171,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    s[0]  s[10]  s[11]  s[12]  s[13]  s[14]  s[15]  s[16]  s[17]  s[18]  ...   \\\n",
      "0      1      0      0      0      0      0      0      1      1      0  ...    \n",
      "0      1      0      0      1      1      1      0      0      1      0  ...    \n",
      "0      1      1      1      0      1      0      0      0      0      0  ...    \n",
      "0      0      1      1      1      1      1      1      0      0      1  ...    \n",
      "0      1      1      0      0      0      1      0      1      1      1  ...    \n",
      "0      0      0      0      1      0      0      0      1      1      0  ...    \n",
      "0      0      1      0      1      0      0      1      0      0      0  ...    \n",
      "0      0      0      0      0      0      0      0      0      1      0  ...    \n",
      "0      1      0      0      1      1      0      0      1      0      0  ...    \n",
      "0      0      0      1      1      0      1      0      0      1      0  ...    \n",
      "0      0      0      1      0      1      0      1      1      1      0  ...    \n",
      "0      0      1      0      1      1      0      0      1      1      0  ...    \n",
      "0      1      1      1      0      1      0      1      0      1      0  ...    \n",
      "0      0      1      0      0      0      1      1      1      1      1  ...    \n",
      "0      1      1      0      1      1      0      0      0      1      0  ...    \n",
      "0      0      0      0      1      1      1      1      1      1      1  ...    \n",
      "0      1      0      1      1      0      0      0      0      1      0  ...    \n",
      "0      0      0      0      1      1      1      0      1      0      1  ...    \n",
      "0      0      0      1      0      0      1      1      0      0      0  ...    \n",
      "0      0      1      0      1      1      1      1      0      1      1  ...    \n",
      "0      1      0      1      1      1      1      1      0      1      1  ...    \n",
      "0      0      0      0      1      1      1      0      1      1      1  ...    \n",
      "0      1      1      0      0      1      0      1      1      1      1  ...    \n",
      "0      0      0      1      1      0      1      1      1      0      0  ...    \n",
      "0      0      1      0      0      0      1      1      1      0      1  ...    \n",
      "0      1      1      0      1      0      0      0      1      0      0  ...    \n",
      "0      0      1      1      1      1      0      0      1      0      1  ...    \n",
      "0      1      1      0      0      0      0      0      0      1      0  ...    \n",
      "0      1      1      1      0      1      0      1      1      0      1  ...    \n",
      "0      1      0      0      1      0      1      1      0      0      0  ...    \n",
      "..   ...    ...    ...    ...    ...    ...    ...    ...    ...    ...  ...    \n",
      "0      1      0      1      1      1      1      0      1      0      1  ...    \n",
      "0      0      1      0      1      1      0      0      0      1      1  ...    \n",
      "0      0      0      0      1      0      0      1      1      1      1  ...    \n",
      "0      0      0      0      0      1      0      1      0      1      0  ...    \n",
      "0      0      0      0      1      0      1      0      0      0      1  ...    \n",
      "0      0      1      0      1      1      1      0      0      0      1  ...    \n",
      "0      0      1      1      1      0      1      0      0      0      0  ...    \n",
      "0      0      1      1      0      0      1      0      1      0      1  ...    \n",
      "0      0      1      1      1      0      0      1      1      1      0  ...    \n",
      "0      0      1      0      0      0      1      0      0      1      1  ...    \n",
      "0      0      1      0      1      0      1      0      0      0      1  ...    \n",
      "0      1      0      0      1      0      1      1      0      1      0  ...    \n",
      "0      1      1      1      0      1      0      1      0      1      0  ...    \n",
      "0      1      0      0      0      1      0      1      1      0      0  ...    \n",
      "0      0      1      1      1      1      0      0      0      0      1  ...    \n",
      "0      0      0      0      0      1      0      1      0      1      0  ...    \n",
      "0      1      0      1      0      1      0      1      0      0      0  ...    \n",
      "0      1      0      1      0      1      1      1      1      1      0  ...    \n",
      "0      1      0      1      0      0      1      0      1      0      1  ...    \n",
      "0      1      0      0      1      1      0      0      1      1      0  ...    \n",
      "0      0      0      0      1      0      1      1      0      0      0  ...    \n",
      "0      0      1      0      1      0      1      0      0      1      0  ...    \n",
      "0      1      1      0      0      0      0      1      1      1      0  ...    \n",
      "0      1      0      0      0      1      1      1      1      1      1  ...    \n",
      "0      1      0      0      0      1      0      0      1      1      0  ...    \n",
      "0      1      1      0      0      1      1      0      1      0      1  ...    \n",
      "0      1      0      0      1      0      1      1      1      0      1  ...    \n",
      "0      0      1      0      0      0      1      1      1      0      0  ...    \n",
      "0      1      0      1      0      0      0      1      0      0      1  ...    \n",
      "0      1      1      1      0      1      0      0      0      1      0  ...    \n",
      "\n",
      "    s[37]  s[38]  s[39]  s[3]  s[4]  s[5]  s[6]  s[7]  s[8]  s[9]  \n",
      "0       1      1      0     1     0     0     1     1     0     0  \n",
      "0       1      1      1     0     0     0     1     1     1     1  \n",
      "0       1      0      1     0     1     0     0     1     1     1  \n",
      "0       1      0      1     0     0     0     0     0     0     1  \n",
      "0       0      1      0     1     0     0     1     1     0     0  \n",
      "0       1      1      0     0     0     0     0     0     1     1  \n",
      "0       1      0      0     1     0     0     1     1     0     0  \n",
      "0       1      0      1     1     0     1     0     1     1     0  \n",
      "0       1      0      1     0     1     0     0     1     1     1  \n",
      "0       0      0      0     0     1     1     0     1     0     0  \n",
      "0       1      0      0     0     0     0     1     1     0     0  \n",
      "0       0      1      0     0     1     1     1     0     0     1  \n",
      "0       0      1      0     1     1     1     1     0     0     0  \n",
      "0       0      0      0     1     1     1     1     1     1     0  \n",
      "0       0      1      1     1     1     0     0     0     0     1  \n",
      "0       1      1      0     0     1     1     0     1     0     0  \n",
      "0       0      1      0     1     1     0     1     0     1     0  \n",
      "0       1      0      1     1     0     0     1     1     1     0  \n",
      "0       1      0      1     1     0     1     0     1     1     0  \n",
      "0       0      0      1     1     0     0     1     0     1     0  \n",
      "0       0      0      1     0     1     1     0     0     0     1  \n",
      "0       1      1      1     1     0     0     1     0     1     0  \n",
      "0       1      0      1     0     1     0     0     1     0     0  \n",
      "0       0      1      0     0     1     1     1     0     1     1  \n",
      "0       1      1      1     0     0     0     1     1     0     1  \n",
      "0       0      1      0     1     1     1     1     0     1     0  \n",
      "0       1      1      0     1     0     0     1     1     0     0  \n",
      "0       0      1      0     1     0     0     0     1     1     0  \n",
      "0       0      0      1     0     1     0     0     0     0     1  \n",
      "0       0      1      0     1     1     1     1     0     1     0  \n",
      "..    ...    ...    ...   ...   ...   ...   ...   ...   ...   ...  \n",
      "0       0      0      0     1     0     1     1     0     0     0  \n",
      "0       0      0      0     0     1     0     0     1     1     0  \n",
      "0       1      0      1     1     0     1     0     0     0     1  \n",
      "0       1      0      0     1     0     1     0     0     1     1  \n",
      "0       0      1      1     1     1     1     0     1     0     0  \n",
      "0       0      0      1     1     1     0     0     1     1     1  \n",
      "0       0      1      1     0     1     1     0     0     1     1  \n",
      "0       0      0      1     0     0     1     0     0     0     0  \n",
      "0       1      1      1     0     0     1     0     0     0     1  \n",
      "0       0      0      0     1     1     0     0     1     0     0  \n",
      "0       0      0      1     1     1     1     1     1     0     0  \n",
      "0       1      1      1     0     1     0     0     0     0     1  \n",
      "0       0      0      0     1     0     0     0     0     1     1  \n",
      "0       1      1      1     1     1     1     0     1     1     0  \n",
      "0       1      0      1     1     0     0     1     1     0     0  \n",
      "0       0      0      0     1     1     0     0     1     1     0  \n",
      "0       1      0      0     1     1     1     1     1     1     0  \n",
      "0       1      0      0     1     0     0     1     0     0     0  \n",
      "0       0      1      0     0     1     0     0     1     0     0  \n",
      "0       0      0      1     0     1     0     1     1     0     0  \n",
      "0       1      0      1     0     0     1     1     0     1     1  \n",
      "0       1      1      0     0     0     0     1     1     0     1  \n",
      "0       0      1      0     0     1     0     0     0     1     1  \n",
      "0       0      1      0     0     1     1     0     0     0     0  \n",
      "0       0      0      1     0     0     0     1     1     0     0  \n",
      "0       1      0      1     0     0     0     0     1     1     0  \n",
      "0       1      1      0     0     1     0     1     0     0     1  \n",
      "0       1      0      1     0     0     1     1     1     1     0  \n",
      "0       1      1      0     1     0     0     0     0     0     1  \n",
      "0       1      0      1     1     1     1     0     1     0     0  \n",
      "\n",
      "[99 rows x 40 columns]\n"
     ]
    }
   ],
   "source": [
    "res=res.drop_duplicates()\n",
    "print(res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualize no. differences by covariate and bin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 172,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAEkBJREFUeJzt3XmQZWV9xvHv4wyCLAo4jaI4DKYMiVguVMc9asQFFR1TZVmgUsSYmorRRKPRaDQat5RZyqDllokiGBE17qHEMFEJrpgGQZZxQUVB0GlEcV9Gf/njnomXZrr7bn3v8Pr9VN3qc895731/feadp892z01VIUm66bvZrAuQJE2GgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXeokeWOSv511HdKoDHTNXJInJFlI8sMk1yQ5K8n9p11HVf1pVb1skLZJzknyJ6P2leTxST6V5MdJzhn1faR+BrpmKsmzgJOBvwduA2wEXg9snnId66bZH3Advd/7lVPuVw0z0DUzSW4FvBR4WlW9t6p+VFW/qKr/rKrndG32TnJykqu7x8lJ9u6WbU9yXN/7rU9ybZKju+f/keRbSa5Pcm6So/ranprkDUk+lORHwB90817eLT8oyZlJFpN8t5s+rFv2CuD3gdd2exWv7eb/TpJtSa5L8sUkj1/ud6+q/66qdwFXT3at6jeZga5Zug+wD/C+Fdq8ALg3cHfgbsA9gRd2y84ATuhr+3Dg2qq6oHt+FnAn4BDgAuD0Je/9BOAVwAHAJ5YsuxnwFuBwensNPwFeC1BVLwA+Djy9qvavqqcn2Q/YBry96+8E4PX9f0SktWaga5ZuTS+Ad67Q5onAS6tqR1UtAi8BTuyWvR14TJJ9u+dP6OYBUFWnVNUPqupnwN8Bd+v2Cnb5QFV9sqp+VVU/7e+0qr5TVe+pqh9X1Q/oBf8DV6jzOOCKqnpLVe3s/qi8B3jcKutAmhgDXbP0HWBDkvUrtLkd8PW+51/v5lFVlwPbgUd3of4YukBPsi7JK5N8Jcn3gSu612/oe68rl+s0yb5J/jXJ17vXnwscuMKx9sOBeyX53q4HvT9Gt13hd5MmykDXLH0a+Cnw2BXaXE0vLHfZyA2PO+867LIZuKwLeehtrW8GHgLcCtjUzU/fa1e61eizgSOBe1XVLYEHLHn90tdeCfxPVR3Y99i/qp66Qh/SRBnompmquh54EfC6JI/ttor3SvKIJP/YNTsDeGGSuSQbuvZv63ubdwAPA55K3+EWesfFf0ZvL2BfelfRDOMAesfNv5fkYODFS5Z/G7hj3/Mzgd9OcmL3O+yV5PeS/O7u3rzbg9gHWA/cLMk+SfYaskbpBgx0zVRVvQp4Fr0TnYv0tnSfDry/a/JyYAH4PHAxvZObL+97/TX0tvTvC7yz763fSu/wzDeBy4DPDFnaycAtgGu71354yfJXA4/rroB5TXec/WHA8fT2IL4F/AOw9zLvfyK9PxhvoHfFzE+AfxuyRukG4hdcSFIb3EKXpEYY6JLUCANdkhphoEtSI1b6QMfEbdiwoTZt2jTNLiXpJu/888+/tqrmVms31UDftGkTCwsL0+xSkm7yknx99VYecpGkZhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEZM9ZOi40iyeqM14j3j2+f4UgvcQpekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEasGuhJTkmyI8klu1n2V0kqyYa1KU+SNKhBttBPBY5dOjPJHYCHAt+YcE2SpBGsGuhVdS5w3W4W/QvwXMAbUUjSHmCkY+hJHgN8s6oumnA9kqQRDX23xST7Ai8AHjZg+y3AFoCNGzcO250kaUCjbKH/FnAEcFGSK4DDgAuS3HZ3jatqa1XNV9X83Nzc6JVKklY09BZ6VV0MHLLreRfq81V17QTrkiQNaZDLFs8APg0cmeSqJE9Z+7IkScNadQu9qk5YZfmmiVUjSRqZnxSVpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGjHIl0SfkmRHkkv65v1Tki8k+XyS9yU5cG3LlCStZpAt9FOBY5fM2wbcparuCnwJeP6E65IkDWnVQK+qc4Hrlsw7u6p2dk8/Axy2BrVJkoYwiWPofwyctdzCJFuSLCRZWFxcnEB3kqTdGSvQk7wA2AmcvlybqtpaVfNVNT83NzdOd5KkFawf9YVJTgKOA46pqppcSZKkUYwU6EmOBf4aeGBV/XiyJUmSRjHIZYtnAJ8GjkxyVZKnAK8FDgC2JbkwyRvXuE5J0ipW3UKvqhN2M/vNa1CLJGkMflJUkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJasQgXxJ9SpIdSS7pm3dwkm1Jvtz9PGhty5QkrWaQLfRTgWOXzHse8JGquhPwke65JGmGVg30qjoXuG7J7M3Aad30acBjJ1yXJGlI60d83W2q6hqAqromySHLNUyyBdgCsHHjxhG7k6TxJZlZ31W15n2s+UnRqtpaVfNVNT83N7fW3UnSb6xRA/3bSQ4F6H7umFxJkqRRjBroHwRO6qZPAj4wmXIkSaMa5LLFM4BPA0cmuSrJU4BXAg9N8mXgod1zSdIMrXpStKpOWGbRMROuRZI0Bj8pKkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWrEqHdb3OOs6X3MVrpD2xTuoKbfZLO7O+Aa/6/SGnALXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGjFWoCf5yySXJrkkyRlJ9plUYZKk4Ywc6EluD/wFMF9VdwHWAcdPqjBJ0nDGPeSyHrhFkvXAvsDV45ckSRrFyIFeVd8E/hn4BnANcH1Vnb20XZItSRaSLCwuLo5eqSRpReMccjkI2AwcAdwO2C/Jk5a2q6qtVTVfVfNzc3OjVypJWtE4h1weAnytqhar6hfAe4H7TqYsSdKwxgn0bwD3TrJvkgDHANsnU5YkaVjjHEM/D3g3cAFwcfdeWydUlyRpSGN9Y1FVvRh48YRqkSSNwU+KSlIjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0YK9CTHJjk3Um+kGR7kvtMqjBJ0nDG+pJo4NXAh6vqcUluDuw7gZokSSMYOdCT3BJ4APBHAFX1c+DnkylLkjSscQ653BFYBN6S5HNJ3pRkv6WNkmxJspBkYXFxcYzu9lSZ4UPTUmv4IFn+IQ1hnEBfDxwNvKGq7gH8CHje0kZVtbWq5qtqfm5ubozuJEkrGSfQrwKuqqrzuufvphfwkqQZGDnQq+pbwJVJjuxmHQNcNpGqJElDG/cqlz8HTu+ucPkq8OTxS5IkjWKsQK+qC4H5CdUiSRqDnxSVpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNWLcT4pK0o3N6k6RtcriVZaPY0+4OaZb6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1IixAz3JuiSfS3LmJAqSJI1mElvozwC2T+B9JEljGCvQkxwGPAp402TKkSSNatwt9JOB5wK/mkAtkqQxjBzoSY4DdlTV+au025JkIcnC4uLiqN1JklYxzhb6/YDHJLkCeAfw4CRvW9qoqrZW1XxVzc/NzY3RnSRpJSMHelU9v6oOq6pNwPHAR6vqSROrTJI0FK9Dl6RGTOQbi6rqHOCcSbyXJGk0bqFLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjRg50JPcIcnHkmxPcmmSZ0yyMEnScMb5kuidwLOr6oIkBwDnJ9lWVZdNqDZJ0hBG3kKvqmuq6oJu+gfAduD2kypMkjSciRxDT7IJuAdw3m6WbUmykGRhcXFxEt1JknZj7EBPsj/wHuCZVfX9pcuramtVzVfV/Nzc3LjdSZKWMVagJ9mLXpifXlXvnUxJkqRRjHOVS4A3A9ur6lWTK0mSNIpxttDvB5wIPDjJhd3jkROqS5I0pJEvW6yqTwCZYC2SpDH4SVFJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEWMFepJjk3wxyeVJnjepoiRJwxs50JOsA14HPAK4M3BCkjtPqjBJ0nDG2UK/J3B5VX21qn4OvAPYPJmyJEnDWj/Ga28PXNn3/CrgXksbJdkCbOme/jDJF0fsbwNw7XILM+Kbji0r17XWna9ghnWt6CZZl+PrRvbMf8eZrq+VJRmntsMHaTROoO/uX7tuNKNqK7B1jH56nSULVTU/7vtMmnUNx7qGY13D2VPrgunUNs4hl6uAO/Q9Pwy4erxyJEmjGifQ/xe4U5IjktwcOB744GTKkiQNa+RDLlW1M8nTgf8C1gGnVNWlE6vsxsY+bLNGrGs41jUc6xrOnloXTKG2VN3osLck6SbIT4pKUiMMdElqxMwDPckpSXYkuWSZ5Unymu72Ap9PcnTfspOSfLl7nDTlup7Y1fP5JJ9Kcre+ZVckuTjJhUkWplzXg5Jc3/V9YZIX9S1bs1s1DFDXc/pquiTJL5Mc3C1by/V1hyQfS7I9yaVJnrGbNlMfYwPWNfUxNmBdUx9jA9Y19TGWZJ8kn01yUVfXS3bTZu8k7+zWyXlJNvUte343/4tJHj52QVU10wfwAOBo4JJllj8SOIvede/3Bs7r5h8MfLX7eVA3fdAU67rvrv7o3f7gvL5lVwAbZrS+HgScuZv564CvAHcEbg5cBNx5WnUtafto4KNTWl+HAkd30wcAX1r6e89ijA1Y19TH2IB1TX2MDVLXLMZYN2b276b3As4D7r2kzZ8Bb+ymjwfe2U3fuVtHewNHdOtu3Tj1zHwLvarOBa5boclm4K3V8xngwCSHAg8HtlXVdVX1XWAbcOy06qqqT3X9AnyG3nX4a26A9bWcNb1Vw5B1nQCcMam+V1JV11TVBd30D4Dt9D7l3G/qY2yQumYxxgZcX8tZszE2Ql1TGWPdmPlh93Sv7rH0SpPNwGnd9LuBY5Kkm/+OqvpZVX0NuJzeOhzZzAN9ALu7xcDtV5g/C0+ht4W3SwFnJzk/vVsfTNt9ul3As5Ic1c3bI9ZXkn3pheJ7+mZPZX11u7r3oLcV1W+mY2yFuvpNfYytUtfMxthq62vaYyzJuiQXAjvobQAsO76qaidwPXBr1mB9jfPR/2lZ7hYDA916YK0l+QN6/9nu3zf7flV1dZJDgG1JvtBtwU7DBcDhVfXDJI8E3g/ciT1kfdHbFf5kVfVvza/5+kqyP73/4M+squ8vXbybl0xljK1S1642Ux9jq9Q1szE2yPpiymOsqn4J3D3JgcD7ktylqvrPJU1tfN0UttCXu8XAzG89kOSuwJuAzVX1nV3zq+rq7ucO4H2MuRs1jKr6/q5dwKr6ELBXejcFmvn66hzPkl3htV5fSfaiFwKnV9V7d9NkJmNsgLpmMsZWq2tWY2yQ9dWZ+hjr3vt7wDnc+LDc/6+XJOuBW9E7PDn59TXJEwSjPoBNLH+S71Hc8ITVZ7v5BwNfo3ey6qBu+uAp1rWR3jGv+y6Zvx9wQN/0p4Bjp1jXbfn1B8buCXyjW3fr6Z3UO4Jfn7A6alp1dct3DeT9prW+ut/9rcDJK7SZ+hgbsK6pj7EB65r6GBukrlmMMWAOOLCbvgXwceC4JW2exg1Pir6rmz6KG54U/SpjnhSd+SGXJGfQO2u+IclVwIvpnVigqt4IfIjeVQiXAz8Gntwtuy7Jy+jdUwbgpXXDXay1rutF9I6Dvb53foOd1buT2m3o7XZBb4C/vao+PMW6Hgc8NclO4CfA8dUbPWt6q4YB6gL4Q+DsqvpR30vXdH0B9wNOBC7ujnMC/A29sJzlGBukrlmMsUHqmsUYG6QumP4YOxQ4Lb0v/LkZvbA+M8lLgYWq+iDwZuDfk1xO74/N8V3NlyZ5F3AZsBN4WvUO34zMj/5LUiNuCsfQJUkDMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSI/4PK5JQtDqmyOQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAE0hJREFUeJzt3XuwZWV95vHvYzeiXCJgH40RsXFKyYilkToqYmK8RG0N2qbKTIGXkMSprnE0YyYZMxozMhpNZRIrQUsN1aUtmmCrUbyEEmNXEkMStc2BAbk0BEZRWzB9kIiAt7T+5o+9WjeHPmdfz97Ny/dTteusvda79/s7q99+ztprrb1WqgpJ0t3fveZdgCRpOgx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOhSJ8m5Sf7XvOuQxmWga+6SvDDJUpLbk9yU5KIkPzvrOqrqv1TV7w/TNsmnk/zncftK8uYk1yW5Lck1SX5l3PeSDjDQNVdJfgs4B/gD4IHACcA7gK0zrmPDLPsD7gCeC9wPOAt4S5LTZlyDGmOga26S3A94A/Dyqrqgqu6oqn+vqr+qqld1bQ5Pck6SG7vHOUkO75btSXJ63/ttTHJzklO653+Z5OtJbk1ycZKT+9qel+TPknwiyR3AU7t5b+yWH5vkwiTLSf6tmz6+W/Ym4OeAt3WfKt7Wzf/pJLuS3JLk2iT/abXfvarOrqprquqHVbUb+AfgiVNdwbrHMdA1T08E7gN8ZI02rwVOBX4GeAzweOD3umU7gTP72j4LuLmqLu2eXwQ8HHgAcClw/or3fiHwJuBo4B9XLLsX8G7gofQ+NXwHeBtAVb2WXgC/oqqOqqpXJDkS2AW8r+vvTOAd/X9EVpPkvsDjgKsGtZXWYqBrnu5PL4D3r9HmRcAbqmpfVS0Drwde0i17H/C8JEd0z1/YzQOgqnZU1W1V9T3gfwOP6T4VHPCxqvqnbiv5u/2dVtU3qurDVfXtqrqNXvD//Bp1ng7cUFXvrqr93R+VDwMvGLAOAM4FLgf+eoi20qo2zrsA3aN9A9iUZOMaof5TwJf7nn+5m0dVXZ9kD/DcJH8FPA94LPxon/ibgF8GFoAfdq/fBNzaTX91tcK6PxJ/CmwBju1mH51kQ1X94CAveSjwhCTf7Ju3Efjz1fro+vlj4FHAU8sr5WlCbqFrnj4LfBd4/hptbqQXlgec0M074MBul63A1VV1fTf/hd28X6B34HFzNz99r10rQH8bOAl4QlX9BPDkFa9f+dqvAn9fVcf0PY6qqpet1kGS1wPPBp5ZVd9aoxZpKAa65qaqbgVeB7w9yfOTHJHksCTPTvJHXbOdwO8lWUiyqWv/F31v837gmcDL6NvdQm+/+PfofQo4gt5ZNKM4mt5+828mOQ44e8XyfwUe1vf8QuARSV7S/Q6HJXlckv94sDdP8hp6f3SeUVXfGLE26aAMdM1VVf0J8Fv0DnQu09vSfQXw0a7JG4El4AvAFfQObr6x7/U30dvSPw34QN9bv5fe7pmvAVcDnxuxtHOA+wI3d6/95IrlbwFe0J0B89ZuP/szgTPofYL4OvB/gMNXef8/oPdp47ruTJnbk/zuiDVKdxJ320lSG9xCl6RGGOiS1AgDXZIaYaBLUiNm+sWiTZs21ebNm2fZpSTd7V1yySU3V9XCoHYzDfTNmzeztLQ0yy4l6W4vyZcHt3KXiyQ1w0CXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJ7ikq6x0gyuNE6mcW9J9xCl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDViYKAn2ZFkX5IrV8z/jSTXJrkqyR+tX4mSpGEMs4V+HrClf0aSpwJbgUdX1cnAm6dfmiRpFAMDvaouBm5ZMftlwB9W1fe6NvvWoTZJ0gjG3Yf+CODnkuxO8vdJHrdawyTbkiwlWVpeXh6zO0nSIOMG+kbgWOBU4FXAB7PKVW+qantVLVbV4sLCwpjdSZIGGTfQ9wIXVM/ngR8Cm6ZXliRpVOMG+keBpwEkeQRwb+DmaRUlSRrdwOuhJ9kJPAXYlGQvcDawA9jRncr4feCsmsXFfiVJqxoY6FV15iqLXjzlWiRJE/CbopLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDViYKAn2ZFkX3czi5XL/keSSuLt5yRpzobZQj8P2LJyZpKHAM8AvjLlmiRJYxgY6FV1MXDLQRb9KfA7gLeek6RDwMBb0B1MkucBX6uqy5MMarsN2AZwwgknjNOdtO4GjeP15O14NS0jHxRNcgTwWuB1w7Svqu1VtVhViwsLC6N2J0ka0jhnufwH4ETg8iQ3AMcDlyb5yWkWJkkazci7XKrqCuABB553ob5YVTdPsS5J0oiGOW1xJ/BZ4KQke5O8dP3LkiSNauAWelWdOWD55qlVI0kam98UlaRGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqRHD3OBiR5J9Sa7sm/fHSa5J8oUkH0lyzPqWKUkaZJgt9POALSvm7QIeVVWPBv4FeM2U65IkjWhgoFfVxcAtK+Z9qqr2d08/R+9G0ZKkOZrGPvRfBy6awvtIkiYwUaAneS2wHzh/jTbbkiwlWVpeXp6kO0nSGsYO9CRnAacDL6qqWq1dVW2vqsWqWlxYWBi3O0nSABvHeVGSLcD/BH6+qr493ZIkSeMY5rTFncBngZOS7E3yUuBtwNHAriSXJTl3neuUJA0wcAu9qs48yOx3rUMtkqQJ+E1RSWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1Ijhrlj0Y4k+5Jc2TfvuCS7klzX/Tx2fcuUJA0yzBb6ecCWFfNeDfxNVT0c+JvuuSRpjgYGelVdDNyyYvZW4D3d9HuA50+5LknSiMbdh/7AqroJoPv5gNUaJtmWZCnJ0vLy8pjdSZIGWfeDolW1vaoWq2pxYWFhvbuTpHuscQP9X5M8CKD7uW96JUmSxjFuoH8cOKubPgv42HTKkSSNa5jTFncCnwVOSrI3yUuBPwSekeQ64Bndc0nSHG0c1KCqzlxl0dOnXIskaQJ+U1SSGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1YuBpi4eKJHPru6rm1rekUaydE+v5X3mOEfUjbqFLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGjFRoCf570muSnJlkp1J7jOtwiRJoxk70JM8GPhvwGJVPQrYAJwxrcIkSaOZdJfLRuC+STYCRwA3Tl6SJGkcYwd6VX0NeDPwFeAm4Naq+tTKdkm2JVlKsrS8vDx+pZLuPpL5PO7hJtnlciywFTgR+CngyCQvXtmuqrZX1WJVLS4sLIxfqSRpTZPscvkF4EtVtVxV/w5cAJw2nbIkSaOaJNC/Apya5Ij0rm37dGDPdMqSJI1qkn3ou4EPAZcCV3TvtX1KdUmSRjTRDS6q6mzg7CnVIkmagN8UlaRGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEROdhy7dU9R6vvlaF5Va144HmWvnGoNb6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGTBToSY5J8qEk1yTZk+SJ0ypMkjSaSb8p+hbgk1X1giT3Bo6YQk2SpDGMHehJfgJ4MvCrAFX1feD70ylLkjSqSXa5PAxYBt6d5P8meWeSI1c2SrItyVKSpeXl5Qm6kyStZZJA3wicAvxZVT0WuAN49cpGVbW9qharanFhYWGC7iRJa5kk0PcCe6tqd/f8Q/QCXpI0B2MHelV9HfhqkpO6WU8Hrp5KVZKkkU16lstvAOd3Z7h8Efi1yUuSJI1jokCvqsuAxSnVIkmagN8UlaRGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY2Y9Dz0Q0at55sna3S8rj1L0tDcQpekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1YuJAT7Khu0n0hdMoSJI0nmlsob8S2DOF95EkTWCiQE9yPPCLwDunU44kaVyTbqGfA/wO8MPVGiTZlmQpydLy8vKE3R2KMseHJP3Y2IGe5HRgX1Vdsla7qtpeVYtVtbiwsDBud5KkASbZQn8S8LwkNwDvB56W5C+mUpUkaWRjB3pVvaaqjq+qzcAZwN9W1YunVpkkaSSehy5JjZjKDS6q6tPAp6fxXpKk8biFLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMmuafoQ5L8XZI9Sa5K8sppFiZJGs0kN7jYD/x2VV2a5GjgkiS7qurqKdUmSRrBJPcUvamqLu2mbwP2AA+eVmGSpNFMZR96ks3AY4HdB1m2LclSkqXl5eVpdCdJOoiJAz3JUcCHgd+sqm+tXF5V26tqsaoWFxYWJu1OkrSKiQI9yWH0wvz8qrpgOiVJksYxyVkuAd4F7KmqP5leSZKkcUyyhf4k4CXA05Jc1j2eM6W6JEkjGvu0xar6RyBTrEWSNAG/KSpJjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJasSk9xTdkuTaJNcnefW0ipIkjW6Se4puAN4OPBt4JHBmkkdOqzBJ0mgm2UJ/PHB9VX2xqr4PvB/YOp2yJEmjGvueosCDga/2Pd8LPGFloyTbgG3d09uTXDtmf5uAm1dbOLebm2btuta78zXMsa413S3rcnzdxaH57zjX9bW2JJPU9tBhGk0S6Af71667zKjaDmyfoJ9eZ8lSVS1O+j7TZl2jsa7RWNdoDtW6YDa1TbLLZS/wkL7nxwM3TlaOJGlckwT6PwMPT3JiknsDZwAfn05ZkqRRjb3Lpar2J3kF8NfABmBHVV01tcruauLdNuvEukZjXaOxrtEcqnXBDGpL1V12e0uS7ob8pqgkNcJAl6RGzD3Qk+xIsi/JlassT5K3dpcX+EKSU/qWnZXkuu5x1ozrelFXzxeSfCbJY/qW3ZDkiiSXJVmacV1PSXJr1/dlSV7Xt2zdLtUwRF2v6qvpyiQ/SHJct2w919dDkvxdkj1JrkryyoO0mfkYG7KumY+xIeua+Rgbsq6Zj7Ek90ny+SSXd3W9/iBtDk/ygW6d7E6yuW/Za7r51yZ51sQFVdVcH8CTgVOAK1dZ/hzgInrnvZ8K7O7mHwd8sft5bDd97AzrOu1Af/Quf7C7b9kNwKY5ra+nABceZP4G4P8BDwPuDVwOPHJWda1o+1zgb2e0vh4EnNJNHw38y8rfex5jbMi6Zj7Ghqxr5mNsmLrmMca6MXNUN30YsBs4dUWb/wqc202fAXygm35kt44OB07s1t2GSeqZ+xZ6VV0M3LJGk63Ae6vnc8AxSR4EPAvYVVW3VNW/AbuALbOqq6o+0/UL8Dl65+GvuyHW12rW9VINI9Z1JrBzWn2vpapuqqpLu+nbgD30vuXcb+ZjbJi65jHGhlxfq1m3MTZGXTMZY92Yub17elj3WHmmyVbgPd30h4CnJ0k3//1V9b2q+hJwPb11OLa5B/oQDnaJgQevMX8eXkpvC++AAj6V5JL0Ln0wa0/sPgJelOTkbt4hsb6SHEEvFD/cN3sm66v7qPtYeltR/eY6xtaoq9/Mx9iAuuY2xgatr1mPsSQbklwG7KO3AbDq+Kqq/cCtwP1Zh/U1yVf/Z2W1SwwMdemB9ZbkqfT+s/1s3+wnVdWNSR4A7EpyTbcFOwuXAg+tqtuTPAf4KPBwDpH1Re+j8D9VVf/W/LqvryRH0fsP/ptV9a2Viw/ykpmMsQF1HWgz8zE2oK65jbFh1hczHmNV9QPgZ5IcA3wkyaOqqv9Y0szG191hC321SwzM/dIDSR4NvBPYWlXfODC/qm7sfu4DPsKEH6NGUVXfOvARsKo+ARyW3kWB5r6+Omew4qPweq+vJIfRC4Hzq+qCgzSZyxgboq65jLFBdc1rjA2zvjozH2Pde38T+DR33S33o/WSZCNwP3q7J6e/vqZ5gGDcB7CZ1Q/y/SJ3PmD1+W7+ccCX6B2sOrabPm6GdZ1Ab5/XaSvmHwkc3Tf9GWDLDOv6SX78hbHHA1/p1t1Gegf1TuTHB6xOnlVd3fIDA/nIWa2v7nd/L3DOGm1mPsaGrGvmY2zIumY+xoapax5jDFgAjumm7wv8A3D6ijYv584HRT/YTZ/MnQ+KfpEJD4rOfZdLkp30jppvSrIXOJvegQWq6lzgE/TOQrge+Dbwa92yW5L8Pr1rygC8oe78EWu963odvf1g7+gd32B/9a6k9kB6H7ugN8DfV1WfnGFdLwBelmQ/8B3gjOqNnnW9VMMQdQH8EvCpqrqj76Xrur6AJwEvAa7o9nMC/C69sJznGBumrnmMsWHqmscYG6YumP0YexDwnvRu+HMvemF9YZI3AEtV9XHgXcCfJ7me3h+bM7qar0ryQeBqYD/w8urtvhmbX/2XpEbcHfahS5KGYKBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRvx/BQWiZ1CHtN8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAEoNJREFUeJzt3XmQZWV9xvHv4wyiLBFw2hWHQWNI1HKhOi5ojIoLKjKaslLgUiQxNYnRRGOi0WA0Gk0ZkzJo3DKluCJoXKKhxDjlEjQqpkGUZURRUUYw04iyuY7+8sc9Ey/NdPfd+t7h9fupujXnnvPe8/768PL0e8+593SqCknSTd/NZl2AJGkyDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6FInyRuT/M2s65BGZaBr5pI8KclCkuuSXJHkzCQPmnYdVfXHVfV3g7RN8skkfzhqX0lemeSyJNck+WaSk0bdl7Sbga6ZSvIc4GTg74HbAhuB1wObp1zHumn2B7wZ+PWq+hXgKOBJSX5nyjWoMQa6ZibJrYCXAs+oqvdX1fVV9dOq+o+qem7XZt8kJye5vHucnGTfbtv2JMf27W99kiuTHNk9/7ck30lydZKzkty9r+1bk7whyYeTXA88tFv3sm77wUnOSLKY5Hvd8qHdtpcDvwW8tntX8dpu/a8n2ZbkqiQXJ/nd5X72qrq4qq7vW/Vz4FcncmD1S8tA1yw9ALgF8IEV2pwE3B+4N3Av4L7AC7ttpwEn9LV9FHBlVZ3bPT8TuCtwG+Bc4NQl+34S8HLgQODTS7bdDHgLcBi9dw0/BF4LUFUnAZ8CnllVB1TVM5PsD2wD3tX1dwLw+v5fIksleX6S64AdwP7da6WRGeiapVvTC+BdK7R5MvDSqtpZVYvAS4CndtveBRyXZL/u+ZPoC8WqOqWqrq2qHwN/C9yre1ew2wer6r+r6udV9aP+Tqvqu1X1vqr6QVVdSy/4f3uFOo8FLq2qt1TVru6XyvuAJy73gqp6Bb1fJkcC7wCuXmH/0qoMdM3Sd4ENSdav0OYOwDf7nn+zW0dVXQJsBx7XhfpxdIGeZF2SVyT5WpJrgEu712/o29dly3WaZL8k/9pdsLwGOAs4aIVz7YcB90vy/d0Per+MbrfCz0b1fIHeO4CXrNRWWo2Brln6LPAj4PErtLmcXljutrFbt9vu0y6bgYu6kIfebH0z8HDgVsCmbn36XrvSrUb/AjgCuF934fLBS16/9LWXAf9VVQf1PQ6oqqev0Ee/9cBdBmwr7ZGBrpmpqquBFwGvS/L4bla8T5JHJ3ll1+w04IVJ5pJs6Nq/s283pwOPBJ7ODc9BHwj8mN67gP3ofYpmGAfSmzV/P8khwIuXbP9f4M59z88Afi3JU7ufYZ8kv5nkN5buOMnNkvxRd+E1Se4LPAP42JA1SjdgoGumqupVwHPoXehcpDfTfSbw712TlwELwJeA8+ld3HxZ3+uvoDfTPwp4d9+u307v9My3gYuAzw1Z2snALYEru9d+ZMn2VwNP7D4B85ruPPsjgePpvYP4DvAPwL7L7P8JwNeAa+n9gvqX7iGNLP6BC0lqgzN0SWqEgS5JjTDQJakRBrokNWKlL3RM3IYNG2rTpk3T7FKSbvLOOeecK6tqbrV2Uw30TZs2sbCwMM0uJekmL8k3V2/lKRdJaoaBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWrEVL8pOo4kqzdaI94zvn2OL7XAGbokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDVi1UBPckqSnUku2MO2v0xSSTasTXmSpEENMkN/K3DM0pVJ7gQ8AvjWhGuSJI1g1UCvqrOAq/aw6Z+B5wHeiEKS9gIjnUNPchzw7ar64oTrkSSNaOi7LSbZDzgJeOSA7bcAWwA2btw4bHeSpAGNMkO/C3A48MUklwKHAucmud2eGlfV1qqar6r5ubm50SuVJK1o6Bl6VZ0P3Gb38y7U56vqygnWJUka0iAfWzwN+CxwRJIdSZ629mVJkoa16gy9qk5YZfumiVUjSRqZ3xSVpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGjHIH4k+JcnOJBf0rfvHJF9O8qUkH0hy0NqWKUlazSAz9LcCxyxZtw24R1XdE/gK8IIJ1yVJGtKqgV5VZwFXLVn30ara1T39HHDoGtQmSRrCJM6h/wFw5nIbk2xJspBkYXFxcQLdSZL2ZKxAT3ISsAs4dbk2VbW1quaran5ubm6c7iRJK1g/6guTnAgcCxxdVTW5kiRJoxgp0JMcA/wV8NtV9YPJliRJGsUgH1s8DfgscESSHUmeBrwWOBDYluS8JG9c4zolSatYdYZeVSfsYfWb16AWSdIY/KaoJDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREj38tlb7OmN5NJVujY29hI2js4Q5ekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUiEH+SPQpSXYmuaBv3SFJtiX5avfvwWtbpiRpNYPM0N8KHLNk3fOBj1XVXYGPdc8lSTO0aqBX1VnAVUtWbwbe1i2/DXj8hOuSJA1p1Lst3raqrgCoqiuS3Ga5hkm2AFsANm7cOGJ3kjS+rHTn1DVWU7gz65pfFK2qrVU1X1Xzc3Nza92dJP3SGjXQ/zfJ7QG6f3dOriRJ0ihGDfQPASd2yycCH5xMOZKkUQ3yscXTgM8CRyTZkeRpwCuARyT5KvCI7rkkaYZWvShaVScss+noCdciSRqD3xSVpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNWLUuy1KmorZ3R0Q1v7ugJosZ+iS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktSIsQI9yZ8nuTDJBUlOS3KLSRUmSRrOyIGe5I7AnwHzVXUPYB1w/KQKkyQNZ9xTLuuBWyZZD+wHXD5+SZKkUYwc6FX1beCfgG8BVwBXV9VHl7ZLsiXJQpKFxcXF0SuVJK1onFMuBwObgcOBOwD7J3nK0nZVtbWq5qtqfm5ubvRKJUkrGueUy8OBb1TVYlX9FHg/cNRkypIkDWucQP8WcP8k+yUJcDSwfTJlSZKGNc459LOB9wLnAud3+9o6obokSUMa6y8WVdWLgRdPqBZJ0hj8pqgkNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEWN9sUiS9i5ZcWvVGva8ctdT4QxdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiPGCvQkByV5b5IvJ9me5AGTKkySNJxx7+XyauAjVfXEJDcH9ptATZKkEYwc6El+BXgw8HsAVfUT4CeTKUuSNKxxTrncGVgE3pLkC0nelGT/pY2SbEmykGRhcXFxjO6k2ak1fJAs/5CGME6grweOBN5QVfcBrgeev7RRVW2tqvmqmp+bmxujO0nSSsYJ9B3Ajqo6u3v+XnoBL0magZEDvaq+A1yW5Ihu1dHARROpSpI0tHE/5fKnwKndJ1y+Dvz++CVJkkYxVqBX1XnA/IRqkSSNwW+KSlIjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqxNiBnmRdki8kOWMSBUmSRjOJGfqzgO0T2I8kaQxjBXqSQ4HHAm+aTDmSpFGNO0M/GXge8PMJ1CJJGsPIgZ7kWGBnVZ2zSrstSRaSLCwuLo7a3V4sM3xI0i+MM0N/IHBckkuB04GHJXnn0kZVtbWq5qtqfm5ubozuJEkrGTnQq+oFVXVoVW0Cjgc+XlVPmVhlkqSh+Dl0SWrE+knspKo+CXxyEvuSJI3GGbokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpERP5YpEk3UBmdPO4mk23ewtn6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaMXKgJ7lTkk8k2Z7kwiTPmmRhkqThjHMvl13AX1TVuUkOBM5Jsq2qLppQbZKkIYw8Q6+qK6rq3G75WmA7cMdJFSZJGs5EzqEn2QTcBzh7D9u2JFlIsrC4uDiJ7iRJezB2oCc5AHgf8Oyqumbp9qraWlXzVTU/Nzc3bneSpGWMFehJ9qEX5qdW1fsnU5IkaRTjfMolwJuB7VX1qsmVJEkaxTgz9AcCTwUeluS87vGYCdUlSRrSyB9brKpPAzP6O1OSpKX8pqgkNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUiLECPckxSS5OckmS50+qKEnS8EYO9CTrgNcBjwbuBpyQ5G6TKkySNJxxZuj3BS6pqq9X1U+A04HNkylLkjSs9WO89o7AZX3PdwD3W9ooyRZgS/f0uiQXj9jfBuDK5TZmxJ2OLSvXtdadr2CGda3oJlmX4+tG9s7/jjM9XitLMk5thw3SaJxA39N/7brRiqqtwNYx+ul1lixU1fy4+5k06xqOdQ3Huoazt9YF06ltnFMuO4A79T0/FLh8vHIkSaMaJ9D/B7hrksOT3Bw4HvjQZMqSJA1r5FMuVbUryTOB/wTWAadU1YUTq+zGxj5ts0asazjWNRzrGs7eWhdMobZU3ei0tyTpJshvikpSIwx0SWrEzAM9ySlJdia5YJntSfKa7vYCX0pyZN+2E5N8tXucOOW6ntzV86Ukn0lyr75tlyY5P8l5SRamXNdDklzd9X1ekhf1bVuzWzUMUNdz+2q6IMnPkhzSbVvL43WnJJ9Isj3JhUmetYc2Ux9jA9Y19TE2YF1TH2MD1jX1MZbkFkk+n+SLXV0v2UObfZO8uzsmZyfZ1LftBd36i5M8auyCqmqmD+DBwJHABctsfwxwJr3Pvd8fOLtbfwjw9e7fg7vlg6dY11G7+6N3+4Oz+7ZdCmyY0fF6CHDGHtavA74G3Bm4OfBF4G7TqmtJ28cBH5/S8bo9cGS3fCDwlaU/9yzG2IB1TX2MDVjX1MfYIHXNYox1Y+aAbnkf4Gzg/kva/Anwxm75eODd3fLdumO0L3B4d+zWjVPPzGfoVXUWcNUKTTYDb6+ezwEHJbk98ChgW1VdVVXfA7YBx0yrrqr6TNcvwOfofQ5/zQ1wvJazprdqGLKuE4DTJtX3Sqrqiqo6t1u+FthO71vO/aY+xgapaxZjbMDjtZw1G2Mj1DWVMdaNmeu6p/t0j6WfNNkMvK1bfi9wdJJ060+vqh9X1TeAS+gdw5HNPNAHsKdbDNxxhfWz8DR6M7zdCvhoknPSu/XBtD2gewt4ZpK7d+v2iuOVZD96ofi+vtVTOV7dW9370JtF9ZvpGFuhrn5TH2Or1DWzMbba8Zr2GEuyLsl5wE56E4Blx1dV7QKuBm7NGhyvcb76Py3L3WJgoFsPrLUkD6X3P9uD+lY/sKouT3IbYFuSL3cz2Gk4Fzisqq5L8hjg34G7spccL3pvhf+7qvpn82t+vJIcQO9/8GdX1TVLN+/hJVMZY6vUtbvN1MfYKnXNbIwNcryY8hirqp8B905yEPCBJPeoqv5rSVMbXzeFGfpytxiY+a0HktwTeBOwuaq+u3t9VV3e/bsT+ABjvo0aRlVds/stYFV9GNgnvZsCzfx4dY5nyVvhtT5eSfahFwKnVtX799BkJmNsgLpmMsZWq2tWY2yQ49WZ+hjr9v194JPc+LTc/x+XJOuBW9E7PTn54zXJCwSjPoBNLH+R77Hc8ILV57v1hwDfoHex6uBu+ZAp1rWR3jmvo5as3x84sG/5M8AxU6zrdvziC2P3Bb7VHbv19C7qHc4vLljdfVp1ddt3D+T9p3W8up/97cDJK7SZ+hgbsK6pj7EB65r6GBukrlmMMWAOOKhbviXwKeDYJW2ewQ0vir6nW747N7wo+nXGvCg681MuSU6jd9V8Q5IdwIvpXVigqt4IfJjepxAuAX4A/H637aokf0fvnjIAL60bvsVa67peRO882Ot71zfYVb07qd2W3tsu6A3wd1XVR6ZY1xOBpyfZBfwQOL56o2dNb9UwQF0ATwA+WlXX9710TY8X8EDgqcD53XlOgL+mF5azHGOD1DWLMTZIXbMYY4PUBdMfY7cH3pbeH/y5Gb2wPiPJS4GFqvoQ8GbgHUkuoffL5viu5guTvAe4CNgFPKN6p29G5lf/JakRN4Vz6JKkARjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqRH/Byu1aLFF+JWUAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAE0BJREFUeJzt3XuQZGddxvHvw24IBCJJ2AFDwmaDBVFCQYgDhIvKRWSBQLAKrYRLRcXaEomCFxQEiSJYeCmIFiC1JUu4hAUM4ZYiyJaKUYHFSQjkskQiBFiSuBMCAQIEFn7+0WehmezM9G26d998P1Vde/qct/v9zdl3njl9zulzUlVIkg5+d5h1AZKkyTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLfZK8IcmfzroOaRQGug4ISZ6RZCHJN5Ncn+SiJI+adh1V9VtV9ReDtE3ykSS/OW6fSY5KspjkP8d9L92+GeiauSS/D5wD/CVwT2Aj8HrgtCnXsW6a/fX5K2DXjPpWQwx0zVSSuwEvB55XVRdU1S1V9b2q+kBVvbBrc2iSc5Jc1z3OSXJot2xXklP73m99khuTnNw9/6ckNyS5OcnFSU7sa3tukn9I8sEktwCP6ea9olt+ZJILu63nr3bTx3bLXgn8HPDa7lPFa7v5P51kR5Kbklyd5FdX+fkfDjwAeNPEVqputwx0zdrDgTsB71mhzUuAU4CTgAcBDwVe2i3bDpzR1/YJwI1VdWn3/CLgvsA9gEuB85a89zOAVwKHA0t3edyBXtAeR+9Tw7eB1wJU1UuA/wDOqqq7VtVZSe4C7ADe3vV3BvD6/j8i/bpPBK8DzgK8BofGZqBr1u5OL4D3rtDmmcDLq2pPVS0Cfw48u1v2duCpSQ7rnj+jmwdAVW2rqm9U1a3AnwEP6j4V7PO+qvqvqvpBVX2nv9Oq+kpVvbuqvlVV36AX/L+wQp2nAtdW1Zuqam/3R+XdwNOXaf+7wM6qumSF95QGtn7WBeh27yvAhiTrVwj1ewFf6Hv+hW4eVXVNkl3AU5J8AHgq8GD44RbwK4FfAeaAH3Sv3wDc3E1/abnCuj8SrwE2A0d2sw9Psq6qvr+flxwHPCzJ1/rmrQfeup/3vhe9QP/Z5fqXhmWga9Y+BnwHeBpw/jJtrqMXlld2zzd28/bZt9vlDsBVVXVNN/8Z9A6s/iJwLXA34KtA+l670q6OPwBOAB5WVTckOQn4ZN/rl772S8C/V9XjV3jPfR4KHA1clQTgzsCdk9wAHLPMHwxpRe5y0UxV1c3Ay4DXJXlaksOSHJLkiUn+umu2HXhpkrkkG7r2b+t7m3cAvwQ8l77dLfT2i99K71PAYfTOohnG4fT2m38tyVHA2UuW/x9wn77nFwL3S/Ls7mc4JMlDkvzMft77ImATveMCJ3U/0yeBkwxzjcpA18xV1auB36d3oHOR3pbuWcB7uyavABaATwOX0zu4+Yq+119Pb0v/EcA7+976LfR2z3wZuAr4+JClnUNvy/nG7rUfWrL874Cnd2fA/H23n/2XgNPpfYK4gd4piYfu52e+tapu2Pegtwvoe920NJJ4gwtJaoNb6JLUCANdkhphoEtSIwx0SWrEVM9D37BhQ23atGmaXUrSQe+SSy65sarmVms31UDftGkTCwsL0+xSkg56Sb6weit3uUhSMwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMOmlvQdbfpmgmvGS/pYOAWuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRqwZ6km1J9iS5Ysn830lydZIrk/z12pUoSRrEIFvo5wKb+2ckeQxwGvDAqjoR+NvJlyZJGsaqgV5VFwM3LZn9XOBVVXVr12bPGtQmSRrCqPvQ7wf8XJKdSf49yUOWa5hkS5KFJAuLi4sjdidJWs2ogb4eOBI4BXgh8K4sc/WsqtpaVfNVNT83Nzdid5Kk1Ywa6LuBC6rnE8APgA2TK0uSNKxRA/29wGMBktwPuCNw46SKkiQNb9XroSfZDjwa2JBkN3A2sA3Y1p3K+F3gzPKi4ZI0U6sGelWdscyiZ024FknSGPymqCQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY1YNdCTbEuyp7uZxdJlf5ikknj7OUmasUG20M8FNi+dmeTewOOBL064JknSCFYN9Kq6GLhpP4teA/wR4K3nJOkAMNI+9CRPBb5cVZ8aoO2WJAtJFhYXF0fpTpI0gKEDPclhwEuAlw3Svqq2VtV8Vc3Pzc0N250kaUCjbKH/FHA88Kkk1wLHApcm+clJFiZJGs76YV9QVZcD99j3vAv1+aq6cYJ1SZKGNMhpi9uBjwEnJNmd5DlrX5YkaVirbqFX1RmrLN80sWokSSPzm6KS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEUN/sUhqUZKZ9V3l9e00GW6hS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhoxyA0utiXZk+SKvnl/k+QzST6d5D1JjljbMiVJqxlkC/1cYPOSeTuAB1TVA4H/AV484bokSUNaNdCr6mLgpiXzPlxVe7unH6d3o2hJ0gxNYh/6bwAXTeB9JEljGCvQk7wE2Auct0KbLUkWkiwsLi6O050kaQUjB3qSM4FTgWfWCpeLq6qtVTVfVfNzc3OjdidJWsVIl89Nshn4Y+AXqupbky1JkjSKQU5b3A58DDghye4kzwFeCxwO7EhyWZI3rHGdkqRVrLqFXlVn7Gf2G9egFknSGPymqCQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjRjpm6IHomWvPTAJyQodr2nPkjQwt9AlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjRjkjkXbkuxJckXfvKOS7Ejy2e7fI9e2TEnSagbZQj8X2Lxk3ouAf6mq+wL/0j2XJM3QqoFeVRcDNy2ZfRrw5m76zcDTJlyXJGlIo17L5Z5VdT1AVV2f5B7LNUyyBdgCsHHjxhG7k6TxZaXrMq2xmsJ1n9b8oGhVba2q+aqan5ubW+vuJOl2a9RA/78kRwN0/+6ZXEmSpFGMGujvB87sps8E3jeZciRJoxrktMXtwMeAE5LsTvIc4FXA45N8Fnh891ySNEOrHhStqjOWWfS4CdciSRqD3xSVpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEWMFepLfS3JlkiuSbE9yp0kVJkkazsiBnuQY4HeB+ap6ALAOOH1ShUmShjPuLpf1wJ2TrAcOA64bvyRJ0ihGDvSq+jLwt8AXgeuBm6vqw0vbJdmSZCHJwuLi4uiVSrdLmeFDB5txdrkcCZwGHA/cC7hLkmctbVdVW6tqvqrm5+bmRq9UkrSicXa5/CLw+aparKrvARcAj5hMWZKkYY0T6F8ETklyWJIAjwN2TaYsSdKwxtmHvhM4H7gUuLx7r60TqkuSNKT147y4qs4Gzp5QLZKkMfhNUUlqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSI8YK9CRHJDk/yWeS7Ery8EkVJkkazlh3LAL+DvhQVT09yR2BwyZQkyRpBCMHepKfAH4e+DWAqvou8N3JlCVJGtY4W+j3ARaBNyV5EHAJ8PyquqW/UZItwBaAjRs3jtHdgSoz7Ltm2Ld0IFr597HW8Fcms4yCzjj70NcDJwP/UFUPBm4BXrS0UVVtrar5qpqfm5sboztJ0krGCfTdwO6q2tk9P59ewEuSZmDkQK+qG4AvJTmhm/U44KqJVCVJGtq4Z7n8DnBed4bL54BfH78kSdIoxgr0qroMmJ9QLZKkMfhNUUlqhIEuSY0w0CWpEQa6JDXCQJekRhjoktSIcc9Dl24X1vSqOStdBMTL9WgIbqFLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGjF2oCdZl+STSS6cREGSpNFMYgv9+cCuCbyPJGkMYwV6kmOBJwP/OJlyJEmjGvdaLucAfwQcvlyDJFuALQAbN24csztJB4WVrk+zlm7n174ZeQs9yanAnqq6ZKV2VbW1quaran5ubm7U7iRJqxhnl8sjgacmuRZ4B/DYJG+bSFWSpKGNHOhV9eKqOraqNgGnA/9aVc+aWGWSpKF4HrokNWIiN7ioqo8AH5nEe0mSRuMWuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0Y556i907yb0l2JbkyyfMnWZgkaTjj3OBiL/AHVXVpksOBS5LsqKqrJlSbJGkI49xT9PqqurSb/gawCzhmUoVJkoYzkX3oSTYBDwZ27mfZliQLSRYWFxcn0Z0kaT/GDvQkdwXeDbygqr6+dHlVba2q+aqan5ubG7c7SdIyxgr0JIfQC/PzquqCyZQkSRrFOGe5BHgjsKuqXj25kiRJoxhnC/2RwLOBxya5rHs8aUJ1SZKGNPJpi1X1n0AmWIskaQx+U1SSGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1Ihx7ym6OcnVSa5J8qJJFSVJGt449xRdB7wOeCJwf+CMJPefVGGSpOGMs4X+UOCaqvpcVX0XeAdw2mTKkiQNa+R7igLHAF/qe74beNjSRkm2AFu6p99McvWI/W0Ablxu4cxubpqV61rrzlcww7pWdFDW5fi6jQPz/3Gm62tlScap7bhBGo0T6Pv7367bzKjaCmwdo59eZ8lCVc2P+z6TZl3Dsa7hWNdwDtS6YDq1jbPLZTdw777nxwLXjVeOJGlU4wT6fwP3TXJ8kjsCpwPvn0xZkqRhjbzLpar2JjkL+GdgHbCtqq6cWGW3NfZumzViXcOxruFY13AO1LpgCrWl6ja7vSVJByG/KSpJjTDQJakRMw/0JNuS7ElyxTLLk+Tvu8sLfDrJyX3Lzkzy2e5x5pTremZXz6eTfDTJg/qWXZvk8iSXJVmYcl2PTnJz1/dlSV7Wt2zNLtUwQF0v7KvpiiTfT3JUt2wt19e9k/xbkl1Jrkzy/P20mfoYG7CuqY+xAeua+hgbsK6pj7Ekd0ryiSSf6ur68/20OTTJO7t1sjPJpr5lL+7mX53kCWMXVFUzfQA/D5wMXLHM8icBF9E77/0UYGc3/yjgc92/R3bTR06xrkfs64/e5Q929i27Ftgwo/X1aODC/cxfB/wvcB/gjsCngPtPq64lbZ8C/OuU1tfRwMnd9OHA/yz9uWcxxgasa+pjbMC6pj7GBqlrFmOsGzN37aYPAXYCpyxp89vAG7rp04F3dtP379bRocDx3bpbN049M99Cr6qLgZtWaHIa8Jbq+ThwRJKjgScAO6rqpqr6KrAD2Dytuqrqo12/AB+ndx7+mhtgfS1nTS/VMGRdZwDbJ9X3Sqrq+qq6tJv+BrCL3rec+019jA1S1yzG2IDrazlrNsZGqGsqY6wbM9/snh7SPZaeaXIa8OZu+nzgcUnSzX9HVd1aVZ8HrqG3Dkc280AfwP4uMXDMCvNn4Tn0tvD2KeDDSS5J79IH0/bw7iPgRUlO7OYdEOsryWH0QvHdfbOnsr66j7oPprcV1W+mY2yFuvpNfYytUtfMxthq62vaYyzJuiSXAXvobQAsO76qai9wM3B31mB9jfPV/2lZ7hIDA116YK0leQy9X7ZH9c1+ZFVdl+QewI4kn+m2YKfhUuC4qvpmkicB7wXuywGyvuh9FP6vqurfml/z9ZXkrvR+wV9QVV9fung/L5nKGFulrn1tpj7GVqlrZmNskPXFlMdYVX0fOCnJEcB7kjygqvqPJU1tfB0MW+jLXWJg5pceSPJA4B+B06rqK/vmV9V13b97gPcw5seoYVTV1/d9BKyqDwKHpHdRoJmvr87pLPkovNbrK8kh9ELgvKq6YD9NZjLGBqhrJmNstbpmNcYGWV+dqY+x7r2/BnyE2+6W++F6SbIeuBu93ZOTX1+TPEAw6gPYxPIH+Z7Mjx+w+kQ3/yjg8/QOVh3ZTR81xbo20tvn9Ygl8+8CHN43/VFg8xTr+kl+9IWxhwJf7NbdenoH9Y7nRwesTpxWXd3yfQP5LtNaX93P/hbgnBXaTH2MDVjX1MfYgHVNfYwNUtcsxhgwBxzRTd8Z+A/g1CVtnsePHxR9Vzd9Ij9+UPRzjHlQdOa7XJJsp3fUfEOS3cDZ9A4sUFVvAD5I7yyEa4BvAb/eLbspyV/Qu6YMwMvrxz9irXVdL6O3H+z1veMb7K3eldTuSe9jF/QG+Nur6kNTrOvpwHOT7AW+DZxevdGzppdqGKAugF8GPlxVt/S9dE3XF/BI4NnA5d1+ToA/oReWsxxjg9Q1izE2SF2zGGOD1AXTH2NHA29O74Y/d6AX1hcmeTmwUFXvB94IvDXJNfT+2Jze1XxlkncBVwF7gedVb/fNyPzqvyQ14mDYhy5JGoCBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhrx/zuDmRsJs3RHAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "group=pd.DataFrame(res.iloc[0,:]).reset_index().drop('index',axis=1).rename({0:'group'},axis=1)\n",
    "\n",
    "colors=['red','yellow']\n",
    "for i in range(n_cov):\n",
    "    cov=pd.DataFrame(bins[:,i]).rename({0:'bins'},axis=1)\n",
    "    df=pd.concat([cov, group], axis=1)\n",
    "    plt.title('Overall data binning')\n",
    "    plt.hist(list(df['bins']),density=False,histtype='bar',color='black')\n",
    "    plt.title('Covariate '+str(i+1))\n",
    "    plt.hist([list(df[df['group']==0]['bins']),list(df[df['group']==1]['bins'])],density=False,histtype='bar',color=colors,bins=no_bins)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate means difference between groups"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 173,
   "metadata": {},
   "outputs": [],
   "source": [
    "##########Variance, mean (Compare with randomization technique)##################################\n",
    "dcov=dict()\n",
    "group=pd.DataFrame(res.iloc[0,:]).reset_index().drop('index',axis=1).rename({0:'group'},axis=1)\n",
    "for i in range(n_cov):\n",
    "    cov=pd.DataFrame(mat[:,i]).rename({0:'bins'},axis=1)\n",
    "    df=pd.concat([cov, group], axis=1)\n",
    "    dcov[i]={'mean':df['bins'].mean(),\n",
    "            'variance':df['bins'].var(),\n",
    "            'mean_control':df[df['group']==0]['bins'].mean(),\n",
    "            'variance_control':df[df['group']==0]['bins'].var(),\n",
    "            'mean_treatment':df[df['group']==1]['bins'].mean(),\n",
    "            'variance_treatment':df[df['group']==1]['bins'].var()}\n",
    "#print(pd.DataFrame.from_dict(dcov))\n",
    "#randomization\n",
    "group=pd.DataFrame(np.random.randint(0,2,(samples,1))).rename({0:'group'},axis=1)\n",
    "dcovrand=dict()\n",
    "for i in range(n_cov):\n",
    "    cov=pd.DataFrame(mat[:,i]).rename({0:'bins'},axis=1)\n",
    "    df=pd.concat([cov, group], axis=1)\n",
    "    dcovrand[i]={'mean':df['bins'].mean(),\n",
    "            'variance':df['bins'].var(),\n",
    "            'mean_control':df[df['group']==0]['bins'].mean(),\n",
    "            'variance_control':df[df['group']==0]['bins'].var(),\n",
    "            'mean_treatment':df[df['group']==1]['bins'].mean(),\n",
    "            'variance_treatment':df[df['group']==1]['bins'].var()}\n",
    "#print(pd.DataFrame.from_dict(dcovrand))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 176,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "different of means-before optimization    1.064594\n",
      "different of means-optimization           0.673086\n",
      "dtype: float64\n"
     ]
    }
   ],
   "source": [
    "df=pd.DataFrame.from_dict(dcov).transpose()\n",
    "df['different of means-optimization']=(df['mean_treatment']-df['mean_control']).abs()\n",
    "df_rand=pd.DataFrame.from_dict(dcovrand).transpose()\n",
    "df_rand['different of means-before optimization']=(df_rand['mean_treatment']-df_rand['mean_control']).abs()\n",
    "result = pd.concat([df_rand['different of means-before optimization'],df['different of means-optimization']], axis=1, sort=False)\n",
    "\n",
    "print(result.mean())"
   ]
  }
 ],
 "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}