Merge pull request #19 from jakekrell/bernoulli_asymptotes

FoKL v3.2.3

docs/_dev/basis_functions/bernoulli_polynomials/BSS-ANOVA_eigendecomposition.ipynb

@@ -0,0 +1,127 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# BSS-ANOVA eigendecomposition at increasing resolution (to later find asymptotic eigenvalue ratios)\n",
+    "\n",
+    "Importing packages."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "from scipy.interpolate import CubicSpline\n",
+    "from scipy.optimize import curve_fit"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Defining the BSS-ANOVA kernel main effect $\\kappa_1$, outlined in equation 8 of [Fast variable selection makes Karhunen-Loève decomposed Gaussian process BSS-ANOVA a speedy and accurate choice for dynamic systems identification](docs/_static/arXiv.2205.13676v2.pdf),\n",
+    "\n",
+    "$\\kappa_1(x,x') = \\mathcal{B}(x)\\mathcal{B}_1(x') + \\mathcal{B}_2(x)\\mathcal{B}_2(x') + \\frac{1}{24}\\mathcal{B}_4(|x-x'|)$\n",
+    "\n",
+    "where\n",
+    "\n",
+    "$\\begin{cases} \\mathcal{B}_1(x) = x - \\frac{1}{2} \\\\ \\mathcal{B}_2(x) = x^2 - x + \\frac{1}{6} \\\\ \\mathcal{B}_4(x) = x^4 - 2x^3 + x^2 - \\frac{1}{30} \\end{cases}$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def b1(x):\n",
+    "    return x - 1/2\n",
+    "\n",
+    "def b2(x):\n",
+    "    return x**2 - x + 1/6\n",
+    "\n",
+    "def b4(x):\n",
+    "    return x**4 - 2*x**3 + x**2 - 1/30\n",
+    "\n",
+    "def k1(xi, xj):\n",
+    "    return b1(xi)*b1(xj) + b2(xi)*b2(xj) - b4(np.abs(xi-xj))/24"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Taking eigenvalues for increasing resolution of covariance matrix (i.e., BSS-ANOVA kernel). Because only 20 Bernoulli polynomials could be computed in MATLAB prior to significant rounding error in plots, only need first 20 eigenvalues."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 22,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/tmp/ipykernel_28833/3318184751.py:18: ComplexWarning: Casting complex values to real discards the imaginary part\n",
+      "  eigvals[res_iter, :] = eigval[:n]  # in future, plot columns which are basis function scales\n"
+     ]
+    }
+   ],
+   "source": [
+    "n = 20  # number of Bernoulli polynomials (i.e., number of eigenvalues to save)\n",
+    "res_n = 5  # number of points to plot\n",
+    "res_lb = 1600  # lower bound (of plot)\n",
+    "res_ub = 2000  # upper bound (of plot)\n",
+    "\n",
+    "eigvals = np.zeros([res_n, n])\n",
+    "res_x = np.linspace(res_lb, res_lb + np.round((res_ub-res_lb)/(res_n-1))*(res_n-1), res_n, dtype=int)\n",
+    "res_iter = 0\n",
+    "for res in res_x:\n",
+    "    x = np.linspace(0, 1, res)\n",
+    "    kernel = np.zeros([res, res])\n",
+    "\n",
+    "    for i in range(res):\n",
+    "        for j in range(res):\n",
+    "            kernel[i, j] = k1(x[i], x[j])\n",
+    "    eigval, eigvec = np.linalg.eig(kernel)\n",
+    "\n",
+    "    eigvals[res_iter, :] = eigval[:n]  # in future, plot columns which are basis function scales\n",
+    "    res_iter += 1\n",
+    "\n",
+    "progress = np.concatenate([res_x[:, np.newaxis], eigvals], axis=1)\n",
+    "np.savetxt(f'current_progress_{res_lb}_{res_ub}.txt', progress)  # res points by basis function order (i.e., 'k' or eigenvalue id)\n",
+    "\n",
+    "# !!! NOTE !!!\n",
+    "# Manually combine multiple 'current_progress_{res_lb}_{res_ub}.txt' files into single 'BSS-ANOVA_eigenvalues_for_20x20_thru_2000x2000.txt'"
+   ]
+  }
+ ],
+ "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.10.12"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}