ms_claa_x2pte.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np, analysis as ana, plottools as pl, prjlib, binning as bins\n",
    "from matplotlib.pyplot import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "mb = bins.multipole_binning(10,spc='p2',lmin=20,lmax=2048)\n",
    "p0 = prjlib.params_init(stype='lcmb',PSA='s14&15_deep56',doreal='True',dearot='True',rlmin='200',ver='v1')\n",
    "p1 = prjlib.params_init(stype='lcmb',PSA='s14&15_deep56',doreal='True',dearot='False',rlmin='200',ver='v1')\n",
    "scb = prjlib.binned_cl_rlz(p0.quad.f['EB'].cl,1,200,mb)\n",
    "ocb0 = prjlib.binned_cl(p0.quad.f['EB'].cl[0],mb)\n",
    "ocb1 = prjlib.binned_cl(p1.quad.f['EB'].cl[0],mb)\n",
    "pa, __ = prjlib.filename_init(rlmin='200',stype='absrot')\n",
    "acb = prjlib.binned_cl_rlz(pa.quad.f['EB'].cl,1,100,mb)\n",
    "dcb = np.mean(acb-scb[:100,:],axis=0)\n",
    "#ocb = prjlib.binned_cl(p.quad.f['EB'].cl[0].replace('dearot','dearot_old'),mb0,mb1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "chi^2: 2.8 , chi^2 (sim): 11.2 , PTE:0.99 (100 rlz)\n",
      "chi^2: 2.8 , chi^2 (sim): 10.5 , PTE:0.99 (200 rlz)\n",
      "chi^2: 2.7 , chi^2 (sim): 10.5 , PTE:0.99 (no derotation)\n",
      "chi^2: 1.5 , chi^2 (sim): 10.5 , PTE:1.0 (diagonal cov)\n",
      "chi^2: 2.9 , chi^2 (sim): 11.2 , PTE:0.99 (absrot error)\n"
     ]
    }
   ],
   "source": [
    "# x2 PTE\n",
    "st = ana.PTEs(ocb0,scb[:100,:],comment='100 rlz')\n",
    "st = ana.PTEs(ocb0,scb,comment='200 rlz')\n",
    "st = ana.PTEs(ocb1,scb,comment='no derotation')\n",
    "st = ana.PTEs(ocb0,scb,diag=True,comment='diagonal cov')\n",
    "st = ana.PTEs(ocb0-dcb*(.06/.08)**2,scb[100:,:],comment='absrot error')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAD8CAYAAACVZ8iyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3X28HVV97/HPN4eQYABJCGKECFi5V21tA0ak1esDKuDDBaxPcFuMijcvXy+p2voE19cVi9KL1op6621NNYoPFSlqSS2WgkDtgygnGsNTNRi0JEUCJMhzIOd87x8zRzeHvc+ec/bsSTLn+85rXmfPmlmz1uycs397Zs1aS7aJiIjZbc7OrkBEROx8CQYREZFgEBERCQYREUGCQUREkGAQEREkGEREVCJptaQtkq7rsV2SPiHpJknrJR3ZsW2FpA3lsqK5WleXYBARUc3ngOOn2P4S4PByWQn8BYCkRcBZwLOAo4CzJC0cak1noHIwkDQi6QeSvlGuf07SzZLWlcuy4VUzImLnsv1tYOsUu5wIfN6Fq4H9JC0BjgMus73V9jbgMqYOKjvFHtPY923AjcC+HWnvsn1R1QOMLFjguQsXTaPImbGGXsSvNFnWnOZ6i8/ZY7yRcvYYaaYcgMeMPNRYWXPnjDVXlpora/85Oxora+367XfYPmCQYxz3ggW+c2v/92ft+u3XAw92JK2yvWqaxR0E3NKxvqlM65W+S6kUDCQdDLwMOAf4o5kWNnfhIpa+ZcbZKxvbq7kPTY80V9b4gub+6B+z8IFGyjlg33sbKQfgGfvf0n+nmjxuz7sbK+sJc+9qrKxT97mjsbJGlmz42aDHuHPrGN+79IlVynrQ9vJBy9udVb1N9DHg3cDkr3HnlA0l50ma1y2jpJWSRiWNjt133yB1jYiYFgPjFf7VZDOwtGP94DKtV/oupW8wkPRyYIvttZM2nQk8BXgmsAh4T7f8tlfZXm57+ciCBYPWNyKiMmMe9ljfpSZrgNeVTxUdDfzC9q3ApcCxkhaWDcfHlmm7lCq3iZ4NnCDppcB8YF9JX7T9++X27ZI+C7xzWJWMiJipur75S/oy8HxgsaRNFE8IzQWw/ZfAJcBLgZuA+4E3lNu2SvoAcE15qLNtT9UQvVP0DQa2z6S4CkDS84F32v59SUts3ypJwElA12dvIyJ2FmPGahqm3/YpfbYbeEuPbauB1bVUZEim8zTRZF+SdADF8zTrgDfXU6WIiPqMkzlbqphWMLB9FXBV+fqYIdQnIqI2BsYSDCoZ5MogImKXlyuDahIMIqK1DDycqX0rSTCIiNYyzm2iihIMIqK9DGOJBZUkGEREaxU9kKOKBIOIaDEx1uhokruvBIOIaK2iATnBoIoEg4horaKfQYJBFQkGEdFq47kyqCTBICJaK1cG1SUYRERrGTGWqd4rSTCIiFbLbaJqEgwiorWMeMgjO7sau4UEg4horaLTWW4TVZFgEBGtlgbkahoNBhaM7TX8gULG5zXXAd0NljVv3+2NlTU23sy3qYXzHmikHIC9R5p7/0YaHBzt1H3uaKysL9yzuLGyYMPAR7DFmHNlUEXepYhotXHUd6lC0vGSfiTpJklndNl+nqR15fJjSXd1bBvr2LamxtOrTeUrA0kjwCiw2fbLJR0GXADsD6wFTrX90HCqGRExfUUD8uA3QMrPv08CLwY2AddIWmP7hl+WZf9hx/5/ABzRcYgHbC8buCJDNJ0rg7cBN3asfwg4z/aTgW3AaXVWLCJiUBMNyP2WCo4CbrK9sfzSewFw4hT7nwJ8efAzaE6ld0HSwcDLgE+X6wKOAS4qdzkfOGkYFYyIGMSY1Xep4CDglo71TWXao0g6BDgMuKIjeb6kUUlXS9olPyurXj99DHg3sE+5vj9wl+0d5fpUb8xKYCXAyMKFM69pRMQ0TaMH8mJJox3rq2yvmmGxJwMX2R7rSDvE9mZJTwKukHSt7Z/M8PhD0TcYSHo5sMX2WknPn24B5Ru6CmDe0qWZcygiGjVe7WmiO2wvn2L7ZmBpx/rBZVo3JwNv6Uywvbn8uVHSVRTtCbtXMACeDZwg6aXAfGBf4OPAfpL2KK8OpnpjIiJ2imKguloemrwGOLx8cGYzxQf+/5i8k6SnAAuB73SkLQTut71d0mKKz9QP11GpOvV9l2yfaftg24dSvAFX2P494ErgVeVuK4CLh1bLiIgZMOJhj/Rd+h6n+NJ7OnApxYM0F9q+XtLZkk7o2PVk4ALbnXdBngqMSvohxefmuZ1PIe0qBnnm6j3ABZI+CPwA+Ew9VYqIqIdNbZ3ObF8CXDIp7X2T1t/fJd+/AU+vpRJDNK1gYPsq4Kry9UaKx60iInZR1TuVzXYZmygiWsvUd2XQdgkGEdFqmdymmgSDiGgto0xuU1GCQUS0loGHaxibaDbIuxQRLabMZ1BRgkFEtJap3AN51kswiIhWy5VBNQkGEdFatnJlUFGCQUS0VtGA3H+4iUgwiIhWyxzIVTUeDDTWf59BjSxqbuLzefMebqwsN/i8tNTMaOP7zn2wkXKa9s5FzY1O/JGtv9ZYWbvb/feiAXn3qvPOkiuDiGi19ECuJsEgIlorPZCrSzCIiFarOOH9rJdgEBGtZcPD4wkGVSQYRERrFbeJEgyqSDCIiFbb3Z6A2lkSDCKitfJoaXV9r58kzZf0PUk/lHS9pD8u0z8n6WZJ68pl2fCrGxExHcVton5LpSNJx0v6kaSbJJ3RZfvrJd3e8Zn4po5tKyRtKJcVNZ5gbapcGWwHjrF9r6S5wL9I+ma57V22Lxpe9SIiBlPHHMiSRoBPAi8GNgHXSFpj+4ZJu37F9umT8i4CzgKWU1ysrC3zbhu4YjXqGxJduLdcnVsuzXRPjYgYQPE00UjfpYKjgJtsb7T9EHABcGLFahwHXGZ7axkALgOOn9EJDVGl6yNJI5LWAVsoTuq75aZzJK2XdJ6keT3yrpQ0Kml0/L77aqp2RER/E53O+i3A4onPqXJZOelQBwG3dKxvKtMme2X5mXiRpKXTzLtTVQoGtsdsLwMOBo6S9BvAmcBTgGcCi4D39Mi7yvZy28vnLFhQU7UjIqoZR30X4I6Jz6lyWTWDov4OONT2b1J8+z+/zvMYtmk9gGv7LuBK4Hjbt5a3kLYDn6W4jIqI2GVMPE1U4cqgn83A0o71g8u0X5Vl31l+HgJ8GnhG1by7gipPEx0gab/y9V4UDSj/LmlJmSbgJOC6YVY0ImImanqa6BrgcEmHSdoTOBlY07nDxGdi6QTgxvL1pcCxkhZKWggcW6btUqo8TbQEOL9sTZ8DXGj7G5KukHQAIGAd8OYh1jMiYtpssaOGHsi2d0g6neJDfARYbft6SWcDo7bXAG+VdAKwA9gKvL7Mu1XSBygCCsDZtrcOXKma9Q0GttcDR3RJP2YoNYqIqFFdnc5sXwJcMintfR2vz6RoS+2WdzWwupaKDEl6IEdEa83WHsiSngZMfGG/okt/iEfJCE4R0Wo1NSDvNiSdCnwNOLBcvibpdf3y5cogIlprlk5u827gdybaJSR9FLgK+PxUmRIMIqLV6hiOYjezo7OB2vY2SeP9MiUYRERr2bBj9k1u8wNJCyfGPiq7Bqzvl6nRYCDDnIeHH6XHb9lr6GVMeGBJc79oCxc2N5zHvQ90HV2kdveP7dlIOQB/fMD1jZV11u2/3lhZTbp3rJnfizrNtttEtt84af0uoO9IqbMuZEbE7DGNsYlaQ9ITJf2tpC3lkNprJB3SL1+CQUS0mq2+S8t8FriIosPw44ELqdDHIcEgIlqt4kB1bbK/7S+WA4yO2f4isH+/TAkGEdFa9uzrZwDcXs66tke5vAG4vV+mBIOIaDExNj6n79IybwD+O/Cf5XJCmTalPFoaEa3WwjaBKdneBLxyuvkSDCKitWbj2ESSVsOjG0JsT3l1kGAQEe3lot1glvlGx+sFwKupMJlOgkFEtFoLnxaaku2vTUr6gqTv9MuXYBARreWyATn4V0kjtsd67ZBgEBGtNttuE0k6q0vyPVMFAkgwiIiWm21PEwH3dLxewCPnY+6pbzCQNB/4NjCv3P8i22dJOgy4gKJn21rgVNsPzaDiERFDYdcXDCQdD3ycYg7kT9s+d9L2PwLeRDEH8u3AG23/rNw2Blxb7voftk+opVJd2P7opHp9CLiiX74qN9O2A8fY/i1gGXC8pKOBDwHn2X4ysA04bdq1jogYsjp6IEsaAT4JvAR4GnBKObVkpx8Ay23/JsXYQB/u2PaA7WXlMrRA0E35Jf1n5Tn01DcYuHBvuTq3XEwxv+ZFZfr5wEkzr25ExHDY/ZcKjgJusr2x/HC9ADjxkeX4Stv3l6tXAwfXeR5VSVoyOc327/VrM6jUzC5pRNI6YAtwGfAT4C7bO8pdNgEH9ci7UtKopNGx+5objz8iwojx8Tl9F2DxxOdUuaycdKiDgFs61nt+5pVOA77ZsT6/PO7Vkob9xfmvZpKpUgNyGVGWlTPmfB14StUCbK8CVgHMP3jpLGvXj4idreKHzh22l9dRnqTfB5YDz+tIPsT2ZklPAq6QdK3tn9RR3mS2Xz6TfNN6ALecMedK4LeB/SRNBJODqdDDLSKiUa5tPoPNwNKO9a6feZJeBLwXOMH29l9Ww95c/txIMTn9ETM/qeHoGwwkHVBeESBpL+DFFI8pXQm8qtxtBXDxsCoZETFjrrD0dw1wuKTDJO0JnAys6dxB0hHApygCwZaO9IWS5pWvFwPPBm4Y7KTqV+U20RLg/LIleg5woe1vSLoBuEDSByla0T8zxHpGRMxIHY+W2t4h6XTgUopHS1fbvl7S2cCo7TXAnwJ7A38jCX71COlTgU9JGqf4DD3X9u4XDGyvp8slTXm5c9QwKhURUQcD4+P19DOwfQlwyaS093W8flGPfP8GPL2WSlRUXoEcXa5ebfuOfnka74HcRGfA8cdv779TTRYvuqf/TjX5xf17NVbWXvOa6T/4N0+6vJFyAF69sevf6lA8ZqS5/pd3Pzy/sbK2bW/ud7AWppkPnV2IpOMoHtr5Z4p34P9JWmn7H6bKl+EoIqLVZtvYRMD/AZ7b0fv5iRRPgSYYRMQsNvuCgSYCAYDt/5DU92GhBIOIaLHKj462yRZJC21vAyifBt312gwiIho1y64MbB83af0uii4BU8qsDxHRXgaPq+/SJmW/ho9JWivp+5I+IWlhv3wJBhHRcqqwtMpqittCrwR+l2I47dX9MuU2UUS02yy7TQQ8yfYrOtY/IOmH/TLlyiAi2q2e4Sh2J/dJev7EiqQXAH2HjM6VQUS01yzsdAaspBhCaHG5vhV4Xb9MCQYR0WqzrdOZ7euAZ0jau1y/t08WIMEgItquZU8L9VPOW/+/geMAS/pH4JyOWdi6SptBRLSa3H9pmY8B+wKnAPOA64H/2y9Trgwior3a2UDcz+/Y/k0ASWO2/1rSW/tlSjCIiBbTbGxAfsQJS3osMLdfptwmioh2m32Plt4saVn5ej/ge8D7+2XKlUFEtNv4zq5As8rZ1SYcTzHj2gP98lWZA3mppCsl3SDpeklvK9PfL2mzpHXl8tKZVz8iYggm+hn0W1pE0iETC/Ag8DhJP5T0+KnyVbky2AG8w/b3Je0DrJV0WbntPNsfGazqERHDU9fTQpKOBz5OMQfyp22fO2n7PODzwDOAO4HX2v5pue1M4DRgDHir7UvrqVVXf0fRbtB55r8G/LOkUdundMtUZQ7kW4Fby9f3SLoROGjw+kZENKCGYCBpBPgkxVDQm4BrJK2ZNLH9acA220+WdDLwIeC1kp4GnAz8OvAE4HJJ/8X22OA1e7SJJ4km1f/7to+UtK5Xvmk1IEs6FDgC+G6ZdLqk9ZJW9xoiVdJKSaOSRsfu6zs8RkTErugo4CbbG20/BFwAnDhpnxOB88vXFwEvlKQy/QLb223fDNxUHq9Jf1v+/GyvHSo3IJddm78KvN323ZL+AvgARdz9APBnwBsn57O9imJyZuYfvLR97fYRsUureJtosaTRjvVV5WfXhIOAWzrWNwHPmnSMX+5je4ekXwD7l+lXT8o7tLsrks7qkvxmSdv5VVB4lErBQNJcikDwJdtfA7B9W8f2vwK+0e84HoEdjx3KldEj6K49h17GhC0P7ddYWQc+4a7Gyrp62UWNlHPE6GsbKQfgge3N/V7sMTL83/MJO8ZGGitrZM5u9miOqTocxR22lw+5Nk25p0vaDoq2jgso7u48St9gUF7mfAa40fZHO9KXlO0JAK8ArptujSMihq6e+xGbgaUd6weXad322SRpD+CxFA3JVfLWpvNzeoKkl9n+E0mv7JWvSpvBs4FTgWMmPUb6YUnXSloPvAD4w5lWPiJiWGoam+ga4HBJh0nak6JBeM2kfdYAK8rXrwKusO0y/WRJ8yQdBhxO0RGsMbZf2FGvrqo8TfQvdJ8X7pIZ1isiojk1XBmUbQCnA5dS3G5Zbft6SWcDo7bXUNxB+YKkmyjmEDi5zHu9pAuBGyhu17xlWE8SAUj6e9sv63EeN/fKlx7IEdFuNT22YvsSJn0Jtv2+jtcPAq/ukfcc4Jx6atLXaTPJlGAQEa3V0iGqp2T75zPJl2AQEe02yya3makEg4hotdl2ZTBTCQYR0W4JBpUkGEREe83CNgNJzwFeBjwV2ItibLl/Bb5qe2uvfAkGEdFusyQYSHomxbBAPwOuBL5JMYT1ARQjqX5D0ndtd+0TlmAQEa2m3WwEjQE8ALzc9t1dtv09cHZ51dBVgkFERAvY7jskUNmJuKvMgRwR7TbL5kCW9Ac90n9L0lt65UswiIj2qjAuUQsbmN/dI30b8PZemXKbKCLarX0f9v0cKOlmijPv7HFn4JBemRIMIqLdZl8w2AIsZ5pnnmAQEa0lZtXTRBMusX3ndDMlGEREe7WzTWBKtlfOJF+CQUS02ywLBgDlBGQn8Ku5ljcDF9v+Zq88CQYR0W6zLBhI+nPgMODzwC1l8lLgdEkvsf3WbvkaDQZ6GOb/fPiTdz/4hB1DL2PCHo9prqymJqkHOHpdz9nxatXkJPV777W9sbIe2D63sbKa/LC7f9tezRVWk9l2mwg43vaTu6R/RdKGXpn69jOQtFTSlZJukHS9pLeV6YskXSZpQ/lz4QCVj4gYjlnW6Qy4T9LRkxMlPQu4r1emKp3OdgDvsP004GjgLZKeBpwBfMv24cC3yvWIiF2Hi6eJ+i2DqvLlWNIySd8pv1Svl/Tajm2fk3SzpHXlsmyA6rwB+Jikf5d0eVmfG4GPA6/vlanvbSLbt1IMgYrte8qDHgScCDy/3O184CrgPQOcQERE/Zr55j/x5fhcSWeU65M/D+8HXmd7g6QnAGslXWr7rnL7u2wPfC/Y9veBoyUdSEcDsu3bpso3rTYDSYcCRwDfBQ4sAwXAz4EDe+RZCawE2OOxuZMUEc1qqM2g75dj2z/ueP2fkrZQDC99FzUrr0xeQEcwKAPPtl55Ko9NJGlv4KvA2ycPkWq7550326tsL7e9fOQxC6oWFxFRj2ptBosljXYs031Wv9KX4wmSjgL2BH7SkXxOefvoPEnzpll+57FfD4xSBKd55fI8YFTSil75Kl0ZSJpLEQi+ZPtrZfJtkpbYvlXSEoou0BERu47qDcR32F4+1Q6SLgce32XTex9RpG2p9/VI+Xn5BWCF7YkWizMpgsiewCqKq4qzK9W8e32OtP2LSeU+FriG4srlUfoGA0kCPgPcaPujHZvWACuAc8ufF8+s3hERwyHqu01k+0U9y5EqfTmWtC/FRDPvtX11x7Enriq2S/os8M5Bqkr3uz5iitBY5crg2cCpwLWS1pVp/4siCFwo6TSKadZeM63qRkQ0oKE2g75fjiXtCXwd+PzkhuKOQCLgJKDvRDVT+CDFLaHLgE1l2kHAscAf98pU5Wmif+GRw6B2euE0KxkR0axmgkHXL8eSlgNvtv2mMu25wP7lfX2A19teB3xJ0gEUn7XrgDfPtCK2Py9pDXAcv2pAvgo4Y/Kto04ZjiIi2q2BYFCOEvqoL8e2R4E3la+/CHyxR/5jaq7PXcBXppMnwSAi2msWjlo6UwkGEdFuCQaVJBhERKvNwsltZiTBICJabbbdJpK0kS4P/dg+bKp8CQYR0V7tHJW0n87OcwuA1wKL+2VKMIiIdptlwcD21o7VrcBHJI32y5dgEBGtVWcP5N3chyXN6Rj+4lESDCKi1TQ+u6KBpNV07yj8EknYfkO3fAkGEdFes7PN4BszyZRgEBGtNgtvE+0NXNw59ISkZ9heO1WmyvMZRETslmbfHMifAP5Z0m90pH26X6ZGrwy8Bzx4wPB7gIw89qGhlzFhw/O6Dg0+FIdf9frGytpxfzO/Go9bUvskTz3dsXWfxsraa8H2xsrafveM50GZtjn3jTRWVl1m4ZXBRuB/Al+RdLbtr1Ah5OU2UUS02+wLBrK9VtJ/owgIz6DCZ31uE0VEe7kYjqLf0jKGX/Y3OBYYB369X6YEg4horYl+Bv2WNrF9ZMdr2z7Ddt/7ewkGEdFudv+lBSS9UtJJkrreEpL0REl/2it/2gwiotXa9s1/CpcB7wDOkfQT4MfAg8DjgCMohqY4t1fmvlcGklZL2iLpuo6090vaLGldubx00LOIiKhdlcdKWxIsbN9t+yzg6RTzIN8A3Ap8E3iZ7eNsX9krf5Urg88Bfw58flL6ebY/MqNaR0Q0pIkGYkmLKKaZPBT4KfAa29u67DcGXFuu/oftE8r0w4ALgP2BtcCptmf0jHw5/tD3yqWyvlcGtr9NcXkREbHbaehpojOAb9k+HPhWud7NA7aXlcsJHekfoviC/WRgG3BaLbWahkEakE+XtL68jbSw106SVkoalTQ6du99AxQXETFNpqkG5BOBiR6o5wMnVc0oScAxwEUzyV+XmQaDvwB+DVhGcU/qz3rtaHuV7eW2l4/svWCGxUVEzEzFR0sXT3xpLZeV0yzmQNu3lq9/DhzYY7/55fGvljTxgb8/cJftHeX6JuCgaZY/sBk9TWT7tonXkv6KGY6SFxExdNW++N9he/lUO0i6HHh8l03vfURxtqWezzAdYnuzpCcBV0i6FvhFj30bNaNgIGlJRxR8BXDdVPtHROwMdU5uY/tFPcuRbpv4XJS0BNjS4xiby58bJV1F8cjnV4H9JO1RXh0cDGyup9bVVXm09MvAd4D/KmmTpNMoZs25VtJ64AXAHw65nhER02ej8f5LDdYAK8rXK4CLJ+8gaaGkeeXrxcCzgRtsG7gSeNVU+Yet75WB7VO6JH9mCHWJiKhfM/0IzgUuLL8s/wx4DYCk5cCbbb8JeCrwKUnjFF/Ez7V9Q5n/PcAFkj4I/ICd8BmbHsgR0WpN9EC2fSfwwi7po8Cbytf/RtEhrFv+jcBRw6xjPwkGEdFeBmbZHMgzlWAQEe2WWFBJgkFEtNosGqhuIAkGEdFqNT0t1HoJBhHRXi0alXTYEgwiorWKTmeJBlU0GwwMIw9q6MVseN75/XeqyeH/tKL/TjUZu3tuY2XNeaiZSfDuuHOfRsoBmPPzeY2VtX1pcxPrantzExZqbPh/v7Vr3xzHQ5Erg4hotVwZVJNgEBHtlTaDyhIMIqLFaht7qPUSDCKi3XKbqJIEg4hoLzczB3IbJBhERLvlyqCSBIOIaLfEgkoSDCKi1TSe+0RVJBhERHuZdDqrqMq0l6slbZF0XUfaIkmXSdpQ/lw43GpGREyfMHL/JSoEA+BzwPGT0s4AvmX7cOBb5XpExK7H7r9E/2Bg+9vA1knJJwITAwCdD5xUc70iIurRQDCocrdE0gskretYHpR0Urntc5Ju7ti2bOBKTdNMR7g60Pat5eufAwf22lHSSkmjkkbH77tvhsVFRMzARJtBv2Vwfe+W2L7S9jLby4BjgPuBf+zY5V0T222vq6VW0zDwcIe2pxz9w/Yq28ttL5+zYMGgxUVETIvGx/suNZju3ZJXAd+0fX8dhddhpsHgNklLAMqfW+qrUkREXSrcIqqnzaDy3ZLSycCXJ6WdI2m9pPMkNTfeemmmwWANMDGQ/wrg4nqqExFRI1M1GCyeuJ1dLisnH0rS5ZKu67Kc+Igi+9wtKb9APx24tCP5TOApwDOBRcB7Bj73aerbz0DSl4HnU7xZm4CzgHOBCyWdBvwMeM0wKxkRMWPV7gLdYXv5VDvYflGvbZJuk7TE9q0V7pa8Bvi67Yc7jj1xVbFd0meBd1aqdY36BgPbp/TY9MKa6xIRUbuG+hFM3C05l/53S06huBL4pY5AIor2huu65hyi5ubLi4jYGZppMzgXeLGkDcCLynUkLZf06YmdJB0KLAX+aVL+L0m6FrgWWAx8sI5KTUeGo4iI9rJhbPjjUdi+ky53S2yPAm/qWP8pcFCX/Y4ZZv2qSDCIiHZLD+NKGg0GT190O9/7vb8cejmHf+nNQy9jwtj85n7R5t/e3F29Bx8/1kg5e2ye30g5AGrwM2F8016NlTVnpLGiGHlAzRVWlwSDSnJlEBHtZSBzIFeSYBARLWZwxrCuIsEgItrLNNKA3AYJBhHRbmkzqCTBICLaLcGgkgSDiGixTF5TVYJBRLSXgXqGqG69BIOIaLdcGVSSYBARLdbMcBRtkGAQEe1lcPoZVJJgEBHtlh7IlSQYRES7pc2gkgSDiGgvO08TVZRgEBHtliuDSgYKBpJ+CtwDjAE7+s0hGhHRLOOxZoZj393VcWXwAtt31HCciIh6ZQjrynKbKCLaLY+WVjLo1FkG/lHSWkkru+0gaaWkUUmjt9+Zy7WIaI4Bj7vvMihJr5Z0vaRxST1vl0s6XtKPJN0k6YyO9MMkfbdM/4qkPQeu1DQNGgyeY/tI4CXAWyQ9d/IOtlfZXm57+QH7Nzg/X0SEy8lt+i2Duw74XeDbvXaQNAJ8kuLz8mnAKZKeVm7+EHCe7ScD24DT6qjUdAwUDGxvLn9uAb4OHFVHpSIi6uKxsb7LwGXYN9r+UZ/djgJusr3R9kPABcCJkgQcA1xU7nc+cNLAlZqmGbcZSFoAzLF9T/n6WODsqfKsXb/93pElG/q9Yf0sBvo0WL+jxmPVpsmymtTW85rKbDznJkx+Xw8Z9IAN1XF5AAADy0lEQVT3sO3Sy33R4gq7zpc02rG+yvaqQcuf5CDglo71TcCzgP2Bu2zv6Eg/qOay+xqkAflA4OtFUGMP4K9t/0OfPD8a9PFTSaN1PcJa57F2pbKa1NbzmspsPOcmDON9tX18XceSdDnw+C6b3mv74rrK2VlmHAxsbwR+q8a6RETssmy/aMBDbAaWdqwfXKbdCewnaY/y6mAivVGDNiBHREQ11wCHl08O7QmcDKyxbeBK4FXlfiuAxq80mg4GddyDq/M+Xt33BHeVsprU1vOaymw85ybstu+rpFdI2gT8NvD3ki4t058g6RKA8lv/6cClwI3AhbavLw/xHuCPJN1E0YbwmcbPwRm3IyJi1sttooiISDCIiIghBQNJSyVdKemGsov228r0RZIuk7Sh/LmwwrF+KulaSesmngOuehxJqyVtkXRdR1rXvCp8ouwOvl7SkcM850HLa9J0/g92p/PqVNfviqQV5f4bJK3YGeeyq6jzbyLvawNs174AS4Ajy9f7AD+m6H79YeCMMv0M4EMVjvVTYPGktErHAZ4LHAlc1y8v8FLgm4CAo4HvDvOcBy2vyWU6/we703nV/bsCLAI2lj8Xlq8X7uxz24nvaS1/E3lfG/r/auiX4mLgxcCPgCUdvyg/qpC32wdR5eMAh076A++aF/gUcEq3/YZxznWXN+T/v8r/B7vTedX9uwKcAnyqI/0R+832ZaZ/E3lfm1mG3mYg6VDgCOC7wIG2by03/ZyiF3M/3UZGnclxJvTK262r+Iy6hFc859rKa8B0/g92p/PqZ7rn2KZzr9WAfxN5Xxsw1PkMJO0NfBV4u+27VQxdAYBtS6ryXOtzbG+W9DjgMkn/3rlxGsd5lEHy9lLTOe9qhvZ/sLuYDec4LC39m2idoV0ZSJpL8QvwJdtfK5Nvk7Sk3L4E2NLvOO4+Muq0j9OhV95eXcUrm+Y5D1xeU6b5f7DbnFcF0z3HNp17LWr6m8j72oBhPU0kih50N9r+aMemNRRdraFCl2tJCyTtM/GaYmTU66Z7nEl65V0DvK58ouFo4Bcdl7J9zeCcByqvKTP4P9gtzqui6Z7jpcCxkhaWT8gcW6bNSjX+TeR9bcIwGiKA51DcZ14PrCuXl1J0s/4WsAG4HFjU5zhPAn5YLtdTjA5I1eMAXwZuBR6muM94Wq+8FE8wfBL4CXAtsHyY5zxoeU0t0/0/2F3Oa1i/K8AbgZvK5Q07+7x28nta299E3tfhLxmOIiIi0gM5IiISDCIiggSDiIggwSAiIkgwiIgIEgwiIoIEg4iIAP4/ZJ5LIMKPeyUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAD8CAYAAACVZ8iyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xu8XWV95/HPNycXaEBJCGKECCjMKK1txIhMdawiQrwMYLUKUzUqTMZXpWqtFxhfIxZlBm0ramtbU43gpUYGtaQWpSBQ2yrIicZwqwYjalIkQkC5SC7nfOePtY5uDvucvc45a68ka3/fea3X2evyrOdZ++Ts317ruck2EREx2Gbt7gJERMTul2AQEREJBhERkWAQEREkGEREBAkGERFBgkFERCWSVkvaKummCfZL0ocl3SZpg6RjOvatkLSxXFY0V+rqEgwiIqq5CFg+yf4XAEeVy0rgrwEkLQTOBZ4BHAucK2lBX0s6DZWDgaQhSd+W9KVy/SJJP5C0vlyW9q+YERG7l+2vAdsmOeQU4JMuXAccIGkxcBJwpe1ttu8BrmTyoLJbzJ7CsW8CbgUe1bHtbbYvrXqCof3me/bChVPIcpqG0qt6pjSrmfdw3uxdjeQDMG+oubxma7TBvEYay+sxQzsay2vdhu132T5oJuc46bnzffe23u/Pug3bbwYe6ti0yvaqKWZ3CPDjjvXN5baJtu9RKgUDSYcCLwLOB94y7cwWLuRxb3vzdJNX5kfv7Hsev6IG82rO7H2aeQ+fePBdjeQD8MT9m8vrwDkPNJbXwXN+1lhef3DA5sbyGlq88YczPcfd20b45hWPr5LXQ7aXzTS/vVnVx0QfBN4OjP+6c35ZUXKhpHndEkpaKWlY0vDI/c39gUREGBit8K8mW4AlHeuHltsm2r5H6RkMJL0Y2Gp73bhd5wBPAp4OLATe0S297VW2l9leNrTf/JmWNyKiMmN2eqTnUpO1wKvLVkXHAT+zfQdwBXCipAVlxfGJ5bY9SpXHRM8ETpb0QmAf4FGSPm37leX+7ZI+Aby1X4WMiJiuur75S/os8BxgkaTNFC2E5gDY/hvgcuCFwG3Ag8Bry33bJL0HuKE81Xm2J6uI3i16BgPb51DcBSDpOcBbbb9S0mLbd0gScCrQte1tRMTuYsxITcP02z69x34Db5hg32pgdS0F6ZOptCYa7zOSDqKoQV0PvL6eIkVE1GeUtC6sYkrBwPa1wLXl6+P7UJ6IiNoYGEkwqGQmdwYREXu83BlUk2AQEa1lYGem9q0kwSAiWss4j4kqSjCIiPYyjCQWVJJgEBGtVfRAjioSDCKixcRIS8cPq1uCQUS0VlGBnGBQRYJBRLRW0c8gwaCKBIOIaLXR3BlUkmAQEa2VO4PqEgwiorWMGMlU75UkGEREq+UxUTUJBhHRWkbs8NDuLsZeIcEgIlqr6HSWx0RVJBhERKulArmaVgYD72rum8CsebXNn9rT0Jzm8po3b1cj+Syc92Aj+QA8avZDjeW1aM59jeX1Bwdsbiyvv7r30Mbygo0zPoMtRpw7gyryLkVEq42inksVkpZL+q6k2ySd3WX/hZLWl8v3JN3bsW+kY9/aGi+vNpXvDCQNAcPAFtsvlnQEsAY4EFgHvMr2jv4UMyJi6ooK5Jk/ACk//z4CPB/YDNwgaa3tW36Zl/1HHcf/IfDUjlP8wvbSGRekj6ZyZ/Am4NaO9fcBF9o+ErgHOKPOgkVEzNRYBXKvpYJjgdtsbyq/9K4BTpnk+NOBz878CppT6V2QdCjwIuBj5bqA44FLy0MuBk7tRwEjImZixOq5VHAI8OOO9c3ltkeQdBhwBHB1x+Z9JA1Luk7SHvlZWfX+6YPA24H9y/UDgXttj9UyTvbGrARWAgwtWDD9kkZETNEUeiAvkjTcsb7K9qppZnsacKntzhYfh9neIukJwNWSbrT9/Wmevy96BgNJLwa22l4n6TlTzaB8Q1cBzHv8ksw5FBGNGq3Wmugu28sm2b8FWNKxfmi5rZvTgDd0brC9pfy5SdK1FPUJe1cwAJ4JnCzphcA+wKOADwEHSJpd3h1M9sZEROwWxUB1tTSavAE4qmw4s4XiA/+/jz9I0pOABcA3OrYtAB60vV3SIorP1PfXUag69XyXbJ9j+1Dbh1O8AVfb/n3gGuBl5WErgMv6VsqIiGkwYqeHei49z1N86T0LuIKiIc0ltm+WdJ6kkzsOPQ1YY7vzKciTgWFJ36H43LygsxXSnmImba7eAayR9F7g28DH6ylSREQ9bGrrdGb7cuDycdveNW793V3SfR14Si2F6KMpBQPb1wLXlq83UTS3iojYQ1XvVDboWjkcRUQElHUGGY6ikgSDiGi1TG5TTYJBRLSWUSa3qSjBICJay8DOGsYmGgR5lyKixZT5DCpKMIiI1jKVeyAPvASDiGi13BlUk2AQEa1lK3cGFSUYRERrFRXIvYebiASDiGi1zIFcVePBoInfy7xHb+9/JqV95zU30+eOnc39umYPjfQ+qAb7z27ud7WrwW+If3jAjxrL6y/ufXxjed21c//eB+1Bigrk1BlUkTuDiGi19ECuJsEgIlorPZCrSzCIiFarOOH9wEswiIjWsmHnaIJBFQkGEdFaxWOiBIMqEgwiotXSA7maBIOIaK00La2u5/2TpH0kfVPSdyTdLOlPyu0XSfqBpPXlsrT/xY2ImIriMVGvpdKZpOWSvivpNklnd9n/Gkk/7fhMPLNj3wpJG8tlRY0XWJsqdwbbgeNt3y9pDvCvkr5c7nub7Uv7V7yIiJmpYw5kSUPAR4DnA5uBGySttX3LuEM/Z/uscWkXAucCyyhuVtaVae+ZccFq1DMkunB/uTqnXNzXUkVE1KBoTTTUc6ngWOA225ts7wDWAKdULMZJwJW2t5UB4Epg+bQuqI8q3R9JGpK0HthKcVHXl7vOl7RB0oWS5k2QdqWkYUnDI/c/UFOxIyJ6G+t01msBFo19TpXLynGnOgT4ccf65nLbeC8tPxMvlbRkiml3q0rBwPaI7aXAocCxkn4DOAd4EvB0YCHwjgnSrrK9zPayof3m11TsiIhqRlHPBbhr7HOqXFZNI6t/AA63/ZsU3/4vrvM6+m1KDXBt3wtcAyy3fUf5CGk78AmK26iIiD3GWGuiCncGvWwBlnSsH1pu+1Ve9t3l5yHAx4CnVU27J6jSmuggSQeUr/elqED5d0mLy20CTgVu6mdBIyKmo6bWRDcAR0k6QtJc4DRgbecBY5+JpZOBW8vXVwAnSlogaQFwYrltj1KlNdFi4OKyNn0WcIntL0m6WtJBgID1wOv7WM6IiCmzxa4aeiDb3iXpLIoP8SFgte2bJZ0HDNteC7xR0snALmAb8Joy7TZJ76EIKADn2d4240LVrGcwsL0BeGqX7cf3pUQRETWqq9OZ7cuBy8dte1fH63Mo6lK7pV0NrK6lIH2SHsgR0VqD2gNZ0tHA2Bf2q7v0h3iEjOAUEa1WUwXyXkPSq4AvAAeXyxckvbpXutwZRERrDejkNm8HfnusXkLSB4BrgU9OlijBICJarY7hKPYyuzorqG3fI2m0V6IEg4hoLRt2Dd7kNt+WtGBs7KOya8CGXomaDQaGWTv6H6V33vFrfc9jjA9u7luH1NyQUL8YmdtIPg80lA/A3xz6jcbyevudj2iA1zezNdJYXj/ftU9jedVl0B4T2X7duPV7gZ4jpQ5cyIyIwTGFsYlaQ9LjJf29pK3lkNprJR3WK12CQUS0mq2eS8t8AriUosPwY4FLqNDHIcEgIlqt4kB1bXKg7U+XA4yO2P40cGCvRAkGEdFa9uD1MwB+Ws66NrtcXgv8tFeiBIOIaDExMjqr59IyrwX+G/Af5XJyuW1SaVoaEa3WwjqBSdneDLx0qukSDCKitQZxbCJJq+GRFSG2J707SDCIiPZyUW8wYL7U8Xo+8HtUmEwnwSAiWq2FrYUmZfsL4zZ9SlLPHpcJBhHRWi4rkIN/kzRke8Lu6gkGEdFqg/aYSNK5XTbfN1kggASDiGi5QWtNBNzX8Xo+D5+PeUI9g4GkfYCvAfPK4y+1fa6kI4A1FD3b1gGvsr1jGgWPiOgLu75gIGk58CGKOZA/ZvuCcfvfApxJMQfyT4HX2f5huW8EuLE89Ee2T66lUF3Y/sC4cr0PuLpXuioP07YDx9v+LWApsFzSccD7gAttHwncA5wx5VJHRPRZHT2QJQ0BHwFeABwNnF5OLdnp28Ay279JMTbQ+zv2/cL20nLpWyDopvyS/sPyGibUMxi4cH+5OqdcTDG/5qXl9ouBU6df3IiI/rB7LxUcC9xme1P54boGOOXh+fga2w+Wq9cBh9Z5HVVJWjx+m+3f71VnUKmaXdKQpPXAVuBK4PvAvbZ3lYdsBg6ZIO1KScOShkceeKBKdhERtTBidHRWzwVYNPY5VS4rx53qEODHHesTfuaVzgC+3LG+T3ne6yT1+4vz304nUaUK5DKiLC1nzPki8KSqGdheBawCmLdkyYDV60fE7lbxQ+cu28vqyE/SK4FlwO90bD7M9hZJTwCulnSj7e/Xkd94tl88nXRTaoBbzphzDfBfgAMkjQWTQ6nQwy0iolGubT6DLcCSjvWun3mSTgDeCZxse/svi2FvKX9uopicvrmp8CrqGQwkHVTeESBpX+D5FM2UrgFeVh62ArisX4WMiJg2V1h6uwE4StIRkuYCpwFrOw+Q9FTgoxSBYGvH9gWS5pWvFwHPBG6Z2UXVr8pjosXAxWVN9CzgEttfknQLsEbSeylq0T/ex3JGRExLHU1Lbe+SdBZwBUXT0tW2b5Z0HjBsey3wp8B+wP+TBL9qQvpk4KOSRik+Qy+wvfcFA9sb6HJLU97uHNuPQkVE1MHA6Gg9/QxsXw5cPm7buzpenzBBuq8DT6mlEBWVdyDHlavX2b6rV5rGeyCrgSrkuY9rrtXSAfN/0VheP71n/8bymjN3V++DavCpw/65kXwAXvXD3+l9UE3mDzXX//K+XfMay2vb9l9rLK9aGBiwHsiSTqJotPMvFO/AX0laafsrk6XLcBQR0WqDNjYR8H+BZ3f0fn48RSvQBIOIGGCDFww0FggAbP9IUs/GQgkGEdFilZuOtslWSQts3wNQtgbd8+oMIiIaNWB3BrZPGrd+L0WXgEll1oeIaC+DR9VzaZOyX8MHJa2T9C1JH5a0oFe6BIOIaDlVWFplNcVjoZcCv0sxnPbqXonymCgi2m3AHhMBT7D9ko7190j6Tq9EuTOIiHarZziKvckDkp4ztiLpuUDPzle5M4iI9hrATmfASoohhBaV69uAV/dKlGAQEa02aJ3ObN8EPE3SfuX6/T2SAAkGEdF2LWst1Es5b/3/Bk4CLOmfgPM7ZmHrKnUGEdFqcu+lZT4IPAo4HZgH3Az8Ra9EuTOIiPZqZwVxL79t+zcBJI3Y/jtJb+yVKMEgIlpMg1iB/LALlvRoYE6vRHlMFBHtNnhNS38gaWn5+gDgm8C7eyXKnUFEtNvo7i5As8rZ1cYsp5hxrefEK1XmQF4i6RpJt0i6WdKbyu3vlrRF0vpyeeH0ix8R0Qdj/Qx6LS0i6bCxBXgIeIyk70h67GTpqtwZ7AL+2Pa3JO0PrJN0ZbnvQtt/NrOiR0T0T12thSQtBz5EMQfyx2xfMG7/POCTwNOAu4FX2L693HcOcAYwArzR9hX1lKqrf6CoN+i88icC/yJp2Pbp3RJVmQP5DuCO8vV9km4FDpl5eSMiGlBDMJA0BHyEYijozcANktaOm9j+DOAe20dKOg14H/AKSUcDpwG/DjwOuErSf7I9MvOSPdJYS6Jx5f+W7WMkrZ8o3ZQqkCUdDjwVuL7cdJakDZJWTzREqqSVkoYlDY8+0NzcxBERNToWuM32Jts7gDXAKeOOOQW4uHx9KfA8SSq3r7G93fYPgNvK8zXp78ufn5jogMoVyGXX5s8Db7b9c0l/DbyHIu6+B/hz4HXj09leRTE5M/OWLGlfvX1E7NEqPiZaJGm4Y31V+dk15hDgxx3rm4FnjDvHL4+xvUvSz4ADy+3XjUvbt6crks7tsvn1krbzq6DwCJWCgaQ5FIHgM7a/AGD7zo79fwt8qeeJZsHI/P5X7Wvj/n3PY8xPFu7bWF6Peux9jeW1/ulrGsnn6K+/spF8AEYbHJZg3312NJbXrpGhxvLavn0va4Boqg5HcZftZX0uTVO6fVDsoqjrWEPxdOcRev5my9ucjwO32v5Ax/bFZX0CwEuAm6Za4oiIvqvnecQWYEnH+qHltm7HbJY0G3g0RUVylbS16fycHiPpRbb/j6SXTpSuSp3BM4FXAcePa0b6fkk3StoAPBf4o+kWPiKiX2oam+gG4ChJR0iaS1EhvHbcMWuBFeXrlwFX23a5/TRJ8yQdARxF0RGsMbaf11Gurqq0JvpXus8Ld/k0yxUR0Zwa7gzKOoCzgCsoHrestn2zpPOAYdtrKZ6gfErSbRRzCJxWpr1Z0iXALRSPa97Qr5ZEAJL+0faLJriOH0yUbi97ABgRMUU1NVuxfTnjvgTbflfH64eA35sg7fnA+fWUpKczppMowSAiWqulQ1RPyvZPppMuwSAi2m3AJreZrgSDiGi1QbszmK4Eg4hotwSDShIMIqK9BrDOQNKzgBcBTwb2pRhb7t+Az9veNlG6BIOIaLcBCQaSnk4xLNAPgWuAL1MMYX0QxUiqX5J0ve2ufcISDCKi1TQ4k9v8Anix7Z932fePwHnlXUNXCQYRES1gu+eQQGUn4q4yB3JEtNuAzYEsaaekM7tsXy7p0xOlSzCIiPaqMC5RCyuYbwdeK+ldnRttfwWYcGTWBIOIaLcBuzMA7geOB5ZKWlWOoDpmwrHVEwwiot0GLxi4nFXtdykCw/WS3izp75hktNRUIEdEa4mBak005gtjL2y/RdJvAy8A/hn42ESJEgwior3aWScwKdvvHbf+deDrvdIlGEREuw1YMAAoJyA7mV/NtbwFuMz2lydKk2AQEe02YMFA0l8CRwCfBH5cbl4CnCXpBbbf2C1do8FAu2DuXf2fvHvH43b2PY8xQ/vuaiyvpiapB1h6w2mN5LNzR3P/BYdm921yqUfYsbO569q5s/9/U2N2PTSnsbzqMmiPiYDlto/ssv1zkjZOlKhnayJJSyRdI+kWSTdLelO5faGkKyVtLH8umEHhIyL6Y/BaEz0g6bjxGyU9A3hgokRVmpbuAv7Y9tHAccAbJB0NnA181fZRwFfL9YiIPYeL1kS9lpmq8uVY0lJJ3yi/VG+Q9IqOfRdJ+oGk9eWydAbFeS3wQUn/Lumqsjy3Ah8CXjNRop73srbvoBgCFdv3lSc9BDgFeE552MXAtcA7ZnABERH1a+ab/9iX4wsknV2uj/88fBB4te2Nkh4HrJN0he17y/1vs33pTAti+1vAcZIOpqMC2fadk6Wb0oNNSYcDTwWuBw4uAwXAT4CDJ0izElgJMPvReZIUEc1qqM6g55dj29/reP0fkrZSDC99LzUr70yeS0cwKAPPPROlqdwDWdJ+wOeBN48fItX2hE/ebK+yvcz2sqH586tmFxFRj2p1BoskDXcsK6eYS6Uvx2MkHQvMBb7fsfn88vHRhZLmTTH/znO/BhimCE7zyuV3gGFJKyZKV+nOQNIcikDwGdtjvdvulLTY9h2SFgNbp1v4iIi+qF5BfJftCQdxA5B0FfDYLrve+bAsbUsT34+Un5efAlbYHquxOIciiMwFVlHcVZxXqeTdy3OM7Z+Ny/fRwA0Udy6P0DMYSBLwceBW2x/o2LUWWAFcUP68bHrljojoD1HfYyLbJ0yYj1Tpy7GkR1FMNPNO29d1nHvsrmK7pE8Ab51JUen+1EdMEhqr3Bk8E3gVcKOk9eW2/0URBC6RdAbFNGsvn1JxIyIa0FCdQc8vx5LmAl8EPjm+orgjkAg4Feg5Uc0k3kvxSOhKYHO57RDgROBPJkpUpTXRv1JElG6eN8VCRkQ0q5lg0PXLsaRlwOttn1luezZwYPlcH+A1ttcDn5F0EMVn7Xrg9dMtiO1PSloLnMSvKpCvBc4e/+ioU4ajiIh2ayAY2L6bLl+ObQ8DZ5avPw10nWnM9vE1l+de4HNTSZNgEBHtNYCjlk5XgkFEtFuCQSUJBhHRagM4uc20JBhERKsN2mMiSZvo0ujH9hGTpUswiIj2aueopL10dp6bD7wCWNQrUYJBRLTbgAUD29s6VrcBfyZpuFe6BIOIaK06eyDv5d4vaVbH8BePkGAQEa2m0cGKBpJW072j8AskYfu13dIlGEREew1mncGXppMowSAiWm0AHxPtB1zWOfSEpKfZXjdZosrzGURE7JUGbw7kDwP/Iuk3OrZ9rFeiRu8MPAtG5jXwzjf4VWDjcy5qLK+jrn1NY3mN3DenkXwe+/htvQ+qyd0/b25ypblzdjWW10MPzm0sr73RAN4ZbAL+B/A5SefZ/hwVQl4eE0VEuw1eMJDtdZL+K0VAeBoVPuvzmCgi2svFcBS9lpYx/LK/wYnAKPDrvRIlGEREa431M+i1tIntYzpe2/bZtod6pUswiIh2s3svLSDppZJOldT1kZCkx0v604nSp84gIlqtbd/8J3El8MfA+ZK+D3wPeAh4DPBUiqEpLpgocc87A0mrJW2VdFPHtndL2iJpfbm8cKZXERFRuyrNSlsSLGz/3Pa5wFMo5kG+BbgD+DLwItsn2b5movRV7gwuAv4S+OS47Rfa/rNplToioiFNVBBLWkgxzeThwO3Ay23f0+W4EeDGcvVHtk8utx8BrAEOBNYBr7K9YzplKccf+ma5VNbzzsD21yhuLyIi9joNtSY6G/iq7aOAr5br3fzC9tJyOblj+/sovmAfCdwDnFFLqaZgJhXIZ0naUD5GWjDRQZJWShqWNDx6/wMzyC4iYopMUxXIpwAXl68vBk6tmlCSgOOBS6eTvi7TDQZ/DTwRWErxTOrPJzrQ9irby2wvm7Vfcz1AIyKgctPSRWNfWstl5RSzOdj2HeXrnwAHT3DcPuX5r5M09oF/IHCv7bFu65uBQ6aY/4xNqzWR7TvHXkv6W6Y5Sl5ERN9V++J/l+1lkx0g6SrgsV12vfNh2dmWJmzDdJjtLZKeAFwt6UbgZxMc26hpBQNJizui4EuAmyY7PiJid6hzchvbJ0yYj3Tn2OeipMXA1gnOsaX8uUnStRRNPj8PHCBpdnl3cCiwpZ5SV1elaelngW8A/1nSZklnUMyac6OkDcBzgT/qczkjIqbORqO9lxqsBVaUr1cAl40/QNICSfPK14uAZwK32DZwDfCyydL3W887A9und9n88T6UJSKifs30I7gAuKT8svxD4OUAkpYBr7d9JvBk4KOSRim+iF9g+5Yy/TuANZLeC3yb3fAZmx7IEdFqTfRAtn038Lwu24eBM8vXX6foENYt/Sbg2H6WsZcEg4hoLwMDNgfydCUYRES7JRZUkmAQEa02QAPVzUiCQUS0Wk2thVovwSAi2qtFo5L2W4JBRLRW0eks0aCKZoOBwHP7/4vZtLy5JrpP+EqDgwtajWU1e1sz/zXuPXDfRvIB2HVnc3n9YnFzH0Cj23vOaFgbzd4LJwzeC4u8O+TOICJaLXcG1SQYRER7pc6gsgSDiGix2sYear0Eg4hotzwmqiTBICLay83MgdwGCQYR0W65M6gkwSAi2i2xoJIEg4hoNY3mOVEVCQYR0V4mnc4qqjLt5WpJWyXd1LFtoaQrJW0sfy7obzEjIqZOGLn3EhWCAXARsHzctrOBr9o+CvhquR4Rseexey/ROxjY/hqwbdzmU4CLy9cXA6fWXK6IiHo0EAyqPC2R9FxJ6zuWhySdWu67SNIPOvYtnXGhpqjKnUE3B9u+o3z9E+DgiQ6UtFLSsKThkfsfmGZ2ERHTMFZn0GuZuZ5PS2xfY3up7aXA8cCDwD91HPK2sf2219dSqimYbjD4JduTjv5he5XtZbaXDe03f6bZRURMiUZHey41mOrTkpcBX7b9YB2Z12G6weBOSYsByp9b6ytSRERdKjwiqqfOoPLTktJpwGfHbTtf0gZJF0qaV0ehpmK6wWAtsKJ8vQK4rJ7iRETUyFQNBovGHmeXy8rxp5J0laSbuiynPCzLHk9Lyi/QTwGu6Nh8DvAk4OnAQuAdM772KerZz0DSZ4HnULxZm4FzgQuASySdAfwQeHk/CxkRMW3VngLdZXvZZAfYPmGifZLulLTY9h0Vnpa8HPii7Z0d5x67q9gu6RPAWyuVukY9g4Ht0yfY9byayxIRUbuG+hGMPS25gN5PS06nuBP4pY5AIor6hpu6puyjGVcgR0Ts0ZqpM7gAeL6kjcAJ5TqSlkn62NhBkg4HlgD/PC79ZyTdCNwILALeW0ehpiLDUUREe9kw0v/xKGzfTZenJbaHgTM71m8HDuly3PH9LF8VCQYR0W7pYVxJo8HgKQt+yjdf/tG+53PkJf+z73mMGdqh5vLa3lxeOxaNNJKPb9+vkXwA1Nzbx46f/Fpjec1qdCC2oSYzq0eCQSW5M4iI9jKQOZArSTCIiBYzOGNYV5FgEBHtZRqpQG6DBIOIaLfUGVSSYBAR7ZZgUEmCQUS0WCavqSrBICLay0A9Q1S3XoJBRLRb7gwqSTCIiBZrZjiKNkgwiIj2Mjj9DCpJMIiIdksP5EoSDCKi3VJnUEmCQUS0l53WRBUlGEREu+XOoJIZBQNJtwP3ASPArl5ziEZENMt4pJnh2Pd2ddwZPNf2XTWcJyKiXhnCurI8JoqIdkvT0kpmzTC9gX+StE7Sym4HSFopaVjS8E/vzu1aRDTHgEfdc5kpSb8n6WZJo5ImfFwuabmk70q6TdLZHduPkHR9uf1zkubOuFBTNNNg8CzbxwAvAN4g6dnjD7C9yvYy28sOOnAvnDIvIvZeLie36bXM3E3A7wJfm+gASUPARyg+L48GTpd0dLn7fcCFto8E7gHOqKNQUzGjYGB7S/lzK/BF4Ng6ChURURePjPRcZpyHfavt7/Y47FjgNtubbO8A1gCnSBJwPHBpedzFwKkzLtQUTbvOQNJ8YJbt+8rXJwLnTZZm3Ybt9w8t3tjrDetlEdCjwvqtNZ6rNk3m1aS2XtdkBvGamzAUzJXWAAAD0klEQVT+fT1spie8j3uuuMqXLqpw6D6ShjvWV9leNdP8xzkE+HHH+mbgGcCBwL22d3VsP6TmvHuaSQXywcAXi6DGbODvbH+lR5rvzrT5qaThupqw1nmuPSmvJrX1uiYziNfchH68r7aX13UuSVcBj+2y6522L6srn91l2sHA9ibgt2osS0TEHsv2CTM8xRZgScf6oeW2u4EDJM0u7w7GtjdqphXIERFRzQ3AUWXLobnAacBa2wauAV5WHrcCaPxOo+lgUMczuDqf49X9THBPyatJbb2uyQziNTdhr31fJb1E0mbgvwD/KOmKcvvjJF0OUH7rPwu4ArgVuMT2zeUp3gG8RdJtFHUIH2/8GpxxOyIiBl4eE0VERIJBRET0KRhIWiLpGkm3lF2031RuXyjpSkkby58LKpzrdkk3Slo/1g646nkkrZa0VdJNHdu6plXhw2V38A2SjunnNc80vyZN5XewN11Xp7r+r0haUR6/UdKK3XEte4o6/ybyvjbAdu0LsBg4pny9P/A9iu7X7wfOLrefDbyvwrluBxaN21bpPMCzgWOAm3qlBV4IfBkQcBxwfT+veab5NblM5XewN11X3f9XgIXApvLngvL1gt19bbvxPa3lbyLva0O/r4b+U1wGPB/4LrC44z/Kdyuk7fZBVPk8wOHj/sC7pgU+Cpze7bh+XHPd+fX591f5d7A3XVfd/1eA04GPdmx/2HGDvkz3byLvazNL3+sMJB0OPBW4HjjY9h3lrp9Q9GLupdvIqNM5z5iJ0nbrKj6tLuEVr7m2/Bowld/B3nRdvUz1Gtt07bWa4d9E3tcG9HU+A0n7AZ8H3mz75yqGrgDAtiVVadf6LNtbJD0GuFLSv3funMJ5HmEmaSdS0zXvafr2O9hbDMI19ktL/yZap293BpLmUPwH+IztL5Sb75S0uNy/GNja6zzuPjLqlM/TYaK0E3UVr2yK1zzj/Joyxd/BXnNdFUz1Gtt07bWo6W8i72sD+tWaSBQ96G61/YGOXWspulpDhS7XkuZL2n/sNcXIqDdN9TzjTJR2LfDqskXDccDPOm5le5rGNc8ov6ZM43ewV1xXRVO9xiuAEyUtKFvInFhuG0g1/k3kfW1CPyoigGdRPGfeAKwvlxdSdLP+KrARuApY2OM8TwC+Uy43U4wOSNXzAJ8F7gB2UjxnPGOitBQtGD4CfB+4EVjWz2ueaX5NLVP9Hewt19Wv/yvA64DbyuW1u/u6dvN7WtvfRN7X/i8ZjiIiItIDOSIiEgwiIoIEg4iIIMEgIiJIMIiICBIMIiKCBIOIiAD+PyrfZvh3zR7FAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# correlation coefficients\n",
    "pl.plot_corr(scb[:100,:],mb.bp,spc='p2',clab='corr. coeff. (100 rlz)')\n",
    "pl.plot_corr(scb,mb.bp,spc='p2',clab='corr. coeff. (200 rlz)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "chi^2: 6.1 , chi^2 (sim): 10.5 , PTE:0.845 (Lmin=10)\n",
      "chi^2: 4.2 , chi^2 (sim): 10.5 , PTE:0.935 (Lmin=30)\n"
     ]
    }
   ],
   "source": [
    "# x2pte dependencies on Lmin\n",
    "for lmin in [10,30]:\n",
    "    mb  = bins.multipole_binning(10,spc='p2',lmin=lmin,lmax=2048)\n",
    "    scb = prjlib.binned_cl_rlz(p0.quad.f['EB'].cl,1,200,mb)\n",
    "    ocb = prjlib.binned_cl(p0.quad.f['EB'].cl[0],mb)\n",
    "    st = ana.PTEs(ocb,scb,decimals=5,comment='Lmin='+str(lmin))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "chi^2: 11.3 , chi^2 (sim): 16.3 , PTE:0.765 (bin num = 15)\n",
      "chi^2: 14.3 , chi^2 (sim): 22.3 , PTE:0.88 (bin num = 20)\n"
     ]
    }
   ],
   "source": [
    "# x2pte dependencies on number of bins\n",
    "for bn in [15,20]:\n",
    "    mb  = bins.multipole_binning(bn,spc='p2',lmin=20,lmax=2048)\n",
    "    scb = prjlib.binned_cl_rlz(p0.quad.f['EB'].cl,1,200,mb)\n",
    "    ocb = prjlib.binned_cl(p0.quad.f['EB'].cl[0],mb)\n",
    "    st = ana.PTEs(ocb,scb,decimals=5,comment='bin num = '+str(bn))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "py4so",
   "language": "python",
   "name": "py4so"
  },
  "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}