JaxPM/notebooks/lensing_demo.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "0ae8de09",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "%pylab is deprecated, use %matplotlib inline and import the required libraries.\n",
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%pylab inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f60695dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "import logging\n",
    "\n",
    "logger = logging.getLogger(\"root\")\n",
    "\n",
    "class CheckTypesFilter(logging.Filter):\n",
    "    def filter(self, record):\n",
    "        return \"check_types\" not in record.getMessage()\n",
    "logger.addFilter(CheckTypesFilter())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "cdb87c17",
   "metadata": {},
   "outputs": [],
   "source": [
    "import jax\n",
    "import jax.numpy as jnp\n",
    "\n",
    "import jax_cosmo as jc\n",
    "\n",
    "from jax.experimental.ode import odeint\n",
    "from jaxpm.painting import cic_paint, cic_paint_2d\n",
    "from jaxpm.pm import linear_field, lpt, make_ode_fn\n",
    "from jaxpm.lensing import density_plane, convergence_Born"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "id": "504049d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Below are a few parameters that generate a low resolution version of the k-TNG simulation\n",
    "\n",
    "box_size = [200.,200.,4000.]    # Transverse comoving size of the simulation volume\n",
    "nc = [64, 64, 320]              # Number of transverse voxels in the simulation volume\n",
    "lensplane_width = 102.5         # Width of each lensplane\n",
    "field_size = 5                  # Size of the lensing field in degrees\n",
    "field_npix = 128                # Number of pixels in the lensing field\n",
    "z_source = jnp.linspace(0,2)    # Source planes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "id": "d5edb051",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Defining the coordinate grid for lensing map\n",
    "import astropy.units as u\n",
    "xgrid, ygrid = np.meshgrid(np.linspace(0, field_size, field_npix, endpoint=False), # range of X coordinates\n",
    "                           np.linspace(0, field_size, field_npix, endpoint=False)) # range of Y coordinates\n",
    "\n",
    "coords = np.stack([xgrid, ygrid], axis=0)*u.deg\n",
    "c = coords.reshape([2, -1]).T.to(u.rad)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "id": "3b6ded99",
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_map(Omega_c=0.2589, sigma8=0.8159):\n",
    "    # Instantiates a cosmology with desired parameters\n",
    "    cosmology = jc.Planck15(Omega_c=Omega_c, sigma8=sigma8)\n",
    "\n",
    "    # Planning out the scale factor stepping to extract desired lensplanes\n",
    "    n_lens = int(box_size[-1] // lensplane_width)\n",
    "    r = jnp.linspace(0., box_size[-1], n_lens+1)\n",
    "    r_center = 0.5*(r[1:] + r[:-1])\n",
    "\n",
    "    # Retrieve the scale factor corresponding to these distances\n",
    "    a = jc.background.a_of_chi(cosmology, r)\n",
    "    a_center = jc.background.a_of_chi(cosmology, r_center)\n",
    "\n",
    "    # Then one step per lens plane\n",
    "    stages = a_center[::-1]\n",
    "\n",
    "    # Create a small function to generate the matter power spectrum\n",
    "    k = jnp.logspace(-4, 1, 128)\n",
    "    pk = jc.power.linear_matter_power(cosmology, k)\n",
    "    pk_fn = lambda x: jc.scipy.interpolate.interp(x.reshape([-1]), k, pk).reshape(x.shape)\n",
    "\n",
    "    # Create initial conditions\n",
    "    initial_conditions = linear_field(nc, box_size, pk_fn, seed=jax.random.PRNGKey(0))\n",
    "\n",
    "    # Create particles\n",
    "    particles = jnp.stack(jnp.meshgrid(*[jnp.arange(s) for s in nc]),axis=-1).reshape([-1,3])\n",
    "\n",
    "    cosmo = jc.Planck15(Omega_c=Omega_c, sigma8=sigma8)\n",
    "\n",
    "    # Initial displacement\n",
    "    dx, p, f = lpt(cosmo, initial_conditions, particles, 0.1)\n",
    "\n",
    "    # Evolve the simulation forward\n",
    "    res = odeint(make_ode_fn(nc), [particles+dx, p], \n",
    "                 jnp.concatenate([jnp.atleast_1d(0.1), stages]), cosmo, rtol=1e-5, atol=1e-5)\n",
    "    \n",
    "    # Extract the lensplanes\n",
    "    lensplanes = []\n",
    "    for i in range(len(a_center)):\n",
    "        plane = density_plane(states[::-1][i],\n",
    "                              nc,\n",
    "                              (i+0.5)*lensplane_width/box_size[-1]*nc[-1],\n",
    "                              width=lensplane_width/box_size[-1]*nc[-1],\n",
    "                              plane_resolution=64\n",
    "                           )\n",
    "        lensplanes.append((r_center[i], stages[::-1][i], plane))\n",
    "    \n",
    "    # Get convergence map\n",
    "    m = convergence_Born(cosmology, \n",
    "                          lensplanes, \n",
    "                          dx=box_size[0]/64, \n",
    "                          dz=lensplane_width,\n",
    "                          coords=jnp.array(c).T.reshape(2,field_npix,field_npix),\n",
    "                          z_source=z_source)\n",
    "    return m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "id": "979726b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = make_map()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "42541aad",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABGYAAAOuCAYAAABFVRGkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAEAAElEQVR4nOz9W49ty7YmBn2tRfTex8ict7X29eyqOlWnipItlwuEDBISlowESNQjD/wFHhA8WELiCSEhvyDxgoR/BI8I2WVuEn7BQgUqwManqlw36tz23muvteYlM8cYvUdE46FdIvrInHMfWz4rt0SGtNacM3Nceo8e0aK1r33tayQieBkv42W8jJfxMl7Gy3gZL+NlvIyX8TJexst4GT/84Oe+gJfxMl7Gy3gZL+NlvIyX8TJexst4GS/jZbyM/38dL8DMy3gZL+NlvIyX8TJexst4GS/jZbyMl/EyXsYzjRdg5mW8jJfxMl7Gy3gZL+NlvIyX8TJexst4GS/jmcYLMPMyXsbLeBkv42W8jJfxMl7Gy3gZL+NlvIyX8UzjBZh5GS/jZbyMl/EyXsbLeBkv42W8jJfxMl7Gy3im8QLMvIyX8TJexst4GS/jZbyMl/EyXsbLeBkv42U808hf+uXf+Uv/IwEz5LgAzEBOABEkkf45JUhmtClBMiF/WsGfzqB1A9YNIH0dctLPIIJkjj+pNND9GVQq0BrQeutuOcz6nbWCmujn1ao/A/TPZH8XAS4r5OFk/xRQYiBngFj/nhLk9qjXY58pdj9g0vctEyQzyusZ5ZgwfSpIDyv4tIFOF71uJr0eGf5jhswTkFjvrwH8cNbrrU3ft0xASmhvjmgTQyaGMCGdCuhSwQ8X0GUFSgVK0Wufsl7jlNGOE7Z3B2yvEh5+mlAnQj0C6QK8+8cb8t2G6bsHYLVrrRVYbA79OlOCJIYcZ8gygU6bfidRf+j2WioVKBVy0jmlwwEA8Hf/6H81vPiHHf+t//q/JSCgzQmSCNvrhLIwPv0+Y3styPeEdAGme0E6A/N9w3RXkS4VfCrgSwHdn/V+7PnLYdLnlxltybh8PaPOhO1W19b8qSGtDfP7DbRVW/MEsbUgmdASIa0NXBqoCHitaJl1nwAgAXit4I+2Ppcc+wcA6FxArYHOq+6DUvV1r28gy4TyekE9JqRTBV9KrFc0ATVBeT3j9OMJkgh1IggDLQN1IZx/AlAFbn4loKJvkwRcviLUGahHgTBw82eE5YPg+JuC+cMGXitoq6CHC+j+tN/LttbbMoFEgNb0tZ8e9POPi97X3QMgAnl1o2tvyQAz6LIBW+lr7OYAOc6orxZsryeI3V4+V+SPF7UT5w3ICe1mRpsY5ZW+7t//u/+TZ1mPf+dv/I8FIvGskBPcVvpaksygJqAqOj9bRbuZUW8yyiGh3CRc3jBOPyNAAN4AXoHlYwNvwHxXgaZrlQRhH7nqn81siBAAAvKpgtcGABCC2SmAtwo6FyAz2pxQbyecfjKjJaBNBC6Cw7cF6dKQP1xAraEdJ0hibK8ntIlQDmr306r3ky4C3hpgJrvcJGyvGA8/Ydz9vmD5jvCjPyzIDw3z9xcAQF0Sym3Ghz+YUBe1X8L6HwSYPwG8CvIJ4GL3LIg/YU96umuYP2xI5wK+v0CYYy+BATCjLno2ba8zWibU2fcF2Vw1cAV40z00fVRb77ZQMkMSox4SylHv7fQ1g4tgehCkFVi+25DOFdO390Cp+Pf+4f/iWdbiX/u3/5fCZ8aP/kNg+VCxfLeprdsaqLW4J7pssfewbZDLqvY9JVDSc4bmWT9U5PG5kOwcH/Zy+AI8nKNE4HMB1QoqDahN9/fNhDbrfPLWkM4VXBr4YQXa1U0x+rNgjucqhLDB/p56TKgzg6qAqyDfF6RP57BZtBa1YU30bE0JcnPoZz8RZM6QRLpvJ8bdLxaUI2H50JAfGqaHAr5U0MXsYlXb5X5Ru5mxfrXg8i7j019mlBvg8qOGfE/46h8I5ruG4y/POi/3Z0AE7fUNZEnYXk1oM6PNfpboHsunCr5UndvUzxyuDbQ2gElt4U3C/c8y6gE4f02QDPzD/9m/+Sxr8b/9X/ufiyRGPWaAoHPWBG1JaKlf0vRxBT9soPNF/RVAz76U1K8jUt8H6GdPYj3ztqKvvayQwW8kf5/9jGY725cJmCe1vccJ+dMF/P5OfbM2LLzWIG7PxX5O7m8O0+nraHgNHRZgniDLDExZP9c+n2oLPw4iQG2grai/Ws3fIwYts97/POnnns6Qoq+TdY2vF79fJn3fPIHmGXRY1McdBjWJuYJIn8952u9v+7ncHCBTAq1Fr3M1e+F7Zd0g26b3Vmufm5TUR8wpPqMdJoAJ/8f/6//02XzGv/NX/03xPR7/+Rh9+HWDtAbKWc9yH+7rJ4ZMQ7jE6kd7LCOZUW4mCNv5LNB9OwyuDWjQvT7puhICeGvgImgzoy6MuhDWW46zkRqQNgFvwPGbLT6HmvqaKA1scZT7H/zppM/cYwiiiJnCngKQlCCLxgV1SXrGCkAi4Iv5oTZ/bVZ7v75JqDPh/BWjHIHlvWC6tzUpeqZyEeT7iulu6zYvJ42xpoRyO6EujMu7jLoQTj8mkADHXzfkiyA/NFAD6kFtoiT1c5YPFflO/VNeq65TjzWZIUvWdWfXgtbA57J7Dv/ef/hvPct6/O/87H9gDpx9ve8fszHk8WzOYes0Lm5q09xWpASaJoiYHWqy34tAj6Wn3Nc8MSinnY0F0O3C9Wiy8/cl6xmLS7dFIILcHiBEet5XUb9/3dSWStPPb6L2bZnDhtDhoHH+lNEOs55xmSFma4QJddE1J0n/vd0ythvC939LUF9XpE8J6UK4/WNdh1x1X5RD9/mENQZKG5DPgvlDCV+BmpifpD46GoDkfgft7H5bUt/Ttkf0FzqH6qdK+PxhWxo0vitN/ZDa7PwQ/N0//V8/uRa/CMyE8R43dCLUVwvazCgHdZCmuwI+2yZprTt1dphKciNkwAxrgEC5IZ0YkkQdSEAvGABtbI6fLQgzjjHYFtZ42NikEfTA0ompEFscdFn7Z6YONMUa29TIJSbQ1pDOxX5W4wCn+pm5YtIAYZl0I13smmNTUHwHi6DaQ5fEwAJIs8M4V1CxIC8nOxAMUEkEXgU331RdvARz4HShoaoTgNYUnGqCRx4vESQl1CUhlQZazeHhYW5b041GBCpTBOSPDvIfeDQ7OMTAwXJgrLeEu39pxY9//hG/+eN3mL5P4D8mzB+bHgwfVwW+LlsEFwAgkxnylAC2TSiCdGmAMNZXCmCUA6FlxnRHoA2g2kBVrwHQAI6zrhcFZhSkSCsCxGyzOWRPOQhia9/XmG1aqQ0kAsoJWQRUZw1Omboj2gSoAl4bpvuGNhMADV64QO8rAW0G7n9Pv69Novf1tkKS4PDLjHwP3HzTsHxfbC9vYWyQWO1AGJkWYB77gZiTXpM7LuO+BCKQq68PaBOD1wlUmgaPpaHNGbLommwTBaDQEukh2wQcgVNCS4w6cwTqzzHEnGxy4NMPznUDtaQ/BzqIC9haSCjHDCEgrQ3TA9C+Y9SZUG70s+tMBqxldb42dfbygzpw/FDAWzNnRe3A9ZwHqFEbUAVUq2IbmW1ugXIknH9E4I3AJWO6r8gfAFQ9SJEMFKoANdrPd9gX+74q4E0w3QmO3xDm93q96aSBrBCBMoM3QT7r/Uz3AAioEwGs11QXUseuIgA6VIAb+rqvovdVms53SnqYEgHVwJ5Jnd82kQFQ+lHLxwpe9VBGE3WQSUEu3MxxP3VhtIWx3TC2W8b6mnD6qSBdCPUDIV0EkIx8Trrvt88dDH/x4/aPEngFDr5/N0tmJLX1qGZjthIJDhEBxM4JEf2zVg26zJECE5AziKg7i/NkoNUQvIiB3ImxvcqoBw7AK58q0qmg+f7O6myBgLQ1SFWAhdy+mM2UlNBuJgUiLmW31hQMZzAs0Elkto9QoYFOLkNgvBXI+dITJckcVLLz13wcEgK3AimMfJk1yCp9nQsRuFa1f77+SwXOF3BrmGy9lsOEugDpwpoguKtIpxbOGlXzR0QgNkd
      "text/plain": [
       "<Figure size 1440x1440 with 39 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# figure(figsize=[20,20])\n",
    "# for i in range(39):\n",
    "#     subplot(7,7,i+1)\n",
    "#     imshow(lensplanes[i][-1])\n",
    "#     axis('off')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "b0e36a9a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "39"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(lensplanes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "id": "3b16fb28",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7f5ae3154760>"
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQEAAAD7CAYAAABqkiE2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABi5klEQVR4nO29bax+y3UX9luzz/nfv50QbDexdbGt2kjmxaDSICvNS9VaGEpIoxg+pHXUVLetq6tKKQSEROzmQ9QPkSwVRUTqm66SELeYpK5JayuiJMYQoUqNiSERxHGMXUydSy6+oaUFhVz/z3n26oeZNbNmZs3s2fvZzzn7+p4lPXqePe/P3nv91m+tmdmbmBkP8iAP8soVd98DeJAHeZD7lQcQeJAHeYXLAwg8yIO8wuUBBB7kQV7h8gACD/Igr3B5AIEHeZBXuFwMBIjoW4nos0T0eSJ636X6eZAHeZDzhC6xToCIJgD/AMAfAfA8gF8A8F3M/Cu7d/YgD/IgZ8nVhdr9BgCfZ+Z/CABE9JMA3g3ABIFH7lX8qqvfFo4oz6SquJ1olisyu2V61bdWBHhV1ZH/3uuseZA3uwb31xgJrn4slFtMHGiukXnJOmvtpnEOx5owSp1hs/8F/tk/ZeavK9MvBQJvBPBr6vh5AP+GLkBEzwJ4FgAeT1+Nb/66f99nOCcFdOGQR+288M2unVcpsyuVzqgbhPVx6USV7e5UtwKQcrwAMKe7Iiq3vun071nKGfmj5fTveSFfpdFs9NMcZyO9Kjc3y3GzvRlFQbtcWVblVey5M44qv9WG1X85BgA8UKY1jr/OH/m/6oKXiwlYNiy3U8zPMfM7mPkdj9yr1I01SwFdOOR1bqTwTeZJUmXKi967kfVNzJyUYkZUFLPd4gJkyrSiLnFhtcv/Fo6zcuZ5wz4AoMdqnTedr65HFwCsa1KWsY4NxUrFOordaq/Ka7d/n0KWISBDjVew10uBwPMA3qyO3wTg14dr7wUExk1Z/dZtlnWM40WF7tTbXrcYawlerTakn3Ls1nlZBRQL+RYgj7Q5CMiV5VV5XQC4C8W2lM/VakZmuTEF3xsILgUCvwDgbUT0ViJ6BOA9AD62WMtC5i5at/O6jEB+d2+YTtsViBT1RllBb3zFccUKVHpV1+jD7Hep3BkA0OxnpM0uYA8ygCXpGYJSenlfAXKRmAAz3xLRfw7gZwBMAH6MmT/druAvIBH5i0OULrZzKS3mhYviKM/TQgSa2ccIJI+K336w6bduV/LKetJ2OI7+vtybAqu63aLPrK5VrxxfVrc+fT1ar/sbZgrNNs9gAD1WV6aXv2P/gzEAs+0OC7gPxgDPBqpxO6rHQ64aEzmq4wNGueo+MuRSgUEw818F8Fc3Vm5TGSuvkxaBoNdOqXgz59SsMx5izgN/M/oKrcpmdXU9q8/WGBas32EYQC+tlw7U1n/J4i8p9RrGcGAZBgIgGDK7nUOtGLQjpgPxAZ1m3eBWjGDp2PJRLZ8a5/r72wOO5ngs374cux4nVgBAqx3FAKogoFVHn9ula7IAAMyMVXGAJQAoy+/lChhxgXbZQZ8fK+IDve5Wlb4DYfMm3AgE6oZqBqmWjpcoa6HQmbSAwDheDSImYBrttVyALUHA1niwcH5H0izlHwCATLbQ+nv2980AYbPwCnVdUfY4IGDNwy4BgVFXNVIlVTfqiJW12u+UNVnBYJ+rAo6tsVjtjJQ7kwGY7Sy12RqbpfxrAcCSLXW2yBmLy6JYFr7VXavsIBAcBwSAYWVefXMVjKA7dbXkHgzS2ItOJVb0GrVlHy2n89co6xIAlGnluSvHCAxb/yEGsOQGjC66KYssuRM9abgETTZwrlvQKZ91s1jiTqQ+sWzeSEZ0eK2VCbKZFdQDbR5fbCqxMaaLMACrnb1cAC1brL8ea0+2Ku5RpwZ3BoKDgAASCiuLsTsQlFbXYgUtK2WxAiuv6Hu1ZW/VLesVaXfCAApG1R3/KAMorX9D+U3rv4UB6P9Y1r1HWcUGgC4QrA0WHgcEAPNCdGcM8oIqf8ySiVQLi0ZZwcplx5n07rkeK5hRAcLu6wAWmNTQQqxefuxvWWG3rrFv9zlY1yi3y47bNbMEsc55Pv9i87u0cglZEyjcSmkLVpBJp2zty57BCgbqSV1L7sQFGGEAZR2LAWT5F7T+ZZ+6/l3JHsHBoX5WugaGHAcE4g2lLuooEFjt6PpLND9I1z0o27b62sIKzggarnYBtgKm9De8X6EPJCaAx+KDyg+0KfwaABhkAbtKJ0C42i0A1rsGZdOLJe5blm683hqCVn2r/BpW0Mrr9VEcn8sKNjEAq9zeswAL42mt/TeVX4+vqtAwAOe2cWTZAAQjcrFlw5uEWS2xndMfC8t45SbJ9hgA/saSPQa+QJHPaZ+BSGsfQTiWm557+wisZce9/QfGcbV0GLCXHZft6Hz9e68YgAWKo+0182vAPlvxyz5G2um11ahzl2/qMvcUiFh7C2JFe8lwZAOnRpMbxri/6P+U3WjqD40ygmb+Cout5CJBw+J4MWjYYS1m22Wbl3QBLBlkAMO0H1gPAD1ZCQBd6exq7MqWAGGsuy8jOAYIAB1/euHmaVmaEUvVUi7DPdg0lTjSp/TRcw/K/6GPjTp3EgTUaSPTgEVatd6/5fP3An89BnAEF2ApONgBgsXlxDsCwXFAAFgGgqU1BEB7HYF1w5Z1q35rMBgtO8wKyj5KpR7Zf3AfQcCyPZ0/6gL0ov09y79V+ddsJ8663MAQdpCzgeBluWwY2JcRVG2st9DlcRcI6gHmxyvcg8Wg4ZILYPWzlgGU+WXa0orB0WfrActK2lP+cixW2z255GzAmVOFZwEBMAQExwkMGoFAn851sFAF3xaDhVUboe5SEE//LvI2Bw3lv0mfC2W7QUOdhnsMAprgYIMwm+Bxhq+/pLw7KP/FWYBzi3GFbqAQSPdvL2DYG0I3966lFQhcYAR5GwOMoJQ1QcM17kGvncG8kZWGdx4E3GA5zZt4ra+v+z8XAI4kA0HCoS3HKxYIaTkIEyiUfA9G4A9qRiDpIxa6VXYPVrCChTQfZYavEAYwovQjsuPzA3ZlANbj70rZgxEAy6zAqjJc8tLSukFGorxLf7i3Rr3nt1d9diz2OVOJvXpFrCDrc1RhW+1LO+cwgIU4RBcARvz8kZt5KZ7wcpG9GAGwihUcBwRGZIAqDgULy7yy3dGgYamklwCC4thcLbiksJ3g3iVXAg67AFk7fBnlH21zrYzO9w8r7xgQDLsHA2BwEHdghSz8+ezkZG/+cWXBdpsr3kzUHc9dbSI5Q9gllyZbaWm5JyJWmqazIT+jr7LSTa9qs56Eu/CWpXwcDYWxgKGlDI22LSVrr+IrxtGi9dbDcM9orxxjd5VhR47PBKzXigH5DbA2ILLquW7jSm4+1XiPMVjVt/isPTDrlRUx96mvAOX4GrkOIFsiFm3Qsq2ZI1/T7rgFdunTbzB9dmhPxrfquYU4EhNwDQW3pAEAQyygZ8l7DMA67o5xx7pbpLfGvJCMDbTG02MD8lvO8zxnVr7JCICaFVSDM8Y1yhbWMIUV7X5FsQQcCQTuQi6omKtYwF3LEr1vlbVEg8uKdk0gANDa9JKNR2RkrhxYBj8NDksxhcF2tdL15/RV373ZgOH/PdZejx0cBAQGaP4CPTVZwCgDKPtq9bcmf2vZS9TvtTsKCPp3CwhEyvgA0GYEQG2xW4ppnYdzmMKavq12B+MJu7KEPQCmrDZcshAiejMR/U0i+gwRfZqIvjekv46IPk5Enwvfr93ah93xgJ+1lLYWABaOKxawFkAOwCLMV7r3frdeAa8BeOlFnC1fXHx6/WmJ9qt75300rjDab9lmR4Z99RHff/j/urH2cF5g8BbAn2Xm3wvgGwF8DxG9HcD7AHyCmd8G4BPheL0ssQDo7DOs8p5xgKPJAcAlijqP1fUaCcyNKuZwsG0wIHhkQPANjgcXW1n9mm1h5heY+e+G3/8CwGcAvBHAuwF8MBT7IIA/PtzohoDgYjBw61SgdWykLbKAuwaRVn8D1n2YDURrP1LG5deC0nVrgsG
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "imshow(m[2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "id": "dd4183e8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "50"
      ]
     },
     "execution_count": 66,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "id": "1839df65",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABGYAAARNCAYAAAAQIkYGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAEAAElEQVR4nOz9u5Ity5I2Cn0ekaPmXGvtS/8XBBCRUVAxQ8XOeQ5EUMAMQwAVAQWJt0DDDDADVOARkLhIx/7+u3v37t57rzWrRmY4gl/CPTIyx6iac9aYzSk3qxp5jbx5eIR/8bkHMTM+5EM+5EM+5EM+5EM+5EM+5EM+5EM+5EM+5P2lPPoGPuRDPuRDPuRDPuRDPuRDPuRDPuRDPuRD/usqH8DMh3zIh3zIh3zIh3zIh3zIh3zIh3zIh3zIg+QDmPmQD/mQD/mQD/mQD/mQD/mQD/mQD/mQD3mQfAAzH/IhH/IhH/IhH/IhH/IhH/IhH/IhH/IhD5IPYOZDPuRDPuRDPuRDPuRDPuRDPuRDPuRDPuRB8gHMfMiHfMiHfMiHfMiHfMiHfMiHfMiHfMiHPEiWs53/xX/zf/T+c2m3b3TJQiAigAgoRf90nQhciy8jLDORwFWl9GUicCGg6H6K2+DbuBCY5Nq+rH9cCFx0meKylQU/nuM5YZudM67bfabzb61jX14/hnfn2DH/r//Z/5S+zQd6vfyX/63/Mfv3BPr3A9K3tX1caH9cPN5k2Md+zLD/1ccgX6sxyNSbWf4AoAFky/G3aXHxWD9n2KbL1Pbb0jnjNQCgtbTO4/EAwG04h/P2Np47HB+3xePG68QyAfDsHkK5/5f2v3uIPv4Xf/c/lBsoBFCwLUBfB0Bme0Y7E+0PkZ/PpbjN2dmjqnalEBCW+3a1GRVhfb+M+MaGV088bJgek/d1ndZl5sm2qOMAtbif5Rp6LDXudUXXyeoDs+y385qty36KejfWq1v1wMqO2/w5D9qlsP3/9P/53zxEF/97/+f/OS+loZaGS9mwlIaFGpayYW0VL63ieVvwvC5YW8F1K2itYG0F29b/eCOw/mIloBFo1b8N8tsgy02XG/Q96+uwNzD8zrYzdd2IepKWEfVnv5+4X9/vpbHfW/7joGPI+haW5fm4l7VpeaYfXhb7fjADG3dbHez3VMYhsbG++zq8P2FtmtV3e5d+/xt73fm//t/+lw/Rxf/j//u/wxdaUcEo1PCEDUUr+xe+4E/b7/Dn7Wf8efsZ/7L+jH9dP+Ov2yf85foZv65P+G294Lf1gi/rgufrgutasV7lj68FuBbQC6G8iF6WTfVxUx3d9FsH/Ul6FISHN5T0Yve71ynbT43lPpLuBT1KOsv7umL3Ntq8zfSriY61lm3aTOK+k/aDl9LbBdOnof1oFWiLrLeq2yvQFt2+AO0CtArwgtSn9e+wAf/P/9X/5GF9xv/7//e/zYUaKhgVjAs1FFWIZ674c/sJf26ij3/ZfsJfts/4tT3ht+2C53bBb9sFVy54ad1VakxoTHhpYk+3VvDbesG1FbysNdnWdS3YrhXtWoErgdYCrISyis6adPuYvy01mtg36nZvQ6gD0HLFFpTNjtdjzE6E+kEc/Y3Q3we83pQNbgPLxqAVKGuwfeMzAN7/6c9F3aeIOrZ0PWuL6FFbRCeb6pvpnOgbi84tAF8Y7akBFwY9bVieNnz6fMVPT1f88vSCn5YrPtcrnurm3/zX9Qn/+//+//Yh+vjf/T/8L9LHLQTU0kDE2FrBunW9aa2gNerV2V4kMUphPD2t+OXTC/7w9IzfP33Bz8sL/rA845flWcpAUT0taCBsTFhbxd9/+T3++ctP+MuXT/jt109YvyzAl4LyXKR9j20swjIPdrb1+p1sodd7Fl2M23jQP/Sy/UFVT3hoH4lDuVu3rf67spff+xz2zva+WPfV5ThexMa1C2G7kOjYpds/0dfw5zrJQW8ZfJE/PDWUpw2XpxWXy4anZcWlNiyloRCjEOPaCv4f/4P/9VQXT4GZbyLfCmh5w3W5QMAZQD98fwfUGFy7w+ydIRPmtE6sTmIdjoleTlxlBoG6QgyvwQ1iOJew3xaLJ9U7UFing/Nvrd9zD5P1hwtrB6WUVJl9n62zNCRcaLcdQO80237qv8QsldqPGc6/+xh9d1oH7gJlhvvbgTIzpzKCMlOAZ9gWl49AmShHoMw9Et95KX49IurXKvQ4O/FGYeZgWxrQAqDi69T1tdZuZ6D2JzY+DUCBAAusdqNq/WsEoIGCNzeYBhUKe5IxysunwMz4oPtnPwZmOHX4+v7eIKdtsXxzRkYHZdDZ3f0BwU7neozC+u70kW2fvcaJ48zV9JU6SBOvMdaPo+3vKJ+WFQsJMCOgTAdnCrF0zFrDVho2JqluAJhJ/lr4OwFl3AFYg2MQnAY2rCB2hHTd25EI0IR2LDrBaVvsKM62RScjOM9+bwGU2R2bygidRgVfysYK/CGBMm5zx09eAGz2fNTbBIRtB+t5QAgK3MoxyWkegNgMKKmz9EBb+vvym4My0Rne9GN/pisutMk+1U+TBnInYmuErWXQEGvQSdXHch0AGVsOduVIyP+FY7nrTNKpoO9lw07f4nHiKFh7bDrXbZtfvwz6YNtv2BIe9Sg5MVapOrBnupQAmYWSgxzBfANhzGkWpwVgc54dnAkOyRLqPPVvUdaHYTIAgM+0uo5VMAoYT9oQXVEcNDRpIFy54soVz23Bc1vw0ipetgVN360BM9dWHZh52SquW8G6VTCgTrU6gUyqW9TtXCPXN5DqSOhsO2gd7KsB4jvbtgZwZkVwXHtd6OAydm2264nqSLTfxEAjRlEjzgygsjZ3+dvuQBg7hNAHrwPoZ+BLBmWoO73BIW4L0J7Yf5MDfNlQLxs+fRKw4penF/zx6Td8rit+qlf8VK9+j5/q+jZF+gZyUN2Tngggo6BMEz1znbDud5Hx6efS8JsCO/triW290DXZ2Oe24LpVXLeK9aliWxWwrdz1UYE/jl3FSb9wB1IHUGYEblJ7y5PytG9moAxXvYzaOmufi+onFbUvVoEKS52Keoegk4WyfhqInGzgAAZGMLoCKGoHS7B7BkyrPvIiQKGBMp+eVvz0JODgp2XFp9rt0ctW5wqBrwFmfgRHitvxPtLWahyZjwDMAL7swRjt2If35w65lwlQ4aBAg84Fn0iAgr5+BMbEbSMYM5Ypx4d74pNjZ9e0m03HWKtwUMaDJDnDBs7IDnVsAwhwBM6EfSAFBMb9CK/XwJeCfo4BIXcdM7y4R4MyMwbLDJQZGTFngM23FipeNhXqrJmwPQE+j5DWwKVIo1QIR+AMq3+Vvq3qidsDLxMduAHAG4AiOs8o2kuS8gUohnaC43sg9Ao7GAo3RGFTPHP2Om8dE8AY389Rr0PDPQNjeChnBGXs9x5VKwC3ocapczwFaAIIgxbqrXUSDITn8f2O7xzYOd3vKJ/qmgCZQoyFhD1TWkVjwloKqDEqMbbwEcWEiJPAWwE26qDMNoAyCtKUFb1jZs5Ds3erzUcAXiLr0iWsJ/2JbAfGvjOnyw7iqH5MARnGjumQwMMIyGxh2RkPtj4wtuKnH/SAaWfxsQNgogTGpXQOuxPtrDhj1AZmQ1t6J9RHLSvtRrHfW35RR8Cc4KrLJgLMrLjQigLOwAyTO7rrVgdQRnSTVkKJoMyK5IQaYwA4sGcY+lvD5+gMsPtAvpEl0x2UsE0eLjskhN6uTTy2W+DMjqk7YemyMjETmBcAmTaCM+Z4mMNsI8YRhDHHRR0TGzUWxky/Z9JvxcsD22gAn1QZagJn4MsVDUUblw2EjYVpsHLFysVBmS/bIo6zKowwD6swEbaKdSvY1LE2m8oMwEEZvWi0YWEb275oExF0zFiKgQVDDGWvIDBmOihTkh7263qfcby+gSil90sYUGaPMCBaBQobKBz7F9nOJ7ZCdIIj8Ge/l8DCGkEZc4yNKfOUWTLl0nC5bLhcVvz8dMXvPz3jj0+/4Y+XL/hlecansuLn+iLflws+lQ7SvLfU0gRogegHpXZYQBgZx1NQphVw6+2zHEhgYqwAXuqC30q2oQD
      "text/plain": [
       "<Figure size 1440x1440 with 49 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "figure(figsize=[20,20])\n",
    "for i in range(49):\n",
    "    subplot(7,7,i+1)\n",
    "    imshow(m[i+1])\n",
    "    axis('off')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bace2afc",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}