diff --git a/denoising-plots/.ipynb_checkpoints/Untitled-checkpoint.ipynb b/denoising-plots/.ipynb_checkpoints/Untitled-checkpoint.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..2d4751c9dd5597e0bc92c06c2ec65ce0bca0e323
--- /dev/null
+++ b/denoising-plots/.ipynb_checkpoints/Untitled-checkpoint.ipynb
@@ -0,0 +1,139 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f5605c04",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import scipy\n",
+    "import matplotlib.pyplot as plt\n",
+    "import quaternion  # load the quaternion module\n",
+    "import bispy as bsp\n",
+    "import torch\n",
+    "import tqdm\n",
+    "import time\n",
+    "from sklearn.preprocessing import normalize\n",
+    "from utils import STIS,optimize_loop,snr_bivariate,param_search,objective_STIS,objective_KReSP,KReSP\n",
+    "from ray import  tune,init\n",
+    "from tempfile import TemporaryFile\n",
+    "import pickle\n",
+    "\n",
+    "init(num_cpus=13)\n",
+    "\n",
+    "## PLOT AN EXAMPLE \n",
+    "\n",
+    "np.random.seed(5)\n",
+    "N = 1024 # length of the signal\n",
+    "t = np.linspace(0, 2*np.pi/4, N) # time vector\n",
+    "dt = t[1]-t[0]\n",
+    "\n",
+    "# ellipse parameters - AM-FM-PM polarized \n",
+    "theta1 = np.pi/4 - 2*t\n",
+    "chi1 = np.pi/16 - t\n",
+    "phi1 = 0 \n",
+    "f0 = 25/N/dt \n",
+    "S0 = bsp.utils.windows.hanning(N)\n",
+    "\n",
+    "x_quad = bsp.signals.bivariateAMFM(S0, theta1, chi1, 2*np.pi*f0*t+ phi1)\n",
+    "x = quaternion.as_float_array(x_quad)[:,:2]\n",
+    "\n",
+    "bsp.utils.visual.plot2D(t,x_quad)\n",
+    "plt.savefig(\"clean_sig.pdf\")\n",
+    "\n",
+    "sigma = 0.05\n",
+    "n = np.zeros([N,4])\n",
+    "noise_complex = np.random.randn(N,2)\n",
+    "y = x + noise_complex\n",
+    "\n",
+    "uH = np.imag(scipy.signal.hilbert(noise_complex[:,0]))\n",
+    "vH = np.imag(scipy.signal.hilbert(noise_complex[:,1]))\n",
+    "n[:,0] = noise_complex[:,0]\n",
+    "n[:,1] = noise_complex[:,1]\n",
+    "n[:,2] = uH\n",
+    "n[:,3] = vH\n",
+    "n = quaternion.from_float_array(n)\n",
+    "y_quad = sigma*n + x_quad # Noisy signal\n",
+    "\n",
+    "bsp.utils.visual.plot2D(t,y_quad)\n",
+    "plt.savefig(\"noisy_sig.pdf\")\n",
+    "n = np.random.randn(N,2)\n",
+    "y = sigma*n + x # Noisy signal\n",
+    "\n",
+    "print(\"sigma: \" + str(sigma) + \" Noise SNR: \"+  str(snr_bivariate(x,y) ) )\n",
+    "search_space = {\"x\":tune.grid_search([x]),\"t\":tune.grid_search([t]),\"y\":tune.grid_search([y]),\"lambdax\": tune.grid_search((0.1)**np.linspace(5,15,7)), \"lambdaS\":  tune.grid_search((0.1)**np.linspace(5,10,7)) , \"beta1\":tune.grid_search((0.10)**np.linspace(-2,2,5)),\"beta2\":tune.grid_search((0.10)**np.linspace(-2,2,5)),\"sigma2\":tune.grid_search([sigma**2])}\n",
+    "config = param_search(objective_STIS,search_space)\n",
+    "model = STIS(t,y,lambdax=config[\"lambdax\"],lambdaS=config[\"lambdaS\"],beta1=config[\"beta1\"],beta2=config[\"beta2\"],sigma2=sigma**2,p=2)\n",
+    "optimizer = torch.optim.Adam(model.parameters(), lr=0.01)\n",
+    "model = optimize_loop(model,optimizer,numit=1000)\n",
+    "\n",
+    "x_stis = quaternion.from_float_array(model.Xquad.detach().numpy())\n",
+    "bsp.utils.visual.plot2D(t,x_stis)\n",
+    "plt.savefig(\"denoised_via_all_terms.pdf\")\n",
+    "\n",
+    "print(\"sigma: \" + str(sigma) + \" STIS SNR: \"+  str(snr_bivariate(x,model.X.detach().numpy())))\n",
+    "\n",
+    "\n",
+    "# NO STOKES REGULARIZATION  \n",
+    "search_space = {\"x\":tune.grid_search([x]),\"t\":tune.grid_search([t]),\"y\":tune.grid_search([y]),\"lambdax\": tune.grid_search((0.1)**np.linspace(5,15,7)), \"lambdaS\":tune.grid_search([0.0]) , \"beta1\":tune.grid_search([0.0]),\"beta2\":tune.grid_search([0.0]),\"sigma2\":tune.grid_search([sigma**2])}\n",
+    "config = param_search(objective_STIS,search_space)\n",
+    "model = STIS(t,y,lambdax=config[\"lambdax\"],lambdaS=config[\"lambdaS\"],beta1=config[\"beta1\"],beta2=config[\"beta2\"],sigma2=sigma**2,p=2)\n",
+    "optimizer = torch.optim.Adam(model.parameters(), lr=0.01)\n",
+    "model = optimize_loop(model,optimizer,numit=1000)\n",
+    "\n",
+    "print(\"sigma: \" + str(sigma) + \" STIS SNR: \"+  str((snr_bivariate(x,model.X.detach().numpy()))))\n",
+    "x_stis = quaternion.from_float_array(model.Xquad.detach().numpy())\n",
+    "bsp.utils.visual.plot2D(t,x_stis)\n",
+    "plt.savefig(\"denoised_via_no_stokes.pdf\")\n",
+    "\n",
+    "# ONLY SMOOTH STOKES\n",
+    "search_space = {\"x\":tune.grid_search([x]),\"t\":tune.grid_search([t]),\"y\":tune.grid_search([y]),\"lambdax\": tune.grid_search([0.0]), \"lambdaS\":  tune.grid_search((0.1)**np.linspace(5,10,7)) , \"beta1\":tune.grid_search((0.1)**np.linspace(-2,2,7)),\"beta2\":tune.grid_search((0.1)**np.linspace(-2,2,5)),\"sigma2\":tune.grid_search([sigma**2])}\n",
+    "config = param_search(objective_STIS,search_space)\n",
+    "model = STIS(t,y,lambdax=config[\"lambdax\"],lambdaS=config[\"lambdaS\"],beta1=config[\"beta1\"],beta2=config[\"beta2\"],sigma2=sigma**2,p=2)\n",
+    "optimizer = torch.optim.Adam(model.parameters(), lr=0.01)\n",
+    "model = optimize_loop(model,optimizer,numit=1000)\n",
+    "x_stis = quaternion.from_float_array(model.Xquad.detach().numpy())\n",
+    "bsp.utils.visual.plot2D(t,x_stis)\n",
+    "plt.savefig(\"denoised_via_no_signal_smoother.pdf\")\n",
+    "\n",
+    "print(\"sigma: \" + str(sigma) +  \" STIS SNR: \"+  str((snr_bivariate(x,model.X.detach().numpy()))))\n",
+    "\n",
+    "# Kernel regression on normalized \n",
+    "search_space = {\"x\":tune.grid_search([x]),\"t\":tune.grid_search([t]),\"y\":tune.grid_search([y]),\"alpha\": tune.grid_search((0.1)**np.linspace(5,15,7)),\"lambda_1\": tune.grid_search((0.1)**np.linspace(5,15,7)), \"lambda_s\":  tune.grid_search((0.1)**np.linspace(5,10,7)) , \"beta\":tune.grid_search((0.10)**np.linspace(-2,2,7)),\"gamma\":tune.grid_search([0.2]),\"sigma2\":tune.grid_search([sigma**2])}\n",
+    "config = param_search(objective_KReSP,search_space)\n",
+    "model = KReSP(t,y,lambda_1=config[\"lambda_1\"],beta=config[\"beta\"],lambda_s=config[\"lambda_s\"],alpha=config[\"alpha\"],gamma=config[\"gamma\"],eps=10**-7,win_width=32,sigma2=sigma**2)\n",
+    "optimizer = torch.optim.Adam(model.parameters(), lr=0.01)\n",
+    "model = optimize_loop(model,optimizer,numit=300)\n",
+    "\n",
+    "\n",
+    "print(\"sigma: \" + str(sigma) + \" KReSP SNR: \"+  str(snr_bivariate(x,model.X.detach().numpy())))\n",
+    "x_kerreg = quaternion.from_float_array(model.Xquad.detach().numpy())\n",
+    "bsp.utils.visual.plot2D(t,x_kerreg)\n",
+    "plt.savefig(\"denoised_via_kerreg.pdf\")"
+   ]
+  }
+ ],
+ "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.3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/denoising-plots/.vscode/settings.json b/denoising-plots/.vscode/settings.json
new file mode 100644
index 0000000000000000000000000000000000000000..9e26dfeeb6e641a33dae4961196235bdb965b21b
--- /dev/null
+++ b/denoising-plots/.vscode/settings.json
@@ -0,0 +1 @@
+{}
\ No newline at end of file
diff --git a/denoising-plots/__pycache__/utils.cpython-310.pyc b/denoising-plots/__pycache__/utils.cpython-310.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..e849215aeb954d30eeab762af367b765ff57cd5a
Binary files /dev/null and b/denoising-plots/__pycache__/utils.cpython-310.pyc differ
diff --git a/denoising-plots/__pycache__/utils.cpython-311.pyc b/denoising-plots/__pycache__/utils.cpython-311.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..114ceef33bc818fe6e5258e6e708523aedb6f567
Binary files /dev/null and b/denoising-plots/__pycache__/utils.cpython-311.pyc differ
diff --git a/denoising-plots/clean_sig.pdf b/denoising-plots/clean_sig.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..c05ffc7df2650b7c9d4547ca4c2342de2961efdf
Binary files /dev/null and b/denoising-plots/clean_sig.pdf differ
diff --git a/denoising-plots/denoised_via_all_terms.pdf b/denoising-plots/denoised_via_all_terms.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..3c4ec18277be4db5d80a782d2ad0ba0f826f28fd
Binary files /dev/null and b/denoising-plots/denoised_via_all_terms.pdf differ
diff --git a/denoising-plots/denoised_via_kerreg.pdf b/denoising-plots/denoised_via_kerreg.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..c129e230aed54d8196152bf497e39825b3c7b3b7
Binary files /dev/null and b/denoising-plots/denoised_via_kerreg.pdf differ
diff --git a/denoising-plots/denoised_via_no_signal_smoother.pdf b/denoising-plots/denoised_via_no_signal_smoother.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..638cbef3ae60660d24815b77182b71118e346dd2
Binary files /dev/null and b/denoising-plots/denoised_via_no_signal_smoother.pdf differ
diff --git a/denoising-plots/denoised_via_no_stokes.pdf b/denoising-plots/denoised_via_no_stokes.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..d392f9ea873628e8ca6f59b8e8a8011482fd8b31
Binary files /dev/null and b/denoising-plots/denoised_via_no_stokes.pdf differ
diff --git a/denoising-plots/denoising-perf.pdf b/denoising-plots/denoising-perf.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..f1b97d0d8805fb8fbafe653447804d2331c1923c
Binary files /dev/null and b/denoising-plots/denoising-perf.pdf differ
diff --git a/denoising-plots/example.pkl b/denoising-plots/example.pkl
new file mode 100644
index 0000000000000000000000000000000000000000..37222505709bb710fd5c3dfe33aa6b2157146ca4
Binary files /dev/null and b/denoising-plots/example.pkl differ
diff --git a/denoising-plots/examples.ipynb b/denoising-plots/examples.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..e9e8b9dd8cc76e39f4ef4255268864b1dab1fe7b
--- /dev/null
+++ b/denoising-plots/examples.ipynb
@@ -0,0 +1,197 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "6db15641",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/html": [],
+      "text/plain": [
+       "<IPython.core.display.HTML object>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "\u001b[36m(pid=3980984)\u001b[0m /usr/lib/python3/dist-packages/scipy/__init__.py:146: UserWarning: A NumPy version >=1.17.3 and <1.25.0 is required for this version of SciPy (detected version 1.26.4\n",
+      "\u001b[36m(pid=3980984)\u001b[0m   warnings.warn(f\"A NumPy version >={np_minversion} and <{np_maxversion}\"\n",
+      "2024-03-14 15:03:03,729\tWARNING tune.py:186 -- Stop signal received (e.g. via SIGINT/Ctrl+C), ending Ray Tune run. This will try to checkpoint the experiment state one last time. Press CTRL+C (or send SIGINT/SIGKILL/SIGTERM) to skip. \n",
+      "2024-03-14 15:03:09,319\tWARNING tune.py:1057 -- Experiment has been interrupted, but the most recent state was saved.\n",
+      "Resume experiment with: Tuner.restore(path=\"/home/ypilavci/ray_results/objective_STIS_2024-03-14_14-56-04\", trainable=...)\n",
+      "2024-03-14 15:03:15,927\tWARNING experiment_analysis.py:193 -- Failed to fetch metrics for 17 trial(s):\n",
+      "- objective_STIS_99690_00970: FileNotFoundError('Could not fetch metrics for objective_STIS_99690_00970: both result.json and progress.csv were not found at /home/ypilavci/ray_results/objective_STIS_2024-03-14_14-56-04/objective_STIS_99690_00970_970_beta1=100.0000,beta2=0.0100,lambdaS=0.0000,lambdax=0.0000,sigma2=0.0400,t=ref_ph_8d1a09f1,x=ref_ph__2024-03-14_15-02-55')\n",
+      "- objective_STIS_99690_00971: FileNotFoundError('Could not fetch metrics for objective_STIS_99690_00971: both result.json and progress.csv were not found at /home/ypilavci/ray_results/objective_STIS_2024-03-14_14-56-04/objective_STIS_99690_00971_971_beta1=10.0000,beta2=0.0100,lambdaS=0.0000,lambdax=0.0000,sigma2=0.0400,t=ref_ph_8d1a09f1,x=ref_ph_a_2024-03-14_15-02-55')\n",
+      "- objective_STIS_99690_00972: FileNotFoundError('Could not fetch metrics for objective_STIS_99690_00972: both result.json and progress.csv were not found at /home/ypilavci/ray_results/objective_STIS_2024-03-14_14-56-04/objective_STIS_99690_00972_972_beta1=1.0000,beta2=0.0100,lambdaS=0.0000,lambdax=0.0000,sigma2=0.0400,t=ref_ph_8d1a09f1,x=ref_ph_a5_2024-03-14_15-02-55')\n",
+      "- objective_STIS_99690_00974: FileNotFoundError('Could not fetch metrics for objective_STIS_99690_00974: both result.json and progress.csv were not found at /home/ypilavci/ray_results/objective_STIS_2024-03-14_14-56-04/objective_STIS_99690_00974_974_beta1=0.0100,beta2=0.0100,lambdaS=0.0000,lambdax=0.0000,sigma2=0.0400,t=ref_ph_8d1a09f1,x=ref_ph_a5_2024-03-14_15-02-55')\n",
+      "- objective_STIS_99690_00976: FileNotFoundError('Could not fetch metrics for objective_STIS_99690_00976: both result.json and progress.csv were not found at /home/ypilavci/ray_results/objective_STIS_2024-03-14_14-56-04/objective_STIS_99690_00976_976_beta1=10.0000,beta2=100.0000,lambdaS=0.0000,lambdax=0.0000,sigma2=0.0400,t=ref_ph_8d1a09f1,x=ref_ph_2024-03-14_15-03-00')\n",
+      "- objective_STIS_99690_00977: FileNotFoundError('Could not fetch metrics for objective_STIS_99690_00977: both result.json and progress.csv were not found at /home/ypilavci/ray_results/objective_STIS_2024-03-14_14-56-04/objective_STIS_99690_00977_977_beta1=1.0000,beta2=100.0000,lambdaS=0.0000,lambdax=0.0000,sigma2=0.0400,t=ref_ph_8d1a09f1,x=ref_ph__2024-03-14_15-03-00')\n",
+      "- objective_STIS_99690_00978: FileNotFoundError('Could not fetch metrics for objective_STIS_99690_00978: both result.json and progress.csv were not found at /home/ypilavci/ray_results/objective_STIS_2024-03-14_14-56-04/objective_STIS_99690_00978_978_beta1=0.1000,beta2=100.0000,lambdaS=0.0000,lambdax=0.0000,sigma2=0.0400,t=ref_ph_8d1a09f1,x=ref_ph__2024-03-14_15-03-00')\n",
+      "- objective_STIS_99690_00980: FileNotFoundError('Could not fetch metrics for objective_STIS_99690_00980: both result.json and progress.csv were not found at /home/ypilavci/ray_results/objective_STIS_2024-03-14_14-56-04/objective_STIS_99690_00980_980_beta1=100.0000,beta2=10.0000,lambdaS=0.0000,lambdax=0.0000,sigma2=0.0400,t=ref_ph_8d1a09f1,x=ref_ph_2024-03-14_15-03-01')\n",
+      "- objective_STIS_99690_00982: FileNotFoundError('Could not fetch metrics for objective_STIS_99690_00982: both result.json and progress.csv were not found at /home/ypilavci/ray_results/objective_STIS_2024-03-14_14-56-04/objective_STIS_99690_00982_982_beta1=1.0000,beta2=10.0000,lambdaS=0.0000,lambdax=0.0000,sigma2=0.0400,t=ref_ph_8d1a09f1,x=ref_ph_a_2024-03-14_15-03-01')\n",
+      "- objective_STIS_99690_00983: FileNotFoundError('Could not fetch metrics for objective_STIS_99690_00983: both result.json and progress.csv were not found at /home/ypilavci/ray_results/objective_STIS_2024-03-14_14-56-04/objective_STIS_99690_00983_983_beta1=0.1000,beta2=10.0000,lambdaS=0.0000,lambdax=0.0000,sigma2=0.0400,t=ref_ph_8d1a09f1,x=ref_ph_a_2024-03-14_15-03-01')\n",
+      "- objective_STIS_99690_00984: FileNotFoundError('Could not fetch metrics for objective_STIS_99690_00984: both result.json and progress.csv were not found at /home/ypilavci/ray_results/objective_STIS_2024-03-14_14-56-04/objective_STIS_99690_00984_984_beta1=0.0100,beta2=10.0000,lambdaS=0.0000,lambdax=0.0000,sigma2=0.0400,t=ref_ph_8d1a09f1,x=ref_ph_a_2024-03-14_15-03-01')\n",
+      "- objective_STIS_99690_00986: FileNotFoundError('Could not fetch metrics for objective_STIS_99690_00986: both result.json and progress.csv were not found at /home/ypilavci/ray_results/objective_STIS_2024-03-14_14-56-04/objective_STIS_99690_00986_986_beta1=10.0000,beta2=1.0000,lambdaS=0.0000,lambdax=0.0000,sigma2=0.0400,t=ref_ph_8d1a09f1,x=ref_ph_a_2024-03-14_15-03-01')\n",
+      "- objective_STIS_99690_00987: FileNotFoundError('Could not fetch metrics for objective_STIS_99690_00987: both result.json and progress.csv were not found at /home/ypilavci/ray_results/objective_STIS_2024-03-14_14-56-04/objective_STIS_99690_00987_987_beta1=1.0000,beta2=1.0000,lambdaS=0.0000,lambdax=0.0000,sigma2=0.0400,t=ref_ph_8d1a09f1,x=ref_ph_a5_2024-03-14_15-03-01')\n",
+      "- objective_STIS_99690_00988: FileNotFoundError('Could not fetch metrics for objective_STIS_99690_00988: both result.json and progress.csv were not found at /home/ypilavci/ray_results/objective_STIS_2024-03-14_14-56-04/objective_STIS_99690_00988_988_beta1=0.1000,beta2=1.0000,lambdaS=0.0000,lambdax=0.0000,sigma2=0.0400,t=ref_ph_8d1a09f1,x=ref_ph_a5_2024-03-14_15-03-01')\n",
+      "- objective_STIS_99690_00989: FileNotFoundError('Could not fetch metrics for objective_STIS_99690_00989: both result.json and progress.csv were not found at /home/ypilavci/ray_results/objective_STIS_2024-03-14_14-56-04/objective_STIS_99690_00989_989_beta1=0.0100,beta2=1.0000,lambdaS=0.0000,lambdax=0.0000,sigma2=0.0400,t=ref_ph_8d1a09f1,x=ref_ph_a5_2024-03-14_15-03-01')\n",
+      "- objective_STIS_99690_00990: FileNotFoundError('Could not fetch metrics for objective_STIS_99690_00990: both result.json and progress.csv were not found at /home/ypilavci/ray_results/objective_STIS_2024-03-14_14-56-04/objective_STIS_99690_00990_990_beta1=100.0000,beta2=0.1000,lambdaS=0.0000,lambdax=0.0000,sigma2=0.0400,t=ref_ph_8d1a09f1,x=ref_ph__2024-03-14_15-03-01')\n",
+      "- objective_STIS_99690_00991: FileNotFoundError('Could not fetch metrics for objective_STIS_99690_00991: both result.json and progress.csv were not found at /home/ypilavci/ray_results/objective_STIS_2024-03-14_14-56-04/objective_STIS_99690_00991_991_beta1=10.0000,beta2=0.1000,lambdaS=0.0000,lambdax=0.0000,sigma2=0.0400,t=ref_ph_8d1a09f1,x=ref_ph_a_2024-03-14_15-03-03')\n"
+     ]
+    }
+   ],
+   "source": [
+    "import numpy as np\n",
+    "import scipy\n",
+    "import matplotlib.pyplot as plt\n",
+    "import quaternion  # load the quaternion module\n",
+    "import bispy as bsp\n",
+    "import torch\n",
+    "from utils import STIS,optimize_loop,snr_bivariate,param_search,objective_STIS,objective_KReSP,KReSP,plot_on_sphere\n",
+    "from ray import  tune,init\n",
+    "import pickle\n",
+    "\n",
+    "\n",
+    "init(num_cpus=16)\n",
+    "\n",
+    "## PLOT AN EXAMPLE \n",
+    "\n",
+    "np.random.seed(5)\n",
+    "N = 1024 # length of the signal\n",
+    "t = np.linspace(0, 2*np.pi/4, N) # time vector\n",
+    "dt = t[1]-t[0]\n",
+    "\n",
+    "# ellipse parameters - AM-FM-PM polarized \n",
+    "theta1 = np.pi/4 - 2*t\n",
+    "chi1 = np.pi/16 - t\n",
+    "phi1 = 0 \n",
+    "f0 = 25/N/dt \n",
+    "S0 = bsp.utils.windows.hanning(N)\n",
+    "\n",
+    "x_quad = bsp.signals.bivariateAMFM(S0, theta1, chi1, 2*np.pi*f0*t+ phi1)\n",
+    "x = quaternion.as_float_array(x_quad)[:,:2]\n",
+    "\n",
+    "\n",
+    "bsp.utils.visual.plot3D(t,x_quad)\n",
+    "plt.savefig(\"clean_sig.pdf\")\n",
+    "plt.close()\n",
+    "\n",
+    "sigma = 0.2\n",
+    "n = np.zeros([N,4])\n",
+    "noise_complex = np.random.randn(N,2)\n",
+    "y = x +sigma*noise_complex\n",
+    "\n",
+    "uH = np.imag(scipy.signal.hilbert(noise_complex[:,0]))\n",
+    "vH = np.imag(scipy.signal.hilbert(noise_complex[:,1]))\n",
+    "n[:,0] = noise_complex[:,0]\n",
+    "n[:,1] = noise_complex[:,1]\n",
+    "n[:,2] = uH\n",
+    "n[:,3] = vH\n",
+    "n = quaternion.from_float_array(n)\n",
+    "y_quad = sigma*n + x_quad # Noisy signal\n",
+    "\n",
+    "ax = plt.figure().add_subplot(projection='3d')\n",
+    "# ax.view_init(15,-135,0)\n",
+    "plt.title(\"Normalized Stokes Parameters\")\n",
+    "plt.tight_layout()\n",
+    "plot_on_sphere(y_quad,ax,label=\"noisy signal\",t=t,scatter=True)\n",
+    "plot_on_sphere(x_quad,ax,label=\"original signal\",scatter=False)\n",
+    "\n",
+    "plt.legend()\n",
+    "plt.tight_layout()\n",
+    "plt.savefig(\"noisy_sig_sphere.pdf\")\n",
+    "plt.close()\n",
+    "\n",
+    "bsp.utils.visual.plot3D(t,y_quad)\n",
+    "plt.savefig(\"noisy_sig.pdf\")\n",
+    "\n",
+    "\n",
+    "print(\"sigma: \" + str(sigma) + \" Noise SNR: \"+  str(snr_bivariate(x,y) ) )\n",
+    "search_space = {\"x\":tune.grid_search([x]),\"t\":tune.grid_search([t]),\"y\":tune.grid_search([y]),\"lambdax\": tune.grid_search((0.1)**np.linspace(5,15,7)), \"lambdaS\":  tune.grid_search((0.1)**np.linspace(5,10,7)) , \"beta1\":tune.grid_search((0.10)**np.linspace(-2,2,5)),\"beta2\":tune.grid_search((0.10)**np.linspace(-2,2,5)),\"sigma2\":tune.grid_search([sigma**2])}\n",
+    "config = param_search(objective_STIS,search_space)\n",
+    "model = STIS(t,y,lambdax=config[\"lambdax\"],lambdaS=config[\"lambdaS\"],beta1=config[\"beta1\"],beta2=config[\"beta2\"],sigma2=sigma**2,p=2)\n",
+    "optimizer = torch.optim.Adam(model.parameters(), lr=0.01)\n",
+    "model = optimize_loop(model,optimizer,numit=1000)\n",
+    "\n",
+    "x_stis = quaternion.from_float_array(model.Xquad.detach().numpy())\n",
+    "bsp.utils.visual.plot2D(t,x_stis)\n",
+    "plt.savefig(\"denoised_via_all_terms.pdf\")\n",
+    "\n",
+    "print(\"sigma: \" + str(sigma) + \" STIS SNR: \"+  str(snr_bivariate(x,model.X.detach().numpy())))\n",
+    "\n",
+    "\n",
+    "# NO STOKES REGULARIZATION  \n",
+    "search_space = {\"x\":tune.grid_search([x]),\"t\":tune.grid_search([t]),\"y\":tune.grid_search([y]),\"lambdax\": tune.grid_search((0.1)**np.linspace(5,15,7)), \"lambdaS\":tune.grid_search([0.0]) , \"beta1\":tune.grid_search([0.0]),\"beta2\":tune.grid_search([0.0]),\"sigma2\":tune.grid_search([sigma**2])}\n",
+    "config = param_search(objective_STIS,search_space)\n",
+    "model = STIS(t,y,lambdax=config[\"lambdax\"],lambdaS=config[\"lambdaS\"],beta1=config[\"beta1\"],beta2=config[\"beta2\"],sigma2=sigma**2,p=2)\n",
+    "optimizer = torch.optim.Adam(model.parameters(), lr=0.01)\n",
+    "model = optimize_loop(model,optimizer,numit=1000)\n",
+    "\n",
+    "print(\"sigma: \" + str(sigma) + \" STIS SNR: \"+  str((snr_bivariate(x,model.X.detach().numpy()))))\n",
+    "x_nostokes = quaternion.from_float_array(model.Xquad.detach().numpy())\n",
+    "bsp.utils.visual.plot2D(t,x_nostokes)\n",
+    "plt.savefig(\"denoised_via_no_stokes.pdf\")\n",
+    "\n",
+    "# ONLY SMOOTH STOKES\n",
+    "search_space = {\"x\":tune.grid_search([x]),\"t\":tune.grid_search([t]),\"y\":tune.grid_search([y]),\"lambdax\": tune.grid_search([0.0]), \"lambdaS\":  tune.grid_search((0.1)**np.linspace(5,10,7)) , \"beta1\":tune.grid_search((0.1)**np.linspace(-2,2,7)),\"beta2\":tune.grid_search((0.1)**np.linspace(-2,2,5)),\"sigma2\":tune.grid_search([sigma**2])}\n",
+    "config = param_search(objective_STIS,search_space)\n",
+    "model = STIS(t,y,lambdax=config[\"lambdax\"],lambdaS=config[\"lambdaS\"],beta1=config[\"beta1\"],beta2=config[\"beta2\"],sigma2=sigma**2,p=2)\n",
+    "optimizer = torch.optim.Adam(model.parameters(), lr=0.01)\n",
+    "model = optimize_loop(model,optimizer,numit=1000)\n",
+    "x_onlystokes = quaternion.from_float_array(model.Xquad.detach().numpy())\n",
+    "bsp.utils.visual.plot2D(t,x_onlystokes)\n",
+    "plt.savefig(\"denoised_via_no_signal_smoother.pdf\")\n",
+    "\n",
+    "print(\"sigma: \" + str(sigma) +  \" STIS SNR: \"+  str((snr_bivariate(x,model.X.detach().numpy()))))\n",
+    "\n",
+    "# Kernel regression on normalized \n",
+    "search_space = {\"x\":tune.grid_search([x]),\"t\":tune.grid_search([t]),\"y\":tune.grid_search([y]),\"alpha\": tune.grid_search((0.1)**np.linspace(5,15,5)),\"lambda_1\": tune.grid_search((0.1)**np.linspace(5,15,5)), \"lambda_s\":  tune.grid_search((0.1)**np.linspace(5,10,5)) , \"beta\":tune.grid_search((0.10)**np.linspace(-2,2,5)),\"gamma\":tune.grid_search((0.1)**np.linspace(0,1,5)),\"sigma2\":tune.grid_search([sigma**2])}\n",
+    "config = param_search(objective_KReSP,search_space)\n",
+    "model = KReSP(t,y,lambda_1=config[\"lambda_1\"],beta=config[\"beta\"],lambda_s=config[\"lambda_s\"],alpha=config[\"alpha\"],gamma=config[\"gamma\"],eps=10**-7,win_width=64,sigma2=sigma**2)\n",
+    "optimizer = torch.optim.Adam(model.parameters(), lr=0.01)\n",
+    "model = optimize_loop(model,optimizer,numit=300)\n",
+    "\n",
+    "\n",
+    "print(\"sigma: \" + str(sigma) + \" KReSP SNR: \"+  str(snr_bivariate(x,model.X.detach().numpy())))\n",
+    "x_kerreg = quaternion.from_float_array(model.Xquad.detach().numpy())\n",
+    "bsp.utils.visual.plot2D(t,x_kerreg)\n",
+    "plt.savefig(\"denoised_via_kerreg.pdf\")\n",
+    "\n",
+    "with open('example.pkl', 'wb') as f:\n",
+    "\tpickle.dump(x,f)\n",
+    "\tpickle.dump(y,f)\n",
+    "\tpickle.dump(x_stis,f)\n",
+    "\tpickle.dump(x_nostokes,f)\n",
+    "\tpickle.dump(x_onlystokes,f)\n",
+    "\tpickle.dump(x_kerreg,f)\n"
+   ]
+  }
+ ],
+ "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.10.12"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/denoising-plots/examples.py b/denoising-plots/examples.py
new file mode 100644
index 0000000000000000000000000000000000000000..9d5044b401a89b7a54436c00af4d7cf0b6f02606
--- /dev/null
+++ b/denoising-plots/examples.py
@@ -0,0 +1,123 @@
+import numpy as np
+import scipy
+import matplotlib.pyplot as plt
+import quaternion  # load the quaternion module
+import bispy as bsp
+import torch
+from utils import STIS,optimize_loop,snr_bivariate,param_search,objective_STIS,objective_KReSP,KReSP,plot_on_sphere
+from ray import  tune,init
+import pickle
+
+
+init(num_cpus=16)
+
+## PLOT AN EXAMPLE 
+
+np.random.seed(5)
+N = 1024 # length of the signal
+t = np.linspace(0, 2*np.pi/4, N) # time vector
+dt = t[1]-t[0]
+
+# ellipse parameters - AM-FM-PM polarized 
+theta1 = np.pi/4 - 2*t
+chi1 = np.pi/16 - t
+phi1 = 0 
+f0 = 25/N/dt 
+S0 = bsp.utils.windows.hanning(N)
+
+x_quad = bsp.signals.bivariateAMFM(S0, theta1, chi1, 2*np.pi*f0*t+ phi1)
+x = quaternion.as_float_array(x_quad)[:,:2]
+
+
+bsp.utils.visual.plot3D(t,x_quad)
+plt.savefig("clean_sig.pdf")
+plt.close()
+
+sigma = 0.1
+n = np.zeros([N,4])
+noise_complex = np.random.randn(N,2)
+y = x +sigma*noise_complex
+
+uH = np.imag(scipy.signal.hilbert(noise_complex[:,0]))
+vH = np.imag(scipy.signal.hilbert(noise_complex[:,1]))
+n[:,0] = noise_complex[:,0]
+n[:,1] = noise_complex[:,1]
+n[:,2] = uH
+n[:,3] = vH
+n = quaternion.from_float_array(n)
+y_quad = sigma*n + x_quad # Noisy signal
+
+ax = plt.figure().add_subplot(projection='3d')
+# ax.view_init(15,-135,0)
+plt.title("Normalized Stokes Parameters")
+plt.tight_layout()
+plot_on_sphere(y_quad,ax,label="noisy signal",t=t,scatter=True)
+plot_on_sphere(x_quad,ax,label="original signal",scatter=False)
+
+plt.legend()
+plt.tight_layout()
+plt.savefig("noisy_sig_sphere.pdf")
+plt.close()
+
+bsp.utils.visual.plot3D(t,y_quad)
+plt.savefig("noisy_sig.pdf")
+
+
+print("sigma: " + str(sigma) + " Noise SNR: "+  str(snr_bivariate(x,y) ) )
+search_space = {"x":tune.grid_search([x]),"t":tune.grid_search([t]),"y":tune.grid_search([y]),"lambdax": tune.grid_search((0.1)**np.linspace(5,15,7)), "lambdaS":  tune.grid_search((0.1)**np.linspace(5,10,7)) , "beta1":tune.grid_search((0.10)**np.linspace(-2,2,5)),"beta2":tune.grid_search((0.10)**np.linspace(-2,2,5)),"sigma2":tune.grid_search([sigma**2])}
+config = param_search(objective_STIS,search_space)
+model = STIS(t,y,lambdax=config["lambdax"],lambdaS=config["lambdaS"],beta1=config["beta1"],beta2=config["beta2"],sigma2=sigma**2,p=2)
+optimizer = torch.optim.Adam(model.parameters(), lr=0.01)
+model = optimize_loop(model,optimizer,numit=1000)
+
+x_stis = quaternion.from_float_array(model.Xquad.detach().numpy())
+bsp.utils.visual.plot2D(t,x_stis)
+plt.savefig("denoised_via_all_terms.pdf")
+
+print("sigma: " + str(sigma) + " STIS SNR: "+  str(snr_bivariate(x,model.X.detach().numpy())))
+
+
+# NO STOKES REGULARIZATION  
+search_space = {"x":tune.grid_search([x]),"t":tune.grid_search([t]),"y":tune.grid_search([y]),"lambdax": tune.grid_search((0.1)**np.linspace(5,15,7)), "lambdaS":tune.grid_search([0.0]) , "beta1":tune.grid_search([0.0]),"beta2":tune.grid_search([0.0]),"sigma2":tune.grid_search([sigma**2])}
+config = param_search(objective_STIS,search_space)
+model = STIS(t,y,lambdax=config["lambdax"],lambdaS=config["lambdaS"],beta1=config["beta1"],beta2=config["beta2"],sigma2=sigma**2,p=2)
+optimizer = torch.optim.Adam(model.parameters(), lr=0.01)
+model = optimize_loop(model,optimizer,numit=1000)
+
+print("sigma: " + str(sigma) + " STIS SNR: "+  str((snr_bivariate(x,model.X.detach().numpy()))))
+x_nostokes = quaternion.from_float_array(model.Xquad.detach().numpy())
+bsp.utils.visual.plot2D(t,x_nostokes)
+plt.savefig("denoised_via_no_stokes.pdf")
+
+# ONLY SMOOTH STOKES
+search_space = {"x":tune.grid_search([x]),"t":tune.grid_search([t]),"y":tune.grid_search([y]),"lambdax": tune.grid_search([0.0]), "lambdaS":  tune.grid_search((0.1)**np.linspace(5,10,7)) , "beta1":tune.grid_search((0.1)**np.linspace(-2,2,7)),"beta2":tune.grid_search((0.1)**np.linspace(-2,2,5)),"sigma2":tune.grid_search([sigma**2])}
+config = param_search(objective_STIS,search_space)
+model = STIS(t,y,lambdax=config["lambdax"],lambdaS=config["lambdaS"],beta1=config["beta1"],beta2=config["beta2"],sigma2=sigma**2,p=2)
+optimizer = torch.optim.Adam(model.parameters(), lr=0.01)
+model = optimize_loop(model,optimizer,numit=1000)
+x_onlystokes = quaternion.from_float_array(model.Xquad.detach().numpy())
+bsp.utils.visual.plot2D(t,x_onlystokes)
+plt.savefig("denoised_via_no_signal_smoother.pdf")
+
+print("sigma: " + str(sigma) +  " STIS SNR: "+  str((snr_bivariate(x,model.X.detach().numpy()))))
+
+# Kernel regression on normalized 
+search_space = {"x":tune.grid_search([x]),"t":tune.grid_search([t]),"y":tune.grid_search([y]),"alpha": tune.grid_search((0.1)**np.linspace(5,15,7)),"lambda_1": tune.grid_search((0.1)**np.linspace(5,15,7)), "lambda_s":  tune.grid_search((0.1)**np.linspace(5,10,7)) , "beta":tune.grid_search((0.10)**np.linspace(-2,2,7)),"gamma":tune.grid_search([0.2]),"sigma2":tune.grid_search([sigma**2])}
+config = param_search(objective_KReSP,search_space)
+model = KReSP(t,y,lambda_1=config["lambda_1"],beta=config["beta"],lambda_s=config["lambda_s"],alpha=config["alpha"],gamma=config["gamma"],eps=10**-7,win_width=32,sigma2=sigma**2)
+optimizer = torch.optim.Adam(model.parameters(), lr=0.01)
+model = optimize_loop(model,optimizer,numit=300)
+
+
+print("sigma: " + str(sigma) + " KReSP SNR: "+  str(snr_bivariate(x,model.X.detach().numpy())))
+x_kerreg = quaternion.from_float_array(model.Xquad.detach().numpy())
+bsp.utils.visual.plot2D(t,x_kerreg)
+plt.savefig("denoised_via_kerreg.pdf")
+
+with open('example.pkl', 'wb') as f:
+	pickle.dump(x,f)
+	pickle.dump(y,f)
+	pickle.dump(x_stis,f)
+	pickle.dump(x_nostokes,f)
+	pickle.dump(x_onlystokes,f)
+	pickle.dump(x_kerreg,f)
diff --git a/denoising-plots/main.py b/denoising-plots/main.py
new file mode 100644
index 0000000000000000000000000000000000000000..4d44da7378e00098ca309a18d904344a7d5fdb9f
--- /dev/null
+++ b/denoising-plots/main.py
@@ -0,0 +1,112 @@
+import numpy as np
+import scipy
+import matplotlib.pyplot as plt
+import quaternion  # load the quaternion module
+import bispy as bsp
+import torch
+import tqdm
+import time
+from sklearn.preprocessing import normalize
+from utils import STIS,optimize_loop,snr_bivariate,param_search,objective_STIS,objective_KReSP,KReSP
+from ray import  tune,init
+import pickle 
+
+np.random.seed(5)
+init(num_cpus=16)
+N = 1024 # length of the signal
+t = np.linspace(0, 2*np.pi/4, N) # time vector
+dt = t[1]-t[0]
+
+# ellipse parameters - AM-FM-PM polarized 
+theta1 = np.pi/4 - 2*t
+chi1 = np.pi/16 - t
+phi1 = 0 
+f0 = 25/N/dt 
+S0 = bsp.utils.windows.hanning(N)
+
+x = bsp.signals.bivariateAMFM(S0, theta1, chi1, 2*np.pi*f0*t+ phi1)
+x = quaternion.as_float_array(x)[:,:2]
+sigmaRange = np.linspace(0.01,0.5,10)
+NREP =  5
+snr_noise = np.zeros([len(sigmaRange),NREP])
+snr_p2 = np.zeros([len(sigmaRange),NREP])
+snr_nostokes = np.zeros([len(sigmaRange),NREP])
+snr_onlystokes = np.zeros([len(sigmaRange),NREP])
+snr_kerreg = np.zeros([len(sigmaRange),NREP])
+
+for nrep in tqdm.tqdm(range(NREP)):
+	s = 0 
+	for sigma in (sigmaRange):
+		n = np.random.randn(N,2)
+		y = sigma*n + x # Noisy signal
+		snr_noise[s,nrep] =  (snr_bivariate(x,y) )
+		print("sigma: " + str(sigma) + " nrep: " + str(nrep) + " Noise SNR: "+  str(snr_bivariate(x,y) ) )
+
+		search_space = {"x":tune.grid_search([x]),"t":tune.grid_search([t]),"y":tune.grid_search([y]),"lambdax": tune.grid_search((0.1)**np.linspace(5,15,7)), "lambdaS":  tune.grid_search((0.1)**np.linspace(5,10,7)) , "beta1":tune.grid_search((0.10)**np.linspace(-2,2,5)),"beta2":tune.grid_search((0.10)**np.linspace(-2,2,5)),"sigma2":tune.grid_search([sigma**2])}
+
+		config = param_search(objective_STIS,search_space)
+		
+		model = STIS(t,y,lambdax=config["lambdax"],lambdaS=config["lambdaS"],beta1=config["beta1"],beta2=config["beta2"],sigma2=sigma**2,p=2)
+#		model = STIS(t,y,lambdax=1e-8,lambdaS=1e-8,beta1=0.1,beta2=0.1,sigma2=sigma**2,p=2)
+
+		optimizer = torch.optim.Adam(model.parameters(), lr=0.01)
+
+		model = optimize_loop(model,optimizer,numit=1000)
+
+
+		snr_p2[s,nrep] = (snr_bivariate(x,model.X.detach().numpy()))
+
+		print("sigma: " + str(sigma) + " nrep: " + str(nrep) + " STIS SNR: "+  str(snr_p2[s,nrep]))
+
+
+		# NO STOKES REGULARIZATION  
+		search_space = {"x":tune.grid_search([x]),"t":tune.grid_search([t]),"y":tune.grid_search([y]),"lambdax": tune.grid_search((0.1)**np.linspace(5,15,7)), "lambdaS":tune.grid_search([0.0]) , "beta1":tune.grid_search([0.0]),"beta2":tune.grid_search([0.0]),"sigma2":tune.grid_search([sigma**2])}
+
+
+		config = param_search(objective_STIS,search_space)
+		
+		model = STIS(t,y,lambdax=config["lambdax"],lambdaS=config["lambdaS"],beta1=config["beta1"],beta2=config["beta2"],sigma2=sigma**2,p=2)
+		optimizer = torch.optim.Adam(model.parameters(), lr=0.01)
+		model = optimize_loop(model,optimizer,numit=1000)
+
+
+		snr_nostokes[s,nrep] = (snr_bivariate(x,model.X.detach().numpy()))
+
+		print("sigma: " + str(sigma) + " nrep: " + str(nrep) + " STIS SNR: "+  str(snr_nostokes[s,nrep]))
+
+
+		# ONLY SMOOTH STOKES
+		search_space = {"x":tune.grid_search([x]),"t":tune.grid_search([t]),"y":tune.grid_search([y]),"lambdax": tune.grid_search([0.0]), "lambdaS":  tune.grid_search((0.1)**np.linspace(5,10,7)) , "beta1":tune.grid_search((0.1)**np.linspace(-2,2,7)),"beta2":tune.grid_search((0.1)**np.linspace(-2,2,5)),"sigma2":tune.grid_search([sigma**2])}
+
+		config = param_search(objective_STIS,search_space)
+		
+		model = STIS(t,y,lambdax=config["lambdax"],lambdaS=config["lambdaS"],beta1=config["beta1"],beta2=config["beta2"],sigma2=sigma**2,p=2)
+		optimizer = torch.optim.Adam(model.parameters(), lr=0.01)
+		model = optimize_loop(model,optimizer,numit=1000)
+
+
+		snr_onlystokes[s,nrep] = (snr_bivariate(x,model.X.detach().numpy()))
+
+		print("sigma: " + str(sigma) + " nrep: " + str(nrep) + " STIS SNR: "+  str(snr_onlystokes[s,nrep]))
+
+		# Kernel regression on normalized 
+		search_space = {"x":tune.grid_search([x]),"t":tune.grid_search([t]),"y":tune.grid_search([y]),"alpha": tune.grid_search((0.1)**np.linspace(5,15,7)),"lambda_1": tune.grid_search((0.1)**np.linspace(5,15,7)), "lambda_s":  tune.grid_search((0.1)**np.linspace(5,10,7)) , "beta":tune.grid_search((0.10)**np.linspace(-2,2,7)),"gamma":tune.grid_search([0.2]),"sigma2":tune.grid_search([sigma**2])}
+		config = param_search(objective_KReSP,search_space)
+
+		model = KReSP(t,y,lambda_1=config["lambda_1"],beta=config["beta"],lambda_s=config["lambda_s"],alpha=config["alpha"],gamma=config["gamma"],eps=10**-7,win_width=32,sigma2=sigma**2)
+
+		optimizer = torch.optim.Adam(model.parameters(), lr=0.01)
+		model = optimize_loop(model,optimizer,numit=300)
+
+		snr_kerreg[s,nrep] = (snr_bivariate(x,model.X.detach().numpy()))
+
+		print("sigma: " + str(sigma) + " nrep: " + str(nrep) + " KReSP SNR: "+  str(snr_kerreg[s,nrep]))
+
+		s+=1 
+
+	with open('results.pkl', 'wb') as f:
+		pickle.dump(snr_noise,f)
+		pickle.dump(snr_p2,f)
+		pickle.dump(snr_nostokes,f)
+		pickle.dump(snr_onlystokes,f)
+		pickle.dump(snr_kerreg,f)
diff --git a/denoising-plots/noisy_sig.pdf b/denoising-plots/noisy_sig.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..3c2a8cd5c97ea9dfbefa5e06b66a4fabc7264933
Binary files /dev/null and b/denoising-plots/noisy_sig.pdf differ
diff --git a/denoising-plots/noisy_sig_sphere.pdf b/denoising-plots/noisy_sig_sphere.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..7338582d5eb66bd405ff09a060e26f273411a7d9
Binary files /dev/null and b/denoising-plots/noisy_sig_sphere.pdf differ
diff --git a/denoising-plots/plot.py b/denoising-plots/plot.py
new file mode 100644
index 0000000000000000000000000000000000000000..b176937b62a4426490a7d0004f502e8311fc04a4
--- /dev/null
+++ b/denoising-plots/plot.py
@@ -0,0 +1,94 @@
+import numpy as np
+import quaternion
+import matplotlib.pyplot as plt
+import bispy as bsp
+from sklearn.preprocessing import normalize
+from utils import plot_on_sphere
+import pickle
+
+## PLOT FOR THE RESULTS
+
+with open('results.pkl', 'rb') as f:
+	snr_noise = pickle.load(f)
+	snr_p2 = pickle.load(f)
+	snr_nostokes = pickle.load(f)
+	snr_onlystokes = pickle.load(f)
+	snr_kerreg = pickle.load(f)
+
+plt.plot(np.mean(snr_noise,axis=1),np.mean(snr_p2,axis=1),label="$\mathbf{x}_{mix}$",marker="d")
+ci = 1.96 * np.std(snr_p2,axis=1)/np.sqrt(snr_p2.shape[1])
+plt.fill_between(np.mean(snr_noise,axis=1), (np.mean(snr_p2,axis=1)-ci), (np.mean(snr_p2,axis=1)+ci), color='b', alpha=.1)
+
+plt.plot(np.mean(snr_noise,axis=1),np.mean(snr_nostokes,axis=1),label="No stokes reg.",marker="s")
+ci = 1.96 * np.std(snr_nostokes,axis=1)/np.sqrt(snr_nostokes.shape[1])
+plt.fill_between(np.mean(snr_noise,axis=1), (np.mean(snr_nostokes,axis=1)-ci), (np.mean(snr_nostokes,axis=1)+ci), color="orange", alpha=.1)
+
+plt.plot(np.mean(snr_noise,axis=1),np.mean(snr_onlystokes,axis=1),label="Only Stokes reg.",marker="o")
+ci = 1.96 * np.std(snr_onlystokes,axis=1)/np.sqrt(snr_onlystokes.shape[1])
+plt.fill_between(np.mean(snr_noise,axis=1), (np.mean(snr_onlystokes,axis=1)-ci), (np.mean(snr_onlystokes,axis=1)+ci), color='g', alpha=.1)
+
+plt.plot(np.mean(snr_noise,axis=1),np.mean(snr_kerreg,axis=1),label="$\mathbf{x}_{ker}$",marker="x")
+ci = 1.96 * np.std(snr_kerreg,axis=1)/np.sqrt(snr_kerreg.shape[1])
+plt.fill_between(np.mean(snr_noise,axis=1), (np.mean(snr_kerreg,axis=1)-ci), (np.mean(snr_kerreg,axis=1)+ci), color='r', alpha=.1)
+
+
+plt.plot(np.mean(snr_noise,axis=1),np.mean(snr_noise,axis=1),label="no denoising",linestyle="--")
+
+plt.ylabel("Output SNR")
+plt.grid()
+plt.xlabel("Input SNR")
+plt.legend()
+plt.tight_layout()
+
+plt.savefig("denoising-perf.pdf")
+plt.close()
+
+
+N = 1024
+t = np.linspace(0, 2*np.pi/4, 1024) # time vector
+dt = t[1] - t[0]
+theta1 = np.pi/4 - 2*t
+chi1 = np.pi/16 - t
+phi1 = 0 
+f0 = 25/N/dt 
+S0 = bsp.utils.windows.hanning(N)
+
+x_quad = bsp.signals.bivariateAMFM(S0, theta1, chi1, 2*np.pi*f0*t+ phi1)
+
+
+with open('example.pkl', 'rb') as f:
+	x=pickle.load(f)
+	y=pickle.load(f)
+	x_stis=pickle.load(f)
+	x_nostokes=pickle.load(f)
+	x_onlystokes=pickle.load(f)
+	x_kerreg=pickle.load(f)
+
+bsp.utils.visual.plot2D(t,x_stis)
+plt.savefig("denoised_via_all_terms.pdf")
+
+bsp.utils.visual.plot2D(t,x_nostokes)
+plt.savefig("denoised_via_no_stokes.pdf")
+
+bsp.utils.visual.plot2D(t,x_onlystokes)
+plt.savefig("denoised_via_no_signal_smoother.pdf")
+
+bsp.utils.visual.plot2D(t,x_kerreg)
+plt.savefig("denoised_via_kerreg.pdf")
+
+ax = plt.figure().add_subplot(projection='3d')
+plot_on_sphere(x_stis,ax,label="$\\mathbf{x}_{mix}$",t=t,scatter=True)
+plot_on_sphere(x_quad,ax,label="True Signal",t=[0],scatter=False)
+
+plt.legend()
+plt.tight_layout()
+plt.savefig("xmix_norm_stokes.pdf")
+plt.close()
+ax = plt.figure().add_subplot(projection='3d')
+
+plot_on_sphere(x_kerreg,ax,label="$\\mathbf{x}_{ker}$",t=t,scatter=True)
+plot_on_sphere(x_quad,ax,label="True Signal",t=[0],scatter=False)
+plt.legend()
+plt.tight_layout()
+plt.savefig("xker_norm_stokes.pdf")
+plt.close()
diff --git a/denoising-plots/results.pkl b/denoising-plots/results.pkl
new file mode 100644
index 0000000000000000000000000000000000000000..8ce7df358b78caa6ef3802275594cf64b9f372df
Binary files /dev/null and b/denoising-plots/results.pkl differ
diff --git a/denoising-plots/utils.py b/denoising-plots/utils.py
new file mode 100644
index 0000000000000000000000000000000000000000..e314279fab344119d893bea67eb15673086500b6
--- /dev/null
+++ b/denoising-plots/utils.py
@@ -0,0 +1,240 @@
+import torch
+import numpy as np
+import quaternion
+from ray import train, tune
+from ray.tune.experimental.output import AirVerbosity, get_air_verbosity 
+from matplotlib import pyplot as plt
+import bispy as bsp
+
+# MODELS
+## Prototyping STIS by pytorch 
+class STIS(torch.nn.Module):
+    def __init__(self,t,y, device='cpu', lambdax=0.0,lambdaS=0.0,beta1=0.0,beta2=0.0,sigma2=1.0,p=2):
+        super().__init__()
+
+        # Set parameters
+        self.dt = torch.tensor(np.diff(t))
+        self.lambdax = lambdax
+        self.lambdaS = lambdaS
+        self.beta1 = beta1
+        self.beta2 = beta2
+        self.sigma2 = sigma2
+        self.p = p
+        self.N = len(y)
+
+        # Generate the quaternion embedding
+        self.y = torch.tensor(y).to(device)
+        ffty = torch.fft.fft(self.y,dim=0)
+        ffty[1:((self.N)//2),:] *= 2.0
+        ffty[(self.N//2):self.N,:] = 0.0 + 0.0j
+
+        y_as = torch.fft.ifft(ffty,dim=0)
+        self.yquad = torch.cat((y_as.real, y_as.imag),1)
+
+        # Initialize the optimization parameter
+        self.X = torch.nn.Parameter(torch.tensor(y).to(device))   
+
+        # Compute the Stokes of the measurements
+        self.Sy = torch.tensor(np.zeros([self.N,4]))
+        self.Sy[:,0] = self.yquad[:,0]**2 + self.yquad[:,1]**2 + self.yquad[:,2]**2 + self.yquad[:,3]**2
+        self.Sy[:,1] = self.yquad[:,0]**2 - self.yquad[:,1]**2 + self.yquad[:,2]**2 - self.yquad[:,3]**2
+        self.Sy[:,2] = 2*( self.yquad[:,0]* self.yquad[:,1] + self.yquad[:,2]*self.yquad[:,3] )
+        self.Sy[:,3] = 2*( self.yquad[:,1]* self.yquad[:,2] - self.yquad[:,0]*self.yquad[:,3] )
+
+    def forward(self):
+        fidelity = torch.norm(self.X.T - self.y.T)**2
+        
+        # Generate the quaternion embedding
+        fftX = torch.fft.fft(self.X,dim=0)
+        fftX[1:self.N//2,:] *= 2.0
+        fftX[(self.N//2):self.N,:] = 0.0 + 0.0j
+
+        Xhilbert = torch.fft.ifft(fftX,dim=0)
+        self.Xquad = torch.cat((Xhilbert.real, Xhilbert.imag),1)
+        
+        ## INSTANTANEOUS STOKES PARAMETERS COMPUTED 
+        self.Sx = torch.tensor(np.zeros([self.N,4]))
+ 
+        self.Sx[:,0] = self.Xquad[:,0]**2 + self.Xquad[:,1]**2 + self.Xquad[:,2]**2 + self.Xquad[:,3]**2
+        self.Sx[:,1] = self.Xquad[:,0]**2 - self.Xquad[:,1]**2 + self.Xquad[:,2]**2 - self.Xquad[:,3]**2
+        self.Sx[:,2] = 2*( self.Xquad[:,0]* self.Xquad[:,1] + self.Xquad[:,2]*self.Xquad[:,3] )
+        self.Sx[:,3] = 2*( self.Xquad[:,1]* self.Xquad[:,2] - self.Xquad[:,0]*self.Xquad[:,3] )       
+
+        fidelity_stokes = self.beta1*(torch.norm( self.Sx[:,1:3] -self.Sy[:,1:3],self.p )**self.p) + self.beta2*(torch.sum(self.Sx[:,0]/(2*self.sigma2) + 0.5*torch.log(self.Sx[:,0])  - torch.log( torch.special.i1e(torch.sqrt(self.Sx[:,0]*self.Sy[:,0])/(self.sigma2)) )  -torch.sqrt(self.Sx[:,0]*self.Sy[:,0])/(self.sigma2)) ) 
+
+        velS = torch.diff(self.Sx,dim=0).T / self.dt 
+        velX = torch.diff(self.X.T) / self.dt
+        reg1 = self.lambdaS*torch.sum(torch.norm(velS,dim=1)**2)
+        reg2 = self.lambdax*(torch.norm(velX)**2)
+	
+        return (fidelity/self.N + reg1) + (fidelity_stokes/self.N + reg2)
+
+class KReSP(torch.nn.Module):
+    def __init__(self,t,y, device='cpu',lambda_1=1,beta=0.1,lambda_s=0.1,alpha=0.1,gamma=1.0,eps=10**-7,win_width=32,sigma2=0.01):
+        super().__init__()
+
+        self.dt = torch.tensor(np.diff(t))
+        self.sigma2 = sigma2
+        self.beta=beta
+        self.lambda_s = lambda_s
+        self.lambda_1 = lambda_1 
+        self.alpha = alpha
+        self.gamma = gamma 
+        self.t = torch.unsqueeze(torch.tensor(t),1)
+        self.N = len(y)
+        self.win_width = win_width
+        # Generate the quaternion embedding
+        self.y = torch.tensor(y).to(device)
+        ffty = torch.fft.fft(self.y,dim=0)
+        ffty[1:((self.N)//2),:] *= 2.0
+        ffty[(self.N//2):self.N,:] = 0.0 + 0.0j
+
+        y_as = torch.fft.ifft(ffty,dim=0)
+        self.yquad = torch.cat((y_as.real, y_as.imag),1)
+
+        # Initialize the optimization parameter
+        self.X = torch.nn.Parameter(torch.tensor(np.random.randn(self.N,2)).to(device))   
+
+        # Compute the Stokes of the data
+        self.Sy = torch.tensor(np.zeros([self.N,4]))
+        self.Sy[:,0] = self.yquad[:,0]**2 + self.yquad[:,1]**2 + self.yquad[:,2]**2 + self.yquad[:,3]**2
+        self.Sy[:,1] = self.yquad[:,0]**2 - self.yquad[:,1]**2 + self.yquad[:,2]**2 - self.yquad[:,3]**2
+        self.Sy[:,2] = 2*( self.yquad[:,0]* self.yquad[:,1] + self.yquad[:,2]*self.yquad[:,3] )
+        self.Sy[:,3] = 2*( self.yquad[:,1]* self.yquad[:,2] - self.yquad[:,0]*self.yquad[:,3] )
+
+        # Compute the normalized parameters 
+        self.SnY = torch.empty(self.N,3)
+        self.SnY[:,0] = self.Sy[:,1]/ self.Sy[:,0]
+        self.SnY[:,1] = self.Sy[:,2]/ self.Sy[:,0]
+        self.SnY[:,2] = self.Sy[:,3]/ self.Sy[:,0]
+        
+        self.X = torch.nn.Parameter(torch.tensor(y).to(device))   
+        # The Gaussian kernel
+        self.K = torch.exp(-torch.cdist(self.t[self.win_width-self.win_width//2:self.win_width+self.win_width//2],self.t[self.win_width-self.win_width//2:self.win_width+self.win_width//2])/(2*(self.gamma**2)))/(self.gamma*np.sqrt(2*torch.pi))
+        self.eps = eps
+    def forward(self):
+        fidelity = (torch.norm(self.X.T - self.y.T)**2)
+        velX = torch.diff(self.X.T)/self.dt
+
+        ## HILBERT TRANSFORM COMPUTED IN SPECTRAL DOMAIN 
+        fftX = torch.fft.fft(self.X,dim=0)
+        fftX[1:(self.N//2),:] *= 2.0
+        fftX[(self.N//2):self.N,:] = 0.0 + 0.0j
+
+        Xhilbert = torch.fft.ifft(fftX,dim=0)
+        self.Xquad = torch.cat((Xhilbert.real, Xhilbert.imag),1)
+        ## INSTANTANEOUS STOKES PARAMETERS COMPUTED 
+        self.Sx = torch.tensor(np.zeros([self.N,4]))
+ 
+        self.Sx[:,0] = self.Xquad[:,0]**2 + self.Xquad[:,1]**2 + self.Xquad[:,2]**2 + self.Xquad[:,3]**2
+        self.Sx[:,1] = self.Xquad[:,0]**2 - self.Xquad[:,1]**2 + self.Xquad[:,2]**2 - self.Xquad[:,3]**2
+        self.Sx[:,2] = 2*( self.Xquad[:,0]* self.Xquad[:,1] + self.Xquad[:,2]*self.Xquad[:,3] )
+        self.Sx[:,3] = 2*( self.Xquad[:,1]* self.Xquad[:,2] - self.Xquad[:,0]*self.Xquad[:,3] )       
+
+        # FIDELITY FUNCTION FOR S0
+        fidelity_S0 = self.beta*(torch.sum(self.Sx[:,0]/(2*self.sigma2) + 0.5*torch.log(self.Sx[:,0]  - torch.log( torch.special.i1e(torch.sqrt(self.Sx[:,0]*self.Sy[:,0])/(self.sigma2)) ) ) -torch.sqrt(self.Sx[:,0]*self.Sy[:,0])/(self.sigma2) )) 
+        velS0 = torch.diff(self.Sx[:,0])/self.dt
+
+        # COMPUTING NORMALIZED STOKES PARAMS 
+        SnX = torch.empty(self.N,3)
+        SnX[:,0] = self.Sx[:,1] / self.Sx[:,0] 
+        SnX[:,1] = self.Sx[:,2] / self.Sx[:,0]  
+        SnX[:,2] = self.Sx[:,3] / self.Sx[:,0]     
+
+        # COMPUTING THE KERNEL COST FUNCTION HERE 
+        D = torch.zeros(1)
+        if(self.alpha >0.0):
+            for i in range(self.win_width//2,self.N-self.win_width//2,self.win_width):
+                D += torch.sum( (self.Sx[i-self.win_width//2:i+self.win_width//2,0]*self.K)*self.K*(torch.acos(torch.clamp(torch.matmul(SnX[i-self.win_width//2:i+self.win_width//2,:], self.SnY[i-self.win_width//2:i+self.win_width//2,:].T),-1+self.eps,1-self.eps))))                     
+        else:
+            D= 0.0
+        
+        kerreg_loss = self.alpha*(D)
+        reg2 = self.lambda_1*(torch.norm(velX)**2)
+        reg3 = self.lambda_s*(torch.norm(velS0)**2)
+
+        return (fidelity/self.N  + fidelity_S0/self.N + kerreg_loss +reg2 + reg3)
+    
+# HELPER FUNCTIONS
+def optimize_loop(model,optimizer,numit=1000): 
+    loss_curve = []
+    for e in range(numit):
+        l = model()
+        l.backward()
+        optimizer.step()
+        optimizer.zero_grad()
+        loss_curve.append(l.item())
+    return model 
+
+def snr_bivariate(s_true,s_noisy):
+    MSE = 10*np.log10(np.mean(np.linalg.norm(s_true-s_noisy))**2)
+    MEAN = 10*np.log10(np.mean(np.linalg.norm(s_true))**2)
+    return MEAN - MSE 
+
+def param_search(objective,search_space):
+
+    tuner = tune.Tuner( 
+        objective,
+        tune_config=tune.TuneConfig(
+            metric="snr",
+            mode="max",
+        ),
+        run_config=train.RunConfig(
+            verbose = 0,
+            stop={"training_iteration": 1},
+        ),
+        param_space=search_space,
+    )
+    results = tuner.fit()
+    return results.get_best_result().config
+
+
+def objective_STIS(config):
+    model = STIS(config["t"],config["y"],lambdax=config["lambdax"],lambdaS=config["lambdaS"],beta1=config["beta1"],beta2=config["beta2"],sigma2=config["sigma2"])
+    optimizer = torch.optim.Adam(model.parameters(), lr=0.01)
+
+    while True:
+        model = optimize_loop(model,optimizer,numit=1000)
+        X_STIS = model.X.detach().numpy()
+        snr = snr_bivariate(config["x"],X_STIS)
+        train.report({"snr": snr})  # Report to Tune
+
+def objective_KReSP(config):
+    model = KReSP(config["t"],config["y"],lambda_1=config["lambda_1"],beta=config["beta"],lambda_s=config["lambda_s"],alpha=config["alpha"],gamma=config["gamma"],eps=10**-7,sigma2=config["sigma2"])
+
+    optimizer = torch.optim.Adam(model.parameters(), lr=0.01)
+
+    while True:
+        model = optimize_loop(model,optimizer,numit=300)
+        X_KReSP = model.X.detach().numpy()
+        snr = snr_bivariate(config["x"],X_KReSP)
+        train.report({"snr": snr})  # Report to Tune
+
+
+
+def plot_on_sphere(x_quad,ax,label="",t=[0],scatter=False):
+
+    if(len(t)==1):
+        t= np.zeros(len(x_quad))
+    S = np.norm(x_quad) + bsp.utils.StokesNorm(x_quad)
+    S0 = quaternion.as_float_array(S)[:,0]
+    S1 = quaternion.as_float_array(S)[:,2]
+    S2 = quaternion.as_float_array(S)[:,3]
+    S3 = quaternion.as_float_array(S)[:,1]
+
+    u = np.linspace(0, 2 * np.pi, 20)
+    v = np.linspace(0, np.pi, 20)
+    a = 1 * np.outer(np.cos(u), np.sin(v))
+    b = 1 * np.outer(np.sin(u), np.sin(v))
+    c = 1 * np.outer(np.ones(np.size(u)), np.cos(v))
+    # Plot the surface
+    ax.plot_surface(a,b,c,alpha=0.051)
+    ax.set_xlabel("$s_1(t)$")
+    ax.set_ylabel("$s_2(t)$")
+    ax.set_zlabel("$s_3(t)$")    
+    if(scatter):
+        ax.scatter(S1/S0,S2/S0,S3/S0,label=label,c=t,s=10,zorder=1)
+    else:
+        ax.plot(S1/S0,S2/S0,S3/S0,label=label,linewidth=6,color="red",alpha=0.5,zorder=100)
+
+
diff --git a/denoising-plots/xker_norm_stokes.pdf b/denoising-plots/xker_norm_stokes.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..fc6d4fae0de7e9c93c59435938500fcfdda8dc98
Binary files /dev/null and b/denoising-plots/xker_norm_stokes.pdf differ
diff --git a/denoising-plots/xmix_norm_stokes.pdf b/denoising-plots/xmix_norm_stokes.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..05ffd10b0468f48806800533d9a2871cd46a60a0
Binary files /dev/null and b/denoising-plots/xmix_norm_stokes.pdf differ
diff --git a/likelihood-plot/main.py b/likelihood-plot/main.py
new file mode 100644
index 0000000000000000000000000000000000000000..d490c4650426016f3414bd4257f976ecb3070bdc
--- /dev/null
+++ b/likelihood-plot/main.py
@@ -0,0 +1,141 @@
+import numpy as np
+import scipy
+import matplotlib.pyplot as plt
+import quaternion  # load the quaternion module
+import bispy as bsp
+import torch
+import tqdm
+from sklearn.preprocessing import normalize
+
+# time vector
+N = 2048 # length of the signal
+t = np.linspace(0, 100, N)
+dt = t[1]-t[0]
+# ellipse parameters 
+theta1 = np.pi/4  -t
+chi1 = np.pi/6  # -t/10
+phi1 = np.pi/6
+f0 = 100/N/dt 
+env = 1.0# bsp.utils.windows.hanning(N)
+x1 = bsp.signals.bivariateAMFM(env, theta1, chi1, 2*np.pi*f0*t+ phi1)
+x1_arr = quaternion.as_float_array(x1)
+
+ts = 105
+
+
+S0 = np.norm(x1)[ts]
+S = bsp.utils.StokesNorm(x1)
+S = quaternion.as_float_array(S)
+S1 = S[ts,2]
+S2 = S[ts,3]
+S3 = S[ts,1]
+
+
+NRep = 100
+NoiseRep = 1000
+Sigma_REP = 100
+Sigma_exp_range =  np.linspace(start=-3, stop=3, num=Sigma_REP)
+
+score_L1 = np.zeros([Sigma_REP,NRep])
+score_L2 = np.zeros([Sigma_REP,NRep])
+score_L3 = np.zeros([Sigma_REP,NRep]) 
+
+score_N1 = np.zeros([Sigma_REP,NRep])
+score_N2 = np.zeros([Sigma_REP,NRep])
+score_N3 = np.zeros([Sigma_REP,NRep])
+
+for exprep in tqdm.tqdm(range(NRep)): # For each experiment trial
+    for (idx,sigma_exp) in enumerate((Sigma_exp_range)): # For each sigma
+        sigma = 10**sigma_exp
+        scale = 2*(sigma**2)
+        std = 2*sigma*np.sqrt(S0)
+
+        S0_emp = np.zeros([NoiseRep]) 
+        S1_emp = np.zeros([NoiseRep]) 
+        S2_emp = np.zeros([NoiseRep]) 
+        S3_emp = np.zeros([NoiseRep]) 
+
+        for i in range(NoiseRep): # For each noise realization 
+
+            ## Hilbert transform of the noise given to j and k components
+            n = np.zeros([N,4])
+            noise_complex = np.random.randn(N,2)*sigma
+            uH = np.imag(scipy.signal.hilbert(noise_complex[:,0]))
+            vH = np.imag(scipy.signal.hilbert(noise_complex[:,1]))
+            n[:,0] = noise_complex[:,0]
+            n[:,1] = noise_complex[:,1]
+            n[:,2] = uH
+            n[:,3] = vH
+
+            n = quaternion.from_float_array(n)
+            x = (x1 + n) 
+
+            S0_emp[i] = np.norm(x[ts])
+
+            S_emp = bsp.utils.StokesNorm(x)
+            S_emp = quaternion.as_float_array(S_emp)
+            S1_emp[i] = S_emp[ts,2]
+            S2_emp[i] = S_emp[ts,3]
+            S3_emp[i] = S_emp[ts,1]
+
+        score_L1[idx,exprep] = scipy.stats.kstest(S1_emp,scipy.stats.laplace(S1,scale).cdf).pvalue
+        score_N1[idx,exprep] = scipy.stats.kstest(S1_emp,scipy.stats.norm(S1,std).cdf).pvalue
+
+        score_L2[idx,exprep] = scipy.stats.kstest(S2_emp,scipy.stats.laplace(S2,scale).cdf).pvalue
+        score_N2[idx,exprep] = scipy.stats.kstest(S2_emp,scipy.stats.norm(S2,std).cdf).pvalue
+
+        score_L3[idx,exprep] = scipy.stats.kstest(S3_emp,scipy.stats.laplace(S3,scale).cdf).pvalue
+        score_N3[idx,exprep] = scipy.stats.kstest(S3_emp,scipy.stats.norm(S3,std).cdf).pvalue
+
+
+plt.figure(figsize=[12,3])
+plt.subplot(1,3,1)
+plt.title("S1")
+plt.plot(10**Sigma_exp_range,np.mean(score_L1,axis=1),label="Laplace")
+ci = 1.96 * np.std(score_L1,axis=1)/np.sqrt(NRep)
+plt.fill_between(10**Sigma_exp_range,(np.mean(score_L1,axis=1)-ci),(np.mean(score_L1,axis=1)+ci), alpha=.5)
+
+plt.plot(10**Sigma_exp_range,np.mean(score_N1,axis=1),label="Gaussian")
+ci = 1.96 * np.std(score_N1,axis=1)/np.sqrt(NRep)
+plt.fill_between(10**Sigma_exp_range,(np.mean(score_N1,axis=1)-ci),(np.mean(score_N1,axis=1)+ci), alpha=.5)
+plt.hlines(0.05,10**Sigma_exp_range[0],10**Sigma_exp_range[-1],'r',linestyles='--',label="p=0.05")
+
+plt.xscale("log")
+plt.ylabel("p-Value")
+plt.xlabel("Noise Variance $\\sigma$")
+plt.legend()
+
+plt.subplot(1,3,2)
+plt.title("S2")
+plt.plot(10**Sigma_exp_range,np.mean(score_L2,axis=1),label="Laplace")
+ci = 1.96 * np.std(score_L2,axis=1)/np.sqrt(NRep)
+plt.fill_between(10**Sigma_exp_range,(np.mean(score_L2,axis=1)-ci),(np.mean(score_L2,axis=1)+ci), alpha=.5)
+
+plt.plot(10**Sigma_exp_range,np.mean(score_N2,axis=1),label="Gaussian")
+ci = 1.96 * np.std(score_N2,axis=1)/np.sqrt(NRep)
+plt.fill_between(10**Sigma_exp_range,(np.mean(score_N2,axis=1)-ci),(np.mean(score_N2,axis=1)+ci), alpha=.5)
+plt.hlines(0.05,10**Sigma_exp_range[0],10**Sigma_exp_range[-1],'r',linestyles='--',label="p=0.05")
+
+plt.xscale("log")
+plt.ylabel("p-Value")
+plt.xlabel("Noise Variance $\\sigma$")
+plt.legend()
+
+plt.subplot(1,3,3)
+plt.title("S3")
+plt.plot(10**Sigma_exp_range,np.mean(score_L3,axis=1),label="Laplace")
+ci = 1.96 * np.std(score_L3,axis=1)/np.sqrt(NRep)
+plt.fill_between(10**Sigma_exp_range,(np.mean(score_L3,axis=1)-ci),(np.mean(score_L3,axis=1)+ci), alpha=.5)
+
+plt.plot(10**Sigma_exp_range,np.mean(score_N3,axis=1),label="Gaussian")
+ci = 1.96 * np.std(score_N3,axis=1)/np.sqrt(NRep)
+plt.fill_between(10**Sigma_exp_range,(np.mean(score_N3,axis=1)-ci),(np.mean(score_N3,axis=1)+ci), alpha=.5)
+plt.hlines(0.05,10**Sigma_exp_range[0],10**Sigma_exp_range[-1],'r',linestyles='--',label="p=0.05")
+
+plt.xscale("log")
+plt.ylabel("p-Value")
+plt.xlabel("Noise Variance $\\sigma$")
+
+plt.legend()
+plt.tight_layout()
+plt.savefig("sigma-vs-pvalue")
\ No newline at end of file
diff --git a/likelihood-plot/sigma-vs-pvalue.png b/likelihood-plot/sigma-vs-pvalue.png
new file mode 100644
index 0000000000000000000000000000000000000000..a3d62b7d5777253c993384c1f0dc72a4cc49ab33
Binary files /dev/null and b/likelihood-plot/sigma-vs-pvalue.png differ