From 1075f995eee7bf432725b57944ce7b0e2c89ae2f Mon Sep 17 00:00:00 2001 From: Michael Pilosov <40366263+mathematicalmichael@users.noreply.github.com> Date: Sat, 6 Feb 2021 20:21:36 -0700 Subject: [PATCH] updates to contour plots (#7) * wip for contour plots * adding notebook * adding updated scripts --- notebooks/Contours.ipynb | 332 +++++++++++++++++++++++++++++++++++ src/mud_examples/helpers.py | 8 + src/mud_examples/plotting.py | 41 ++++- 3 files changed, 374 insertions(+), 7 deletions(-) create mode 100644 notebooks/Contours.ipynb diff --git a/notebooks/Contours.ipynb b/notebooks/Contours.ipynb new file mode 100644 index 0000000..1fe4e1a --- /dev/null +++ b/notebooks/Contours.ipynb @@ -0,0 +1,332 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "presentation = False\n", + "save = True" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "import ipywidgets as wd\n", + "from mud_examples.plotting import plot_2d_contour_example as plot_full\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import matplotlib\n", + "if not presentation:\n", + " matplotlib.rcParams['mathtext.fontset'] = 'stix'\n", + " matplotlib.rcParams['font.family'] = 'STIXGeneral'\n", + " fdir = 'contours'\n", + "else:\n", + " fdir = '../presentation/figures'\n", + "\n", + "matplotlib.rcParams['font.size'] = 24\n", + "matplotlib.backend = 'Agg'" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "plt.rcParams['figure.figsize'] = 10,10" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "lam_true = np.array([0.7, 0.3])\n", + "initial_mean = np.array([0.25, 0.25])\n", + "A = np.array([[1, 1]])\n", + "b = np.zeros((1,1))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Projection is to the nullspace of $A\\Sigma_{in}$, not $A$. If $\\Sigma_{in}$ is non-diagonal, then performing $\\Sigma_{in} = U \\Sigma W^T$ yields a transformation matrix $U$ that \n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Slider objects" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "63f3e6e6ef95429dacf80e99c531278f", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "VBox(children=(interactive(children=(FloatSlider(value=0.5, continuous_update=False, description='$\\\\sigma_{11…" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "out = wd.Output()\n", + "cov_11_slider = wd.FloatSlider(value=1, min=0.05, max=1, step=0.05, description='$\\sigma_{11}$', continuous_update=False)\n", + "cov_01_slider = wd.FloatSlider(value=0, min=-0.95-1E-4, max=0.95+1E-4, step=0.01, description='$\\sigma_{01}=\\sigma_{10}$', continuous_update=False)\n", + "\n", + "def mv_slider(val):\n", + " pr_slide.max = tk_slide.value + 1E-4\n", + "\n", + "def ensure_positive_determinant(val):\n", + " \"\"\"\n", + " Bound limits of covariance to ensure det(C)>0\n", + " \"\"\"\n", + " rt = np.sqrt(cov_11_slider.value)\n", + " # need to account for slider numerical rounding\n", + " cov_01_slider.max = rt-1E-4\n", + " cov_01_slider.min = -rt+1E-4\n", + "\n", + "cov_11_slider.observe(ensure_positive_determinant)\n", + "\n", + "tk_slide = wd.FloatSlider(value=1, min=0, max=1, step=0.02, description='Tikonov', continuous_update=False)\n", + "pr_slide = wd.FloatSlider(value=1, min=0, max=1, step=0.02, description='Unreg', continuous_update=False)\n", + "tk_slide.observe(mv_slider)\n", + "full_check = wd.Checkbox(value=True, description='Show Full')\n", + "data_check = wd.Checkbox(value=True, description='Show Data')\n", + "numr_check = wd.Checkbox(value=False, description='Numerical')\n", + "comparison = wd.Checkbox(value=False, description='Compare MAP')\n", + "\n", + "fig_title = wd.Text(value='latest_figure.png')\n", + "obs_std_slider = wd.FloatSlider(value=1, min=0.01, max=1, step=0.01, description='Obs. Std', continuous_update=False)\n", + "# pr_slide.observe(mv_slider)\n", + "cov_11_slider.value = 0.5\n", + "cov_01_slider.value = -0.25\n", + "obs_std_slider.value = 0.5\n", + "\n", + "I = wd.interactive(plot_full, A=wd.fixed(A), b=wd.fixed(b), save=wd.fixed(save),\n", + " param_ref=wd.fixed(lam_true), compare=comparison,\n", + " cov_01=cov_01_slider, cov_11=cov_11_slider, \n", + " initial_mean=wd.fixed(initial_mean),\n", + " alpha=tk_slide,\n", + " omega=pr_slide,\n", + " show_full=full_check,\n", + " show_data=data_check,\n", + " show_est=numr_check,\n", + " obs_std = obs_std_slider,\n", + " figname=fig_title)\n", + "\n", + "display(wd.VBox([I, out]))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Figures\n", + "---\n" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "assert save is True, \"No saving requested. Stopping notebook.\"" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Data Mismatch" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "fig_title.value = f'{fdir}/data_mismatch_contour.png'\n", + "data_check.value = True\n", + "full_check.value = False\n", + "tk_slide.value = 0\n", + "pr_slide.value = 0" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Tikonov Regularization" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "fig_title.value = f'{fdir}/tikonov_contour.png'\n", + "tk_slide.value = 1\n", + "pr_slide.value = 0\n", + "data_check.value = False\n", + "full_check.value = False" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Modified Regularization" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "fig_title.value = f'{fdir}/consistent_contour.png'\n", + "tk_slide.value = 1\n", + "pr_slide.value = 1\n", + "data_check.value = False\n", + "full_check.value = False" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Classical Solution" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "fig_title.value = f'{fdir}/classical_solution.png'\n", + "tk_slide.value = 1\n", + "pr_slide.value = 0\n", + "data_check.value = True\n", + "full_check.value = True" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Consistent Solution" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [], + "source": [ + "fig_title.value = f'{fdir}/consistent_solution.png'\n", + "tk_slide.value = 1\n", + "pr_slide.value = 1\n", + "data_check.value = True\n", + "full_check.value = True" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Comparison" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [], + "source": [ + "fig_title.value = f'{fdir}/map_compare_contour.png'\n", + "data_check.value = True\n", + "full_check.value = True\n", + "tk_slide.value = 1\n", + "pr_slide.value = 0\n", + "comparison.value = True\n", + "cov_01_slider.value = -0.5" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "---\n", + "---\n", + "---" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Cleanup" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [], + "source": [ + "!rm latest_figure.png" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "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.8.6" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/src/mud_examples/helpers.py b/src/mud_examples/helpers.py index db10680..121b904 100644 --- a/src/mud_examples/helpers.py +++ b/src/mud_examples/helpers.py @@ -1,6 +1,14 @@ +import os + import numpy as np from mud.funs import mud_sol, map_sol + +def check_dir(directory): + if not os.path.exists(directory): + os.makedirs(directory) + + def fit_log_linear_regression(input_values, output_values): x, y = np.log10(input_values), np.log10(output_values) X, Y = np.vander(x, 2), np.array(y).reshape(-1, 1) diff --git a/src/mud_examples/plotting.py b/src/mud_examples/plotting.py index a63b277..d98f2ac 100644 --- a/src/mud_examples/plotting.py +++ b/src/mud_examples/plotting.py @@ -13,6 +13,7 @@ from mud.funs import mud_sol, map_sol from mud.norm import full_functional, norm_input, norm_data, norm_predicted from mud.plot import make_2d_unit_mesh +from mud_examples.helpers import check_dir def fit_log_linear_regression(input_values, output_values): @@ -27,8 +28,10 @@ def plot_2d_contour_example(A=np.array([[1, 1]]), b=np.zeros([1, 1]), # noqa: C cov_11=0.5, cov_01=-0.25, initial_mean=np.array([0.25, 0.25]), alpha=1, omega=1, obs_std=1, - show_full=True, show_data=True, show_est=False, - fsize=42, figname='latest_figure.png', save=False): + show_full=True, show_data=True, + show_est=False, param_ref=None, compare=False, + fsize=42, figname='latest_figure.png', save=False, + ): """ alpha: float in [0, 1], weight of Tikhonov regularization omega: float in [0, 1], weight of Modified regularization @@ -94,6 +97,7 @@ def plot_2d_contour_example(A=np.array([[1, 1]]), b=np.zeros([1, 1]), # noqa: C cmap=cm.viridis, alpha=0.25) plt.axis('equal') + if alpha + omega > 0: plt.scatter(initial_mean[0], initial_mean[1], label='Initial Mean', @@ -102,6 +106,17 @@ def plot_2d_contour_example(A=np.array([[1, 1]]), b=np.zeros([1, 1]), # noqa: C plt.annotate('Initial Mean', (initial_mean[0] + 0.001 * fsize, initial_mean[1] - 0.001 * fsize), fontsize=fsize, backgroundcolor="w") + else: + if param_ref is not None: + plt.scatter(param_ref[0], param_ref[1], + label='$\\lambda^\dagger$', + color='k', s=msize, marker='*') + if compare: + plt.annotate('Truth', + (param_ref[0] + 0.0005 * fsize, param_ref[1] + 0.0005 * fsize), + fontsize=fsize, backgroundcolor="w") + + show_mud = omega > 0 or compare if show_full: # scatter and line from origin to least squares @@ -125,6 +140,7 @@ def plot_2d_contour_example(A=np.array([[1, 1]]), b=np.zeros([1, 1]), # noqa: C plt.scatter(map_pt[0], map_pt[1], label='min: Tk', color='xkcd:blue', marker='o', s=3 * msize, zorder=10) + if (alpha > 0 and omega != 1): # analytical MAP point map_pt_eq = map_sol(A, b, observed_data_mean, initial_mean, initial_cov, @@ -132,10 +148,18 @@ def plot_2d_contour_example(A=np.array([[1, 1]]), b=np.zeros([1, 1]), # noqa: C plt.scatter(map_pt_eq[0], map_pt_eq[1], label='MAP', color='xkcd:orange', marker='x', s=msize, lw=10, zorder=10) - plt.annotate('MAP', - (map_pt_eq[0] + 0.001 * fsize, map_pt_eq[1] - 0.001 * fsize), - fontsize=fsize, backgroundcolor="w") - if omega > 0: # analytical MUD point + + if compare: # second map point has half the regularization strength + plt.annotate('MAP$_{\\alpha}$', + (map_pt_eq[0] - 0.004 * fsize, map_pt_eq[1] - 0.002 * fsize), + fontsize=fsize, backgroundcolor="w") + + else: + plt.annotate('MAP$_{\\alpha}$', + (map_pt_eq[0] + 0.0001 * fsize, map_pt_eq[1] - 0.002 * fsize), + fontsize=fsize, backgroundcolor="w") + + if show_mud: # analytical MUD point mud_pt_eq = mud_sol(A, b, observed_data_mean, initial_mean, initial_cov) plt.scatter(mud_pt_eq[0], mud_pt_eq[1], @@ -151,7 +175,7 @@ def plot_2d_contour_example(A=np.array([[1, 1]]), b=np.zeros([1, 1]), # noqa: C v = v[::-1] # in 2D, we can just swap entries and put a negative sign in front of one v[0] = - v[0] - if show_full and omega > 0: + if show_full and show_mud: # grid search to find upper/lower bounds of line being drawn. # importance is the direction, image is nicer with a proper origin/termination s = np.linspace(-1, 1, 1000) @@ -173,6 +197,9 @@ def plot_2d_contour_example(A=np.array([[1, 1]]), b=np.zeros([1, 1]), # noqa: C plt.yticks(fontsize=0.75 * fsize) plt.tight_layout() if save: + if '/' in figname: + fdir = ''.join(figname.split('/')[:-1]) + check_dir(fdir) plt.savefig(figname, dpi=300) # plt.title('Predicted Covariance: {}'.format((A@initial_cov@A.T).ravel() ))