cDMDf.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import all packages\n",
    " \n",
    "import numpy as np\n",
    "import time\n",
    "from sklearn.preprocessing import MinMaxScaler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_data():\n",
    "    # load data\n",
    "    data = np.loadtxt('DMDdata200.csv', delimiter=',')\n",
    "    data = data[:, 1:]\n",
    "    avg = np.mean(data)\n",
    "    data = data - avg\n",
    "    scaler = MinMaxScaler()\n",
    "    X_one_column = data.reshape([-1,1])\n",
    "    result_one_column = scaler.fit_transform(X_one_column)\n",
    "    result = result_one_column.reshape(data.shape)\n",
    "    x1 = result[:,:-1]\n",
    "    x2 = result[:,1:]\n",
    "    \n",
    "    return x1, x2, data\t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "x1, x2, data = load_data()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_rows = 10000\n",
    "\n",
    "rand_int = []\n",
    "for i in range(p_rows):\n",
    "    rand_int.append(np.random.randint(0, x1.shape[0])) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "y1 = []\n",
    "for i in range(len(rand_int)):\n",
    "    y1.append(x1[i, :])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "y2 = []\n",
    "for i in range(len(rand_int)):\n",
    "    y2.append(x2[i, :])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "y1 = np.array(y1)\n",
    "y2 = np.array(y2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_shape = (90824, p_rows)\n",
    "percent_Trunc = 0.99\n",
    "dt = 2.5e-05"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compression matrix generation\n",
    "\n",
    "C_gaussian = np.random.normal(0, 1, size=data_shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compress data\n",
    "\n",
    "y1 = np.dot(C_gaussian.T, x1)\n",
    "y2 = np.dot(C_gaussian.T, x2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time taken for cDMD:  9.832293272018433\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\sumuk\\AppData\\Local\\Temp\\ipykernel_14328\\3659512206.py:27: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  bc = np.linalg.lstsq(phic, x1[:, 0])\n"
     ]
    }
   ],
   "source": [
    "#cDMD\n",
    "\n",
    "startc = 0\n",
    "endc = 0\n",
    "\n",
    "startc = time.time()\n",
    "\n",
    "[Uc, Sc, Vc] = np.linalg.svd(y1, full_matrices=False, compute_uv=True, hermitian=False)\n",
    "\n",
    "menergyc = 0.0\n",
    "rc = 0\n",
    "\n",
    "for i in range(0, np.size(Sc)):\n",
    "    menergyc += Sc[i]/np.sum(Sc)\n",
    "    if menergyc > percent_Trunc:\n",
    "        break\n",
    "    rc = i\n",
    "    \n",
    "menergyc, rc\n",
    "\n",
    "Atildec = Uc[:, :rc].T @ y2 @ Vc[:rc, :].T * np.reciprocal(Sc[: rc])\n",
    "[eivValc,eivVecc] = np.linalg.eig(Atildec)\n",
    "\n",
    "phic = x2 @ Vc[:rc, :].T @ np.diag(np.reciprocal(Sc[: rc])) @ eivVecc\n",
    "phic = phic / np.linalg.norm(phic, axis=0)\n",
    "\n",
    "bc = np.linalg.lstsq(phic, x1[:, 0])\n",
    "S = np.diag(Sc[: rc])\n",
    "num_modesc = phic.shape[1]\n",
    "omegac = np.log(eivValc) / dt\n",
    "\n",
    "endc = time.time()\n",
    "\n",
    "print(\"Time taken for cDMD: \", endc - startc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time taken for DMD:  105.60917139053345\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\sumuk\\AppData\\Local\\Temp\\ipykernel_14328\\1636692992.py:26: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  b = np.linalg.lstsq(phi, x1[:, 0])\n"
     ]
    }
   ],
   "source": [
    "# DMD\n",
    "\n",
    "start = 0\n",
    "end = 0\n",
    "\n",
    "start = time.time()\n",
    "\n",
    "[U, S, V] = np.linalg.svd(x1, full_matrices=False, compute_uv=True, hermitian=False)\n",
    "\n",
    "menergy = 0.0\n",
    "r = 0\n",
    "\n",
    "for i in range(0, np.size(S)):\n",
    "    menergy += S[i]/np.sum(S)\n",
    "    if menergy > percent_Trunc:\n",
    "        break\n",
    "    r = i\n",
    "    \n",
    "menergy, r\n",
    "\n",
    "Atilde = U[:, :r].T @ x2 @ V[:, :r] * np.reciprocal(S[: r])\n",
    "[eivVal,eivVec] = np.linalg.eig(Atilde)\n",
    "\n",
    "phi = x2 @ V[:r, :].T @ np.diag(np.reciprocal(S[: r])) @ eivVec\n",
    "phi = phi / np.linalg.norm(phi, axis=0)\n",
    "b = np.linalg.lstsq(phi, x1[:, 0])\n",
    "S = np.diag(S[: r])\n",
    "num_modes = phi.shape[1]\n",
    "omega = np.log(eivVal) / dt\n",
    "\n",
    "end = time.time()\n",
    "\n",
    "print(\"Time taken for DMD: \", end - start)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Reconstruction of the data\n",
    "\n",
    "def recon(phi, b, omega, num_steps, num_modes, dt):\n",
    "\n",
    "    time = np.arange(num_steps) * dt\n",
    "    time_dynamics = np.zeros((90824, num_steps), dtype=np.complex128)\n",
    "    for i in range(num_steps):\n",
    "            for j in range(num_modes):\n",
    "                time_dynamics[:, i] = b[j] * phi[:, j] * np.exp(omega[j].imag * time[i])\n",
    "\n",
    "    return time_dynamics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\sumuk\\AppData\\Local\\Temp\\ipykernel_14328\\152581265.py:9: RuntimeWarning: overflow encountered in exp\n",
      "  time_dynamics[:, i] = b[j] * phi[:, j] * np.exp(omega[j].imag * time[i])\n"
     ]
    }
   ],
   "source": [
    "reconstructed_fieldc = recon(phic, bc[0], omegac, 2000, num_modesc, dt)\n",
    "reconstructed_field = recon(phi, b[0], omega, 2000, num_modes, dt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.5069136383017739, 0.49039389883163015)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.metrics import mean_squared_error\n",
    "from math import sqrt\n",
    "\n",
    "rmsc = sqrt(mean_squared_error(reconstructed_fieldc[:, :1999].real, x2[:, :]))\n",
    "rms = sqrt(mean_squared_error(x2[:, :], reconstructed_field[:, :1999].real))\n",
    "rmsc, rms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# from matplotlib import pyplot as plt\n",
    "\n",
    "# x = [5000, 10000, 20000, 30000, 40000, 50000]\n",
    "# t_p = [9.74, 11.73, 14.93, 16.83, 20.81, 24.2]\n",
    "# t_g = [5.7, 7.5, 12.8, 19.5, 26.03, 31.3]\n",
    "\n",
    "# err_g = [0.44, 0.44, 0.57, 0.5, 0.57, 0.57]\n",
    "\n",
    "# err_p = [0.58, 0.57, 0.56, 0.51, 0.49, 0.49]\n",
    "\n",
    "# plt.title(\"Time taken for cDMD - Acoustics\")\n",
    "# plt.xlabel(\"Number of rows\")\n",
    "# plt.ylabel(\"Time (s)\")\n",
    "\n",
    "# plt.plot(x, t_p, label = \"p random rows\")\n",
    "# plt.plot(x, t_g, label = \"Gaussian\")\n",
    "# plt.grid(color = 'green', linestyle = '--', linewidth = 0.5)\n",
    "# plt.legend()\n",
    "# plt.savefig(\"time.png\", dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plt.title(\"Error for cDMD - Acoustics\")\n",
    "# plt.xlabel(\"Number of rows\")\n",
    "# plt.ylabel(\"Error\")\n",
    "\n",
    "# plt.plot(x, err_p, label = \"p random rows\")\n",
    "# plt.plot(x, err_g, label = \"Gaussian\")\n",
    "# plt.grid(color = 'green', linestyle = '--', linewidth = 0.5)\n",
    "# plt.legend()\n",
    "\n",
    "# plt.savefig(\"error.png\", dpi=300)"
   ]
  }
 ],
 "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.11.3"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}