Linear interpolation with only xarray native methods

vertical_linear_interpolation_xarray_native.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray as xr\n",
    "from dask import array as da\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from dask.distributed import Client"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table style=\"border: 2px solid white;\">\n",
       "<tr>\n",
       "<td style=\"vertical-align: top; border: 0px solid white\">\n",
       "<h3>Client</h3>\n",
       "<ul>\n",
       "  <li><b>Scheduler: </b>tcp://127.0.0.1:38086\n",
       "  <li><b>Dashboard: </b><a href='http://127.0.0.1:8787/status' target='_blank'>http://127.0.0.1:8787/status</a>\n",
       "</ul>\n",
       "</td>\n",
       "<td style=\"vertical-align: top; border: 0px solid white\">\n",
       "<h3>Cluster</h3>\n",
       "<ul>\n",
       "  <li><b>Workers: </b>1</li>\n",
       "  <li><b>Cores: </b>6</li>\n",
       "  <li><b>Memory: </b>6.00 GB</li>\n",
       "</ul>\n",
       "</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<Client: scheduler='tcp://127.0.0.1:38086' processes=1 cores=6>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "client = Client(n_workers=1, threads_per_worker=6, memory_limit=6e9)\n",
    "client"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# array dimensions\n",
    "L, l = 12, 1    # time\n",
    "K, k = 100, 50   # vertical\n",
    "J, j = 720, 180  # lat\n",
    "I, i = 360, 10   # lon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.DataArray 'add-fc977cbe09d05845638bd549eb3a2b63' (time: 12, s: 100, lat: 720, lon: 360)>\n",
       "dask.array<shape=(12, 100, 720, 360), dtype=float64, chunksize=(1, 50, 180, 10)>\n",
       "Coordinates:\n",
       "  * s        (s) float64 0.0 0.0101 0.0202 0.0303 ... 0.9697 0.9798 0.9899 1.0\n",
       "Dimensions without coordinates: time, lat, lon"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "darr = xr.DataArray(\n",
    "    (\n",
    "        da.ones(shape=(L, K, J, I), chunks=(l, k, j, i)).cumsum(axis=1) - 1.0\n",
    "        + da.random.uniform(0, 0.5, size=(L, K, J, I), chunks=(l, k, j, i))\n",
    "    ),\n",
    "    dims=('time', 's', 'lat', 'lon'),\n",
    "    coords={'s': np.linspace(0, 1, K)}\n",
    ")\n",
    "darr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2.48832, 'GB')"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "darr.nbytes / 1e9, \"GB\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.DataArray 'uniform-d733d72a09d4a647a1d6b2d00bb56032' (time: 12, lat: 720, lon: 360)>\n",
       "dask.array<shape=(12, 720, 360), dtype=float64, chunksize=(1, 180, 10)>\n",
       "Dimensions without coordinates: time, lat, lon"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "s_random_levels = xr.DataArray(\n",
    "    da.random.uniform(0, 1, size=(L, J, I), chunks=(l, j, i)),\n",
    "    dims=('time', 'lat','lon')\n",
    ")\n",
    "s_random_levels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = darr.coords['s']\n",
    "delta_s_l = s.diff('s', label='lower')\n",
    "delta_s_u = s.diff('s', label='upper')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "s_in_interval_l = (s <= s_random_levels) & (s_random_levels < s + delta_s_l)\n",
    "s_in_interval_u = (s - delta_s_u <= s_random_levels) & (s_random_levels < s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# linear interpolation\n",
    "darr_int = (\n",
    "    darr.where(s_in_interval_l).mean('s')\n",
    "    + (\n",
    "        (s_random_levels - s.where(s_in_interval_l).mean('s'))\n",
    "        / (s.where(s_in_interval_u).mean('s') - s.where(s_in_interval_l).mean('s'))\n",
    "        * (darr.where(s_in_interval_u).mean('s') - darr.where(s_in_interval_l).mean('s'))\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.DataArray (time: 12, lat: 720, lon: 360)>\n",
       "dask.array<shape=(12, 720, 360), dtype=float64, chunksize=(1, 180, 10)>\n",
       "Dimensions without coordinates: time, lat, lon"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "darr_int"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "112324"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(darr_int.__dask_graph__())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1min 1s, sys: 2.47 s, total: 1min 4s\n",
      "Wall time: 1min 27s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(array([118121., 260575., 348070., 401735., 427422., 426774., 401618.,\n",
       "        347858., 260482., 117745.]),\n",
       " array([7.99377563e-06, 5.00010185e-02, 9.99940432e-02, 1.49987068e-01,\n",
       "        1.99980093e-01, 2.49973117e-01, 2.99966142e-01, 3.49959167e-01,\n",
       "        3.99952191e-01, 4.49945216e-01, 4.99938241e-01]),\n",
       " <a list of 10 Patch objects>)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEICAYAAACj2qi6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAGj5JREFUeJzt3XuwXeV93vHvYzA2ic3FoDBUkiNsK5ORmQTbKlaaTupAAgInFplgj4gTZI9qJQHS3NpYJGlx7NBA2piasU1DimrhphEqiYNiyyEq4PG4DRc5xmCBXU64DFIwKBIXU2Jcwa9/7Ffx5vjonFeXc7ak8/3M7Nlr/da71vsuncN5WJe9V6oKSZJ6vGzUA5AkHToMDUlSN0NDktTN0JAkdTM0JEndDA1JUjdDQxonyZYkbxv1OKSDkaGhWSfJw0l+bFztPUm+AFBVb6yqz02xjQVJKsmR0zhU6aBjaEgHIcNIBytDQxpn+EgkyelJNid5JsnjST7cmn2+vT+V5NkkP5TkZUl+O8kjSZ5Icn2SY4e2e2FbtiPJvx3XzweS3JjkvyV5BnhP6/uvkzyV5LEkH01y1ND2KslFSR5I8o0kH0ry+iT/u413/XB76UAwNKTJfQT4SFUdA7weWN/qP9Lej6uqV1XVXwPvaa8fBV4HvAr4KECSRcDHgXcDJwPHAnPH9bUMuBE4Dvhj4AXgV4ETgR8CzgQuGrfO2cBbgCXAbwDXAj8LzAdOBS7Yj32XvoOhodnqz9v/wT+V5CkGf9An8v+ANyQ5saqerarbJ9nmu4EPV9WDVfUscCmwvJ1qOh/4i6r6QlV9C/h3wPgvfvvrqvrzqnqxqv6hqr5YVbdX1a6qehj4Q+BfjFvn96vqmaraAnwF+KvW/9PAZ4E39f+TSFMzNDRbnVdVx+1+8Z3/B7/bSuD7gK8muSvJT0yyzX8CPDI0/whwJHBSW/bo7gVV9RywY9z6jw7PJPm+JJ9O8vV2yurfMzjqGPb40PQ/TDD/qknGK+01Q0OaRFU9UFUXAN8DXAncmOS7+c6jBIC/A753aP61wC4Gf8gfA+btXpDkaOCE8d2Nm78G+CqwsJ0e+00g+7430v4zNKRJJPnZJHOq6kXgqVZ+Edje3l831PxPgF9NckqSVzE4MrihqnYxuFbxk0n+Wbs4/QGmDoBXA88Azyb5fuAXD9R+SfvK0JAmtxTYkuRZBhfFl7frDc8BlwP/q10XWQKsAT7J4M6qh4BvAr8E0K45/BKwjsFRx7PAE8Dzk/T9r4GfAb4B/BFww4HfPWnvxIcwSTOvHYk8xeDU00OjHo/UyyMNaYYk+ckk39WuifxH4F7g4dGOSto7hoY0c5YxuFj+d8BCBqe6PNTXIcXTU5Kkbh5pSJK6HXZfinbiiSfWggULRj0MSTqkfPGLX/z7qpozVbvDLjQWLFjA5s2bRz0MSTqkJHlk6laenpIk7QVDQ5LUzdCQJHUzNCRJ3QwNSVI3Q0OS1M3QkCR1MzQkSd0MDUlSt8PuE+HSVBas/syohzDjHr7i7aMegg4THmlIkroZGpKkboaGJKmboSFJ6mZoSJK6GRqSpG6GhiSpm6EhSepmaEiSunWHRpIjknwpyafb/ClJ7kgyluSGJEe1+iva/FhbvmBoG5e2+teSnD1UX9pqY0lWD9Un7EOSNBp7c6Txy8D9Q/NXAldV1RuAJ4GVrb4SeLLVr2rtSLIIWA68EVgKfLwF0RHAx4BzgEXABa3tZH1Ikkag67unkswD3g5cDvxakgBnAD/TmqwFPgBcAyxr0wA3Ah9t7ZcB66rqeeChJGPA6a3dWFU92PpaByxLcv8kfegQNxu//2mURvXv7XdeHX56jzT+E/AbwItt/gTgqara1ea3AnPb9FzgUYC2/OnW/h/r49bZU32yPl4iyaokm5Ns3r59e+cuSZL21pShkeQngCeq6oszMJ59UlXXVtXiqlo8Z86cUQ9Hkg5bPaenfhh4R5JzgVcCxwAfAY5LcmQ7EpgHbGvttwHzga1JjgSOBXYM1XcbXmei+o5J+pAkjcCURxpVdWlVzauqBQwuZN9aVe8GbgPOb81WADe16Q1tnrb81qqqVl/e7q46BVgI3AncBSxsd0od1frY0NbZUx+SpBHYn89pvJ/BRfExBtcfrmv164ATWv3XgNUAVbUFWA/cB/wlcHFVvdCOIi4BbmZwd9b61nayPiRJI7BXT+6rqs8Bn2vTD/Ltu5+G23wTeOce1r+cwR1Y4+sbgY0T1CfsQ5I0Gn4iXJLUzdCQJHUzNCRJ3QwNSVI3Q0OS1M3QkCR1MzQkSd0MDUlSN0NDktTN0JAkdTM0JEndDA1JUjdDQ5LUzdCQJHUzNCRJ3XqeEf7KJHcm+XKSLUl+p9U/keShJHe312mtniRXJxlLck+SNw9ta0WSB9prxVD9LUnubetcnSSt/pokm1r7TUmOP/D/BJKkXj1HGs8DZ1TVDwKnAUuTLGnL/k1VndZed7faOQwe5boQWAVcA4MAAC4D3srgwUqXDYXANcD7htZb2uqrgVuqaiFwS5uXJI1IzzPCq6qebbMvb6+aZJVlwPVtvduB45KcDJwNbKqqnVX1JLCJQQCdDBxTVbe354JfD5w3tK21bXrtUF2SNAJd1zSSHJHkbuAJBn/472iLLm+noK5K8opWmws8OrT61labrL51gjrASVX1WJv+OnDSHsa3KsnmJJu3b9/es0uSpH3Q9YzwqnoBOC3JccCnkpwKXMrgD/lRwLXA+4EPTtdAq6qSTHiEU1XXtjGwePHiyY6CNM6C1Z8Z9RB0GBvl79fDV7x9ZH0fzvbq7qmqegq4DVhaVY+1U1DPA/+VwXUKgG3A/KHV5rXaZPV5E9QBHm+nr2jvT+zNeCVJB1bP3VNz2hEGSY4Gfhz46tAf8zC41vCVtsoG4MJ2F9US4Ol2iulm4Kwkx7cL4GcBN7dlzyRZ0rZ1IXDT0LZ232W1YqguSRqBntNTJwNrkxzBIGTWV9Wnk9yaZA4Q4G7gF1r7jcC5wBjwHPBegKrameRDwF2t3Qeramebvgj4BHA08Nn2ArgCWJ9kJfAI8K593VFJ0v6bMjSq6h7gTRPUz9hD+wIu3sOyNcCaCeqbgVMnqO8AzpxqjJKkmeEnwiVJ3QwNSVI3Q0OS1M3QkCR1MzQkSd0MDUlSN0NDktTN0JAkdTM0JEndDA1JUjdDQ5LUzdCQJHUzNCRJ3QwNSVI3Q0OS1K3nyX2vTHJnki8n2ZLkd1r9lCR3JBlLckOSo1r9FW1+rC1fMLStS1v9a0nOHqovbbWxJKuH6hP2IUkajZ4jjeeBM6rqB4HTgKXtMa5XAldV1RuAJ4GVrf1K4MlWv6q1I8kiYDnwRmAp8PEkR7QnAn4MOAdYBFzQ2jJJH5KkEZgyNGrg2Tb78vYq4AzgxlZfy+A54QDL2jxt+Znt2d/LgHVV9XxVPcTgcbCnt9dYVT1YVd8C1gHL2jp76kOSNAJd1zTaEcHdwBPAJuBvgaeqaldrshWY26bnAo8CtOVPAycM18ets6f6CZP0MX58q5JsTrJ5+/btPbskSdoHXaFRVS9U1WnAPAZHBt8/raPaS1V1bVUtrqrFc+bMGfVwJOmwtVd3T1XVU8BtwA8BxyU5si2aB2xr09uA+QBt+bHAjuH6uHX2VN8xSR+SpBHouXtqTpLj2vTRwI8D9zMIj/NbsxXATW16Q5unLb+1qqrVl7e7q04BFgJ3AncBC9udUkcxuFi+oa2zpz4kSSNw5NRNOBlY2+5yehmwvqo+neQ+YF2S3wW+BFzX2l8HfDLJGLCTQQhQVVuSrAfuA3YBF1fVCwBJLgFuBo4A1lTVlrat9++hD0nSCEwZGlV1D/CmCeoPMri+Mb7+TeCde9jW5cDlE9Q3Aht7+5AkjYafCJckdTM0JEndDA1JUjdDQ5LUzdCQJHUzNCRJ3QwNSVI3Q0OS1M3QkCR1MzQkSd0MDUlSN0NDktTN0JAkdTM0JEndep6noRmwYPVnRj0E6bAyqv+mHr7i7SPpd6Z4pCFJ6tbzuNf5SW5Lcl+SLUl+udU/kGRbkrvb69yhdS5NMpbka0nOHqovbbWxJKuH6qckuaPVb2iPfaU9GvaGVr8jyYIDufOSpL3Tc6SxC/j1qloELAEuTrKoLbuqqk5rr40Abdly4I3AUuDjSY5oj4v9GHAOsAi4YGg7V7ZtvQF4EljZ6iuBJ1v9qtZOkjQiU4ZGVT1WVX/Tpr8B3A/MnWSVZcC6qnq+qh4Cxhg8svV0YKyqHqyqbwHrgGVJApwB3NjWXwucN7SttW36RuDM1l6SNAJ7dU2jnR56E3BHK12S5J4ka5Ic32pzgUeHVtvaanuqnwA8VVW7xtVfsq22/OnWfvy4ViXZnGTz9u3b92aXJEl7oTs0krwK+FPgV6rqGeAa4PXAacBjwB9Mywg7VNW1VbW4qhbPmTNnVMOQpMNeV2gkeTmDwPjjqvozgKp6vKpeqKoXgT9icPoJYBswf2j1ea22p/oO4LgkR46rv2Rbbfmxrb0kaQR67p4KcB1wf1V9eKh+8lCznwK+0qY3AMvbnU+nAAuBO4G7gIXtTqmjGFws31BVBdwGnN/WXwHcNLStFW36fODW1l6SNAI9H+77YeDngHuT3N1qv8ng7qfTgAIeBn4eoKq2JFkP3MfgzquLq+oFgCSXADcDRwBrqmpL2977gXVJfhf4EoOQor1/MskYsJNB0EiSRmTK0KiqLwAT3bG0cZJ1Lgcun6C+caL1qupBvn16a7j+TeCdU41RkjQz/ES4JKmboSFJ6mZoSJK6GRqSpG6GhiSpm6EhSepmaEiSuhkakqRuhoYkqZuhIUnqZmhIkroZGpKkboaGJKmboSFJ6mZoSJK69Ty5b36S25Lcl2RLkl9u9dck2ZTkgfZ+fKsnydVJxpLck+TNQ9ta0do/kGTFUP0tSe5t61zdnha4xz4kSaPRc6SxC/j1qloELAEuTrIIWA3cUlULgVvaPMA5DB7xuhBYBVwDgwAALgPeyuCBS5cNhcA1wPuG1lva6nvqQ5I0AlOGRlU9VlV/06a/AdwPzAWWAWtbs7XAeW16GXB9DdwOHNeeJ342sKmqdlbVk8AmYGlbdkxV3d6e/339uG1N1IckaQT26ppGkgXAm4A7gJOq6rG26OvASW16LvDo0GpbW22y+tYJ6kzSx/hxrUqyOcnm7du3780uSZL2QndoJHkV8KfAr1TVM8PL2hFCHeCxvcRkfVTVtVW1uKoWz5kzZzqHIUmzWldoJHk5g8D446r6s1Z+vJ1aor0/0erbgPlDq89rtcnq8yaoT9aHJGkEeu6eCnAdcH9VfXho0QZg9x1QK4CbhuoXtruolgBPt1NMNwNnJTm+XQA/C7i5LXsmyZLW14XjtjVRH5KkETiyo80PAz8H3Jvk7lb7TeAKYH2SlcAjwLvaso3AucAY8BzwXoCq2pnkQ8Bdrd0Hq2pnm74I+ARwNPDZ9mKSPiRJIzBlaFTVF4DsYfGZE7Qv4OI9bGsNsGaC+mbg1AnqOybqQ5I0Gn4iXJLUzdCQJHUzNCRJ3QwNSVI3Q0OS1M3QkCR1MzQkSd0MDUlSN0NDktTN0JAkdTM0JEndDA1JUjdDQ5LUzdCQJHUzNCRJ3Xqe3LcmyRNJvjJU+0CSbUnubq9zh5ZdmmQsydeSnD1UX9pqY0lWD9VPSXJHq9+Q5KhWf0WbH2vLFxyonZYk7ZueI41PAEsnqF9VVae110aAJIuA5cAb2zofT3JEkiOAjwHnAIuAC1pbgCvbtt4APAmsbPWVwJOtflVrJ0kaoSlDo6o+D+ycql2zDFhXVc9X1UMMHvl6enuNVdWDVfUtYB2wrD0T/Azgxrb+WuC8oW2tbdM3Ame29pKkEdmfaxqXJLmnnb46vtXmAo8OtdnaanuqnwA8VVW7xtVfsq22/OnW/jskWZVkc5LN27dv349dkiRNZl9D4xrg9cBpwGPAHxywEe2Dqrq2qhZX1eI5c+aMciiSdFjbp9Coqser6oWqehH4IwannwC2AfOHms5rtT3VdwDHJTlyXP0l22rLj23tJUkjsk+hkeTkodmfAnbfWbUBWN7ufDoFWAjcCdwFLGx3Sh3F4GL5hqoq4Dbg/Lb+CuCmoW2taNPnA7e29pKkETlyqgZJ/gR4G3Bikq3AZcDbkpwGFPAw8PMAVbUlyXrgPmAXcHFVvdC2cwlwM3AEsKaqtrQu3g+sS/K7wJeA61r9OuCTScYYXIhfvt97K0naL1OGRlVdMEH5uglqu9tfDlw+QX0jsHGC+oN8+/TWcP2bwDunGp8kaeb4iXBJUjdDQ5LUzdCQJHUzNCRJ3QwNSVI3Q0OS1M3QkCR1MzQkSd0MDUlSN0NDktRtyq8RmU0WrP7MqIcg6RA3yr8jD1/x9mnvwyMNSVI3Q0OS1M3QkCR1MzQkSd2mDI0ka5I8keQrQ7XXJNmU5IH2fnyrJ8nVScaS3JPkzUPrrGjtH0iyYqj+liT3tnWuTpLJ+pAkjU7PkcYngKXjaquBW6pqIXBLmwc4h8EjXhcCq4BrYBAADJ7491YGD1y6bCgErgHeN7Te0in6kCSNyJShUVWfZ/C41WHLgLVtei1w3lD9+hq4HTiuPU/8bGBTVe2sqieBTcDStuyYqrq9Pf/7+nHbmqgPSdKI7Os1jZOq6rE2/XXgpDY9F3h0qN3WVpusvnWC+mR9fIckq5JsTrJ5+/bt+7A7kqQe+30hvB0h1AEYyz73UVXXVtXiqlo8Z86c6RyKJM1q+xoaj7dTS7T3J1p9GzB/qN28VpusPm+C+mR9SJJGZF9DYwOw+w6oFcBNQ/UL211US4Cn2ymmm4GzkhzfLoCfBdzclj2TZEm7a+rCcduaqA9J0ohM+d1TSf4EeBtwYpKtDO6CugJYn2Ql8AjwrtZ8I3AuMAY8B7wXoKp2JvkQcFdr98Gq2n1x/SIGd2gdDXy2vZikD0nSiEwZGlV1wR4WnTlB2wIu3sN21gBrJqhvBk6doL5joj4kSaPjJ8IlSd0MDUlSN0NDktTN0JAkdTM0JEndDA1JUjdDQ5LUzdCQJHUzNCRJ3QwNSVI3Q0OS1M3QkCR1MzQkSd0MDUlSN0NDktRtv0IjycNJ7k1yd5LNrfaaJJuSPNDej2/1JLk6yViSe5K8eWg7K1r7B5KsGKq/pW1/rK2b/RmvJGn/HIgjjR+tqtOqanGbXw3cUlULgVvaPMA5wML2WgVcA4OQYfA0wLcCpwOX7Q6a1uZ9Q+stPQDjlSTto+k4PbUMWNum1wLnDdWvr4HbgeOSnAycDWyqqp1V9SSwCVjalh1TVbe3JwJeP7QtSdII7G9oFPBXSb6YZFWrnVRVj7XprwMntem5wKND625ttcnqWyeoS5JGZMpnhE/hn1fVtiTfA2xK8tXhhVVVSWo/+5hSC6xVAK997WunuztJmrX260ijqra19yeATzG4JvF4O7VEe3+iNd8GzB9afV6rTVafN0F9onFcW1WLq2rxnDlz9meXJEmT2OfQSPLdSV69exo4C/gKsAHYfQfUCuCmNr0BuLDdRbUEeLqdxroZOCvJ8e0C+FnAzW3ZM0mWtLumLhzaliRpBPbn9NRJwKfaXbBHAv+9qv4yyV3A+iQrgUeAd7X2G4FzgTHgOeC9AFW1M8mHgLtauw9W1c42fRHwCeBo4LPtJUkakX0Ojap6EPjBCeo7gDMnqBdw8R62tQZYM0F9M3Dqvo5RknRg+YlwSVI3Q0OS1M3QkCR1MzQkSd0MDUlSN0NDktTN0JAkdTM0JEndDA1JUjdDQ5LUzdCQJHUzNCRJ3QwNSVI3Q0OS1M3QkCR1MzQkSd0O+tBIsjTJ15KMJVk96vFI0mx2UIdGkiOAjwHnAIuAC5IsGu2oJGn2OqhDAzgdGKuqB6vqW8A6YNmIxyRJs9Y+PyN8hswFHh2a3wq8dXyjJKuAVW322SRf28f+TgT+fh/XPVS5z7OD+zwL5Mr92ufv7Wl0sIdGl6q6Frh2f7eTZHNVLT4AQzpkuM+zg/s8O8zEPh/sp6e2AfOH5ue1miRpBA720LgLWJjklCRHAcuBDSMekyTNWgf16amq2pXkEuBm4AhgTVVtmcYu9/sU1yHIfZ4d3OfZYdr3OVU13X1Ikg4TB/vpKUnSQcTQkCR1m5WhMdVXkyR5RZIb2vI7kiyY+VEeWB37/CNJ/ibJriTnj2KMB1rHPv9akvuS3JPkliRd96kfzDr2+ReS3Jvk7iRfONS/YaH3a4aS/HSSSnLI34Lb8TN+T5Lt7Wd8d5J/eUAHUFWz6sXggvrfAq8DjgK+DCwa1+Yi4D+36eXADaMe9wzs8wLgB4DrgfNHPeYZ2ucfBb6rTf/iLPk5HzM0/Q7gL0c97unc39bu1cDngduBxaMe9wz8jN8DfHS6xjAbjzR6vppkGbC2Td8InJkkMzjGA23Kfa6qh6vqHuDFUQxwGvTs821V9VybvZ3B54AOZT37/MzQ7HcDh/KdML1fM/Qh4ErgmzM5uGky8q9Wmo2hMdFXk8zdU5uq2gU8DZwwI6ObHj37fLjZ231eCXx2Wkc0/br2OcnFSf4W+H3gX83Q2KbDlPub5M3A/Kr6zEwObBr1/l7/dDvtemOS+RMs32ezMTSkl0jys8Bi4D+Meiwzoao+VlWvB94P/PaoxzNdkrwM+DDw66Meywz7C2BBVf0AsIlvnzU5IGZjaPR8Nck/tklyJHAssGNGRjc9ZuPXsXTtc5IfA34LeEdVPT9DY5sue/tzXgecN60jml5T7e+rgVOBzyV5GFgCbDjEL4ZP+TOuqh1Dv8v/BXjLgRzAbAyNnq8m2QCsaNPnA7dWu8J0iJqNX8cy5T4neRPwhwwC44kRjPFA69nnhUOzbwcemMHxHWiT7m9VPV1VJ1bVgqpawOC61TuqavNohntA9PyMTx6afQdw/wEdwajvBhjRHQjnAv+HwV0Iv9VqH2TwCwXwSuB/AGPAncDrRj3mGdjnf8rg/Oj/ZXBUtWXUY56Bff6fwOPA3e21YdRjnoF9/giwpe3vbcAbRz3m6dzfcW0/xyF+91Tnz/j32s/4y+1n/P0Hsn+/RkSS1G02np6SJO0jQ0OS1M3QkCR1MzQkSd0MDUlSN0NDktTN0JAkdfv/avnJWI/DQA0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "(darr_int - (K-1) * s_random_levels).plot.hist();  # should be in [0, 0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}