Inno4Scale.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "45734234-7021-48e2-b3e0-3fab28ec47d7",
   "metadata": {},
   "source": [
    "## Vlasiator - A Global Hybrid-Vlasov Simulation Model \n",
    "Vlasiator [@palmroth2018] is an open-source simulation software used to model the behavior of plasma in the Earth's magnetosphere, a region of space where the solar wind interacts with the Earth’s magnetic field. Vlasiator models collisionless space plasma dynamics by solving the 6-dimensional Vlasov equation, using a hybrid-Vlasov approach. It uses a 3D Cartesian grid in real space, with each cell storing another 3D Cartesian grid in velocity space. The velocity mesh contained in each spatial cell in the simulation domain has been represented so far by a sparse grid approach, fundamentally based on an associative container such as a key-value hashtable. Storing a 3D VDF at every spatial cells increases the memory requirements exponentially both during runtime and for storing purposes. Our proposal revolves around developing innovative solutions to compressing the VDFs during runtmime."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "796e2f0c-f6d8-4600-9e57-80a57358bc7d",
   "metadata": {},
   "source": [
    "![title](images/egi.png)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ac9891a-860b-4411-8326-1dd98e4b5717",
   "metadata": {},
   "source": [
    "## VDF Compression\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ffaa767-6fdc-42f4-a9c1-93009df7cfeb",
   "metadata": {},
   "source": [
    "### Let's read in a vdf from a sample file and see what that looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aae7d622-c42c-470b-8d9d-509efb1719e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys,os\n",
    "# sys.path.append('/home/mjalho/analysator')\n",
    "import tools as project_tools\n",
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "mpl.rcParams.update(mpl.rcParamsDefault)\n",
    "import matplotlib.colors as colors\n",
    "# plt.rcParams['figure.figsize'] = [7, 7]\n",
    "import ctypes\n",
    "import pyzfp,zlib\n",
    "import mlp_compress\n",
    "\n",
    "import pytools"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be5f3337-4d10-4115-9927-2d456b737d2f",
   "metadata": {},
   "outputs": [],
   "source": [
    "file=\"/home/kstppd/Desktop/bulk1.0001280.vlsv\";cid=356780649;\n",
    "#Read the VDF into a 3D uniform mesh and plot it\n",
    "vdf=project_tools.extract_vdf(file,cid,25)\n",
    "# np.save(\"sample_vdf.bin\",np.array(vdf,dtype=np.double));\n",
    "np.array(vdf,dtype=np.double).tofile(\"sample_vdf.bin\")\n",
    "nx,ny,nz=np.shape(vdf)\n",
    "print(f\"VDF shape = {np.shape(vdf)}\")\n",
    "fig,(ax1, ax2) = plt.subplots(1, 2)\n",
    "cax = fig.add_axes([0.95,0.25,0.05,0.5])\n",
    "im1 = ax1.imshow(vdf[:,:,nz//2],norm=colors.LogNorm(vmin=1e-15,vmax=3e-13))\n",
    "im2 = ax2.imshow(vdf[:,ny//2,:],norm=colors.LogNorm(vmin=1e-15,vmax=3e-13))\n",
    "\n",
    "\n",
    "fig.colorbar(im1, cax=cax)\n",
    "fig.suptitle(\"Original VDF\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f75b96f-3ba2-49d0-b3bf-c30966de5abe",
   "metadata": {},
   "source": [
    "### The vdf shown above is sampled on a uniform 3D velocity mesh and contains 64bit floating point numbers that represent the phase space density. We can calculate the total size of this VDF is bytes using ```sys.getsizeof(vdf)```."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e456c54-ab54-4dd4-8112-ebe75da96c96",
   "metadata": {},
   "outputs": [],
   "source": [
    "vdf_mem=sys.getsizeof(vdf)\n",
    "num_stored_elements=len(vdf[vdf>1e-15])\n",
    "print(f\"VDF takes {vdf_mem} B.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74ae8266-9333-445b-9a07-2702795fa0fd",
   "metadata": {},
   "source": [
    "### Now in Vlasiator we have countlesss VDFs since there is one per spatial cell. It would be great if we could compress them efficiently. We can try to do so by using zlib which is a form of lossless compression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed5662d9-8f88-43f7-b1f6-563d8fde0f81",
   "metadata": {},
   "outputs": [],
   "source": [
    "compressed_vdf = zlib.compress(vdf)\n",
    "compressed_vdf_mem=len(compressed_vdf)\n",
    "compression_ratio=vdf_mem/compressed_vdf_mem\n",
    "print(f\"Achieved compression ratio using zlib= {round(compression_ratio,2)}.\")\n",
    "decompressed_vdf = zlib.decompress(compressed_vdf)\n",
    "recon = np.frombuffer(decompressed_vdf, dtype=vdf.dtype).reshape(vdf.shape)\n",
    "project_tools.plot_vdfs(vdf,recon)\n",
    "project_tools.print_comparison_stats(vdf,recon)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2e4811c-6d2b-482f-bdeb-6c60e6453e1e",
   "metadata": {},
   "source": [
    "### We can use a lossy compression method like zfp[@zfp] to get even higher compression ratios."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d4c6ba0b-b56a-4ef5-8c97-7b289762e570",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Compresses a VDF  using ZFP (Zstandard Compressed FP)\n",
    "Input:VDF - numpy array \n",
    "Output: recon (Reconstructed VDF) - numpy array \n",
    "\"\"\"\n",
    "tolerance = 1e-13\n",
    "compressed_vdf = pyzfp.compress(vdf, tolerance=tolerance)\n",
    "compressed_vdf_mem=len(compressed_vdf)\n",
    "compression_ratio=vdf_mem/compressed_vdf_mem\n",
    "print(f\"Achieved compression ratio using zfp= {round(compression_ratio,2)}.\")\n",
    "recon = pyzfp.decompress(compressed_vdf,vdf.shape,vdf.dtype,tolerance)\n",
    "project_tools.plot_vdfs(vdf,recon)\n",
    "project_tools.print_comparison_stats(vdf,recon)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e59226f-a344-47f7-aa2b-d016b9e31ce3",
   "metadata": {},
   "source": [
    "### We will compress the VDF using an MLP. [@park2019]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "169f5b61-cda9-48d0-a550-b103f6355e5c",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Compresses a VDF using an MLP (Multilayer Perceptron).\n",
    " Input: \"sample_vdf.bin\" - Binary file containing the VDF data\n",
    "        order - Order of the fourier features\n",
    "        epochs - Number of training epochs for the MLP model\n",
    "        n_layers - Number of layers in the MLP model\n",
    "        n_neurons - Number of neurons in each layer of the MLP model\n",
    " Output: recon (Reconstructed VDF) - NumPy array representing the reconstructed volume data\n",
    "\"\"\"\n",
    "order=0\n",
    "epochs=10\n",
    "n_layers=4\n",
    "n_neurons=25\n",
    "recon=mlp_compress.compress_mlp(\"sample_vdf.bin\",order,epochs,n_layers,n_neurons)\n",
    "recon=np.array(recon,dtype=np.double)\n",
    "recon= np.reshape(recon,np.shape(vdf),order='C')\n",
    "project_tools.plot_vdfs(vdf,recon)\n",
    "project_tools.print_comparison_stats(vdf,recon)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c753cd6c-eb1c-4f88-bb5b-2781a3e4d351",
   "metadata": {},
   "source": [
    "### We will compress the VDF using an MLP with Fourier Features. [@2020fourier]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "afc424ea-8b00-4dd3-bb3c-af0a20daeba2",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Compresses a VDF using an MLP (Multilayer Perceptron).\n",
    " Input: \"sample_vdf.bin\" - Binary file containing the VDF data\n",
    "        order - Order of the fourier features\n",
    "        epochs - Number of training epochs for the MLP model\n",
    "        n_layers - Number of layers in the MLP model\n",
    "        n_neurons - Number of neurons in each layer of the MLP model\n",
    " Output: recon (Reconstructed VDF) - NumPy array representing the reconstructed volume data\n",
    "\"\"\"\n",
    "order=16\n",
    "epochs=12\n",
    "n_layers=4\n",
    "n_neurons=25\n",
    "recon=mlp_compress.compress_mlp(\"sample_vdf.bin\",order,epochs,n_layers,n_neurons)\n",
    "recon=np.array(recon,dtype=np.double)\n",
    "recon= np.reshape(recon,np.shape(vdf),order='C')\n",
    "project_tools.plot_vdfs(vdf,recon)\n",
    "project_tools.print_comparison_stats(vdf,recon)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e8837a4-50cd-422d-8c70-6df8771f514c",
   "metadata": {},
   "source": [
    "### Now we use a Spherical Harmonic Decomposition to perform the compression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3815300e-e181-4705-8e21-58e7d5b533a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Compresses a VDF using a spherical harmonic decompostion\n",
    "Input: \"sample_vdf.bin\" - Binary file containing the VDF data\n",
    "       degree - Degree of the spherical harmonic decomposition (l)\n",
    "Output: recon (Reconstructed VDF) - NumPy array representing the reconstructed volume data\n",
    "\"\"\"\n",
    "\n",
    "degree=10\n",
    "recon=mlp_compress.compress_sph(\"sample_vdf.bin\",degree)\n",
    "recon=np.array(recon,dtype=np.double)\n",
    "recon= np.reshape(recon,np.shape(vdf),order='C')\n",
    "project_tools.plot_vdfs(vdf,recon)\n",
    "project_tools.print_comparison_stats(vdf,recon)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "136ceeae-4283-42f0-84e5-c60017b17e80",
   "metadata": {},
   "source": [
    "### Now we use a CNN to perform the compression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad5fb16b-01ba-450b-81cc-0f9a65717323",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Function: train_and_reconstruct\n",
    "\n",
    "Description:\n",
    "This function takes an input array and trains a Convolutional Neural Network (CNN) model to reconstruct the input array.\n",
    "It uses Mean Squared Error (MSE) loss and the Adam optimizer for training.\n",
    "\n",
    "Inputs:\n",
    "- input_array (numpy array): The input array to be reconstructed.\n",
    "- num_epochs (int, optional): The number of training epochs.\n",
    "- learning_rate (float, optional): The learning rate for the Adam optimize.\n",
    "Outputs:\n",
    "    Reconstructed vdf array \n",
    "    Size of model in bytes\n",
    "\"\"\"\n",
    "\n",
    "import torch\n",
    "import torch.nn as nn\n",
    "import torch.optim as optim\n",
    "import numpy as np\n",
    "\n",
    "class CNN(nn.Module):\n",
    "    def __init__(self):\n",
    "        super(CNN, self).__init__()\n",
    "        self.conv1 = nn.Conv3d(1, 16, kernel_size=3, padding=1)\n",
    "        self.conv2 = nn.Conv3d(16, 32, kernel_size=3, padding=1)\n",
    "        self.conv3 = nn.Conv3d(32, 64, kernel_size=3, padding=1)\n",
    "        self.conv4 = nn.Conv3d(64, 1, kernel_size=3, padding=1)\n",
    "        self.relu = nn.ReLU()\n",
    "        \n",
    "    def forward(self, x):\n",
    "        x = self.relu(self.conv1(x))\n",
    "        x = self.relu(self.conv2(x))\n",
    "        x = self.relu(self.conv3(x))\n",
    "        x = self.conv4(x)\n",
    "        return x\n",
    "\n",
    "def train_and_reconstruct(input_array, num_epochs=30, learning_rate=0.001):\n",
    "    device = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\n",
    "    input_tensor = torch.tensor(input_array, dtype=torch.float32).unsqueeze(0).unsqueeze(0).to(device)  # Add batch and channel dimensions, move to device\n",
    "    model = CNN().to(device)  # Move model to device\n",
    "    criterion = nn.MSELoss()\n",
    "    optimizer = optim.Adam(model.parameters(), lr=learning_rate)\n",
    "    for epoch in range(num_epochs):\n",
    "        optimizer.zero_grad()\n",
    "        output_tensor = model(input_tensor)\n",
    "        loss = criterion(output_tensor, input_tensor)\n",
    "        loss.backward()\n",
    "        optimizer.step()\n",
    "        \n",
    "        if (epoch+1) % 100== 0:\n",
    "            print(f'Epoch [{epoch+1}/{num_epochs}], Loss: {loss.item():.4f}')\n",
    "    \n",
    "    with torch.no_grad():\n",
    "        output_tensor = model(input_tensor)\n",
    "    reconstructed_array = output_tensor.squeeze(0).squeeze(0).cpu().numpy()\n",
    "    \n",
    "    param_size = 0\n",
    "    for param in model.parameters():\n",
    "        param_size += param.nelement() * param.element_size()\n",
    "    buffer_size = 0\n",
    "    for buffer in model.buffers():\n",
    "        buffer_size += buffer.nelement() * buffer.element_size()\n",
    "    size = (param_size + buffer_size)   \n",
    "    return reconstructed_array, size\n",
    "\n",
    "\n",
    "vdf_temp=vdf.copy()\n",
    "vdf_temp[vdf_temp<1e-16]=1e-16\n",
    "vdf_temp = np.log10(vdf_temp)\n",
    "input_array=vdf_temp\n",
    "recon,total_size= train_and_reconstruct(input_array,100)\n",
    "recon = 10 ** recon\n",
    "recon[recon <= 1e-16] = 0\n",
    "vdf_size=nx*ny*nz*8\n",
    "print(f\"Compresion achieved using a CNN = {round(vdf_size/total_size,2)}\")\n",
    "project_tools.plot_vdfs(vdf,recon)\n",
    "project_tools.print_comparison_stats(vdf,recon)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "805d6d35-924d-4f20-b337-6efc0af0746a",
   "metadata": {},
   "source": [
    "### Here we still use a CNN but this time we use minibatch training and batch normalization layers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a3dcc16f-ddbd-4e2d-afcf-2a4172ce37e8",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Function: train_and_reconstruct\n",
    "\n",
    "Description:\n",
    "This function takes an input array and trains a Convolutional Neural Network (CNN) model to reconstruct the input array.\n",
    "It uses Mean Squared Error (MSE) loss and the Adam optimizer for training.\n",
    "\n",
    "Inputs:\n",
    "- input_array (numpy array): The input array to be reconstructed.\n",
    "- num_epochs (int, optional): The number of training epochs.\n",
    "- learning_rate (float, optional): The learning rate for the Adam optimizer\n",
    "Outputs:\n",
    "    Reconstructed vdf array \n",
    "    Size of model in bytes\n",
    "\"\"\"\n",
    "import torch\n",
    "import torch.nn as nn\n",
    "import torch.optim as optim\n",
    "import numpy as np\n",
    "\n",
    "class CNN(nn.Module):\n",
    "    def __init__(self):\n",
    "        super(CNN, self).__init__()\n",
    "        self.conv1 = nn.Conv3d(1, 16, kernel_size=3, padding=1)\n",
    "        self.bn1 = nn.BatchNorm3d(16)\n",
    "        self.conv2 = nn.Conv3d(16, 32, kernel_size=3, padding=1)\n",
    "        self.bn2 = nn.BatchNorm3d(32)\n",
    "        self.conv3 = nn.Conv3d(32, 64, kernel_size=3, padding=1)\n",
    "        self.bn3 = nn.BatchNorm3d(64)\n",
    "        self.conv4 = nn.Conv3d(64, 1, kernel_size=3, padding=1)\n",
    "        self.relu = nn.ReLU()\n",
    "        \n",
    "    def forward(self, x):\n",
    "        x = self.relu(self.bn1(self.conv1(x)))\n",
    "        x = self.relu(self.bn2(self.conv2(x)))\n",
    "        x = self.relu(self.bn3(self.conv3(x)))\n",
    "        x = self.conv4(x)\n",
    "        return x\n",
    "\n",
    "def train_and_reconstruct(input_array, num_epochs=30, learning_rate=0.001, batch_size=32):\n",
    "    device = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\n",
    "    input_tensor = torch.tensor(input_array, dtype=torch.float32).unsqueeze(0).unsqueeze(0).to(device)  # Move input tensor to device\n",
    "    model = CNN().to(device) \n",
    "    criterion = nn.MSELoss()\n",
    "    optimizer = optim.Adam(model.parameters(), lr=learning_rate)\n",
    "    \n",
    "    for epoch in range(num_epochs):\n",
    "        for i in range(0, input_tensor.size(0), batch_size):\n",
    "            optimizer.zero_grad()\n",
    "            batch_input = input_tensor[i:i+batch_size]\n",
    "            output_tensor = model(batch_input)\n",
    "            loss = criterion(output_tensor, batch_input)\n",
    "            loss.backward()\n",
    "            optimizer.step()\n",
    "    \n",
    "    with torch.no_grad():\n",
    "        output_tensor = model(input_tensor)\n",
    "    reconstructed_array = output_tensor.squeeze(0).squeeze(0).cpu().numpy()\n",
    "    \n",
    "    param_size = 0\n",
    "    for param in model.parameters():\n",
    "        param_size += param.nelement() * param.element_size()\n",
    "    buffer_size = 0\n",
    "    for buffer in model.buffers():\n",
    "        buffer_size += buffer.nelement() * buffer.element_size()\n",
    "    size = (param_size + buffer_size)   \n",
    "    return reconstructed_array, size\n",
    "\n",
    "vdf_temp = vdf.copy()\n",
    "vdf_temp[vdf_temp < 1e-16] = 1e-16\n",
    "vdf_temp = np.log10(vdf_temp)\n",
    "input_array = vdf_temp\n",
    "recon, total_size = train_and_reconstruct(input_array, 100)\n",
    "\n",
    "recon = 10 ** recon\n",
    "recon[recon <= 1e-16] = 0\n",
    "vdf_size = nx * ny * nz * 8\n",
    "print(f\"Compression achieved using a CNN = {round(vdf_size / total_size, 2)}\")\n",
    "\n",
    "project_tools.plot_vdfs(vdf, recon)\n",
    "project_tools.print_comparison_stats(vdf, recon)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a69fd4a-b5d2-4561-b6f7-b0475729714b",
   "metadata": {},
   "source": [
    "### Now we use Hermite Decomposition to perform the compresion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e681c831-3cd0-425f-8ddb-ad5a0057fc32",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Loads the original 3D VDF and fits it to a Maxwellian distribution.\n",
    "  Input: vdf - numpy array representing the original 3D VDF\n",
    "  Output: vdf_herm_3d Reconstructed VDF using Hermite Decomposition\n",
    "\"\"\"\n",
    "### load original 3d vdf and fit Maxwellian\n",
    "vdf_3d=vdf.copy()\n",
    "print('loading done')\n",
    "vdf_size=nx*ny*nz*8\n",
    "\n",
    "#### Fit Maxwellian\n",
    "v_min,v_max,n_bins=0,nx,nx ### define limits and size of velocity axes\n",
    "\n",
    "amp,ux,uy,uz,vthx,vthy,vthz=1e-14,nx,nx,nx,10,10,10 ### initial guess for scipy curve fit\n",
    "guess=amp,ux,uy,uz,vthx,vthy,vthz ### initial guess for scipy curve fit\n",
    "\n",
    "max_fit_3d,params=project_tools.max_fit(vdf_3d,v_min,v_max,n_bins,guess) ### fitting\n",
    "print('Maxwell fit done')\n",
    "\n",
    "\n",
    "#### forward transform ####\n",
    "mm=15 ### PUT THE NUMBER OF HARMONICS\n",
    "norm_amp,u,vth=params[0],params[1:4],params[4:7] ### get the maxwellin fit parameters of thermal and bulk velocity\n",
    "\n",
    "vdf_3d_norm=vdf_3d/norm_amp ### normalize data\n",
    "vdf_3d_flat= vdf_3d_norm.flatten() ### flatten data\n",
    "\n",
    "v_xyz=project_tools.get_flat_mesh(v_min,v_max,n_bins) ### flattening the mesh nodes coordinates\n",
    "herm_array=np.array(project_tools.herm_mpl_arr(m_pol=mm,v_ax=v_xyz,u=params[1:4],vth=params[4:7])) ### create array of hermite polynomials\n",
    "\n",
    "hermite_matrix=project_tools.coefficient_matrix(vdf_3d_flat,mm,herm_array,v_xyz) ### calculate the coefficients of the Hermite transform\n",
    "print('Forward transform done')\n",
    "total_size =5*8+8*np.prod(np.shape(hermite_matrix))\n",
    "\n",
    "#### inverse transform ####\n",
    "inv_herm_flat=project_tools.inv_herm_trans(hermite_matrix, herm_array, v_xyz) ### inverse Hermite transform\n",
    "vdf_herm_3d = (np.reshape(inv_herm_flat,(n_bins,n_bins,n_bins)))*norm_amp ### reshaping back to 3d array and renormalization\n",
    "print('Inverse transform done')\n",
    "print(f\"Compresion achieved using Hermite = {round(vdf_size/total_size,2)}\")\n",
    "project_tools.plot_vdfs(vdf,vdf_herm_3d)\n",
    "project_tools.print_comparison_stats(vdf,vdf_herm_3d)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "66f07c18-105f-4f3e-8489-c90bf18cc962",
   "metadata": {},
   "source": [
    "### Now we use a Gausian Mixture Morel to perform the compresion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2728768-2967-4669-ba52-b201c7284027",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    " Loads the original 3D VDF and performs Gaussian Mixture Model (GMM) decomposition.\n",
    " Input: vdf - NumPy array representing the original 3D VDF\n",
    " Output: vdf_rec Reconstructed VDF using GMM\n",
    "\"\"\"\n",
    "#### load original 3d vdf\n",
    "vdf_3d=vdf.copy()\n",
    "\n",
    "### define number of populations and normalization parameter\n",
    "n_pop=15\n",
    "norm_range=300\n",
    "\n",
    "### RUN GMM\n",
    "means,weights,covs,norm_unit=project_tools.run_gmm(vdf_3d,n_pop,norm_range)\n",
    "### reconstruction resolution and limits of v_space axes\n",
    "n_bins=nx\n",
    "v_min,v_max=0,nx\n",
    "\n",
    "### reconstruction of the vdf \n",
    "vdf_rec=project_tools.reconstruct_vdf(n_pop,means,covs,weights,n_bins,v_min,v_max)\n",
    "vdf_rec=vdf_rec*norm_unit*norm_range\n",
    "total_size =5*8+8*np.prod(np.shape(np.array(covs)))+8*np.prod(np.shape(np.array(weights)))+8*np.prod(np.shape(np.array(means)))\n",
    "print(f\"Compresion achieved using GMM = {round(vdf_size/total_size,2)}\")\n",
    "\n",
    "project_tools.plot_vdfs(vdf,vdf_rec)\n",
    "project_tools.print_comparison_stats(vdf,vdf_rec)"
   ]
  }
 ],
 "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.11.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}