PCA Comparison (Scikit-allel vs Scikit-allel vs Dask-ML)

pca_comparison.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PCA Comparisons\n",
    "\n",
    "Dask-ML, Scikit-learn, and Scikit-allel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import allel\n",
    "import numpy as np\n",
    "import xarray as xr\n",
    "import matplotlib.pyplot as plt\n",
    "from sgkit.testing import simulate_genotype_call_dataset\n",
    "from sklearn.decomposition import PCA as skPCA\n",
    "from dask_ml.decomposition import PCA as daPCA\n",
    "xr.set_options(display_style='text');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1000, 100)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Count non-reference alleles from simulated data\n",
    "ds = simulate_genotype_call_dataset(n_variant=1000, n_sample=100)\n",
    "gn = ds.call_genotype.sum(dim='ploidy')\n",
    "gn.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre>&lt;xarray.DataArray &#x27;call_genotype&#x27; (variants: 5, samples: 5)&gt;\n",
       "array([[0, 1, 2, 2, 1],\n",
       "       [1, 1, 1, 0, 2],\n",
       "       [1, 1, 1, 1, 2],\n",
       "       [1, 1, 1, 1, 0],\n",
       "       [1, 0, 1, 1, 2]])\n",
       "Dimensions without coordinates: variants, samples</pre>"
      ],
      "text/plain": [
       "<xarray.DataArray 'call_genotype' (variants: 5, samples: 5)>\n",
       "array([[0, 1, 2, 2, 1],\n",
       "       [1, 1, 1, 0, 2],\n",
       "       [1, 1, 1, 1, 2],\n",
       "       [1, 1, 1, 1, 0],\n",
       "       [1, 0, 1, 1, 2]])\n",
       "Dimensions without coordinates: variants, samples"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gn[:5, :5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1000, 100)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Apply Patterson scaling to counts\n",
    "def rescale(x, ploidy=2):\n",
    "    mean = x.mean(dim='samples')\n",
    "    p = mean / ploidy\n",
    "    std = np.sqrt(p * (1 - p))\n",
    "    return (x - mean) / std\n",
    "gnr = rescale(gn).chunk(chunks=gn.shape)\n",
    "gnr.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Deterministic PCA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(100, 2)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Run Dask-ML PCA\n",
    "pca = daPCA(n_components=2, svd_solver='full')\n",
    "pcs1 = pca.fit(gnr.data.T).transform(gnr.data.T).compute()\n",
    "pcs1.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(100, 2)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Run scikit-learn PCA\n",
    "pca = skPCA(n_components=2, svd_solver='full')\n",
    "pcs2 = pca.fit(np.asarray(gnr.data.T)).transform(np.asarray(gnr.data.T))\n",
    "pcs2.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(100, 2)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Run scikit-allel PCA\n",
    "pcs3 = allel.pca(gn, n_components=2, scaler='patterson', ploidy=2)[0]\n",
    "pcs3.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAADQCAYAAAAK/RswAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAj9klEQVR4nO3de5RkdXnu8e8zTaM9RG2QAWFkHIKEKBJuLZAQPUbF8RLkJgEC0RNzMnEtNQvxzArIRDDKgeV4SzSJGRKiMYgilxKCx1GIEUMAHeyBZpSJgATp4cAYGRDokKbnPX/UrrGmpi67qnvX3rXr+azVq7tuXe8Uzbve/bu8P0UEZmZmZmWyKO8AzMzMzBaaCxwzMzMrHRc4ZmZmVjoucMzMzKx0XOCYmZlZ6bjAMTMzs9JxgWNmZmal4wLH+kLSA5JmJD0p6RFJfy/pl5LHVki6WdLPJW2R9G1Jb00e20fSdZI2SwpJy3P9h5hZYc0jz7xF0r9K2irp/0m6VNLz8v3X2Hy5wLF+Oj4ifgk4AnglsFrS24CvAP8AvBjYG/ggcHzymm3A14FT+h+umQ2gXvLMC4CPAPsCL0ues6bPcdsCkzsZWz9IegD4XxFxY3J7DfBy4BDg0xHRNplI2gWYBfaPiAeyjdbMBtF880zd7zkZ+FBEHJJVrJY9j+BY30naD3gz8DSwH3BVvhGZWdnMM8+8GtiYRVzWP7vkHYANlYqkZ4HHgRuAy4G3AQ/nGpWZlcm88oyk44B3AEdnFqH1hQsc66cTa0PHAJJ+NflxH+DH+YRkZiXTc56RdAzwReBtEfHv2YVo/eApKsvTJuAneAGxmWUnVZ6RdDhwHfDOiLipH4FZtlzgWG6iusL9HOBPJf2+pOdLWiTpNyWtrT1P0nOB5yQ3n5PcNjPrKE2ekfQKqrs13xsR1+cZry0cFziWq4i4CjgNeCewGXiE6nbNr9Y9bQZ4Mvn5nuS2mVkqKfLM+4ElwN8lPXSelORFxgPO28TNzMysdDyCY2ZmZqXjAsfMzMxKxwWOmZmZlY4LHDMzMyudUjT623PPPWP58uV5h2Fmde64446fRsSSvONYSM41ZsXTKteUosBZvnw569evzzsMM6sj6T/yjmGhOdeYFU+rXOMpKjMzMysdFzhmZmZWOqWYojKz7FUmp1mzbhObt86w7/gYq1YcxImHL807LDMrkYXMMy5wzKyjyuQ0510zxczsHADTW2c475opABc5ZrYgVlemuPy2B6mdrzDfPOMpKjPraM26TduLm5qZ2TnWrNuUU0RmViaVyekdipua+eQZj+CY2U4ah4mntzY/33Rzi/vNzDqpzzOLpJ2Km5pe84wLHDPbQbNhYkHT5LPv+FgfIzOzsmic9p5rc/B3r3nGU1Rmtl2rYeIA1HDf2OgIq1Yc1KfIzKxMmk17NyPoOc+4wDGz7das29RymDiApeNjKPl+8cmHDMQCY0mXSXpU0t119+0h6ZuSfpR83z3PGM2GTZppJwFnHrPMu6jMbP7aJZ2l42Pccu5r+xjNgvkc8BngH+ruOxe4KSIukXRucvtPcojNbCg0rusbXzzKY0/P7vS8EYltEQvSisIFjtmQatZvotWC4vkME+ctIm6WtLzh7hOA1yQ/fx74F1zgmGWiWZuJ0UVidETMzv1izHhsdGRBR4Zd4JgNoVZ9bU45cilX3zG9w9z4fIeJC2rviHgYICIelrRX3gGZlUnjDqnGRcSz24LxsVF2e84umTUPzbXAkXQZ8NvAoxHxiuS+PYAvA8uBB4DfiYjH8orRrIxa9bX51j1buPjkQ9yxuI6klcBKgGXLluUcjVnxpd0h9fjMLBsueENmceQ9gvM5PDdulrlu+tqcePjSYShoHpG0TzJ6sw/waKsnRsRaYC3AxMRE672sZkOulmda5ZdGWbeZyLXA8dy4WXbqk019Hxv3tQHgOuAdwCXJ96/mG47ZYGvsn9VJP9pM5D2C00yquXEPG5u11phsWvW1qb+/rH1tJF1B9aJpT0kPARdQLWyulPQHwIPAqflFaDa4KpPTXHjdRrbO7LwjqtFC7pBKo4gFTioeNjZrbnVlin+87cGOz6v1tSn7WpuIOKPFQ6/rayBmJdO41qadhd4hlUYRC5zUc+NmtqO0xQ0MdF8bMyuAtN2Il+Z0AVXEAsdz42Y9qB2zkEZZp6PMrH86dSMW8MnTDsttZDjvbeKeGzfrUePOqKeeeTbVAr+8rqbMrFza7cgsQv+svHdReW7crAfNTvzu5KxjlvGREw/JNjAzGxqrVhzUdA3O7otHueD4g3O/iCriFJWZtVCZnOZD129seoZLOy5uzKxbzY5zqS9aaj8XtTGoCxyzAVGZnGbVVXfucHZLJ7VhYhc3ZpZGu/5Z510zBbBTkVOUgqaRCxyzAZB2d9Tui0dZvGt2Z7uYWXl16p81MzvHmnWbBianuMAxK7gzL72VW+77WcfnCQox721mg6e2C7PT+HCnnVNF4gLHrMAqk9Opi5u8dyyY2eDppncWDNZxLi5wzApszbpNHZ+Td68JMxtMaUeHawatf5YLHLMC6zQcPLpIrDn1UBc3ZtaVtKPDNYPYP8sFjllBNNuS2a6R1tjoIi4++dcGKuGYWTGkHR0e5F2YLnDMctbsNN7alsxTjlzK1XdM79RI69gD9uDyP/z1fodqZgOq8QIqTXPQQZ/6XpR3AGbDrHYab31xUzMzO8e37tnCxScfwtLxMUR1mPhTpx3m4sbMUqvlmemtMwRs73HTzlkl2LTgERyzHHU6jXfz1plCN9Iys+JrlmcCdmjkV68snc9d4JjlqNMi4kHakmlmxdQqzwTVUeGyNgZ1gWPWJ90vIh6sLZmDSNIDwM+BOeDZiJjINyKz+ekmzywdH+OWc1+bQ5T94QLHLEOtDsfstIi4KKfxDonfioif5h2E2Xysrkxx+e0PEnVzTu3yzDBcQBV2kbGkByRNSdogaX3e8Zh1q3Y4ZquTv9stIp784Btc3JhZKrVuxNFkQU2rPHPxyYeUPsf0NIIj6ZUR8b2FDqYJX1nZwFqzblPHk7+9iLi9PuSaAL4hKYC/iYi1TWJYCawEWLZsWYahmPXm8g5HLQxrnkld4Eh6OXA6cAbwOOC5arM20hxK50XEO+tzrjk2IjZL2gv4pqR7IuLm+ickRc9agImJiU5nEZr1Rf1am05/lMOaZ9oWOJJeQjXJnAE8C7wEmIiIB7IPrfOVlVmRdWqmNQxz4GnllWsiYnPy/VFJ1wJHATe3f5VZfro9HHOY80zLNTiS/g34GjAKvC0ijgR+3qfiBqpXVkcAbwLeLenVDfGtlLRe0votW7b0KSSz9FatOIjRkebttHZfPDoUc+Bp5JVrJO0m6Xm1n4E3AHdn+Z5m83HcJ/6lq+IGGOo8024EZwvwYmBvYAnwI5r3BMpEpysrDxtb0dWSSv0uqvGxUS58q3dHNcgr1+wNXCsJqrnwixHx9T68r1nXzrz0Vn706FNdvaYM3Yjno2WBExEnSHoBcArwIUkvBcYlHRUR380yqORqalFE/LzuyurPsnxPs26srkxxxe0/YS6CEYkzjt6vaefPYVzY1628ck1E3A8cmtXvN5uv2jqbNOdG1WuXk4ZJ2zU4EfE4cBlwWbII7zTgU5L2i4j9MozLV1ZWSJXJaVZ9ZQOz235x31zE9mHjYU8ovcox15gVUjXX3MnstvSDmZ8a8MMxF1rLAkfSc4HnRcQWqE4VAZ+WdCWwR5ZB+crKiui4T/xL2yHiK27/iQucHuSZa8yKqNuFxADHHrCHi5sG7Rr9/QXwqib3vx7442zCMSumNPPfc826bFkazjVmVPPM8nNv6Lq4OeuYZVz+h7+eUVSDq90U1W9GxMrGOyPickkfyDAms8K55b6fdXzOiJrvmLKOnGtsqFUmpzn7yxu6ft2xB+zhwqaNdgVOu2xd2CMezBZC4yLiNM442ktFeuRcY0Orl+mopSU8+TsL7QqcR5vtYpD0SqrbOs1KqTHhpJl6OvaAPbz+pnfONTZUWh3C28noiFjztkNd2KTUrsBZBVwp6XPAHcl9E8DbqbZRNyuNbtqe19tlkfjYqU448+RcY0Ojdghvp3PqGu226wgXnTS8Tft60a4PznclHQW8G/ifyd0bgaOTXQ5mpXDmpbemWmPTyPPfC8O5xoZJmkN4G511zDKPEPeg02GbS4Ep4EsR8cM+xGPWV6srU6mLmxGJ+y5+c8YRDS3nGhsKaQ7hrdn7ebty+/nHZRhNubXrg/NB4CyqQ8YflXRxRFzat8jMMtRLh1AvIs6Gc42VVX2eGZG2b1pIs67Pozbz124E5zTgsIh4WtILga8DTjo28FZXprj8tgdTr7Vx2/PMOddY6TR2Iq4VNe2Km90Xj3LB8T6rbqG0K3D+KyKeBoiI/5Tk7Zo28LrZkum2533jXGOlkibPLBLUTmHwIbzZaFfgHCDpuuRnNdwmIt6aaWRmC8RtzwvPucZKI22+iYAHLnlLHyIaXu0KnBMabn8sy0DMstDtDilPR+XCucZK44rbf5LqefuOj2UcibXbJv7tfgZitpAqk9NceN1Gts6ka6Ql4JOeksqFc42VSZoFxGOjI6xacVAfohlunbaJmw2U6sK+Dcxu6+51Zx6zzMWNmc1bp11SXkjcPy5wrDR6WWsjqsWNp6SGk6Q3An8OjAB/GxGX5BySDYD6zuf7NpwLdcbR+zXNQ2Oji7j45F9zYdNHhS1wnHisG5XJ6a6Lm11HxEd9rsvQkjQC/CVwHPAQ8D1J10XED/KNzIqscV3f9NYZzrtmCoATD1+6/WKp/rBer+vLR7tGf9dD61YhWe5scOKxbq1Zt6mr5/uYheLIMdccBdwbEfcncXyJ6oJn5xnbSbv+WTOzc6xZt2n7xdJHTjzEBU0BtBvByXMngxOPdSVt+3P3timkvHLNUqB+y8tDwNGNT5K0ElgJsGzZsv5EZoWSZjdmN0cwWH+k2kUlaQxYFhHdXSb3LlXiseG0ujK10/DvvuNjHY9dOMsLiQspx1yjZuHsdEfEWmAtwMTERHenJNrAqkxOc941dzGTcseCt30XT8eOoZKOBzZQbZ+OpMPqm3BlpGPikbRS0npJ67ds2ZJxOFYElclpXv6n/5d/vO3BHdqe/+NtD7L8hWOMLmr2Z1Plc12KL4dc8xBQf8DYi4HNGb6fDYjK5DRnf3lD6uIG8LbvAkrTEv1CqlNGWwEiYgOwPKuAEh0TT0SsjYiJiJhYsmRJxuFY3lZXpjj7yxt4ukXCue3+x1hz6qGMj41uv2/3xaN86rTDeOCSt7i4GQwX0t9c8z3gQEn7S9oVOB3I+uLNBsAHrrmrq+e783kxpdlF9WxEPC61vjrOwPbEA0xTTTy/288ArDjSbP+ei+DEw5c6yQy2vuaaiHhW0nuAdVR3a14WERv78uZWaK0upJrx6HBxpSlw7pb0u8CIpAOBPwb+LcugnHgMqoXNF29/cPuBdO2M9LcAt2zkkWu+Bnwty/ewYqtMTnP+tVM89d9zQPP1Ec14J2bxpSlw3gucDzwDXEG16PhwlkGBE88wq0xO84Fr7urqKuqMo/fr/CQrulxyjQ2vZrujOl1PuWHf4OhY4ETE01STzvnZh2PDrtvDMcFDxGXhXGP91EuuOXCv3fjmOa/JJiBbcIVs9GfDaXVlqquEs0jwid9xX5tB51xj/dTNQby1c6XcjXgwFbXRnw2R+r42ae226wgXnXSIi5tycK6xvqj2tpliZnau43NHJO67+M19iMqykqrRn1kW2rU+b8Xz3+XjXGNZqx2O2akZaD2v6xt87aaopmg/bPxrmURkpdfLqd/gtTZl5VxjWem2G3GNc005tJui+u2+RWFDoTI5zaqvbKDLXOO1NuXnXGMLrppv7mQ2TZ+JhLd+l0u7Kar/6GcgVm6VyWnO+fIGuqxtvNZmCDjX2EJaXZni8tsfpIslfey+eJQLjj/YeaZkOm4Tl3QM8GngZcCuVBvvPRURz884NiuBXua+wUPEw8i5xuar2+nvpeNjrFpxkAubkkrT6O8zVI9K+AowAbwdeGmWQVk59LKIGFzcDDHnGutaL1PfAj55mqe9yy5NgUNE3CtpJCLmgL+XlGn7dBt8lcnprosbFzbmXGPd6HXDwpnHLHNxMwTSFDhPJyftbpD0UeBhYLdsw7JB1c35UfU+5aspc66xLlQmp7subrzWZrikKXB+D1gEvAd4H7AfcEqWQdlg6uVqyn1trI5zjbXV65o+78QcTmnOoqrtcPgvSddHxPczjskGTC9nunhxnzVyrrF2uulCXG/x6CL+jy+ihlKqNTh1/hY4IotAbPB0u/VbVOe+vc7GUnCuse0qk9O8/8o7uzrOxT1trNsCR5lEYQOn2+JmROLjv3Oor6IsrUxzjaQLgT8EtiR3fSAivpble1pvej3128WNdVvgfCiTKOo48RRfZXKa9125IfUOqUXg4sa6lXmuAT4ZET7os8BWV6a6Km68iNjqpWn0d1NEvA4gIiqN92XEiaeguu1t40XEllZOucYKqNvFxGOjI1x8sjue247aHbb5XGAxsKek3fnFkPHzgX37EJsVRC87Fzz/bWnlmGveI+ntwHrg/RHxWIv4VgIrAZYtW5ZhOFaZnOZD12/ksadnU7/GGxaslXYjOH8EnE01wdTvZngC+MsMY4KUicey10s3Ys9/W5cyyTWSbgRe1OSh84G/Bj5M9RTzDwMfB97Z7PdExFpgLcDExES3jbktpW5zjXtnWSftDtv8c+DPJb03Ij69kG+6EInHV1XZczdi64esck1EvD7N8yRdCvzTQr2vda/bXHPsAXu4uLGO2k1RvTYi/hmYlnRy4+MRcU2vb7oQicdXVdnpZUpqfGyUC9/qxX3WvSxzTZv33CciHk5ungTcvdDvYemtWbcpVXEzInHG0fv5IspSaTdF9T+AfwaOb/JYAAuedMCJJ2/dNtNybxtbAHnkmo9KOiz5/Q9QnSaznGzucDHlRcTWi3ZTVBckP74rIp6pf0zSHhnG5MSTozXrNqUubry4zxZCHrkmIn4vi99rvdl3fKzliLFHh61XafrgXCPphIh4FkDSi4AbgCOzCMiJp38qk9N84Jq7eHq22q5Pgk6NQj1iYxnqa66x/qpNfW/eOsO+DRdHq1YctNPIsXONzVeaAqcCXCXpFKqH310H/O8sg7LsVSanOefKDTuc+t2puPGIjWWsgnNNKTVOfU9vneG8a6YAOPHwpdtzSqsCyKwXaQ7bvFTSrlSTz3LgjyLi3zKOyzK2Zt2mHYqbdjz/bf3gXFNezaa+Z2bnWLNu0/a8Ul/omC2Edruozqm/SfWKagNwjKRjIuITGcdmGeq0qG/p+JivpKwvnGvKp3E6qtX6mk55yGw+2o3gPK/h9rUt7rcB1C7pjEjccu5r+xyRDTHnmhJpNh0laLoNfN/xsb7GZsOl3S6qfhx2ZxlrtbBv1YqDdlqDU3PG0fv1P1AbWs415dJsOipgpyJnbHSEVSsO6mdoNmTaTVF9KiLOlnQ9TYrviHhrppHZvHVa2AfstIvqzKO9a8H6y7mmXFpNOwWe+rb+ajdF9YXku0/1HlCdFvZ5UZ8VhHNNibSa/l46Puapb+urdlNUdyQ/rgdmImIbgKQR4Dl9iM3mqdWVlBf2WZE41wymdtPfjT1tPB1leUjTB+cm4PXAk8ntMeAbwG9kFZR1rzI5zYeu38hjT88C1e6fLxgbZevM7E7P9cI+KyjnmgGwujLF5bc/uEPfrGbT3+5pY3lLU+A8NyJqCYeIeFLS4gxjsi5VJqdZddWdzM79IuNsnZllkWB0kZitW0nsKykrMOeagjvz0lu55b6fNX3M099WNGkKnKckHRER3weQdCTgOY4C6HTq97aAF4ztwuJdd/GVlA0C55oCq0xOtyxuajz9bUWSpsA5G/iKpM3J7X2A0zKLyFJJe+r31qdnmfzgG/oUldm8nI1zTaHUr7NZJHV8vqe/rUjSHNXwPUm/ChxEtZXBPRGx88IO66u0p3474digcK4plsaLqLkOh9V5+tuKZlGrByS9MjnNlyTJHAF8BPi4pD36FJ+1kGYoeHSRnHCs8JxriqcyOc37r7wz1UUUwOLRRT6vzgqnZYED/A3w3wCSXg1cAvwD8DiwNvvQrJ1OIzPjY6OsOfVQJxwbBM41BVIbuek0YlNz7AF78IMPv8m5xgqn3RTVSETUVpSdBqyNiKuBqyVtmM+bSjoVuBB4GXBURKyve+w84A+AOeCPI2LdfN5rkLXqMwG07DXhqygbQM41OWtca9OuuBmR2BbhTQtWeG0LHEm7RMSzwOuAlSlfl8bdwMlUr9y2k/Ry4HTgYGBf4EZJvxIR6cZJS6Kxpw3s3GfCvSasRJxrctTNWhtfRNkgaZc8rgC+LemnVLdqfgdA0kupDh33LCJ+mPyuxodOAL4UEc8AP5Z0L3AUcOt83m+QtNsdVd9nAnCvCSsL55ocpd2wMCK5uLGB0u6ohosk3UR1q+Y3IraX9YuA92YUz1LgtrrbDyX37UTSSpIrvWXLlmUUTv91SjbuM2Fl41yTrzQ5xSM3NojaDv9GxG1N7vv3NL9Y0o3Ai5o8dH5EfLXVy5qF0SK2tSQLECcmJtKthhsAnZKNt31bGTnX9EezdX2tDsf0WhsbdPOd324pIl7fw8seAvaru/1iYHOL55ZSq2QD7jNh1oxzTTqN09+1dX2nHLmUq++Y9oYFK51228TzcB1wuqTnSNofOBD4bs4x9dWqFQcxNjqy0/3jY6NOOGYLZ+hyTbPp75nZOb51zxYuPvkQlo6PIWDp+JhzjZVCZiM47Ug6Cfg0sAS4QdKGiFgRERslXQn8AHgWePew7Wrw7iizheNc8wutpr83b53xhgUrpVwKnIi4Fri2xWMXARf1N6JicbIxWxjDmGta9c9qNf3tdX1WVrkUOGZmtrAqk9NceN1Gts4075/Vqjmo1/VZWbnAMTMbcKsrU1x+24NNt4HV+mfdcu5rAU9/2/BwgWNmNsAqk9Mti5ua2vobT3/bMCnaLiozM+vCmnWb2hY34HU2Npxc4JiZDbBOzUG9zsaGlQscM7MB1m50ZvfF7p9lw8sFjpnZAGvWHFTAWccsY/KDb3BxY0PLi4zNzAaYm4OaNecCx8xswHl3lNnOPEVlZmZmpeMCx8zMzErHBY6ZmZmVjgscMzMzKx0vMm6i1Wm8ZmZmNhhc4NSpTE7zoes38tjTzU/jdZFjZmY2GHKZopJ0qqSNkrZJmqi7f7mkGUkbkq/P9iumyuQ0510ztUNxU1M7jdfMzMwGQ15rcO4GTgZubvLYfRFxWPL1rn4FtGbdJmZm51o+3um8FzMrniJeTJlZf+QyRRURPwSQlMfbN9WpgPFpvGYDqXYx9TdNHrsvIg7rbzhm1i9F3EW1v6RJSd+W9KpWT5K0UtJ6Seu3bNky7zdtV8D4NF6zwRQRP4wIzy+bDaHMChxJN0q6u8nXCW1e9jCwLCIOB84Bvijp+c2eGBFrI2IiIiaWLFky73ibHVgHMD7m03jNSiqXiykz64/Mpqgi4vU9vOYZ4Jnk5zsk3Qf8CrB+vvF02vrtA+vMBpOkG4EXNXno/Ij4aouX1S6m/lPSkUBF0sER8UTjEyNiLbAWYGJiItrF4hYTZsVRqG3ikpYAP4uIOUm/DBwI3D/f31vbIVVbRNxq67cPrDMbPEW5mEqbZ8ysP/LaJn6SpIeAXwdukLQueejVwF2S7gSuAt4VET+b7/s12yHlrd9mw0vSEkkjyc8LcjHlPGNWLHntoroWuLbJ/VcDVy/0+7XaIeWt32blJukk4NPAEqoXUxsiYgXVi6k/k/QsMMcCXEw5z5gVS6GmqLKy7/gY002SjLd+m5VbPy+mnGfMiqWI28QXXLMdUt76bWYLyXnGrFiGYgTHO6TMLGvOM2bFMhQFDniHlJllz3nGrDiGYorKzMzMhosLHDMzMysdRbRtzDkQJG0B/mMBftWewE8X4PcspCLGBMWMyzGl06+YXhIR8z9HpUAWKNcU8W+ipqixOa7uFTW2LOJqmmtKUeAsFEnrI2Ii7zjqFTEmKGZcjimdIsY0TIr8+Rc1NsfVvaLG1s+4PEVlZmZmpeMCx8zMzErHBc6O1uYdQBNFjAmKGZdjSqeIMQ2TIn/+RY3NcXWvqLH1LS6vwTEzM7PS8QiOmZmZlY4LHDMzMysdFziApFMlbZS0TdJE3f3LJc1I2pB8fTbvmJLHzpN0r6RNklb0K6aGGC6UNF332bw5jziSWN6YfBb3Sjo3rzgaSXpA0lTy+azPKYbLJD0q6e66+/aQ9E1JP0q+755HbMOmiHmmXVzJY7nnmrpYCpNzknicdzrHkmv+cYFTdTdwMnBzk8fui4jDkq935R2TpJcDpwMHA28E/krSyM4v74tP1n02X8sjgOTf/pfAm4CXA2ckn1FR/Fby+eTVj+JzVP9O6p0L3BQRBwI3Jbcte0XMMy3jKliuqck954DzThc+R475xwUOEBE/jIhNecdRr01MJwBfiohnIuLHwL3AUf2NrlCOAu6NiPsj4r+BL1H9jAyIiJuBnzXcfQLw+eTnzwMn9jOmYVXEPAPONT1y3kkh7/zjAqez/SVNSvq2pFflHQywFPhJ3e2Hkvvy8B5JdyXDkHlNcxTp82gUwDck3SFpZd7B1Nk7Ih4GSL7vlXM8Vrw8A8X8f6sIOQeK+dnUFDXv1PQt/+yS1S8uGkk3Ai9q8tD5EfHVFi97GFgWEf8p6UigIungiHgix5jU5L5M9vq3iw/4a+DDyXt/GPg48M4s4uigb59HD46NiM2S9gK+Keme5IrGSqqIeWYecfX9/60ByTngvDMQhqbAiYjX9/CaZ4Bnkp/vkHQf8CvAgizc6iUmqlcK+9XdfjGweSHiaZQ2PkmXAv+URQwp9O3z6FZEbE6+PyrpWqrD2kVINI9I2iciHpa0D/Bo3gGVRRHzTK9xkcP/WwOSc8B5Zz76ln88RdWGpCW1RXWSfhk4ELg/36i4Djhd0nMk7Z/E9N1+B5H8YdacRHWhYh6+BxwoaX9Ju1JdFHldTrFsJ2k3Sc+r/Qy8gfw+o0bXAe9Ifn4H0OoK3vqgoHkGCpJragqUc8B5Zz76l38iYui/qP7P8hDVq6hHgHXJ/acAG4E7ge8Dx+cdU/LY+cB9wCbgTTl9Zl8ApoC7kj/YfXL87/dm4N+Tz+T8vP+ekph+Ofm7uTP5G8olLuAKqlMgs8nf0x8AL6S6e+FHyfc98v68huGriHmmXVzJY7nnmrpYCpNzknicdzrHk2v+8VENZmZmVjqeojIzM7PScYFjZmZmpeMCx8zMzErHBY6ZmZmVjgscMzMzKx0XOLYDSXPJKbQbJd0p6RxJPf2dSHoyxXMulBSSXlp33/uS+yaS2w9I2rPN73iNpMeTVvc/lHRB3WNHSbo5OfX3Hkl/K2lx3eNflXRrL/8+M+udc41lzQWONZqJ6im0BwPHUe31cEGH18zXFNVGWTVvA37Q5e/4TkQcDkwAZ0k6UtLewFeAP4mIg4CXAV8Hao2wxoEjgPGkkZmZ9Y9zjWXKBY61FBGPAiupHnAnScslfUfS95Ov34Bqh9HkymWDpLvVcFigpD0l3SrpLS3eqkJyEm/SyfVxYEuPMT8F3AEcALwb+HxE3Jo8FhFxVUQ8kjz9FOB6qicBn97s95lZ9pxrLAsucKytiLif6t/JXlTPDDkuIo4ATgP+Inna71LtfnoYcCiwofb65MrmBuCDEXFDi7d5AviJpFcAZwBf7jVeSS8EjqHaxfMVVBNQK2dQ7bR5RfKzmeXEucYW2tActmnzUjs5dxT4jKTDgDmqBwJC9VyWyySNApWI2FD3/JuAd0fEtzu8R+3KZgXwOuD3u4zxVZImgW3AJRGxUWp24G/yD6omw5cC/xoRIelZSa+IiKKd22I2TJxrbMF4BMfaSoZx56heUb2P6lk1h1Kdf94VICJuBl4NTANfkPT25OXPUr2qWVH3+y5Khpc3NLzV9cDvAQ9GxBM9hPqdiDg8Io6MiM8m920Ejmzx/NOA3YEfS3oAWI6Hjs1y41xjC80FjrUkaQnwWeAzUT207AXAwxGxjWqCqJ2A/BLg0Yi4FPg7qovpAAJ4J/Crks4FiIjzk4WFh9W/V0TMAH8CXLSA/4TPAO+QdHTdv+ksSS+iOkz8xohYHhHLqSYnJx2zHDjXWBY8RWWNxpIrnlGqV0VfAD6RPPZXwNWSTgW+BTyV3P8aYJWkWeBJoHZVRUTMSToduF7SExHxV63eOCK+1CauuyRtS36+MiLO6fQPiYhHkvf+mKS9qA4p30z1xOZlwG11z/2xpCckHR0Rt3f63WY2b841zjWZ8mniZmZmVjqeojIzM7PScYFjZmZmpeMCx8zMzErHBY6ZmZmVjgscMzMzKx0XOGZmZlY6LnDMzMysdP4/gHlHdmVr5rUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots(1, 2)\n",
    "fig.set_size_inches((8, 3))\n",
    "for i in range(2):\n",
    "    axs[i].set_title(f'PC{i+1}')\n",
    "    axs[i].scatter(pcs1[:, i], pcs3[:, i])\n",
    "    axs[i].set_xlabel('Dask-ML PCA')\n",
    "    axs[i].set_ylabel('Scikit-allel PCA')\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.testing.assert_allclose(pcs1[:, 0], pcs3[:, 0], atol=.1)\n",
    "np.testing.assert_allclose(pcs1[:, 1], pcs3[:, 1], atol=.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAADQCAYAAAAK/RswAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAo1UlEQVR4nO3de5SddX3v8fcn46iDpQaWUclIGEQaDzhlRofbSesBK0ShyBbFkEJrT1tSV72cBJt1kooIGJuspgpLqx7jOlZbIg5iiNCgXESkRaMMZiBERIiEwMCB2BAFMoVh8j1/7GfHnZl9efbM7Pvntdas7P3s59n7m1k73/ye3+X7U0RgZmZm1kpm1TsAMzMzs5nmBo6ZmZm1HDdwzMzMrOW4gWNmZmYtxw0cMzMzazlu4JiZmVnLcQPHzMzMWo4bOFYTknZIGpX0rKQnJf2zpN9JXlso6Q5Jz0jaJekHkt6VvHaYpOslPS4pJPXU9S9iZg1rGnnmTEn/IWmPpP8n6cuSDq7v38amyw0cq6WzIuJ3gDcDxwMXS3ov8E3gX4DXAa8BLgHOSq7ZB3wXeE/twzWzJjSVPPNKYBUwF/hvyTlraxy3zTC5krHVgqQdwF9FxK3J87XAMUAv8LmIKJlMJL0EGAOOjIgd1Y3WzJrRdPNM3vucA1wWEb3VitWqzz04VnOSDgfOAPYChwPX1jciM2s108wzbwW2VSMuq52X1DsAaysbJb0I/BrYBKwH3gs8UdeozKyVTCvPSDoNeD9wYtUitJpwA8dqKZPrOgaQ9Mbk4WHAw/UJycxazJTzjKSTgK8D742IX1QvRKsFD1FZPT0APIonEJtZ9aTKM5L6geuBv4iI79UiMKsuN3CsbiI7w/0i4OOS/qek35U0S9IfSFqXO0/Sy4GXJU9fljw3MysrTZ6R9CayqzU/HBE31DNemzlu4FhdRcS1wCLgL4DHgSfJLtf8dt5po8CzyeOfJ8/NzFJJkWc+CswB/m9SQ+dZSZ5k3OS8TNzMzMxajntwzMzMrOW4gWNmZmYtxw0cMzMzazlu4JiZmVnLaYlCf6961auip6en3mGYWZ677777VxExp95xzCTnGrPGUyzXtEQDp6enh6GhoXqHYWZ5JD1S7xhmmnONWeMplms8RGVmZmYtxw0cMzMzazktMUSVxsYtI6y96QEe3zPK3NldLF84n0x/d73DMrMW4jxj1jjaogdn45YRVm7YysieUQIY2TPKyg1b2bhlpN6hmVmLKJRnlg0Oc/HGrfUOzawttUUDZ+1NDzA6Nn7AsdGxcdbe9ECdIjKzVlMozwSwfvNO30yZ1UFbDFE9vqfw3owje0ZZsOY2dyeb2bQVyzMBfPSae1g2OOw8Y1ZDbdGDM3d2V8HjAncnm9mMKJZnAMYjPDxuVmNt0cBZvnA+XZ0dBxwT2TurfO5ONrOpWr5wPkpxnofHzWqjLRo4mf5uVp/TS/fsLgR0z+6a1LjJCXDyMWshkr4i6SlJ9+UdO1TSLZIeTP48ZLqfk+nv5vyT5qVq5BQbzjKzmdMWc3Agm3zyx70XrLmNkSJJxsnHrKV8Ffgn4F/yjq0AvhcRayStSJ7/7+l+0KpMLwNHHLp/qfgsifGYfDs1d3aXl5SbVVnbNHAmWr5wPssGhwv25Dj5mLWOiLhDUs+Ew2cDpySPvwbczgw0cODAm6nc0vH81VVdnR2c+sY5BxzPzc3JXW9m09cWQ1SFFOtOzk8+rptj1rJeExFPACR/vroaH1JoeHz1Ob18/+e7Cpau+Og193Dkik0sWHOb843ZNNW1B0fSV4A/Bp6KiDclxw4FBoEeYAfwvoh4uhqfP7E7OddTU6puju+uzNqLpCXAEoB58+ZVfP3E4XGAZYPDBc/NDWe5R8ds+urdg/NV4B0TjuXGxo8Gvpc8r5pMfzd3rngbD685kztXvI1Mf3fJujm+uzJrCU9KOgwg+fOpYidGxLqIGIiIgTlz5szIh5daUp6T69FxrjGbmro2cCLiDmD3hMNnkx0TJ/kzU8uYoHTyya+Z0+PGjlmzuh54f/L4/cC3a/nhhUpXFDIewdLBYfovv9l5xqxC9e7BKSTV2LikJZKGJA3t2rVrRgNIk3xyk5NH9oyydHCYvsucgMwakaSrgR8B8yU9JukvgTXAaZIeBE5LntfMxLk5HSq9uPzpvWPOM2YVatpVVBGxDlgHMDAwUKyszZTkxrxzc3PSvPme0TGWX3vPAdebWf1FxOIiL/1RTQOZoNxqq0L2jI55bo5ZSo3Yg5N6bLya8ufmdKcYLwcYGw+WXTPsOywzq0iuR6dcTw64ErJZWo3YwKnr2HghacfLASLwknIzq1imv5tPv++4VLnGxUjNyqtrA6cRx8YLyR8vB8qWYnc9CzObilyumd3VWfK8NKuwzNqdokAZ8WYzMDAQQ0NDNfu8jVtGWHbNMJX86g45qJNPnHWsx82tbUi6OyIG6h3HTKplrtm4ZYTLbtjG03vHDjje1dnB6nN6nUvMEsVyTdNOMq6nXGJZ/s17GNuXrpXz9N4xlg0OM/TIblZleqsZnpm1gNwk5DTbxnhrGbPJ3MCZolzyKHSHVUwAV23eCeBGjpmlUqgScr6JK7BypSsuvX4bl77LvcbWvhpxknHTyPR3s+WS07lyUV/qehaQbeQc8/HveG6OmU1boa1lILukfNngMBdv3FqHqMzqzz04M2Aq9Sz2ju3jomuG919vZjYVpVZUBbB+804GjjjUecbajntwZljaVRAA+yI7xGVmNlXlVlQFsDTZWubIFZvco2Ntww2cKsj0dzP8idO54KTyOw+nnb9jZlZIRXW6yA6Rn//lH1U3KLMG4AZOFa3K9KZq5JiZTVWu1/iQg8r3GufcuX235wBay3MDp8rKNXJyQ1kbt4ywYM1tLg5oZhXLLXi44KR5ZQuR5ni7B2t1buDUQLFGTucscem7jt0/MXkk2dhzZM+oVz+YWcVWZXq5YlFfqtWcI3tGfUNlLc0NnBpZlek9YDl59+wu1p57HJn+7oLLPHOrH5x4zKwSuT2tOmeVb+Tkbqi8f561Ii8Tr6FiBbuKLfMMst3IXt5pZpXI5YxLr9/GntHyCxlyO5Q711grcQOnAcyd3cVIkUaOdw02s6mYWJ8rt5VDsc1lnGus1biB0wCWL5zPssHhgoknV+PCe82YzTxJO4BngHHgxVbbHDQnv7GzYM1tBW+o5s7ucp6xluI5OA0g09/N+QVWP3R1drB84fyCk5CXDg5z7CXf9bi52fSdGhF9rdq4mahQ3Zyuzg5OfeOcgnnmyJUuDmjNqWEbOJJ2SNoqaVjSUL3jqbbc6of8Scirz+ktOgkZ4LkXxll+7T1u5JhZarm6ORNzzfd/vqtgnonIFgc87TO31zxWs+lQRLER2RIXScdHxF1ViCf/M3YAAxHxq3LnDgwMxNBQ67aBjlyxqei4OWQT1J0r3lazeMzSkHT3dHtFqp1rJD0MPE12Tv+XImJdgXOWAEsA5s2b95ZHHnmkWuHUVbk8A3DBSfNYlemtSTxmaRXLNal7cCQdI+lySQ8CX5zR6KykcnvNeHKgtZIa55oFEfFm4J3AByW9deIJEbEuIgYiYmDOnDlVDqd+yuUZyPbkuHaONYuSDRxJR0haIeke4F+BvwFOq9FYdQA3S7o7uYNqW8sXzi9ZnTRNYjJrZPXKNRHxePLnU8B1wAnV/LxGlnZPq/z5OT3evNMaWNEGjqQfAjcCncB7I+ItwDMRsaNGsZW8s5K0RNKQpKFdu3bVKKT6yE1CLqSzQyxfOL/GEZnNnHrlGkmvkHRw7jFwOnBfNT+zkeXm5hzUWdnUTG/eaY2q1Dd5F3Aw8Bog1y9b+YSdKSp3Z9Uu3cY5uUrIub2rAA45qJO17z3Oyzit2dUr17wG+I+k1+gnwKaI+G4NPrdhZfq7+dkn38nRr35FRdfduX23e3Ks4ZScZCzplcB7gMXAG4DZwMKI+ElVg8reTc2KiGeSx7cAlxdLPq0+yXgqLt64lat//CjjEXRILD7xcE8OtJqqZJJxvXJNpdop1+TnkLRmd3Vy6buO9U2X1VSxXJN6FZWkVwOLyCagwyPi8JkN8YDPej3ZXhvIFiP8ekR8qtj57ZR00rh441au2rxz0nEB53sVhNXIVFdR1TLXVKpdc83GLSMsHRxOda7zjNVaxauoJL1c0v6xn4h4KiI+B7yb7Fh11UTELyPiuOTn2FKNG5vs6h8/WvB4kB0vP2rlJq+AsIZRz1xj6WT6u1lw1KGpzs3lGc/LsXorNQfns8AfFjj+duAj1QnHZkK5LuXxgI9+0wUCrWE41zSB9ReenLqRA9l5OT0rNrniutVNqQbOH0TEhokHI2I9MKlWhDWODpVaVJ41vi9Ye9MDNYjGrCznmiax/sKTuXJRX0UrrZ57YZylg8Pu0bGaK/UtLfW/ZMNu8WCw+MR0UxZG9oxy1Mob6VmxiaNW3uhVEFYvzjVNJLfS6spFfRxyUGf5CxJ3bt/tRo7VVKnk8ZSkSUWvJB1PdlmnNahVmV4uKFI3Z6LccNZ4BFdt3ulGjtWDc00TyvR3s+WS01PnGvjtsFX/5Td72MqqrugqqiThXAN8Fbg7OTwA/BlwXkT8uBYBptGuKxvSOO0zt/PgU89VdE337C6WL5zvpZ42LWlXUTnXNL/zv/wj7ty+u6JrOjvkOl42IypeRZXUnziBbPfxnyc/Ak5spIRjpd1y0SlccNK8A8YAXvHS0uXYc2XY3Z1steBc0/zWX3hyRT05AGPjngdo1VWu0F8/cBSwLSLur1lUFfJdVeWOWnljqgJe3j3YpqrCQn/ONS2ikt4cAQ+vObO6AVnLm0odnEuAQbLVRTdJurCK8VmNpZ2I7N2Drdqca1pLbqVVmoVWr+zqpP/ym+lZsYmeFZvou8xzc2zmlPoKLgL6ImIxcDzQ1jt6t5rcROQ0S8pzuwcvd+0cqw7nmhaT6e/mwb8/s+RKq1mCZ55/kaf3ju0/tmd0jKWDw17sYDOiVAPnvyJiL0BE/GeZc60Jrcr0sn31GVy5qC/V+WP7gpUb7q1uUNaOnGtaVG6lVaGNgl/Z1cn4vsLD5Fdt3ukCgTZtLynx2lGSrk8ea8JzIuJdVY3MaibT3803h3amGjcfHdtHz4pNACw46lDWX3hytcOz1udc0+Iy/d2TVksdmeSRYp57YZzl196z/3qzSpVq4Jw94fk/VjMQq6/1F55c8e7Bd27fzWmfuZ1bLjqlusFZq3OuaUNzZ3cxsme05Dm5lVZu4NhUFG3gRMQPahmI1d+qTO/+FVP9l998wNh4MQ8+9Rx9l93Mpe861knIpsS5pj0tXzif5dfew9h46Ruqx8s0gsyK8Vi3FfSJs45NfW5uYmDPik2unWNmqWT6u1n73uPK1uWaO7urRhFZq3EDxwrK9HdPKhCYxp3bd3Pip26pSkxmM03SOyQ9IOkhSSvqHU+7yfR3s+3ydxQtEtjZIZYvnA/Axi0jLFhzm8tWWGoN28Bx4qm/VZlerqhw52CAJ595wcs8reFJ6gA+D7wTOAZYLOmY+kbVnlZleguutMpt5bBxywgrN2xlZM/o/rIVSweHef3KTc41VlSpvahuIFsCpaBqrmxIEs8vgNOAx4C7gMUR8bNC57u6aPVdvHErV23emfr8Dontq8+oYkTW6CrYi6ouuUbSycClEbEweb4y+bzVxa5xrqmPBWtuKzkhuXMWrD23z/MA21SxXFNqFVU9VzKcADwUEb8EkPQNsistCjZwrPpyE5A3bhnhshu2lZ2AnHYllhn1yzXdwKN5zx8DTpx4kqQlJMUH582rbL8lmxnlJhqP7YOLBocBLym330q1ikpSFzAvImq1M1qqxGO1l6tnUW6/mQ7pgGXnHRKLTzzc+1rZJHXMNYWmmE1qmUfEOmAdZHtwqh2UTZZmSfk+YOngMEsHh12jy4AUc3AknQUMA99NnvflF+GqkrKJR9ISSUOShnbt2lXlcGyicrsHv37OQVy1eef+npzxCK7avJMeTxC0IuqQax4D8jdlex3weBU/z6Zo+cL5dHWWXm2Vz4sdDNJNMr6U7JDRHoCIGAZ6qhVQomziiYh1ETEQEQNz5sypcjhWyKpMLzvWnHnAnlYdEhecNI9f7tpb9LqRPaMs834zNtml1DbX3AUcLelISS8FzgOqffNmU5Dp72b1Ob10VbDgwYsdLM235cWI+HXVIzmQE08Tye1ptWPNmWxffQarMr1l5+AE2f1m3JNjeWqaayLiReBDwE3A/cA1EbGtVp9vlcn0d3P/J9/JgqMOTX3N1T9+tPxJ1rJKTTLOuU/SnwAdko4GPgL8sJpBRcSLknKJpwP4ihNPc+mQUk00Xjo4zNAjuz03x6A+ueZG4MZqfobNrNy2MmlWdY5H7N87b5bgT06c51zTRtL04HwYOBZ4Hrga+A2wtIoxAdnEExG/FxFHRcSnqv15NrMWn3h4+ZMSV23eyTEf/457c6wuucaaT65uTncFVY73RTbXnPaZ26sXmDWUonVwmolrUzSmSmvnABz96ld4884WkbYOTjNxrmlMJ37qFp585oXU5x/UOYu/P+f3vaS8RRTLNQ1Z6K9STjqNK1eBdHRsPPU1XuLZGhq90N9UONc0rvyyFGnN7ur0RsEtoNkK/VmLyCWONMUBc+7cvpueFZtcO6d9ONfYtOWKkUL5ysc5e0bHWLkhu9LKjZzW4yEqq5lsb869jI7tq+g6Aeef5MmBzcZDVFYvG7eMcNHgMGkzTYfEp993nBs5TWoqQ1RbKd1t/PszF970OOk0l4s3bmX9j3dSadva83OaSwVDVM41NuMqbeTkeNiq+UylgXNEqTeMiEdmKLZpc9JpTqd95nYefOq5iq+7wL05TaGCBo5zjVXNVG6oOmeJtee6R6dZVNzAaSZOOs2r3J5WxbiR0/g8RGWNJO1GwfmcZ5pDsVyTZi+qkyTdJelZSS9IGpf0m+qEae1m/YUns2PNmfu3fEjrqs07vadVi3GusWrK9Hez5ZLTuXJR3/6tZcq5avNOb/fQxNIU+vsnYDHwINAF/BXwuWoGZe1pVaa3okaO97RqOc41VnWZ/m4+/b7jUm/emdsk+MgVm5xrmkyqncsi4iGgIyLGI+KfgVOrG5a1q9wGnke/+hWpzg9gvfe0ahnONVYLuc07Z3d1pr4mt3+eKyE3jzQNnL3JhpfDkv5B0jIg3f8+ZlN0y0WnpN5UL8juafX6lb7DanLONVYzmf5uhj+RHbJKOWIFwINPPec80yTSNHD+NDnvQ8BzwOHAe6oZlBlk5+dUMmSV22vGyadpOddYzWX6u7nifX3MqqCRkxu26r/8ZvceN7Cyu4nnLdH8L0k3RMRPqxyT2X6rMr0MHHEoa296IFVlUsgmn6s27/SWD03GucbqJbcc/O823MveCgqRPr13jOXX3nPAe1jjqGiZuKSfRsSbqxjPlHjpZvu4eONW1m/eWbwqXAFe6lkf01km7lxj9VTpRsGuhFxfU14mPvF9ZigesylZlenligqWeUK2R+f8L/+oilFZFVQ110i6VNKIpOHk54xqfp41l0pXdI5HsHRw2CutGkylDZzLqhJFHiceKye3zLOSL++d23d7vLy5VD3XAFdERF/yc2MNPs+aSG5F5wUnzUt9Q5VbaeVGTmNIU+jve7nHEbFx4rEqceKxkjL93XxmUR9dnembOU/vHXPdnAZWp1xjVtKqTC/bV5/BlYv66Ew5E9mFSBtD0f8dJL1c0qHAqyQdIunQ5KcHmFuzCM2KyPR3c/8n38mONWdWtKT8qs073ZvTQOqYaz4k6V5JX5F0SIn4lkgakjS0a9euKoZjjSzT383ac49LXTvHhUjrr9Tt718DdwNvBH6aPL4b+Dbw+SrHlSrxmOWsv/Dk1I0cyPbmLB0cpu8yN3QaQFVyjaRbJd1X4Ods4IvAUUAf8ATw6WLvExHrImIgIgbmzJkz1XCsBeTXzklTCdlDVvVVdhWVpA9HxIyWS5d0K/DaAi99DNgM/Irsd+OTwGER8RcF3mMJsARg3rx5b3nkkYbZcNjqrNINPLs6O1h9Tq9XQMywSldRVSPXpPzcHuDfIuJN5c71KirL2bhlhI9dt5XnXhgve66AKxb1OcdUScW7iUt6W0TcJumcQq9HxIYZjrFQDD2kSDxOOjbRxi0jXHr9NvaMpts52Ms8Z17aBk49co2kwyLiieTxMuDEiDiv3HXONTbRxRu3cvWPH2W8TGdB9+wu7lzxthpF1V6K5ZpShf7+B3AbcFaB1wKoSgMnP/EA7wbuq8bnWGvL9HeT6e9OXTdnPIKVG7buv9Zqqh655h8k9SXvv4PsMJlZxVZlelmV6WXjlhGWDQ4XzTWPpyxUajOnaAMnIj6RPPxARDyf/1oyIbBanHhsxuQqIafpzRkdG2ftTQ+4gVNj9cg1EfGn1Xhfa1+Z/m6GHtldtEDg3NldNY7Iym7VAGyQdHZEvAgg6bXAJuAt1QjIicdmWq43Z+OWES67YRtP7y3e0Hl8z+ikOTze8qFmapprzGZarmL6xF7jrs4Oli+cv//5xi0jrL3pAR7fM8rc2V0sXzjfN1ZVkKaIyEbgWkkdyZyYm4GV1QzKrBoy/d1suSS7AqJY4a6XvmTWpAnKd27f7UrItbER5xprcrlq692zuxDZuTf5ixg2bhlh5YatjOwZJcguJ1+5YatXc1ZBms02vyzppWSTTw/w1xHxwyrHZVY1uUSzcsNWRsd+uwKiq7PjgOf5KlmVZVPjXGOtItdrXMjamx6YlGc8PF4dRRs4ki7KfwocDgwDJ0k6KSI+U+XYzKoml0gmdhMvHRyub2BtyLnG2kmxycYje0ZZsOY2D1vNoFI9OAdPeH5dkeNmTanQXZYbOHXhXGNtY+7sLkYKNHIE+4/nhq3Aqzqno9QqqlpsdmfWUBYcdWjB4aj8KsmeIDiznGusnSxfOH/S8Lhg0vJyD1tNX6khqisjYqmkG5j8uyci3lXVyMzqYP2FJ5dcRZWbIJhLTr7Tmj7nGmsnhYbHC/XogGvnTFepIap/Tf78x1oEYtYoSi0J9wTBqnCusbYycXh8wZrbCjZyXDtnekoNUd2dPBwCRiNiH4CkDuBlNYjNrOEUu6PyndbUOddYuys0bOXaOdOXpg7O94CD8p53AbdWJxyzxlbsjurlnbM4auWN9KzYxFErb/TuwVPjXGNtKdPfzepzeiuqnbN0cJgjV25yrikhTSXjl0fEs7knEfGspINKXWDWqgrdac0CRsf27X8+HrG/XHuusqml4lxjbavS2jkAEXDV5p08vOtZV1svIE0PznOS3px7IuktgPvjrS0VutOicFFkrtq8kyNXbGLBmttcpTQd5xqzAsoNgd+5fbdzTAFpenCWAt+U9Hjy/DBgUdUiMmtwE++0elZsKnpufin23LVW1FKca8wmKbXSKmfp4DBrb3rAc3PypNmq4S5JbwTmk71X/XlElN6W2ayNdEiMx6TVzQfwSqvynGvMCis0NF6Ib6YOVHSIStLxyW6+JEnmzcAq4NOSDi12nVm7WXzi4anO80qrwpxrzErLDY0f1Fl+Vsno2DhLB4c9NE7pOThfAl4AkPRWYA3wL8CvgXXVD82sOazK9HLBSfOK7lCe45oWRTnXmJWR6e/mZ5985wFV1UvxLuWlGzgdEZEr57oIWBcR34qIjwNvmM6HSjpX0jZJ+yQNTHhtpaSHJD0gaeF0PsesVlZletm++gx2rDmTKxf10dXZccDrhWpaLFhzmychZznXmKW0/sKTuXJRX3aBQxnt3ptTsoEjKTdH54+A2/JeSzM5uZT7gHOAO/IPSjoGOA84FngH8IWk2JdZ05hKTYtlg8PtXM/CucasApn+bu5c8baCN1OFtGtvTqnkcTXwA0m/IrtU898BJL2BbNfxlEXE/cl7TXzpbOAbEfE88LCkh4ATgB9N5/PMaq3SmhYBrN+8k4EjDm3HyYHONWZTkL+vVblVVu240KFoD05EfAr4KPBV4A8i9i8TmQV8uErxdAOP5j1/LDk2iaQlkoYkDe3atatK4ZjNvGKTjYNsomo3zjVmU1dJb067LXQo2f0bEZsLHPtFmjeWdCvw2gIvfSwivl3sskJhFIltHckExIGBgdJrdM0aiHcPnsy5xmx60vTmzJ3d1VZ7Wk13fLuoiHj7FC57DMhfc/s64PEi55o1peUL57NscLjg/6ZeaVU55xqzrNzQeG6e38TNO09945wDjrd63Zw0WzXU0vXAeZJeJulI4GjgJ3WOyWxGZfq7Of+keZO6ECautLKqcq6xllVsocP3f75r0vy/3NycVlS1HpxSJL0b+BwwB9gkaTgiFkbENknXAD8DXgQ+GBGlSzeaNaFVmV4Gjji0bbqK68W5xtpVoYUOywaHC57bqkPjdWngRMR1wHVFXvsU8KnaRmRWe6VWWtnMcK4x+61i8/9yQ+OtNj+n0YaozMzMrAqWL5xftAhpofpcSweH6bvs5qatn+MGjpmZWRsoVYS0UH0ugD2jY01biLQuQ1RmZmZWe8WGxkvNw2nWQqTuwTEzM2tz5UpUNGMhUjdwzMzM2lyh+TkTNdtqKw9RmZmZtbnc0NNlN2zj6b1jBc9ptkKk7sExMzMzMv3dbLnkdC5okUKkbuCYmZnZfqsyvVyxqK/gaqtm4iEqMzMzO0ArFCJ1D46ZmZm1HDdwzMzMrOW4gWNmZmYtxw0cMzMzazmeZGzWplpt52Azs3zuwTFrQ8V2Du6/vHl3DjYzy1eXBo6kcyVtk7RP0kDe8R5Jo5KGk5//U4/4zFpdsZ2Dn947xsoNW93IMbOmV68enPuAc4A7Cry2PSL6kp8P1Dgus7ZQak+Z0bHxpttUrxjfTJm1r7rMwYmI+wGkicWgzawW5s7uYqREI6fZNtUrIXcz9aUCr22PiL7ahmNmtdKIc3COlLRF0g8k/WGxkyQtkTQkaWjXrl21jM+s6ZXbObjZNtUrJiLuj4jW6I4ys4pUrYEj6VZJ9xX4ObvEZU8A8yKiH7gI+Lqk3y10YkSsi4iBiBiYM2dONf4KZi0r09/N6nN6md3VOem1ZtxUb4p8M2XWwqo2RBURb5/CNc8DzyeP75a0Hfg9YGiGwzNre7m9ZtIuF2/UZeWSbgVeW+Clj0XEt4tclruZ+k9JbwE2Sjo2In4z8cSIWAesAxgYGIiZitvMJpvJPNNQdXAkzQF2R8S4pNcDRwO/rHNYZi0tzaZ6uWXluZVXI3tGWblh6/7r68k3U2atYabzTL2Wib9b0mPAycAmSTclL70VuFfSPcC1wAciYnc9YjSz3yq0rLyZV1tJmiOpI3nsmymzBjDTeaZeq6iuA64rcPxbwLdqH5GZlVJsVVWjr7aS9G7gc8AcsjdTwxGxkOzN1OWSXgTG8c2UWd3NdJ5pqCEqM2tMxZaVN/pqK99MmTWPmc4zjbhM3MwaTKFl5W202srMamCm84x7cMysrNwEv0ZcRWVmrWGm84wbOGaWSprVVmZm0zGTecZDVGZmZtZy3MAxMzOzlqOI5i/MKWkX8EgNPupVwK9q8DmVcEzpNWJcrRzTERHRUvuoVCnXNOJ3ABzXVDRqbK0eV8Fc0xINnFqRNBQRA/WOI59jSq8R43JM1qi/b8dVuUaNrV3j8hCVmZmZtRw3cMzMzKzluIFTmXX1DqAAx5ReI8blmKxRf9+Oq3KNGltbxuU5OGZmZtZy3INjZmZmLccNHDMzM2s5buCUIelcSdsk7ZM0MOG1lZIekvSApIV1jPFSSSOShpOfM+oYyzuS38dDklbUK458knZI2pr8bobqGMdXJD0l6b68Y4dKukXSg8mfhzRATA3zfWpVzZBXklga6rvQiPkFGifHJLE0XJ4pEVdVv19u4JR3H3AOcEf+QUnHAOcBxwLvAL4gqWPy5TVzRUT0JT831iOA5O//eeCdwDHA4uT31AhOTX439awF8VWy35V8K4DvRcTRwPeS5/WOCRrg+9TimiWvQIN8Fxo8v0Bj5BhozDwDdcg1buCUERH3R8QDBV46G/hGRDwfEQ8DDwEn1Da6hnMC8FBE/DIiXgC+Qfb3ZEBE3AHsnnD4bOBryeOvAZkGiMmqzHllSpxfUmjEPAP1yTVu4ExdN/Bo3vPHkmP18iFJ9ybdgDXvfkw02u8kJ4CbJd0taUm9g5ngNRHxBEDy56vrHE9OI3yf2lEj/htqlO9CI/5ucho5x0Dj5hmo4vfLDRxA0q2S7ivwU+ruQAWOVW3NfZkYvwgcBfQBTwCfrlYc5cIscKwR6hAsiIg3k+3a/qCkt9Y7oAbXKN+nptYMeQWaJrdA4+YXcI6Zqqp+v14yk2/WrCLi7VO47DHg8LznrwMen5mIJksbo6QvA/9WrTjKqOnvJK2IeDz58ylJ15Ht6r6j9FU186SkwyLiCUmHAU/VO6CIeDL3uM7fp6bWDHkFmia3QIPmF2j4HAMNmGeg+rnGPThTdz1wnqSXSToSOBr4ST0CSb6wOe8mO4GxHu4CjpZ0pKSXkp0seX2dYgFA0iskHZx7DJxO/X4/hVwPvD95/H7g23WMBWio71M7api8Ag33XWi4/AJNkWOgAfMMVP/75R6cMiS9G/gcMAfYJGk4IhZGxDZJ1wA/A14EPhgR43UK8x8k9ZHtrt0B/HU9goiIFyV9CLgJ6AC+EhHb6hFLntcA10mC7Pf96xHx3XoEIulq4BTgVZIeAz4BrAGukfSXwE7g3AaI6ZRG+D61sibJK9AguQUaNr9AA+UYaMw8UyKuquYab9VgZmZmLcdDVGZmZtZy3MAxMzOzluMGjpmZmbUcN3DMzMys5biBY2ZmZi3HDZw2Jeljyu5mfG+yi+uJJc4dkPTZ5PGlkv62wDmXS3p78nippIOKvNcpkmpaLEzSVyU9nPw9fyrp5LzX/lbSz5PKrfdI+rO81+ZIGpPkZdJmU+A8s/8155k6cAOnDSX/8P4YeHNE/D7wdg7c4+UAETEUER8p9Z4RcUlE3Jo8XQoUTDwzRZXvsLw8IvrI7qL7peQ9PgCcBpwQEW8C3sqB5eDPBTYDi6cdsFmbcZ5xnqk3N3Da02HAryLieYCI+FWu1Lik4yX9MLnL+Imkg4vdDUm6UNJ3JHUldy/vlfQRYC7wfUnfLxVEUgH0K5LukrRFyR49knok/XtyF/RTSf89OX6KpO9L+jqwNXl+u6Rrk7uj9UqqbZVwB/CG5PHfAX8TEb9Jfg+/joiv5Z27GPgo8DpJjbKpn1mzcJ7Jcp6pEzdw2tPNwOGSfiHpC5L+B4Cy5c8Hgf8VEceRveMaLfQGylYUPQvIRMT+cyLis2T3hzk1Ik4tE8fHgNsi4njgVGCtsqXOnwJOSzavWwR8Nu+aE4CPRcQxyfN+sndyxwCvBxaU+cyzyCatg4GDI2J7kb/f4cBrI+InwDVJHGaWnvOM80xduYHThiLiWeAtwBJgFzAo6c+B+cATEXFXct5vIuLFAm/xp2R3zX1P7u5sik4HVkgaBm4HXg7MAzqBL0vaCnyTbFLJ+UlEPDzh+WMRsQ8YBnqKfNba5HOWAH9Jtou4VBnv88gmHIBv4O5js4o4zzjP1Jv3ompTyf42twO3J//A3w/8lNL/GHPuI7u9/euAh0udqOyeO59Inv7VxJfJJq8HJlxzKfAkcBzZRvh/5b383IT3yE984xT/Ti+PiGsnfM5zkl4fEb8scP5i4DWSzk+ez5V0dEQ8WOT9zWwC5xnnmXpyD04bkjRf0tF5h/qAR4Cfk/0Hdnxy3sGSCv1D3kJ2U7TrJc0t8PozwMEAEXFdRPQlP0MTzrsJ+HBuPFtSf3L8lWTv8PaRvYurdKJfWquBz0v63eTzf1fSEknzgVdERHdE9ERET3LueVWKw6zlOM/s5zxTJ27gtKffAb4m6WeS7iXbNXtpRLxAdgz4c5LuAW4h2507SUT8B/C3ZHdCftWEl9cB3yk3+Q/4JNlu4nsl3Zc8B/gC8H5Jm4HfY/Ld1Ez5IvB94K7k838A7CV7V3XdhHO/hbuPzSrhPJPlPFMn3k3czMzMWo57cMzMzKzluIFjZmZmLccNHDMzM2s5buCYmZlZy3EDx8zMzFqOGzhmZmbWctzAMTMzs5bz/wEqa2TVdfu2PQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots(1, 2)\n",
    "fig.set_size_inches((8, 3))\n",
    "for i in range(2):\n",
    "    axs[i].set_title(f'PC{i+1}')\n",
    "    axs[i].scatter(pcs2[:, i], pcs3[:, i])\n",
    "    axs[i].set_xlabel('Scikit-learn PCA')\n",
    "    axs[i].set_ylabel('Scikit-allel PCA')\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.testing.assert_allclose(pcs2[:, 0], -pcs3[:, 0], atol=.1)\n",
    "np.testing.assert_allclose(pcs2[:, 1], -pcs3[:, 1], atol=.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Randomized PCA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(100, 2)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Run Dask-ML Randomized PCA\n",
    "pca = daPCA(n_components=2, svd_solver='randomized', random_state=0, iterated_power=25)\n",
    "pcs1 = pca.fit(gnr.data.T).transform(gnr.data.T).compute()\n",
    "pcs1.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(100, 2)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Run Scikit-allel Randomized PCA\n",
    "pcs2 = allel.randomized_pca(gn, n_components=2, scaler='patterson', ploidy=2, random_state=0, iterated_power=25)[0]\n",
    "pcs2.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAADQCAYAAAAK/RswAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAoDElEQVR4nO3df5xcdX3v8dcn64IbSt1QApKVkBRplBjZ1S2JN63FHxDQImuABoTWW3uNfTzUPgLcrcmFmoCh5JoWaLVaQ6VaSSGAcU0ar+FHWrGU0G7cDSHVrSAQmFCIDQtCtrBsPvePOSeZ3cyPM7t75pw5834+HvvIzJkzO5+dzH72c74/zd0RERERyZIpSQcgIiIiMtlU4IiIiEjmqMARERGRzFGBIyIiIpmjAkdEREQyRwWOiIiIZI4KHBEREckcFThSE2b2pJkNmdnLZvacmf2tmf1S8NgiM3vAzH5hZvvM7Adm9pHgsZPMbJOZ7TUzN7NZif4gIpJaE8gzHzazfzazQTP7TzO7xcyOTfankYlSgSO1dL67/xLwLuDXgWvM7CLgLuDvgLcAJwKfB84PnnMQ+D5wYe3DFZE6NJ488yZgNTADeHtwztoaxy2TzLSSsdSCmT0J/C93vy+4vxY4HZgHfMndyyYTM3sDMAzMdvcn441WROrRRPNMwfdZDFzr7vPiilXipxYcqTkzOxn4EHAAOBm4O9mIRCRrJphn3gvsjiMuqZ03JB2ANJQeM3sdeBHYAqwHLgKeTTQqEcmSCeUZMzsb+DgwP7YIpSZU4EgtdYVNxwBm9rbg5knAE8mEJCIZM+48Y2YLgL8HLnL3/4gvRKkFdVFJkgaAp9EAYhGJT6Q8Y2YdwCbgE+5+fy0Ck3ipwJHEeH6E+5XAn5jZ75vZL5vZFDP7DTNbF55nZm8Ejg7uHh3cFxGpKEqeMbN3kJ+t+Vl335xkvDJ5VOBIotz9bmAJ8AlgL/Ac+ema3y04bQh4Obj9k+C+iEgkEfLMVcB04OvBGjovm5kGGdc5TRMXERGRzFELjoiIiGSOChwRERHJHBU4IiIikjkqcERERCRzMrHQ3/HHH++zZs1KOgwRKbBjx46fu/v0pOOYTMo1IulTKtdkosCZNWsWvb29SYchIgXM7KmkY5hsyjUi6VMq16iLSkRERDJHBY6IiIhkTia6qEQkfj19OdZuHWDv4BAzWlvoXjSHro62pMNKFb1HIumhAkdEKurpy7Fi4y6GhkcAyA0OsWLjLgD9AQ/oPRJJF3VRiUhFa7cOHPrDHRoaHmHt1oGEIkofvUci6aICR0Qq2jtYfH/TUscbkd4jkXRRgSMiFc1obanqeCPSeySSLokWOGZ2q5k9b2aPFhxbZWY5M+sPvj6UZIwiAt2L5tDS3DTqWEtzE92L5iQUUXS1yjP1/B6JZFHSLTjfAM4tcvwmd28Pvr5X45hEZIyujjZuWDyPttYWDGhrbeGGxfPqZfDsN6hBnqnz90gkcxKdReXuD5jZrCRjEGl0Uac2d3W01eUf61rmmXp9j0SyKOkWnFI+Y2aPBE3L04qdYGZLzazXzHr37dtX6/hE6l5PX472a+9h2YZ+coNDOIenNvf05ZIOrxYq5hlQrhGpV2kscL4KnAq0A88Cf17sJHdf5+6d7t45fXqm9vMTiV24Zsvg0PARjzXI1OZIeQaUa0TqVeoKHHd/zt1H3P0gcAtwZtIxiWRNsTVbCmV9arPyjEj2pa7AMbOTCu5+FHi01LkiMj6VCpisT21WnhHJvkQHGZvZ7cBZwPFm9gywEjjLzNoBB54EPpVUfCJZNaO1hVyJIidrU5uVZ0QaU9KzqC4tcvjrNQ9EpMF0L5ozat+k0LSpzaw8f26mZgIpz4g0Jm22KdKAwgJGO1+LSFapwBFpUFqzRUSyLHWDjEVEREQmSgWOiIiIZI4KHBEREckcFTgiIiKSORpkLCJS56JumCrSSNSCIyJSx8J9xQo3TL1iQz/X9OxKOjSRRKnAERGpY8X2FXPgtu176LjunkbZGV7kCCpwRETqWLl9xV44MMyKjbtU5EhDUoEjIlLHKm2MOjQ8wtqtAzWKRiQ9VOCIiNSx7kVzsArnVNo9XiSLNItKRKSOdXW00fvUftZv34OXOCds5dFsK2kkasEREalzq7vmcdOSdlpbmo94rKW5ie5Fc4rOttL4HMkyFTgiCenpy7FwzTZmL9/CwjXb9IdGJqSro43+ledw85J22lpbMKCttYUbFs+jq6Ot6Gwrjc+RLEu0i8rMbgV+G3je3d8RHDsO2ADMAp4EfsfdX0gqRpE4XNOza1SXQng1DajLYJI1Wp4ptUt8qXE4eweH1HUlmZR0C843gHPHHFsO3O/upwH3B/dFMqOnL1d0vISupmPzDZRnSs62ap3arK4ryaRECxx3fwDYP+bwBcA3g9vfBLpqGZNI3NZuHSg5GFSzXSaf8kxe96I5tDQ3jTrW0tyEO+q6kkxKugWnmBPd/VmA4N8Tip1kZkvNrNfMevft21fTAEUmolwRU2lNE5k0kfIMZCfXdHW0ccPieUeMz3lxaLjo+Sq2pd7V7TRxd18HrAPo7OwsdUEskjozWlvIFfnjYeSvsiVdspRrio3PWbt1oOjncYoZs5ZvocmMEXfaNDZH6sy4WnDM7NcnO5ACz5nZScHrnAQ8H+NricSi3AypYl0FBly2YKb+eIwRY65RngkU+zwCjLiP+ldjc6TeRC5wzOx0M7vOzH4KfDXGmDYBHw9ufxz4boyvJTLpKq03Uqyr4KYl7azumpdo3GlRo1yjPBMY+3lsstLrImtsjtQTcy/d4mpmpwCXBl+vA6cAne7+5KS8uNntwFnA8cBzwEqgB7gTmAnsAS5297EDBEfp7Oz03t7eyQhJZFwKp9lOCZr0x2prbeHB5e9PILpkmNkOd++MeG5suWay8gw0Rq6ZvXxLyUHwY02b2szK8+eq5VESVSrXlByDY2b/ArwJuAO4yN1/amZPTFZxA+Dul5Z46AOT9Roicevpy9F9906GR0Y36Y+lQZvFxZ1rlGeqU2qMWDEvHBjmqrt2Alq/SdKnXBfVPuBY4ERgenCsrgfYicTh2s27DxU35WiGVEnKNSlSakxOKSMHnWs3744xIpHxKVnguPsFwDzgR8C1ZvYEMM3MzqxVcCL14IUDxafZFgr3A5IjKdekS+GYHCg/Jif0woFhDT6W1Ck7TdzdXwRuBW41sxOAJcDNZnayu59ciwBF0mbssvblGGjp+wiUa9Kl2HTyhWu2le260lYjkjblxuC8ETjW3fcBuPvzwJfM7E7guBrFJ5Iq4QypcOXXcgm/taWZ/pXn1Cq0uqVcUx+6F81h2Yb+ko+HM6y6Otq0t5WkQrkxOH8J/GaR4x8E/iiecETSrdiOzMU0TzFWfWRuDSLKBOWaOtDV0cblC2aWPSfcuFN7W0kalCtwfsPdN4496O7rgffGF5JIepWbCVW4rs3ai8/QFWt0yjV1YnXXPG5e0l5yXM6M1paiFwFaP0eSUG4MTrmRZWncw0okdqWm0DbaGjeTTLmmjoSFe2FXLRweSH9FiW4sLZMgtVYueTxfbBZDsHR6/e44J1JGT1+OuZ//PrOWb2HW8i3MXr6Fa3p2HXq81I7MmiE1Ico1dabUxp1dHW0lB95rmQSptXItON3AnWb2DWBHcKwT+D3gkpjjEqm5nr4cV921k5GDh5dgceC27XuAfPN8ePWqAZSTSrmmDhWbaQX5i4BirTvve9t0Fq7Zpt8bqZmSBY67/2twVfVp4H8Gh3cD84NZDiKZsnbrwKjiptDtDz99aK+oUoldxke5JluKXQS8723T+faO3KjZh5pWLnEruw4O0AbsAu5w9x/XIB6RxJQbI1Bq+wWZNMo1GTL2ImDhmm0lBx6rwJG4lFsH5/PA5eSbjL9oZje4+y01i0wkZmPX6nhTSzODQ8VXJY6ymquMj3JN9pW6eNDAY4lTuRacJUC7ux8ws18Bvg8o6UhdC4ua3OAQxuENj3KDQzQ32ahjhS6dr8V0Y6Rck3GlZh9q4LHEqVyB89/ufgDA3f/LzDRdU+raNT27WL99z6ECZmwhMzziTJvazGuvH+SV1/LN6QZctmDmofE3EgvlmowrNfC4cPahVj+WyVauwDnVzDYFt23Mfdz9I7FGJjKJevpyo4qbUgYPDPPEmg/XJCY5RLkm4yrNPiy2BcqyDf0s29BPm4odGadyBc4FY+7/WZyBjGVmTwK/AEaA1929s5avL9nR05fjqjt3VixuQE3mCUk010htlJt9WG4LFM24kvEqN038B7UMpIT3ufvPkw5C6tfYbqlytGBfMpLMNbqQSodKg42HhkdYtqGftVsH1JojkVWaJi5Sd67p2cX6h/dQzcxuNYM3NF1IJazUIOSxcoNDdN+1E1BrjlSW5gLHgXvMzIGvufu6wgfNbCmwFGDmzPI73EpjuKZn16FVh6PQAGKRdAj3sIpyTTJ80Fmx8REVOFJRmmcrLHT3dwHnAZ82s1G7Crv7OnfvdPfO6dOnJxOhpEa1xU2TGTctaVdxI+GF1I7goukIZrbUzHrNrHffPm2NFYeujjYuWzCz7K6rhYaGD9LTl4s1Jql/5iXa8c1sM8WXBAFqO7PBzFYBL7t70cGHnZ2d3tvbW6twJEWu6dnF7Q8/XdVKwwbctKRdV4AxM7MdUca0JJlrzGyGu+81sxOAe4HPuvsDpc5XrolX4TpVTWaRfq/VvSylck25LqrEZjKY2THAFHf/RXD7HOC6pOKRdLrslod48PH9VT0n7JZSMkyVxHKNu+8N/n3ezL4DnAmULHAkXmNnWnVcdw8vHCi+ungoHJdz7ebdDB4Y1ho6ckikWVRm1gLMdPeBmkQFJwLfsfzy+G8A/t7dv1+j15Y60NOXq7q4mTa1mZXnz1XiS5mkco0upNJv5flzWbahv+J5wwf9UCGUGxyi+24NRJYIg4zN7HzyV1hHAbPNrB24Ls5mY3f/GXBGXN9f6ldPX45rN++ueFVXaGrzFP508TuV7FIugVyjC6mU6+poo/ep/ZGXeggNjzjXbt6t3/kGF2UW1Sryzbb/BODu/WY2K76QRIrr6cvRffdOhkeipbopBjf+jsba1JFV1DDX6EKqPqzumkfnKcdx1Z07qxprV81FkGRTlALndXd/0bSbsiRkPGNtNPCwLinXSFHh7/HY/awqmbV8i1pwG1iUAudRM/sY0GRmpwF/BPxLvGGJVD/1G2Dhqcex/pPviSkiiZlyjZRUuJ9VlEUBQweGD3Llnf2jvoc0hijr4HwWmAu8CtwOvAQsizEmEc6+8Z+qXtfm5iXtKm7qm3KNlNXV0caDy9/Pk2s+zM1L2iMv5HbQYdmGfjquu0fr5zSQii047n4AuDr4EonNeAYQAzRPMdZefIauzuqcco1UI/x9776rn+GD0Z7zwoFhlm3op/ep/VrkswGULHDStNCfZF+1A4hDrS3NrPqIpn7XM+UaGa9w3ZzCBQKjWL99D52nHKe8kXGpXOhPGsd4ViIGjbXJGOUamZCw0Ik6bs9Bu5M3gEgL/YnEYTyzowBOO+EYFTcZolwjkyXsdlr/8B6iXDPlBodYtqGfZRv61RqcQeW6qHZRvtn4nbFEJJk33rE2U4AbtYdU5ijXyGRa3TWP1V3z8t3ed+1k+GC01uHBIY3PyZpyXVS/XbMopGGMZ6zNFIOPzZ+ppJNdyjUy6Q6vnfMIQ1FHIQO3bd/DE/teVitxBpTronqqloFI9lW7ro2ajBuDco3EZbyDkB98fD/X9OzSRVWdi7IX1QLgS8Dbye8R0wS84u6/HHNskhE9fTn++O6dvBax1eaYo5rYfd25MUclaaNcI3EpLHSiroa8fvse/vEn+9g7OKQdyutUlJWMvwxcAtwFdAK/B7w1zqAkG3r6cqzatJvBoehjbaYYXP9RXTU1KOUaiVVYoETJSw6HWnzCwcgan1NfIi0E6e6PAU3uPuLufwu8L96wpN6FA/yqKW6mTW3W5pgNTrlG4tbV0Ub/ynO4fMHMqp972/Y9Wgm5jkRpwTlgZkcB/Wb2ReBZ4Jh4w5J6t2rT7sizF7QSsQSUa6Rmwl3KP/ftR3j19eiDkJdt6OeKO/tx16a+aRelBed3g/M+A7wCnAxcGGdQAGZ2rpkNmNljZrY87teTyRW15WZq8xQVNxJKJNdI4+rqaGNg9XlcvmAmTcEu9k0RdrMP19jJDQ6xYuMuteqkVJS9qMIZDv9tZpvd/Ucxx4SZNQF/BZwNPAP8m5ltcvd/j/u1pXqFqxE3mXHp/JMjPU+rEUuhJHKNCBxeOyf0qyu2ELEBmqHhEdZuHdBFWgpF3Yw19DexRHGkM4HH3P1n7v4acAdwQY1eW6ow//p7uW37nkNbLYy4c9v2PRz9htIfrWlTm7Xzt1RSk1yjlmIp5mPzqxufE3X6udRWtQVO5ba7ydEGPF1w/5ng2OFAzJaaWa+Z9e7bt69GYUmopy/H7OVbeO4XrxV9/LXXD9LcNPrj0txk3Lyknb7Pn6OrHakk9lxT0FJ8HnA6cKmZnR7360r6re6ax8JTj4t8voG6qVKo2gLn2liiOFKx5DaqwdDd17l7p7t3Tp8+vUZhSU9fjo7r7mHZhv7Sa+uT/89ae9EZtLW2YOQH4629SGNtJLJa5Bq1FEtJ6z/5Hm5e0s60qc0Vz3Vg7daB+IOSqkRZ6O9+d/8AgLv3jD0Wk2fIDzAMvQXYG+PrSQRn3/hP/PT5VyKd22R2aHEtkSgSyDXFWornF4lrKbAUYObM6qcWS/2qZpfyveqmSp1ym22+EZgKHG9m0zjcqvLLwIyY4/o34DQzmw3kyC/+9bGYX1PKmH/9vSW7o4qJOtBYJMFcU7GlGPKtxcA6gM7OzuibqElmhFPKy20SPKO15dDtcGsIrYKcrHItOJ8ClpFPMIWzGV4i328dG3d/3cw+A2wlv1z7re6+O87XlCONd9fvE489Sqt9SjWSyjVqKZbIym330NLcRPeiOQBH7GKeGxyi+66dh76H1E65zTb/AvgLM/usu3+phjGFr/894Hu1fl3Jq3ZjzNCJxx7Fw1efHUNEklUJ5hq1FEvVwiKlVAtNsUVOhw86V2zoH/V8iV+5Lqr3u/s2IGdmi8c+7u4bY41MEtPTl2N9lcWNGVw2f6ZabqRqSeUatRTLeJUbX1hqkVMnvwrytZt3s/L8uSp0aqBcF9VvAduA84s85oAKnIwpXLAvqtNOOIZ7rzwrvqCkESSWa9RSLLX2woFhVmzcBag1J27luqhWBjf/0N1fLXzMzKIvECB14bJbHuLBx/dX9ZwTjz1KxY1MmHKNZMm0qc0Vxy1q9ePaiLIOzkYzO1QImdmbgXvjC0lqracvV1VxYwaXL5ipsTYy2ZRrpO6tPH/uEYucFpMbHGL28i0sXLNNiwTGJMpu4j3A3WZ2IfkZB5uA/x1nUFIb1XZJTW2ewp8ufqeuOiQuPSjXSJ0L82O443g5Tr7QuWJDP71P7dcYxkkWZbPNW8zsKPLJZxbwKXf/l5jjkphVM0uqTes4SA0o10hWhLmycLp4OQ7ctn0Pt23fQ2tLM6s+okHIk6HcLKorC++Sv6LqBxaY2QJ3vzHm2CRGtz/8dOWTyHdF6apC4qRcI1lUOJ08NzhEk1mk1vLBoWGWbejnc99+hP97oVrMJ6JcC86xY+5/p8RxSbmevhyrNu0+NH1x2tTmSL9oC089TsWN1IJyjWRSsenkC9dsi7T7+KuvH9TaORNUbhZVrTbWlBgV64qqNMJfXVJSS8o10ki6F83higqbFYccVORMQLkuqpvdfZmZbab4/iwfiTUymbCevlzVqxGrS0pqTblGGklXRxu9T+1n/fY9kYscLRA4PuW6qL4V/PtntQhEJlfUQcRhv3CTGZfOP1nFjSRBuUYaSrh5Z+HQgUpeODBM993a06oa5hXGYpjZMcCQux8M7jcBR7v7gRrEF0lnZ6f39vYmHUZqRC1u2lpbeHD5+2sQkTQiM9vh7p1VnK9cIw2npy/H5779CK++frCq5xlwmVrcgdK5JspCf/cDUwvutwD3TVZgMvmizpAKd78VSQnlGmk4XR1tDKw+j9NOOKaq54VTy6/p2RVPYBkQpcB5o7u/HN4Jbk8tc74kLMoMqcsXzFQzp6SNco00rHuvPIvLF8ys+nm3bd+jlZBLiFLgvGJm7wrvmNm7gcpz3CQxTVZ+mXANJJaUUq6Rhra6ax43L2mPtNVDoWUb+tWSU0SUrRqWAXeZ2d7g/knAkrgCMrNVwCeBfcGh/xPs+CsRXTr/5KJjcI45qonrPzpPLTeSVsuoYa4RSaMwP1czABnyLTmdpxyn/F6g4iBjADNrBuaQH9f0E3eP/q5XG1C+wHnZ3SPPqGikgX89fTnWbh1g7+AQM8qsV1O4z5RmSEkSqh1kHDynZrlmPBop10jywnwfZWHAsRaeehzrP/meGKJKn1K5ptw6OL8OPO3u/+nuw0HT8YXAU2a2yt2jbz8tk+Kanl2j1k7IDQ6xYmO+WXJskbO6a54KGqkLSeQatRRLPShcCbma/QMBHnx8P5fd8lDDFDnFlBuD8zXgNQAzey+wBvg74EVgXcxxfcbMHjGzW81sWrETzGypmfWaWe++ffuKnZIp4Yd7bHvb0PAIa7cOJBKTyCRJKtfc5O7twZeKG0m11V3zqh6E/ODj+xt6AHK5Aqep4MppCbDO3b/t7n8CvHUiL2pm95nZo0W+LgC+CpwKtAPPAn9e7Hu4+zp373T3zunTp08knNTr6cuxvkzlvncczZciKRJbrhHJkvEUOSs27mrYIqdsgWNmYRfWB4BtBY9FGZxckrt/0N3fUeTru+7+nLuPBIt93QKcOZHXyoK1WwfKLuk9o7WlZrGIxCC2XFNBxZZiaLzWYkm3cKZVa0tzpPMbuZW/XPK4HfiBmf2c/FTNHwKY2VvJNx3HwsxOcvdng7sfBR6N67XqRbkWGkML9kndiyXXmNl9wJuLPHQ1+ZbiL5BfL+0L5FuKP1Hs+7j7OoKuss7OzijbB4nEauwu5b929fd4baT0R3Pv4FDkCSpZUm438evN7H7yUzXv8cPTraYAn40xpi+aWTv5xPMk8KkYXytVSn0AZ7S2lBxFf5kW7JM6F1eucfcPRjnPzG4B/mG8ryOStC9edEbZHcpbpzazYuMuhoZHgPITVLKkbPOvu28vcuw/4gsH3P134/z+aVVuhlT3ojmjPpygfUgkW2qda9RSLFlSbu2cluYm3Bn19wMOd101bIEjtREOIi41QyrcELPRmhdFYtSwLcWSTWG3VbGegCs29Bd9Tm5wiIVrtmX274oKnBQoN4g4HH8zts9VRMavUVuKJfuK/a0otVigwaHjucEhrtjQT+9T+zPTMxBlLyqJWblBxJohJSIiE9G9aA4tzU2jjhkccWEd7lCelWnlasGpkZ6+HFd/ZxevvJbvBy0cQ1NqELFmSImIyESFLTqFXVfltn/447t3ZqLHQAVODfT05bjqrp2MHDxcL4eVMpQfRJyFD5mIiCRrbNfVwjXbShY5r404s5ZvYdrUZlaeP7du/w6piypmPX05rryzf1RxU+j2h5+mq6ONGxbPo621BQPaWlu4aUl7ZvpBRUQkXaL0DrxwYJjuu3fWbZeVWnBiFGVztJFgyQ8NIhYRkVrp6mhjxcZHGBo+WPa84RFn2YZ+1m4dqLtZVmrBiUFPX472a++JtPNrk1kNIhIRERnthsXvjHxuuDZbPbXmqMCZZD19OVZs3HXEYkulXDr/5JgjEhEROVJXRxsLTz0u8vn1tq+VCpxJtnbrwBErRpZyuVYiFhGRBK3/5Hu4fMFMonYmlFvWJG00BmeCevpyXLt5Ny8ciNZiA9A8xVh78Rl11ZcpIiLZtLprHqu75tHTlyu63UOhelqbTQXOBEQZRDxWS/MUblj8ThU3IiKSKoWTXcLhFoU9Ei3NTXW1NpsKnHEYT6tNva8nICIijaPY4oDFZlEV2/sqLX/nzL3ULkj1o7Oz03t7e2vyWsWq2lIMUvcfLlIrZrbD3TuTjmMy1TLXiKRdsb+H4RYQbTX821cq1yTSgmNmFwOrgLcDZ7p7b8FjK4A/AEaAP3L3rUnEOFZYpZZb3rpQW2vLoV3ARUREsqbYpJqwySScVg4kdoGf1CyqR4HFwAOFB83sdOASYC5wLvAVM2s68um1FVapUYub5ilWV/2UIiIi1ao0oyrpaeWJtOC4+48B7Mh5aRcAd7j7q8ATZvYYcCbwUG0jHN2vOMXs0IrDlWgQsYiINIJKm3ZCviVn4ZptiQzVSNs6OG3A0wX3nwmO1VRhi41DpOKmtaWZm5e08+MvnKfiRkREMq970Rxamit3siS1CnJsBY6Z3Wdmjxb5uqDc04ocK1pdmNlSM+s1s959+/ZNTtCBahbra2tt4eYl7fSvPEeFjUjKmNnFZrbbzA6aWeeYx1aY2WNmNmBmi5KKUaRehRtFR9lyaGh4hKvu3Mns5VtYuGZbTYqd2Lqo3P2D43jaM0Dh3gVvAfaW+P7rgHWQn9kwjtcqKcpKjS3NTdyweJ6KGpF0C8f7fa3w4JjxfjOA+8zs19w92pWNiACHBxBHmV0c9obUagBy2rqoNgGXmNnRZjYbOA3411oHUWqlxiYzjHyrjYobkfRz9x+7e7FRjofG+7n7E0A43k9EqhS25LS1tmBE20S6FgOQk5om/lHgS8B0YIuZ9bv7InffbWZ3Av8OvA58Os4rqrHLUoeL8XUvmlN0BUcVNSKZ0QZsL7hfcryfmS0FlgLMnDkz/shE6lClVZCLCQcgx7VIYCItOO7+HXd/i7sf7e4nuvuigseud/dT3X2Ou/+/uGLo6cvRfdfOUXtuvHBgmO67dwKMqkbVYiOSXnGP93P3de7e6e6d06dPn5ygRTIsaouOwaHJPHEMRG7YrRrWbh1g+OCR+Wx4xFm7dYAHl79fBY1IHYh7vJ+IVK9Si0644nGhsNtqsv72NkyBM3a/jHJz9+tpO3gRGZdNwN+b2Y3kBxknMt5PpBEU29eq1N/g3OAQs5ZvASa+h2NDFDhjq8fc4FDR6jFUT9vBi0hpaRnvJ9LoClt0ABau2VZxkcDCYSPjKXLSNosqFuX2yxiruUnbLIhkRRrG+4nIkaIuEhgOGxmPhihwynU5tbY0H7o9bWozay86Q2NvREREYjR2IHI54x020hBdVKX6+7Tjt4iISDIKu63KdVmNd9hIQ7TgFGsKa2luUleUiIhICnQvmkPzlCPbciYybKQhWnCKjeBOYmdTEREROVL497jY4ruaRVXB2BHcIiIikh6T/Xe6IbqoREREpLGowBEREZHMMfdSK8LUDzPbBzxVg5c6Hvh5DV6nGoopujTGleWYTnH3TG3eNEm5Jo3/56G0xqa4qpfW2OKIq2iuyUSBUytm1uvunUnHUUgxRZfGuBRT40nz+5vW2BRX9dIaWy3jUheViIiIZI4KHBEREckcFTjVWZd0AEUopujSGJdiajxpfn/TGpviql5aY6tZXBqDIyIiIpmjFhwRERHJHBU4IiIikjkqcCows4vNbLeZHTSzzjGPrTCzx8xswMwWJRjjKjPLmVl/8PWhBGM5N3g/HjOz5UnFUcjMnjSzXcF705tgHLea2fNm9mjBsePM7F4z+2nw77QUxJSaz1OWlMolZjbLzIYK3u+/TkNcwWOpyHFBLKn6XKYx10F68l0QS6I5TwVOZY8Ci4EHCg+a2enAJcBc4FzgK2bWdOTTa+Ymd28Pvr6XRADBz/9XwHnA6cClwfuUBu8L3psk14X4BvnPSqHlwP3ufhpwf3A/6ZggBZ+nDCqaSwKPF7zff5iGuFKY4yAln8uU5zpIR76DhHOeCpwK3P3H7j5Q5KELgDvc/VV3fwJ4DDizttGlzpnAY+7+M3d/DbiD/PskgLs/AOwfc/gC4JvB7W8CXSmISWJQJpckSjluXJTrIkg656nAGb824OmC+88Ex5LyGTN7JGgSrGk3R4G0vSchB+4xsx1mtjTpYMY40d2fBQj+PSHheEJp+Dw1ktlm1mdmPzCz30w6mEAaf5/T8rlM43sTSnO+gxrmvDfE9Y3riZndB7y5yENXu/t3Sz2tyLHY5tyXixH4KvCF4PW/APw58Im4Yimjpu9JFRa6+14zOwG418x+ElxZSHFp+TzVnXHmkmeBme7+X2b2bqDHzOa6+0sJx1Xz3+c6yXOQ3lwHyneHqMAB3P2D43jaM8DJBfffAuydnIiOFDVGM7sF+Ie44qigpu9JVO6+N/j3eTP7Dvnm5bT8wj9nZie5+7NmdhLwfNIBuftz4e2EP091Zzy5xN1fBV4Nbu8ws8eBXwMmbYBoPeQ4qJs8BynNdZD6fAc1zHnqohq/TcAlZna0mc0GTgP+NYlAgg9J6KPkBw0m4d+A08xstpkdRX6A4qaEYgHAzI4xs2PD28A5JPf+FLMJ+Hhw++NAqavpmknR56khmNn0cPCumf0q+Vzys2SjAlKU4yB1n8vU5Tqoi3wHtcx57q6vMl/kf5GeIX+F9RywteCxq4HHgQHgvARj/BawC3gk+PCclGAsHwL+I3hfrk7B/9+vAjuDr91JxgTcTr47Yjj4TP0B8CvkZxL8NPj3uBTElJrPU5a+SuUS4MLgs7kT+BFwfhriCh5LRY4LYknV5zJtuS6IKTX5Logn0ZynrRpEREQkc9RFJSIiIpmjAkdEREQyRwWOiIiIZI4KHBEREckcFTgiIiKSOSpwZBQzGwl2od1tZjvN7EozG9fnxMxejnDOKjNzM3trwbErgmOdwf0nzez4Mt/jLDN7MVjq/sdmtrLgsTPN7IFg19+fmNnfmNnUgse/a2YPjefnE5HxU66RuKnAkbGGPL8L7VzgbPJrPays8JyJ2kV+oazQRcC/V/k9fujuHUAncLmZvdvMTgTuAj7n7nOAtwPfB8KFsFqBdwGtwUJmIlI7yjUSKxU4UpK7Pw8sJb/BnZnZLDP7oZn9KPj6H5BfYTS4cuk3s0dtzGaBZna8mT1kZh8u8VI9BDvxBiu5vgjsG2fMrwA7gFOBTwPfdPeHgsfc3e/2w1sRXAhsJr8T8CXFvp+IxE+5RuKgAkfKcvefkf+cnEB+z5Cz3f1dwBLgL4PTPkZ+9dN24AygP3x+cGWzBfi8u28p8TIvAU+b2TuAS4EN443XzH4FWEB+Fc93kE9ApVxKfqXN24PbIpIQ5RqZbNpsU6IId85tBr5sZu3ACPkNASG/L8utZtYM9Lh7f8H59wOfdvcfVHiN8MpmEfAB4PerjPE3zawPOAiscffdZsU2/A1+oHwyfCvwz+7uZva6mb3D3dO2b4tII1GukUmjFhwpK2jGHSF/RXUF+b1qziDf/3wUgLs/ALwXyAHfMrPfC57+OvmrmkUF3+/6oHm5f8xLbQZ+F9jj7i+NI9QfunuHu7/b3f86OLYbeHeJ85cA04AnzOxJYBZqOhZJjHKNTDYVOFKSmU0H/hr4suc3LXsT8Ky7HySfIMIdkE8Bnnf3W4Cvkx9MB+DAJ4C3mdlyAHe/OhhY2F74Wu4+BHwOuH4Sf4QvAx83s/kFP9PlZvZm8s3E57r7LHefRT45KemIJEC5RuKgLioZqyW44mkmf1X0LeDG4LGvAN82s4uBfwReCY6fBXSb2TDwMhBeVeHuI2Z2CbDZzF5y96+UemF3v6NMXI+Y2cHg9p3ufmWlH8Tdnwte+8/M7ATyTcoPkN+xeSawveDcJ8zsJTOb7+4PV/reIjJhyjXKNbHSbuIiIiKSOeqiEhERkcxRgSMiIiKZowJHREREMkcFjoiIiGSOChwRERHJHBU4IiIikjkqcERERCRz/j9iYGU6P3bDzgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots(1, 2)\n",
    "fig.set_size_inches((8, 3))\n",
    "for i in range(2):\n",
    "    axs[i].set_title(f'PC{i+1}')\n",
    "    axs[i].scatter(pcs1[:, i], pcs2[:, i])\n",
    "    axs[i].set_xlabel('Dask-ML PCA')\n",
    "    axs[i].set_ylabel('Scikit-allel PCA')\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.testing.assert_allclose(pcs1[:, 0], pcs2[:, 0], atol=.5)\n",
    "np.testing.assert_allclose(pcs1[:, 1], -pcs2[:, 1], atol=.5)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}