data_processing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6d8143c9",
   "metadata": {},
   "source": [
    "The \"Experiment.csv\" file contains simulation data output with the following columns:timestep, atom id, atom type, vx (velocity along the x-axis), vy (velocity along the y-axis), and vz (velocity along the z-axis), x, y, z. The data is organized in ascending order based on the time-step as the primary sorting condition, and within each time-step, it is further sorted by atom id as the secondary sorting condition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "af20055e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "script_dir = os.path.dirname('__file__') #<-- absolute dir the script is in\n",
    "rel_path = \"output/simulation_hydro1/hydro_test001/vtk/experiment.csv\"\n",
    "abs_file_path = os.path.join(script_dir, rel_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "34eb8ea7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "8b5a1c55",
   "metadata": {},
   "outputs": [],
   "source": [
    "array=np.loadtxt(open(abs_file_path, \"rb\"), delimiter=\",\", skiprows=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "8f37de21",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(258964, 9)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "array.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e5c64631",
   "metadata": {},
   "source": [
    "We extra the velocities data from the file. Arrange the data in the following way:\n",
    "$$V=\\left[\\begin{array}{c}VX\\\\\n",
    "VY\\\\\n",
    "VZ\\end{array}\\right]=\\left[\\begin{array}{cccc} v_x\\left(x_{1}, t_{0}\\right) & v_x\\left(x_{1}, t_{1}\\right) & \\cdots & v_x\\left(x_{1}, t_{100}\\right) \\\\ v_x\\left(x_{2}, t_{0}\\right) & v_x\\left(x_{2}, t_{1}\\right) & \\cdots & v_x\\left(x_{2}, t_{100}\\right) \\\\ \\vdots & \\vdots & \\ddots & \\vdots \\\\ v_x\\left(x_{2563}, t_{0}\\right) & u\\left(x_{2563}, t_{1}\\right) & \\cdots & v_x\\left(x_{2563}, t_{100}\\right) \\\\ v_y\\left(x_{1}, t_{0}\\right) & v_y\\left(x_{1}, t_{1}\\right) & \\cdots & v_y\\left(x_{1}, t_{100}\\right) \\\\ v_y\\left(x_{2}, t_{0}\\right) & v_y\\left(x_{2}, t_{1}\\right) & \\cdots & v_y\\left(x_{2}, t_{100}\\right) \\\\ \\vdots & \\vdots & \\ddots & \\vdots \\\\ v_y\\left(x_{2563}, t_{0}\\right) & v_y\\left(x_{2563}, t_{1}\\right) & \\cdots & v_y\\left(x_{2563}, t_{100}\\right)\\\\ v_z\\left(x_{1}, t_{0}\\right) & v_z\\left(x_{1}, t_{1}\\right) & \\cdots & v_z\\left(x_{1}, t_{100}\\right) \\\\ v_z\\left(x_{2}, t_{0}\\right) & v_z\\left(x_{2}, t_{1}\\right) & \\cdots & v_z\\left(x_{2}, t_{100}\\right) \\\\ \\vdots & \\vdots & \\ddots & \\vdots \\\\ v_z\\left(x_{2563}, t_{0}\\right) & v_z\\left(x_{2563}, t_{1}\\right) & \\cdots & v_z\\left(x_{2563}, t_{100}\\right)\\end{array}\\right]$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "d9b7d195",
   "metadata": {},
   "outputs": [],
   "source": [
    "VX=np.reshape(array[0*2564:(0+1)*2564-1,3],(2563,1))\n",
    "for i in range(1,101):\n",
    "    VX=np.append(VX,np.reshape(array[i*2564:(i+1)*2564-1,3],(2563,1)),axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "95c3b358",
   "metadata": {},
   "outputs": [],
   "source": [
    "VY=np.reshape(array[0*2564:(0+1)*2564-1,4],(2563,1))\n",
    "for i in range(1,101):\n",
    "    VY=np.append(VY,np.reshape(array[i*2564:(i+1)*2564-1,4],(2563,1)),axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "fdc9fbf3",
   "metadata": {},
   "outputs": [],
   "source": [
    "VZ=np.reshape(array[0*2564:(0+1)*2564-1,5],(2563,1))\n",
    "for i in range(1,101):\n",
    "    VZ=np.append(VZ,np.reshape(array[i*2564:(i+1)*2564-1,5],(2563,1)),axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "01ac882c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2563, 101)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "VZ.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "73cb1af8",
   "metadata": {},
   "outputs": [],
   "source": [
    "V=np.concatenate((VX,VY),axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "346e0a41",
   "metadata": {},
   "outputs": [],
   "source": [
    "V=np.concatenate((V,VZ),axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "d4b9008d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(7689, 101)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "V.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44aa7607",
   "metadata": {},
   "source": [
    "Utilizing Singular Value Decomposition (SVD) on matrix V and visualizing the singular values through plotting to observe their trends."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "848a8ffa",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(7689, 101)\n",
      "(101,)\n",
      "(101, 101)\n"
     ]
    }
   ],
   "source": [
    "U, S, Vh = np.linalg.svd(V,full_matrices=False)\n",
    "print(np.shape(U))\n",
    "print(np.shape(S))\n",
    "print(np.shape(Vh))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "8ea440f1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'singlar value')"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEMCAYAAADal/HVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAmTklEQVR4nO3deXhV9bn28e+TORAI8xjmMAoiEhmcoI6gIrYqjtW2Ktq39rW1nh61dvDUvu2pdtDWqqhYe7TOHmelalW0IgKKzPMYpjAmISEhw/P+sUOMlIQMe2ft7H1/rovLrLWz135cGu781m8yd0dERKQ2CUEXICIi0U1BISIidVJQiIhInRQUIiJSJwWFiIjUSUEhIiJ1UlCIiEidFBQiIlKnqA8KM7vAzB42s5fN7Kyg6xERiTeBBIWZzTSzPDNbctj5SWa20szWmNmtAO7+krtfB3wLuCSAckVE4lpQLYq/ApNqnjCzROB+YDIwDLjMzIbV+JY7ql4XEZFmlBTEh7r7bDPre9jpMcAad18HYGZPA1PNbDnwG+BNd/+stmua2XRgOkDr1q1HDxkyJCK1i4jEogULFuxy985Hei2QoKhFT2BzjeNcYCzwfeAMINPMst39wSO92d1nADMAcnJyfP78+REuV0QkdpjZxtpei6agsCOcc3e/D7ivuYsREZGQaBr1lAv0qnGcBWxtyAXMbIqZzcjPzw9rYSIi8SyagmIeMNDM+plZCnAp8EpDLuDur7r79MzMzIgUKCISj4IaHvsUMAcYbGa5ZnaNu5cDNwKzgOXAs+6+NIj6RETkS0GNerqslvNvAG809rpmNgWYkp2d3dhLiIjIYaLp0VOT6dGTiEj4xVRQiIhI+MVUUGjUk4hI+MVUUOjRk4hI+MVUUIiISPjFVFDo0ZOISPjFVFDo0ZOISPjFVFCIiEj4KShERKROCgoREalTTAWFOrNFRMIvpoJCndkiIuEXU0EhIiLhp6AQEZE6KShERKROMRUU6swWEQm/mAoKdWaLiIRfTAWFiIiEn4JCRETqpKAQEZE6KShERKROCgoREalTTAWFhseKiIRfTAWFhseKiIRfTAXFIbv3lwZdgohIzIjJoCgoKQ+6BBGRmBGTQXHgYAXuHnQZIiIxISaDosKd9buKgi5DRCQmxGRQACzcvC/oEkREYkJMBkWCmYJCRCRMYjIo0pMTFRQiImESU0FxaMJdopezfFsBJWUVQZckItLixVRQHJpw175NK8oqnGXbCoIuSUSkxUsKuoBISE9JYkx2JzRCVkSk6WIyKJITjSeuHRt0GSIiMSGmHj0d7sBB9VGIiDRVzAbFM/M2MfwXs7Tuk4hIE8VsUByb1Y7vfS0bdVOIiDRNTPZRAAzt3pah3dsGXYaISIsXsy0KgMKSMpZu1SZGIiJNEdNB8bt/rOIbf/mYolItOy4i0lgxHRSThnejtLyS91bmBV2KiEiLFdNBcULfDnTKSOHNJduDLkVEpMWK6aBITDDOPqYb763I05wKEZFGivqgMLP+ZvaomT3fmPdPHt6d4oMVfLBqZ7hLExGJC4EEhZnNNLM8M1ty2PlJZrbSzNaY2a0A7r7O3a9p7GeN7d+B9q2SeWvJtqaWLSISl4JqUfwVmFTzhJklAvcDk4FhwGVmNqypH5ScmMBZw7rxzvI8Ssv1+ElEpKECCQp3nw3sOez0GGBNVQviIPA0MDUcnzdpRDf2l5bz0epd4biciEhciaY+ip7A5hrHuUBPM+toZg8Co8zsttrebGbTzWy+mc3fufOr/REnDehE27Qk7XonItII0bSEhx3hnLv7buCGo73Z3WcAMwBycnK+ssRTSlIC7//H1+jQOiUshYqIxJNoalHkAr1qHGcBWxtygUNboebn//uyHYdC4sm5G1m1o7AJZYqIxJdoCop5wEAz62dmKcClwCsNucChrVAzMzOP+Pr+0nL++M5qHvvXhiYXKyISLwJ59GRmTwETgU5mlgv83N0fNbMbgVlAIjDT3ZeG83MzUpN466ZTSE6KpnwUEYlugQSFu19Wy/k3gDcae10zmwJMyc7OrvV7OmakArBg4x5+9vJSHvv2CXRpk9bYjxQRiXkx9av10R491dSxdSrLtxXoMZSIyFHEVFA0RN9OrTlnRHeemLORgpKyoMsREYlaMRUUdY16OpIbJgygsLScJz/ZFOHKRERarpgKioY8egIY3jOTUwd15tGP1lNSpuU9RESOJKaCojFumNCfXftLeX5BbtCliIhEpbgPivH9OzKyVzse/GAt+7VlqojIv4mpoGhoH0XVe7h98hC27jvAfz6/CHc/+ptEROJITAVFQ/soDhnbvyM/njSEVTsK2VesEVAiIjVF06KAgbr+1P5cPb4v6SmJQZciIhJVYqpF0RRmRnpKIsUHy/n5y0vYsu9A0CWJiESFmAqKxvRRHG7L3gO8+NkW5m84fF8lEZH4ZLHYeZuTk+Pz589v9Pt37y+tXhMqr7BEa0GJSMwzswXunnOk12KqRREuh0JicW4+E377Pn94exUHDmpCnojEJwVFHbK7ZHDWMV25993VTLj7PZ74ZCNlFZVBlyUi0qwUFHVIT0nk3ktH8dwN4+ndoRV3vLSEs/4wm3nqvxCROBJTQRGOzuwjOaFvB567YTyPXp1DRaUz7aE5/ObNFZSW63GUiMQ+dWY30P7Scu56bRlPz9vMkG5t+Ns1Y+jSJo3ig+WkJCaQlBhT2SsicUKd2WGUkZrEby48lkeuymFP0UESzAB4ePZ6zvzDbPYVHwy4QhGR8DrqzGwzGwQ8AHR19+FmdixwvrvfFfHqotgZw7py2pAuVOUEI3tlUulOZnpysIWJiIRZfVoUDwO3AWUA7r4IuDSSRbUUCQmGVSXFxMFd+OGZgzAzNuwqolC75olIjKhPULRy908PO6f1uGtx4GAFl874hGsfn0/xQd0mEWn56hMUu8xsAOAAZnYRsC2iVbVg6SmJ3HbOED7dsIdpD81he35J0CWJiDRJfYLie8BDwBAz2wL8APhuJItqrEgNj22oqcf15JGrcli/s4ip93/Eki3B1iMi0hT1Hh5rZq2BBHcvjGxJTRfJ4bENsXxbAdc+Pp/dRaX81/nDuWh0FgkJFnRZIiL/pq7hsUcNCjP72ZHOu/t/haG2iIiWoADYWVjK/3lyAfM27GVo97bcfs4QThnYOeiyRES+oqnzKIpq/KkAJgN9w1ZdjOvcJpVnpo/n3kuPo7CkjG89No/cvcVBlyUiUm9HnUfh7r+reWxm9wCvRKyiGJSQYEw9rieThndj3vq9ZLVvBcCv31jOuAEd+drgLgFXKCJSu8bMzG4F9A93IfEgNSmRkwd2AkJLgby1dDtLqzq684vL+O1bK/h8094gSxQR+Tf1mZm9mKqhsUAi0BmI2v6JliIjNYl3b55AeWXo1q7cUciM2ev4y/trmXRMN348aTD9O2cEXKWISP06s/vUOCwHdrh7VM8ki6bO7IbYX1rOYx+t58EP1lJaXsnFOb2YOLgzI7Pa0S1Tu+yJSOQ0atSTmXWo66LuHrWbMrTUoDhkZ2Ep9727mqfnbaKsIvTfp0ubVM47tgc/PW9o9bIhIiLhUldQ1PXoaQGhR05H+lvJUT9FxHRuk8ovLxjOT84dytKtBSzK3cf8jXvZU1RaHRJb9x2gR7v0gCsVkXgQU/tRmNkUYEp2dvZ1q1evDrqcsHN3zIxVOwqZfO+H/OGS4zh/ZI+gyxKRGNDk/SjMrL2ZjTGzUw/9CW+J4eHur7r79MzMzKBLiYhDrYmubdK4YUJ/JlRN3Fu1o5B5G/awvzSqu45EpIWqz6ina4GbgCxgITAOmAOcFtHKpFaZrZL5j7OHAKFWxs3PLmTJlgLMoG/H1ozq1Y5xAzpy4oCO1XM2REQa66hBQSgkTgA+cfevmdkQ4M7IliX1ZWbMvPoEFm/JZ9nWApZszeeDVTt58fMtAGS1T2dkr3ZMGNSZaTm9Aq5WRFqi+gRFibuXmBlmluruK8xscMQrk3rr0jaN09umcfrQrkColbFqx37mrN3Fpxv2sHDTPhLNmJbTC3fnV68vZ/KIbozuU+fANhERoH5BkWtm7YCXgLfNbC+wNZJFSdOYGYO7tWFwtzZ866R+AJRVVAKQV1jKcwty6de5NaP7dGBf8UGWbi1gXP+OJGplWxE5ggaNejKzCUAm8Ja7H4xYVU3U0udRRFpZRSUVlU5aciJPfLKRO15aQte2qXzj+Cym5fSiX6fWQZcoIs2ssfMoDr35XuAZd//Y3T8Ie3XS7JITE0hODH190egs2rdK4cXPcpkxex0PvL+WMX07cPWJfZk8vJv2zxCRei3hcTVwCTAI+F9CoRHVv66rRdE4OwpKeOGzXJ6dt5kNu4sZ2CWDG0/L5rxje+ixlEiMa9LGRTUu0gG4ELgU6O3uA8NXYngpKJqmotJ5Y/E2/vTP1azasZ9Lcnrx3xcdG3RZIhJBTXr0VEM2MITQpkXLwlCXRKnEBGPKyB6cO6I7s5Zup2NGKgC5e4u5e9ZKbj5zEH06qh9DJF4cdWa2mf23ma0mtLT4EmC0u0+JeGUSuIQEY/KI7ozpFxpGu3rHfj5YtbP6MdSGXUXsKYraMQ0iEib1aVGsB8a7+65IFyPR7WtDuvDp7WeQkhT6/eIXry7lw9W7OKFve84a1o2zjumqmeAiMSjqFwU0s9bAX4CDwPvu/uTR3qM+iuaxfFsBry/axj+WbWfVjv2YwVnDuvLdidkc16td0OWJSAOEpTM7nMxsJnAekOfuw2ucnwTcS2gnvUfc/Tdm9k1gn7u/ambPuPslR7u+gqL5rd9VxAsLcvnbnA0UlJQzvn9HLhydxYRBnencJjXo8kTkKJq8emwE/BWYVPOEmSUC9wOTgWHAZWY2jNBihJurvq2iGWuUBujXqTW3nD2Yj287nZ+cM5QNu4u45bkveG9FHgCbdhfz2L/Ws2t/acCVikhD1dlHYWYJwKKav/WHg7vPNrO+h50eA6xx93VVn/00MBXI5cuVa2sNNjObDkwH6N27dzjLlQbISE3iulP7c83J/Vi5o5DuVVu4Lti0hztfXcaEQZ3plJHKR6t3sWpHIVNG9lCLQyTK1RkU7l5pZl+YWW933xThWnryZcsBQgExFrgP+LOZnQu8Wtub3X0GMANCj54iWKfUQ0KCMbR72+rjC47ryUnZnejYOhQK7yzfwV8/3sBdry9jbL+OnD60C2cM7UpfLR8iEnXqM+qpO7DUzD4Fig6ddPfzw1zLEbdcdfci4Nth/ixpZmZGlzZp1ce/OP8Yrhjbm5cXbuUfy7Zz1+vLuev15QzsksEFo3pywaie9NRWryJRoT5B0Vx7T+QCNTdMyKKBq9TW2Ao1nHVJhAzs2oZbzh7MLWcPZtPuYv65YgdvLN7O3bNWcveslVw9vg93Th2Ou/Pc/FxG9mrH4G5tgi5bJO4ENjy2qo/itUP9H2aWBKwCTge2APOAy919aUOvrVFPLdvmPcW89PkW+nZqzZSRPcg/UMbIO//BrZOHcMOEAUGXJxKTmjQ81szGAX8ChgIphIauFrl72zrfWPc1nwImAp2AHcDP3f1RMzsH+GPVZ8x091818LqHWhTXrV69urHlSZSpqHS27D1Al7appCUn8vLCLby/ciff+9oAsruohSESDk0NivmEFgJ8DsgBrgIGuvvt4S40XNSiiG2PfLiOe/6xktLySk3wEwmTJi8K6O5rzCzR3SuAx8zs47BWKNIA157Sn6+P6snjczby+McbmLX0X3TKSCWrfTpZ7dMZ1qMtl4/pTbtWKUGXKhIT6hMUxWaWAiw0s98C2wCNYZRAdcxI5eYzB3H9qf154bNclm0tIHfvARZvyee1RdtYm1fE76aNDLpMkZhQn6D4JqE+gxuBHxIamXRhJItqLI16ij+tU5O4anzfr5xbsb2A9Kot/JZvK+Dvczdx0xkD6ZShiX0ijXHUJTzcfaO7H3D3Ane/091vdvc1zVFcQ7n7q+4+PTMzM+hSJEBDurWt3i9jcW4+ry7aSnJC6H/1Rbn7WLWjkNJyrQYjUl+1tijMbDFQa0+3u2vLM4l6007oxZSRPUhPCbUwfvryUr7YvI8Eg+wuGdx42kCmHNsdM231KlKbWkc9mVmfut7o7hsjUlETaHisHM2yrQWszitk7c4i3l62g+XbChjZqx0/OWdo9QZNIvEo6pYZjzQNj5X6qKx0Xvx8C/fMWsn2ghJG92nPRaOzOO/Y7rRJSw66PJFm1aRlxs2s0MwKDvuz2cz+18z6h79ckeaRkGBcNDqL926ZyE/OGUr+gTJue3Ex63eFljRbu3M/y7YWEIu/TIk0RH1GPf2e0JpLfye0cN+lQDdgJTCT0AxrkRYrPSWR607tz7Wn9GPZtgKGVa16+8iH63lt0VYW/uwsEg2enLuR/ANlnDuie3VnuUg8qM/M7LnuPvawc5+4+zgz+8Ldo2awuvooJJy27jvA+l1FnJTdCYCrZn7K7FU7adcqmce/PYaRmg0uMaSpO9xVmtk0M0uo+jOtxmtR1SbX8FgJpx7t0qtDAuBv3xnD+7dMpE1aEpc//Alz1u4OsDqR5lOfoLiC0KS7PEIL+H0TuNLM0glNwhOJG307tea560+ke7t0vvXYp8xaul19GBLzNOpJpBH2FB3k6pmfsnhLPlnt0zlzWFcuPD6L4T3VmpWWqUmLAppZZ+A6oG/N73f374SrQJGWpkPrFJ65fhwvL9zK28t28OTcTWS1b8XwnpnsKCjh3ndXc9kJvRmRpeCQlq8+o55eBj4E3gG07oFIlVYpSVw2pjeXjelNUWl5dYfdhl1FvPrFVq4cG5qzumXfAVKTErTWlLRY9Rn1tNDdj2uecppGo54kWlRWOgkJoWVBfvD057y5ZDv/cfZgvn1SPxITtFyIRJ+mjnp6rWrnuainUU8SLRJqhMGNp2VzysBO3PX6ci55aE71hD6RlqI+LYpCQvtPlAJlhCbdeVO2Qo00dWZLtHF3XvxsC794dSllFZV8+6R+XDmuDz3bpQddmgigtZ5Eosb2/BJ++doy3lyyDYAzhnbl4pxeTBzcmeTE+jTwRSKjUaOezGyIu68ws+OP9Lq7fxauAkXiRbfMNO6/4nhy9xbz5NxNPP3pJj5YtZOld54NwN/nbqKsopKrT+wbbKEiNdQ16ulmYDrwuyO85sBpEalIJA5ktW/Ff04awk2nD2RN3n6SqloTH6/dxa79pdVB8bc5Gzg2qx0jszK1Z4YERo+eRKLMgYMVpKcksr+0nLG/eoeigxX079Sa80b24LhemQzt3pZubdMUHBJWTeqjMLOLgbfcvdDM7gCOB37p7p+Hv9Sm0fBYiTUFJWW8tXg7L3yWy9z1e6rPt2+VzJSRPbj25P707tgqwAolVjQ1KBa5+7FmdjLwa+Ae4PbDV5SNJmpRSCwqLCljxfZClm8rYP6Gvby5ZBsVlc7j3xnDKQM7B12etHBNWsKDL2djnws84O4vm9kvwlWciNRPm7RkTujbgRP6duCq8X3Znj+UJ+du5IS+oS1cX/liKyVlFVw8OkuPpSSs6jMeb4uZPQRMA94ws9R6vk9EIqhbZho/OmswacmJALyycCvPzttcHRKx2P8owajPX/jTgFnAJHffB3QA/iOSRYlIwz181Wgevir05CCvsIQz/zCbFxbkUlGpwJCmOeqjJ3cvBl6scbwN2BbJokSk4cyM9q1TAMgvLqNVSiI/eu4LfjtrBWcf041Jx3RjTL8O1UNxRepLw2NFYlRlpTNr6XZeXriV91flUVJWSXKi0SYtmYzUJNqmJzHzWyfQpU0aB8srSUlSgMSzpnZmi0gLlJBgTB7RnckjunPgYAUfrMrji9x89peUU1hSRklZJZ1ah5Y+/9nLS1i3s4i/XzdWLQ75NwoKkTiQnpLIpOHdmTS8+xFfH92nPT3bpVeHRElZRXUnuUhMBUWNCXdBlyLSolyc06v66zlrd3PDEwsY268Do/u0Z3Sf9hzfu/1Xlk6X+BJTbUztRyHSdO1aJXPG0K6sztvPr99cwUUPzuEHzyykUqOn4lZMtShEpOmGdm/L76aNBGDX/lL+9vEG7vvnGnq2T+c/Jw0JuDoJgoJCRGrVKSOVH545iN1FB3ng/bVktU/niqq9wCV+KChEpE5mxp3nH8O2/BJ++tISurVN4/ShXYMuS5pRTPVRiEhkJCUm8KfLRnFMj0w27ykGYNPuYqY9OIf5G/Yc5d3S0qlFISL10jo1iWevH8+BstA6oUUHyyktr6BNWjIA/1yxg7V5RVx3av8gy5QIUFCISL2lpySSnhKaXzG0e1tevvHk6teWbS3gnn+soqCkjJvPHKQVbGOIgkJEwuK7E7PJ3XuAP/1zDQfLK7l18hCFRYxQUIhIWCQmGP/v6yNITkzgodnr2FdcxtUn9mVo9zYKjBZOQSEiYZOQYPzX1GNolZLIQ7PX8cz8zXRuk8op2Z04eWAnvnF8FgDlFZVaU6oF0X8pEQkrM+O2c4Yy9/bTufuiYxnbrwMfrNrJu8vzgNA+4Gf9cTYvfb4l4EqlvtSiEJGI6No2jYtzelWvI3VoCZCy8koGdM6gT8dWQZYnDRD1QWFm/YGfAJnuflHQ9YhI4xxaVLBjRmr1TnwAd89aQZ+OrbXXdxSL6KMnM5tpZnlmtuSw85PMbKWZrTGzW+u6hruvc/drIlmniASjtLyCeev38uPnF3HFI3PZsKso6JLkCCLdR/FXYFLNE2aWCNwPTAaGAZeZ2TAzG2Fmrx32p0uE6xORAKUmJfL09HH86uvDWZybz9l/nM3ds1awe39p0KVJDRF99OTus82s72GnxwBr3H0dgJk9DUx1918D5zX2s8xsOjAdoHfv3o29jIg0s4QE44qxfTh9SFd+9cZy/vL+WmZ+tIHLx/Zm+qn96do2LegS414Qo556AptrHOdWnTsiM+toZg8Co8zsttq+z91nuHuOu+d07tw5fNWKSLPolpnGny4bxds/PJXJI7rx1483cPJ//5M7XlocdGlxL4jO7CP1VtW6I4q77wZuiFw5IhJNsru04ffTjuOHZwzikQ/XfWW+xeY9xfTqoNFSzS2IFkUu0KvGcRawNRwXNrMpZjYjPz8/HJcTkQD16tCKO6cO56fnDQNCW7ROuPs93l2+I+DK4k8QQTEPGGhm/cwsBbgUeCUcF9ZWqCKxa1iPtnz/tIGclN0JgHkb9lQveS6RFdFHT2b2FDAR6GRmucDP3f1RM7sRmAUkAjPdfWkk6xCRli8zPZkfnjkIAHfn1hcWsXZnEWP7deCi0VlMPa4nKUlabCISzD12Nkw3synAlOzs7OtWr14ddDkiEkG5e4t56fMtvPDZFtbvKqJHZhrfnTiAi3N6kZacGHR5LY6ZLXD3nCO+FktBcUhOTo7Pnz8/6DJEpBm4O7NX7+K+d1ezYONeurZN5abTB3HJCb1ITNBM7/qqKyjUThORFs3MmDCoM8/fMJ6/XzuWXu1bcfv/Lubc+z5kR0FJ0OXFhKhf66khajx6CroUEWlmZsaJ2Z0YP6Ajby7ZzisLt9I5IzXosmJCTLUoNOpJRMyMc0Z058FvjiYhwcgrLGHag3NYub0w6NJarJgKChGRw+UVlJJXWIJXzevdX1pOLPbNRpIePYlITBveM5N3bp5QPcP76/f/i017iumWmcbXBnfhxtOy6aRHVHXSqCcRiStPzt3Ixt3FbNxdxDvL80hLSuC6U/tz7Sn9yUiNqd+dG6SuUU/xe1dEJC5dMbZP9ddrd+7nnlkr+eM7q3lm3maeum4cfTu1DrC66KQ+ChGJWwM6Z/DAlaN54bvjKSmr4IpH5pK7V8uCHE5BISJxb3SfDvzPNWMpLCnj92+vCrqcqBNTj57UmS0ijTW8ZybP3jCe3lXLmFdUumZ2V4mpFoXmUYhIUwzp1pZWKUkcLK/k/D9/xN/nbgq6pKgQU0EhIhIOxQfL6d2hFd0zQ9uw7i06SGFJWcBVBSemHj2JiIRDu1YpPHDl6Orjhz9cx8MfrmNc/46cMbQrk0d0o0ub+NnLW0EhInIUk4Z3o7zSeWf5Dn7+ylJ+8epSxvXryJSRPThvZHfapiUHXWJExdSEO+1HISKRtiavkFe/2Mari7aybmcRGalJXDamF98+qR892qUHXV6jaT8KEZEwc3cWb8nn0Y/W89qibVw8OovfXHhs0GU1moJCRCSCcvcWk2BGj3bp5BeX0SYtiYQWNrRWGxeJiERQVvtW9GiXTmFJGV9/4F/c9fryoEsKK3Vmi4iESUZqElNH9uTE7I5BlxJWCgoRkTAxM246Y2D18X3vrianT3tOzO4UYFVNp0dPIiIRsL+0nGfmbebyR+Zy+cOfsHxbQdAlNVpMBYWZTTGzGfn5+UGXIiJxLiM1iXd/NIGfnTeMldsLufjBOSzYuCfosholpoJCaz2JSDRJS07kOyf347X/ezKd26Ry1aOfMnfd7qDLarCYCgoRkWjUPTOdZ6aPo1tmGt96bB7/WrMr6JIaREEhItIMurRN4+npoWXMr3x0LpfOmEN+cctYaFBBISLSTDq3SeWZ68dx0+kDSU9OpG16aOBpXkFJwJXVTcNjRUSaUbtWKfzgjEHVx8u3FXDhAx/zy6nDuXB0VoCV1U4tChGRAA3onMHlY3pz6qDOQZdSKwWFiEiAUpISuOO8YXRuk0pZRSUvL9xCeUVl0GV9hYJCRCRKPDNvMzc9vZBvPPAxy7ZGzwS9mAoKTbgTkZbsirG9+fPlo9i67wDn//kjfvvWCsqioHURU0GhCXci0pKZGecd24N3bp7ABaN68pf313L9/yzgwMGKQOuKqaAQEYkF7VqlcM/FI7nrguG8tzKPKx75hL1FBwOrR0EhIhKlrhzXh79cfjxLthRw8UNzWJNXGEgdCgoRkSg2eUR3Hv/OGHbkl/DEJ5uA0DasJWXN9zhKQSEiEuXGD+jIu7dM4LsTBwAwd/0eJt/7YbOFhWZmi4i0AF3apFV/XVhSzoDOrSk4UEZacmLEP1tBISLSwpw5rCtnDuvabJ+nR08iIi3Uhl1FzbIZkoJCRKSFuvGpz/jx84uorPSIfo6CQkSkhbr+1AGs3VnEW0u3R/RzFBQiIi3UOSO6079Ta+5/bw3ukWtVqDNbRKSFSkwwbpg4gB8/v4jTfvcBg7pmMKhrG84a1o0RWeFbyijqg8LMLgDOBboA97v7P4KtSEQkenxjVE/yi8tYsHEvq/IKeWd5Hn95fy0/PXco3zqpX1g+I6JBYWYzgfOAPHcfXuP8JOBeIBF4xN1/U9s13P0l4CUzaw/cAygoRESqJCUmcN2p/bmu6jj/QBl3vLSEXh1ahe8zwnalI/sr8Gfgb4dOmFkicD9wJpALzDOzVwiFxq8Pe/933D2v6us7qt4nIiK1yExP5k+Xjao+fnLuRs4Y2pWubdPqeFfdIhoU7j7bzPoednoMsMbd1wGY2dPAVHf/NaHWx1eYmQG/Ad50988iWa+ISCw5cLCCX7+xgo27i7n9nKGNvk4QfRQ9gc01jnOBsXV8//eBM4BMM8t29weP9E1mNh2YXnVYamZLwlFsDOgE7Aq6iCii+/FVuh9fitl78ZOqP0fRp7YXgggKO8K5Wsd1uft9wH1Hu6i7zwBmAJjZfHfPaXSFMUT34qt0P75K9+NLuhe1C2IeRS7Qq8ZxFrA1gDpERKQeggiKecBAM+tnZinApcArAdQhIiL1ENGgMLOngDnAYDPLNbNr3L0cuBGYBSwHnnX3pWH+6Blhvl5LpnvxVbofX6X78SXdi1pYJKd9i4hIy6e1nkREpE4KChERqVNMBYWZTTKzlWa2xsxuDbqe5mZmvczsPTNbbmZLzeymqvMdzOxtM1td9c/2QdfaXMws0cw+N7PXqo7j+V60M7PnzWxF1f8j4+P1fpjZD6t+RpaY2VNmlhav96I+YiYoaiwNMhkYBlxmZsOCrarZlQM/cvehwDjge1X34FbgXXcfCLxbdRwvbiI0aOKQeL4X9wJvufsQYCSh+xJ398PMegL/F8ipWoMukdDoy7i7F/UVM0FBjaVB3P0g8DQwNeCampW7bzu0zIm7FxL6i6AnofvweNW3PQ5cEEiBzczMsgitPPxIjdPxei/aAqcCjwK4+0F330ec3g9Ck43TzSwJaEVoLle83oujiqWgONLSID0DqiVwVWtsjQLmAl3dfRuEwoTQku3x4I/Aj4HKGufi9V70B3YCj1U9invEzFoTh/fD3bcQWol6E7ANyK/aviDu7kV9xVJQNGhpkFhmZhnAC8AP3L0g6HqCYGaHlrdfEHQtUSIJOB54wN1HAUXE6aOVqr6HqUA/oAfQ2syuDLaq6BZLQaGlQQAzSyYUEk+6+4tVp3eYWfeq17sDebW9P4acBJxvZhsIPYY8zcyeID7vBYR+PnLdfW7V8fOEgiMe78cZwHp33+nuZcCLwInE572ol1gKirhfGqRqSfZHgeXu/vsaL70CXF319dXAy81dW3Nz99vcPcvd+xL6f+Gf7n4lcXgvANx9O7DZzAZXnTodWEZ83o9NwDgza1X1M3M6of68eLwX9RJTM7PN7BxCz6UTgZnu/qtgK2peZnYy8CGwmC+fy99OqJ/iWaA3oR+Si919TyBFBsDMJgK3uPt5ZtaROL0XZnYcoY79FGAd8G1CvyzG3f0wszuBSwiNFPwcuBbIIA7vRX3EVFCIiEj4xdKjJxERiQAFhYiI1ElBISIidVJQiIhInRQUIiJSJwWFSCOZ2f4Gfv/EQ6vYirQkCgoREamTgkKkiapaCu/X2OvhyaoZv4f2SFlhZh8B36jxntZmNtPM5lUt0je16vx9Zvazqq/PNrPZZqafUwlUUtAFiMSIUcAxhNYX+xdwkpnNBx4GTgPWAM/U+P6fEFpW5Dtm1g741MzeIbRQ3zwz+xC4DzjH3WuufivS7PSbikh4fOruuVV/qS8E+gJDCC0+t9pDSyA8UeP7zwJuNbOFwPtAGtDb3YuB64C3gT+7+9pm+zcQqYVaFCLhUVrj6wq+/NmqbY0cAy5095VHeG0EsJvQEtgigVOLQiRyVgD9zGxA1fFlNV6bBXy/Rl/GqKp/9gF+ROhR1mQzG9uM9YockYJCJELcvQSYDrxe1Zm9scbLvwSSgUVmtgT4ZY1l4m9x963ANcAjZpbWzKWLfIVWjxURkTqpRSEiInVSUIiISJ0UFCIiUicFhYiI1ElBISIidVJQiIhInRQUIiJSp/8PrERFvc4dEEoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "plt.semilogy(S,'-.')\n",
    "plt.xlim(0,99)\n",
    "plt.ylim(10**(-2),10**2)\n",
    "plt.xlabel('Index',fontsize=10)\n",
    "plt.ylabel('singlar value',fontsize=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c5845d7",
   "metadata": {},
   "source": [
    "The first 40 singular values account for a significant proportion of the total sum of all singular values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "68dd02e3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9229389181761992"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(S[0:40])/sum(S)"
   ]
  }
 ],
 "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}