doc/src/examples/notebooks/dataset_representation/meshes_with_scatter_plots.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "56cc26a2",
   "metadata": {},
   "source": [
    "# Représentation de scores sur maillages\n",
    "\n",
    "Par défaut la représentation des scores sur maillages dans valjean se fait centrée sur la maille. Si le score Tripoli4 a été tourné avec l'option ``MESH_INFO`` les coordonnées du centre de la maille sont en plus disponibles. Il est alors possible de faire un graphique davantage à l'échelle, même si les axes des mailles ne correspondent pas aux axes par défaut.\n",
    "\n",
    "Différents exemples sont présentés pour illustrer ces représentations. Elles sont faites grâce à *matplotlib*, hors *valjean*.\n",
    "\n",
    "## Les imports par défaut"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ab2c1a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from valjean.eponine.tripoli4.parse import Parser\n",
    "from valjean.eponine.dataset import Dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2eaae31d",
   "metadata": {},
   "source": [
    "## Quelques fonctions\n",
    "\n",
    "Certaines fonctionnalités ont été encapulées dans des fonctions pour plus de simplicité. Ce sont des exemples.\n",
    "\n",
    "### Sélection des intervalles\n",
    "\n",
    "Cette fonction permet de sélectionner les intervalles que l'on veut représenter et permet de *squeezer* le ``Dataset`` et les coordonnées en même temps (suppresion des dimensions triviales). La sélection se fait sur un dictionnaire dont les clefs correspondent au nom des coordonnées à réduire et dont les valeurs correspondent aux indices des intervalles à sélectionner. Deux cas sont possibles :\n",
    "\n",
    "- seul le deuxième intervalle est gardé sur l'axe `'u'` : `{'u': 1}`\n",
    "- les intervalles d'indices 3 à 8 sont gardés sur l'axe `'v'` (qui en a au moins 9) : `{'v': (3, 9)}`\n",
    "\n",
    "Dans le dernier cas les indices sont donnés à une `slice` ([slice doc](https://docs.python.org/3/library/functions.html#slice)) dans la fonction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be77f9bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "def select_space_bins(dset, coords, selection):\n",
    "    \"\"\"Select slice to be represented and squeeze dataset and coordinates accordingly.\n",
    "    \n",
    "    The slice can be done on one value or on a range of values if a tuple is given.\n",
    "    \n",
    "    :param Dataset dset: dataset\n",
    "    :param np.ndarray coords: coordinates (structured array which names are the axis names)\n",
    "    :param dict selection: selection to be applied, keys correspond to the axis names\n",
    "    :rtype: Dataset, np.ndarray\n",
    "    :returns: sliced and squeezed dataset and coordinates\n",
    "    \"\"\"\n",
    "    islices = tuple(slice(None) for _ in range(dset.ndim))\n",
    "    oslices = list(islices)\n",
    "    for key, val in selection.items():\n",
    "        icoord = list(dset.bins.keys()).index(key)\n",
    "        if isinstance(val, int):\n",
    "            oslices[icoord] = slice(val, val+1)\n",
    "        elif isinstance(val, (list, tuple)):\n",
    "            oslices[icoord] = slice(val[0], val[1])\n",
    "    oslices = tuple(oslices)\n",
    "    return dset[oslices].squeeze(), coords[oslices[:3]].squeeze()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b16efe05",
   "metadata": {},
   "source": [
    "### Représentation des mailles individuelles\n",
    "\n",
    "Le choix fait ici est de donner à la maille sa forme \"réelle\". Pour cela le marqueur utilisé est défini comme un chemin entre les différents vertex donnés ([marker doc](https://matplotlib.org/stable/api/markers_api.html#module-matplotlib.markers) et [marker examples](https://matplotlib.org/stable/gallery/lines_bars_and_markers/scatter_star_poly.html)).\n",
    "\n",
    "Dans les cas présentés ici le maillage est supposé régulier (une petite vérification est faite, cela suppose en général une précision dégradée sur l'estimation de la distance entre deux points consécutifs, ce qui apparaît dans l'argument `decimals` demandé par la fonction). Le chemin proposé est calculé à partir des 4 premiers points (2 sur chacune des 2 dimensions utilisées). Un parallélogramme est ainsi défini. Le sens de parcours des sommets est le sens direct pour s'assurer d'avoir un quadrilatère convexe.\n",
    "\n",
    "Les dimensions fournies sont celles qui sont gardées."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bbb13ca5",
   "metadata": {},
   "outputs": [],
   "source": [
    "def vertices(coords, decimals, dims=('u', 'v')):\n",
    "    \"\"\"Compute the vertices of the elementary cell from regular coordinates.\n",
    "    \n",
    "    :param np.ndarray coords: coordinates (structured array which names are axis names)\n",
    "    :param int decimals: number of decimals to consider to identify unique coordinates values\n",
    "    :param list, tuple dims: dimension names of the corrdinated considered\n",
    "    :rtype: list(tuple(float))\n",
    "    :returns: verticies of the path matching the cell (4 points here)\n",
    "    \"\"\"\n",
    "    uu = np.unique(np.diff(coords[dims[0]]).round(decimals))\n",
    "    uv = np.unique(np.diff(coords[dims[1]]).round(decimals))\n",
    "    if len(uu) != len(uv):\n",
    "        print(\"Not a rectangle or parallelogram, returning arbitrary triangle\")\n",
    "        print(f'u: {uu}, v: {uv}')\n",
    "        return [[-1, -1], [1, -1], [1, 1], [-1, -1]]\n",
    "    if len(uu) > 2:\n",
    "        print('Not a parallelogram, returning arbitrary triangle')\n",
    "        print(uu, np.ediff1d(coords[dims[0]]).round(decimals))\n",
    "        return [[-1, -1], [1, -1], [1, 1], [-1, -1]]\n",
    "    f4pts = coords[:2, :2]\n",
    "    meanu = np.mean([np.min(f4pts[dims[0]]), np.max(f4pts[dims[0]])])\n",
    "    meanv = np.mean([np.min(f4pts[dims[1]]), np.max(f4pts[dims[1]])])\n",
    "    # vertices: points rotated in direct order\n",
    "    return [(f4pts[dims[0]][0, 0]-meanu, f4pts[dims[1]][0, 0]-meanv),\n",
    "            (f4pts[dims[0]][1, 0]-meanu, f4pts[dims[1]][1, 0]-meanv),\n",
    "            (f4pts[dims[0]][1, 1]-meanu, f4pts[dims[1]][1, 1]-meanv),\n",
    "            (f4pts[dims[0]][0, 1]-meanu, f4pts[dims[1]][0, 1]-meanv)]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07f0c741",
   "metadata": {},
   "source": [
    "## Cas d'une maille cartésienne alignée sur les axes\n",
    "\n",
    "Lecture du jeu de données"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef4efa35",
   "metadata": {},
   "outputs": [],
   "source": [
    "t4p_bd = Parser(\"box_dyn.res\")\n",
    "t4b_bd = t4p_bd.parse_from_index().to_browser()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c742f0a",
   "metadata": {},
   "source": [
    "Sélection du score (mesh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "99887459",
   "metadata": {},
   "outputs": [],
   "source": [
    "neut_flux = t4b_bd.select_by(score_name=\"neutron_flux_mesh_score\")\n",
    "print(neut_flux['results']['score_eintegrated'].shape)\n",
    "nfm = neut_flux['results']['score_eintegrated'].squeeze()\n",
    "print(nfm)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad0ac183",
   "metadata": {},
   "source": [
    "Sélection / accès aux coordonnées associées"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7acd1afe",
   "metadata": {},
   "outputs": [],
   "source": [
    "nfc = neut_flux['results']['coordinates']\n",
    "print(nfc.shape, nfc.dtype.names)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad181144",
   "metadata": {},
   "source": [
    "Les représentations proposées ici sont en 2 dimensions. Pour représenter complètement le score neuf graphiques sont nécessaires : 3 pour chaque couple de coordonnées. Pour simplifier seul un graphique pour chaque couple de données sera présenté."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61fde55f",
   "metadata": {},
   "outputs": [],
   "source": [
    "nfm_u0, nfc_u0 = select_space_bins(nfm, nfc, {'u': 0})\n",
    "nfm_v0, nfc_v0 = select_space_bins(nfm, nfc, {'v': 0})\n",
    "nfm_w0, nfc_w0 = select_space_bins(nfm, nfc, {'w': 0})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "951c2d83",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(nfm_u0)\n",
    "print(nfc_u0.shape, nfc.dtype.names)\n",
    "print('coordonnées restantes pour u=0')\n",
    "print(nfc_u0)\n",
    "print('coordonnées restantes pour w=0')\n",
    "print(nfc_w0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e8ad88c",
   "metadata": {},
   "source": [
    "Calcul des sommets (vertex) du parallélogramme représentant la maille (ici un rectangle en fait)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "159b5a8a",
   "metadata": {},
   "outputs": [],
   "source": [
    "nfv_u0 = vertices(nfc_u0, 3, ['v', 'w'])\n",
    "nfv_v0 = vertices(nfc_v0, 3, ['u', 'w'])\n",
    "nfv_w0 = vertices(nfc_w0, 3, ['u', 'v'])\n",
    "print(nfv_u0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20681e67",
   "metadata": {},
   "source": [
    "Représentation graphique grâce à *matplotlib*.\n",
    "\n",
    "Le graphique est représenté 3 fois :\n",
    "\n",
    "- marqueur correspondant au chemin défini plus haut en grisé transparent pour pouvoir ajuster sa taille (argument `s`), la transparence permettant de mieux visualiser les recouvrements\n",
    "- marqueur correspondant au chemin défini plus haut avec le contenu attendu (argument `c`=couleur contenant la valeur du `Dataset`)\n",
    "- marqueur + rouge pour visualiser le centre de la maille"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ff55ee8",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(1, 3, figsize=(18, 6), subplot_kw={'aspect': \"equal\"}, constrained_layout=True)\n",
    "axs[0].scatter(nfc_u0['v'].flatten(), nfc_u0['w'].flatten(),\n",
    "                  c=\"gray\", s=30000, marker=nfv_u0, alpha=0.3)\n",
    "im0 = axs[0].scatter(nfc_u0['v'].flatten(), nfc_u0['w'].flatten(),\n",
    "                  c=nfm_u0.value.flatten(), s=28000, marker=nfv_u0)\n",
    "axs[0].scatter(nfc_u0['v'].flatten(), nfc_u0['w'].flatten(), c=\"red\", marker=\"+\")\n",
    "axs[0].set(xlabel='v', ylabel='w')\n",
    "cbar0 = fig.colorbar(im0, ax=axs[0], label='neutron flux')\n",
    "axs[1].scatter(nfc_v0['u'].flatten(), nfc_v0['w'].flatten(),\n",
    "                  c=\"gray\", s=30000, marker=nfv_v0, alpha=0.2)\n",
    "im1 = axs[1].scatter(nfc_v0['u'].flatten(), nfc_v0['w'].flatten(),\n",
    "                  c=nfm_v0.value.flatten(), s=27000, marker=nfv_v0, alpha=0.9)\n",
    "axs[1].scatter(nfc_v0['u'].flatten(), nfc_v0['w'].flatten(), c=\"red\", marker=\"+\")\n",
    "axs[1].set(xlabel='u', ylabel='v')\n",
    "cbar1 = fig.colorbar(im1, ax=axs[1], label='neutron flux')\n",
    "axs[2].scatter(nfc_w0['u'].flatten(), nfc_w0['v'].flatten(),\n",
    "               c=\"gray\", s=24000, marker=nfv_w0, alpha=0.2,)\n",
    "im2 = axs[2].scatter(nfc_w0['u'].flatten(), nfc_w0['v'].flatten(),\n",
    "                  c=nfm_w0.value.flatten(), s=20000, marker=nfv_w0)\n",
    "axs[2].scatter(nfc_w0['u'].flatten(), nfc_w0['v'].flatten(), c=\"red\", marker=\"+\")\n",
    "axs[2].set(xlabel='u', ylabel='v')\n",
    "cbar2 = fig.colorbar(im2, ax=axs[2], label='neutron flux')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4fa4263",
   "metadata": {},
   "source": [
    "Un des arguments les plus imporants ici est : `subplot_kw={'aspect': \"equal\"}` sans lui les axes sont optimisés et ne sont plus à l'échelle. L'effet sera plus flagrant plus loin, sur des mailles non rectangulaires."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9bf4803",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(1, 3, figsize=(18, 6))\n",
    "axs[0].scatter(nfc_u0['v'].flatten(), nfc_u0['w'].flatten(),\n",
    "                  c=\"gray\", s=22000, marker=nfv_u0, alpha=0.2)\n",
    "im0 = axs[0].scatter(nfc_u0['v'].flatten(), nfc_u0['w'].flatten(),\n",
    "                  c=nfm_u0.value.flatten(), s=20000, marker=nfv_u0)\n",
    "axs[0].scatter(nfc_u0['v'].flatten(), nfc_u0['w'].flatten(), c=\"red\", marker=\"+\")\n",
    "axs[0].set(xlabel='v', ylabel='w')\n",
    "cbar0 = fig.colorbar(im0, ax=axs[0], label='neutron flux')\n",
    "im1 = axs[1].scatter(nfc_v0['u'].flatten(), nfc_v0['w'].flatten(),\n",
    "                  c=nfm_v0.value.flatten(), s=21000, marker=nfv_v0, alpha=0.9)\n",
    "axs[1].scatter(nfc_v0['u'].flatten(), nfc_v0['w'].flatten(), c=\"red\", marker=\"+\")\n",
    "axs[1].set(xlabel='u', ylabel='v')\n",
    "cbar1 = fig.colorbar(im1, ax=axs[1], label='neutron flux')\n",
    "axs[2].scatter(nfc_w0['u'].flatten(), nfc_w0['v'].flatten(),\n",
    "               c=\"gray\", s=11000, marker=nfv_w0, alpha=0.2,)\n",
    "im2 = axs[2].scatter(nfc_w0['u'].flatten(), nfc_w0['v'].flatten(),\n",
    "                  c=nfm_w0.value.flatten(), s=10000, marker=nfv_w0)\n",
    "axs[2].scatter(nfc_w0['u'].flatten(), nfc_w0['v'].flatten(), c=\"red\", marker=\"+\")\n",
    "axs[2].set(xlabel='u', ylabel='v')\n",
    "cbar2 = fig.colorbar(im2, ax=axs[2], label='neutron flux')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4c6753e",
   "metadata": {},
   "source": [
    "## Cas d'une maille cartésienne dont les axes ont été tournés\n",
    "\n",
    "Lecture du jeu de données"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e45e6cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "t4p_em = Parser(\"extended_mesh_cartesian_info.res\")\n",
    "t4b_em = t4p_em.parse_from_index().to_browser()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45b99891",
   "metadata": {},
   "source": [
    "Sélection du score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "afab1a60",
   "metadata": {},
   "outputs": [],
   "source": [
    "phot_flux_m8 = t4b_em.select_by(score_name=\"mesh8_reg\")\n",
    "pfm8 = phot_flux_m8['results']['score_eintegrated']\n",
    "cm8 = phot_flux_m8['results']['coordinates']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d3084623",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(pfm8.shape, cm8.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab367463",
   "metadata": {},
   "source": [
    "La représentation sera également faite dans le premier intervalle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98f888f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "pfm8u0, cm8u0 = select_space_bins(pfm8, cm8, {'u': 0})\n",
    "pfm8v0, cm8v0 = select_space_bins(pfm8, cm8, {'v': 0})\n",
    "pfm8w0, cm8w0 = select_space_bins(pfm8, cm8, {'w': 0})\n",
    "print(pfm8u0.shape, cm8u0.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6aca9b9",
   "metadata": {},
   "source": [
    "Calcul des sommets des parallélogrammes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d9f5440",
   "metadata": {},
   "outputs": [],
   "source": [
    "verts_u0 = vertices(cm8u0, 3, ['v', 'w'])\n",
    "verts_v0 = vertices(cm8v0, 3, ['u', 'w'])\n",
    "verts_w0 = vertices(cm8w0, 3, ['u', 'v'])\n",
    "print(verts_u0)\n",
    "print(verts_v0)\n",
    "print(verts_w0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae5a0e2b",
   "metadata": {},
   "source": [
    "Représentation **sans** `aspect='equal'`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5042ea02",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(1, 3, figsize=(18, 6))\n",
    "axs[0].scatter(cm8u0['v'].flatten(), cm8u0['w'].flatten(),\n",
    "                  c=pfm8u0.value.flatten(), s=35000, marker=verts_u0)\n",
    "axs[0].scatter(cm8u0['v'].flatten(), cm8u0['w'].flatten(), c=\"red\", marker=\"+\")\n",
    "axs[0].set(xlabel='v', ylabel='w')\n",
    "axs[1].scatter(cm8v0['u'].flatten(), cm8v0['w'].flatten(),\n",
    "                  c=pfm8v0.value.flatten(), s=18000, marker=verts_v0)\n",
    "axs[1].scatter(cm8v0['u'].flatten(), cm8v0['w'].flatten(), c=\"red\", marker=\"+\")\n",
    "axs[1].set(xlabel='u', ylabel='w')\n",
    "axs[2].scatter(cm8w0['u'].flatten(), cm8w0['v'].flatten(),\n",
    "                  c=pfm8w0.value.flatten(), s=30000, marker=verts_w0)\n",
    "axs[2].scatter(cm8w0['u'].flatten(), cm8w0['v'].flatten(), c=\"red\", marker=\"+\")\n",
    "axs[2].set(xlabel='u', ylabel='v')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dfda85ea",
   "metadata": {},
   "source": [
    "Les parallélogrammes se chevauchent, cela ne correspond pas à ce que l'on attend -> `aspect='equal'`..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f312f55",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(1, 3, figsize=(18, 6), subplot_kw={'aspect': \"equal\"})\n",
    "axs[0].scatter(cm8u0['v'].flatten(), cm8u0['w'].flatten(),\n",
    "                  c=\"gray\", alpha=0.2, s=40000, marker=verts_u0)\n",
    "axs[0].scatter(cm8u0['v'].flatten(), cm8u0['w'].flatten(),\n",
    "                  c=pfm8u0.value.flatten(), s=30000, marker=verts_u0)\n",
    "axs[0].scatter(cm8u0['v'].flatten(), cm8u0['w'].flatten(), c=\"red\", marker=\"+\")\n",
    "axs[0].set(xlabel='v', ylabel='w')\n",
    "axs[1].scatter(cm8v0['u'].flatten(), cm8v0['w'].flatten(),\n",
    "                  c=\"gray\", alpha=0.2, s=13000, marker=verts_v0)\n",
    "axs[1].scatter(cm8v0['u'].flatten(), cm8v0['w'].flatten(),\n",
    "                  c=pfm8v0.value.flatten(), s=10000, marker=verts_v0)\n",
    "axs[1].scatter(cm8v0['u'].flatten(), cm8v0['w'].flatten(), c=\"red\", marker=\"+\")\n",
    "axs[1].set(xlabel='u', ylabel='w')\n",
    "axs[2].scatter(cm8w0['u'].flatten(), cm8w0['v'].flatten(),\n",
    "                  c=\"gray\", alpha=0.2, s=31000, marker=verts_w0)\n",
    "axs[2].scatter(cm8w0['u'].flatten(), cm8w0['v'].flatten(),\n",
    "                  c=pfm8w0.value.flatten(), s=25000, marker=verts_w0)\n",
    "axs[2].scatter(cm8w0['u'].flatten(), cm8w0['v'].flatten(), c=\"red\", marker=\"+\")\n",
    "axs[2].set(xlabel='u', ylabel='v')"
   ]
  }
 ],
 "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.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}