From 59b667581f26cc55f5c3ff4fb4dffae9406e1087 Mon Sep 17 00:00:00 2001 From: "Michael Pilosov, PhD" <40366263+mathematicalmichael@users.noreply.github.com> Date: Mon, 22 Jan 2024 12:07:18 -0700 Subject: [PATCH] 2024 update (#56) * reflect update in mud * cleanup + updates + linting * rtd * readme update * typo * try supporting missing tex * typo * temporarily remove 3.7 with pip * bring back in 3.7 on 0.1.1rc3 * use latest mud version --- .github/workflows/build.yml | 12 +- .github/workflows/docker.yml | 4 +- .github/workflows/examples.yml | 38 +- .github/workflows/main.yml | 13 +- .github/workflows/publish.yml | 8 +- .readthedocs.yml | 6 +- README.md | 20 + docs/conf.py | 20 +- notebooks/Contours.ipynb | 107 ++-- notebooks/Iteration.ipynb | 182 ++++--- notebooks/MUD-Discretized.ipynb | 887 ++++++++++++++++++++------------ notebooks/MUD_SVD.ipynb | 75 ++- scripts/process_data.py | 4 +- setup.cfg | 2 +- src/mud_examples/linear/lin.py | 8 +- src/mud_examples/monomial.py | 4 +- src/mud_examples/plotting.py | 4 +- 17 files changed, 888 insertions(+), 506 deletions(-) diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index 27d9fbe..e5956c9 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -2,21 +2,21 @@ name: builds on: push: - branches: [ main ] + branches: [main] paths: - "src/**.py" - "setup.py" - "setup.cfg" - ".github/workflows/build.yml" pull_request: - branches-ignore: '**docker**' + branches-ignore: ["**docker**"] paths: - "src/**.py" - "setup.py" - "setup.cfg" - ".github/workflows/build.yml" schedule: - - cron: "0 0 1 * *" + - cron: "0 0 1 * *" jobs: build: @@ -27,12 +27,12 @@ jobs: runs-on: ubuntu-latest steps: - name: Checkout - uses: actions/checkout@v2 + uses: actions/checkout@v3 with: fetch-depth: 0 - name: Set up Python ${{ matrix.python-version }} - uses: actions/setup-python@v1 + uses: actions/setup-python@v2 with: python-version: ${{ matrix.python-version }} @@ -42,7 +42,7 @@ jobs: pip install --upgrade wheel setuptools setuptools_scm - name: Inspect version info - run: | + run: | python setup.py --version git describe --dirty --tags --long --match "*[0-9]*" diff --git a/.github/workflows/docker.yml b/.github/workflows/docker.yml index 6b3f49a..49586f3 100644 --- a/.github/workflows/docker.yml +++ b/.github/workflows/docker.yml @@ -4,7 +4,7 @@ on: release: types: [published, created, edited] push: - branches: '**docker**' + branches: ["**docker**"] jobs: build_image: @@ -12,7 +12,7 @@ jobs: runs-on: ubuntu-latest steps: - name: Check out the repo - uses: actions/checkout@v2 + uses: actions/checkout@v3 with: fetch-depth: 0 diff --git a/.github/workflows/examples.yml b/.github/workflows/examples.yml index 5ea5122..9023f0e 100644 --- a/.github/workflows/examples.yml +++ b/.github/workflows/examples.yml @@ -2,27 +2,27 @@ name: examples on: push: - branches: [ main ] + branches: [main] paths: - "**.py" - ".github/workflows/examples.yml" pull_request: - paths: + paths: - "**.py" - ".github/workflows/examples.yml" schedule: - - cron: "0 0 1 * *" + - cron: "0 0 1 * *" jobs: pip-build: name: Default ${{ matrix.python-version }} strategy: matrix: - python-version: ["3.7", "3.8", "3.9", "3.10"] + python-version: ["3.7", "3.8", "3.9", "3.10", "3.11", "3.12"] runs-on: ubuntu-latest steps: - name: Checkout - uses: actions/checkout@v2 + uses: actions/checkout@v3 - name: Set up Python ${{ matrix.python-version }} uses: actions/setup-python@v2 @@ -32,12 +32,33 @@ jobs: - name: Install Dependencies run: | pip install . + # pip install mud==0.1.1 + pip install -U mud - - name: Test CLI + - name: Test CLI (no TeX) + continue-on-error: true run: | cd /tmp mud_run_all -v - + + - name: Install apt dependencies + run: | + sudo apt-get install -yqq \ + texlive-base \ + texlive-latex-base \ + texlive-latex-extra \ + texlive-fonts-recommended \ + texlive-fonts-extra \ + texlive-science \ + latexmk \ + dvipng \ + cm-super + + - name: Test CLI (w TeX) + run: | + cd /tmp + mud_run_all -v + - name: Generate figures using makefile run: | cd scripts @@ -72,7 +93,8 @@ jobs: - name: Install Dependencies run: | pip install . - pip install mud==0.1rc1 + # pip install mud==0.1.1 + pip install -U mud - name: Conda information run: conda list diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 7f5ea1a..bee492f 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -2,37 +2,37 @@ name: tests on: push: - branches: [ main ] + branches: [main] paths: - "src/**.py" - "setup.py" - "setup.cfg" - ".github/workflows/main.yml" pull_request: - branches-ignore: "**docker**" + branches-ignore: ["**docker**"] paths: - "src/**.py" - "setup.py" - "setup.cfg" - ".github/workflows/main.yml" schedule: - - cron: "0 0 */7 * *" + - cron: "0 0 */7 * *" jobs: test: name: Run unit tests strategy: matrix: - python-version: ["3.7", "3.10"] + python-version: ["3.7", "3.10", "3.12"] runs-on: ubuntu-latest steps: - name: Checkout - uses: actions/checkout@v2 + uses: actions/checkout@v3 with: fetch-depth: 0 - name: Set up Python ${{ matrix.python-version }} - uses: actions/setup-python@v1 + uses: actions/setup-python@v2 with: python-version: ${{ matrix.python-version }} @@ -43,4 +43,3 @@ jobs: - name: Run unit tests run: pytest - diff --git a/.github/workflows/publish.yml b/.github/workflows/publish.yml index 747e49e..09ba472 100644 --- a/.github/workflows/publish.yml +++ b/.github/workflows/publish.yml @@ -4,7 +4,7 @@ on: # release: [published, created, edited] push: tags: - - 'v*' + - "v*" jobs: pypi: @@ -15,12 +15,12 @@ jobs: runs-on: ubuntu-latest steps: - name: Checkout - uses: actions/checkout@v2 + uses: actions/checkout@v3 with: fetch-depth: 0 - name: Set up Python ${{ matrix.python-version }} - uses: actions/setup-python@v1 + uses: actions/setup-python@v2 with: python-version: ${{ matrix.python-version }} @@ -45,7 +45,7 @@ jobs: runs-on: ubuntu-latest steps: - name: Check out the repo - uses: actions/checkout@v2 + uses: actions/checkout@v3 - name: Login to DockerHub uses: docker/login-action@v1 diff --git a/.readthedocs.yml b/.readthedocs.yml index 299280c..3b79cad 100644 --- a/.readthedocs.yml +++ b/.readthedocs.yml @@ -4,6 +4,11 @@ # Required version: 2 +build: + os: "ubuntu-20.04" + tools: + python: "3.8" + # Build documentation in the docs/ directory with Sphinx sphinx: configuration: docs/conf.py @@ -18,7 +23,6 @@ formats: - epub python: - version: 3.8 install: - requirements: docs/requirements.txt - method: pip diff --git a/README.md b/README.md index 4ed1c62..a9cfd53 100644 --- a/README.md +++ b/README.md @@ -11,11 +11,31 @@ Authors: Troy Butler & Michael Pilosov # Installation +For Python 3.7-3.12: ```sh pip install mud-examples ``` +To reproduce the results in Michael's thesis, use `mud-examples==0.1`. However, this comes with `mud==0.0.28`. +Newer versions should still produce the same figures. + +TeX is recommended (but not required): + +``` +apt-get install -yqq \ + texlive-base \ + texlive-latex-base \ + texlive-latex-extra \ + texlive-fonts-recommended \ + texlive-fonts-extra \ + texlive-science \ + latexmk \ + dvipng \ + cm-super +``` + + # Quickstart Generate all of the figures the way they are referenced in the paper: diff --git a/docs/conf.py b/docs/conf.py index 7beab81..a083da0 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -86,8 +86,8 @@ templates_path = ["_templates"] # The suffix of source filenames. -#source_suffix = ".rst" -source_suffix = ['.rst', '.md'] +# source_suffix = ".rst" +source_suffix = [".rst", ".md"] # The encoding of source files. # source_encoding = 'utf-8-sig' @@ -149,9 +149,9 @@ # -- Options for HTML output ------------------------------------------------- -#html_theme = "alabaster" +# html_theme = "alabaster" -#html_theme_options = {"sidebar_width": "300px", "page_width": "1200px"} +# html_theme_options = {"sidebar_width": "300px", "page_width": "1200px"} html_theme = "furo" @@ -205,8 +205,8 @@ html_css_files = [ - "custom.css", - ] + "custom.css", +] # If not '', a 'Last updated on:' timestamp is inserted at every page bottom, # using the given strftime format. @@ -267,7 +267,13 @@ # Grouping the document tree into LaTeX files. List of tuples # (source start file, target name, title, author, documentclass [howto/manual]). latex_documents = [ - ("index", "user_guide.tex", "mud-examples Documentation", "Mathematical Michael", "manual") + ( + "index", + "user_guide.tex", + "mud-examples Documentation", + "Mathematical Michael", + "manual", + ) ] # The name of an image file (relative to this directory) to place at the top of diff --git a/notebooks/Contours.ipynb b/notebooks/Contours.ipynb index 9824ac2..d8b6a35 100644 --- a/notebooks/Contours.ipynb +++ b/notebooks/Contours.ipynb @@ -21,15 +21,16 @@ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import matplotlib\n", + "\n", "if not presentation:\n", - " matplotlib.rcParams['mathtext.fontset'] = 'stix'\n", - " matplotlib.rcParams['font.family'] = 'STIXGeneral'\n", - " fdir = 'contours'\n", + " matplotlib.rcParams[\"mathtext.fontset\"] = \"stix\"\n", + " matplotlib.rcParams[\"font.family\"] = \"STIXGeneral\"\n", + " fdir = \"contours\"\n", "else:\n", - " fdir = '../presentation/figures'\n", + " fdir = \"../presentation/figures\"\n", "\n", - "matplotlib.rcParams['font.size'] = 24\n", - "matplotlib.backend = 'Agg'" + "matplotlib.rcParams[\"font.size\"] = 24\n", + "matplotlib.backend = \"Agg\"" ] }, { @@ -38,7 +39,7 @@ "metadata": {}, "outputs": [], "source": [ - "plt.rcParams['figure.figsize'] = 10,10" + "plt.rcParams[\"figure.figsize\"] = 10, 10" ] }, { @@ -50,7 +51,7 @@ "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))" + "b = np.zeros((1, 1))" ] }, { @@ -74,11 +75,27 @@ "outputs": [], "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", + "cov_11_slider = wd.FloatSlider(\n", + " value=1,\n", + " min=0.05,\n", + " max=1,\n", + " step=0.05,\n", + " description=\"$\\sigma_{11}$\",\n", + " continuous_update=False,\n", + ")\n", + "cov_01_slider = wd.FloatSlider(\n", + " value=0,\n", + " min=-0.95 - 1e-4,\n", + " max=0.95 + 1e-4,\n", + " step=0.01,\n", + " description=\"$\\sigma_{01}=\\sigma_{10}$\",\n", + " continuous_update=False,\n", + ")\n", + "\n", "\n", "def mv_slider(val):\n", - " pr_slide.max = tk_slide.value + 1E-4\n", + " pr_slide.max = tk_slide.value + 1e-4\n", + "\n", "\n", "def ensure_positive_determinant(val):\n", " \"\"\"\n", @@ -86,37 +103,51 @@ " \"\"\"\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", + " cov_01_slider.max = rt - 1e-4\n", + " cov_01_slider.min = -rt + 1e-4\n", + "\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='Tikhonov', 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 = wd.FloatSlider(\n", + " value=1, min=0, max=1, step=0.02, description=\"Tikhonov\", continuous_update=False\n", + ")\n", + "pr_slide = wd.FloatSlider(\n", + " value=1, min=0, max=1, step=0.02, description=\"Unreg\", continuous_update=False\n", + ")\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", + "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", + "fig_title = wd.Text(value=\"latest_figure.png\")\n", + "obs_std_slider = wd.FloatSlider(\n", + " value=1, min=0.01, max=1, step=0.01, description=\"Obs. Std\", continuous_update=False\n", + ")\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(contour_example, 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", + "I = wd.interactive(\n", + " contour_example,\n", + " A=wd.fixed(A),\n", + " b=wd.fixed(b),\n", + " save=wd.fixed(save),\n", + " param_ref=wd.fixed(lam_true),\n", + " compare=comparison,\n", + " cov_01=cov_01_slider,\n", + " 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", "\n", "display(wd.VBox([I, out]))" ] @@ -151,7 +182,7 @@ "metadata": {}, "outputs": [], "source": [ - "fig_title.value = f'{fdir}/data_mismatch_contour.png'\n", + "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", @@ -171,7 +202,7 @@ "metadata": {}, "outputs": [], "source": [ - "fig_title.value = f'{fdir}/tikonov_contour.png'\n", + "fig_title.value = f\"{fdir}/tikonov_contour.png\"\n", "tk_slide.value = 1\n", "pr_slide.value = 0\n", "data_check.value = False\n", @@ -191,7 +222,7 @@ "metadata": {}, "outputs": [], "source": [ - "fig_title.value = f'{fdir}/consistent_contour.png'\n", + "fig_title.value = f\"{fdir}/consistent_contour.png\"\n", "tk_slide.value = 1\n", "pr_slide.value = 1\n", "data_check.value = False\n", @@ -211,7 +242,7 @@ "metadata": {}, "outputs": [], "source": [ - "fig_title.value = f'{fdir}/classical_solution.png'\n", + "fig_title.value = f\"{fdir}/classical_solution.png\"\n", "tk_slide.value = 1\n", "pr_slide.value = 0\n", "data_check.value = True\n", @@ -231,7 +262,7 @@ "metadata": {}, "outputs": [], "source": [ - "fig_title.value = f'{fdir}/consistent_solution.png'\n", + "fig_title.value = f\"{fdir}/consistent_solution.png\"\n", "tk_slide.value = 1\n", "pr_slide.value = 1\n", "data_check.value = True\n", @@ -251,7 +282,7 @@ "metadata": {}, "outputs": [], "source": [ - "fig_title.value = f'{fdir}/map_compare_contour.png'\n", + "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", diff --git a/notebooks/Iteration.ipynb b/notebooks/Iteration.ipynb index d1f5287..68f5b67 100644 --- a/notebooks/Iteration.ipynb +++ b/notebooks/Iteration.ipynb @@ -18,14 +18,15 @@ "import numpy as np\n", "from matplotlib import pyplot as plt\n", "import matplotlib\n", + "\n", "if not presentation:\n", - " matplotlib.rcParams['mathtext.fontset'] = 'stix'\n", - " matplotlib.rcParams['font.family'] = 'STIXGeneral'\n", - "matplotlib.rcParams['font.size'] = 24\n", - "matplotlib.rcParams['figure.figsize'] = 10, 10\n", + " matplotlib.rcParams[\"mathtext.fontset\"] = \"stix\"\n", + " matplotlib.rcParams[\"font.family\"] = \"STIXGeneral\"\n", + "matplotlib.rcParams[\"font.size\"] = 24\n", + "matplotlib.rcParams[\"figure.figsize\"] = 10, 10\n", "# matplotlib.rcParams['xtick.bottom'] = False\n", "# matplotlib.rcParams['ytick.left'] = False\n", - "matplotlib.backend = 'Agg'" + "matplotlib.backend = \"Agg\"" ] }, { @@ -58,14 +59,14 @@ "np.random.seed(24)\n", "numQoI = 10\n", "A = utils.rotationMap(numQoI, orth=False)\n", - "A = A.reshape(-1,2)\n", - "b = np.zeros((numQoI,1))\n", + "A = A.reshape(-1, 2)\n", + "b = np.zeros((numQoI, 1))\n", "# A, b = A[np.array([2,6]),:].reshape(-1,2), b[np.array([2,6])].reshape(-1,1)\n", - "ref_param = np.array([0.5, 0.5]).reshape(-1,1)\n", - "y = A@ref_param + b\n", - "initial_mean = np.random.randn(2).reshape(-1,1)\n", + "ref_param = np.array([0.5, 0.5]).reshape(-1, 1)\n", + "y = A @ ref_param + b\n", + "initial_mean = np.random.randn(2).reshape(-1, 1)\n", "tol = 0.1\n", - "initial_cov = np.eye(2)*utils.std_from_equipment(tol)" + "initial_cov = np.eye(2) * utils.std_from_equipment(tol)" ] }, { @@ -84,7 +85,7 @@ "outputs": [], "source": [ "def applystyle(plotobj):\n", - " plotobj.axis('equal')\n", + " plotobj.axis(\"equal\")\n", " plotobj.ylim([-0.8, 1.5])\n", " plotobj.xlim([-0.5, 1.5])\n", " # plt.xticks([])\n", @@ -97,14 +98,16 @@ "metadata": {}, "outputs": [], "source": [ - "plt.figure(figsize=(10,10))\n", - "mud_chain = funs.iterate(A, b, y, initial_mean, initial_cov, data_cov=None, num_epochs=1, idx=None)\n", + "plt.figure(figsize=(10, 10))\n", + "mud_chain = funs.iterate(\n", + " A, b, y, initial_mean, initial_cov, data_cov=None, num_epochs=1, idx=None\n", + ")\n", "plot.plotChain(mud_chain, ref_param)\n", "plot.plot_contours(A, ref_param)\n", "\n", "applystyle(plt)\n", "\n", - "plt.savefig(f'iterative/{numQoI}D-firstepoch.png')\n", + "plt.savefig(f\"iterative/{numQoI}D-firstepoch.png\")\n", "plt.show()" ] }, @@ -114,16 +117,20 @@ "metadata": {}, "outputs": [], "source": [ - "plt.figure(figsize=(10,10))\n", - "mud_chain = funs.iterate(A, b, y, initial_mean, initial_cov, data_cov=None, num_epochs=1)\n", + "plt.figure(figsize=(10, 10))\n", + "mud_chain = funs.iterate(\n", + " A, b, y, initial_mean, initial_cov, data_cov=None, num_epochs=1\n", + ")\n", "plot.plotChain(mud_chain, ref_param)\n", "for _ in range(3):\n", - " mud_chain = funs.iterate(A, b, y, mud_chain[-1], initial_cov, data_cov=None, num_epochs=1)\n", + " mud_chain = funs.iterate(\n", + " A, b, y, mud_chain[-1], initial_cov, data_cov=None, num_epochs=1\n", + " )\n", " plot.plotChain(mud_chain, ref_param)\n", - " \n", + "\n", "plot.plot_contours(A, ref_param)\n", "applystyle(plt)\n", - "plt.savefig(f'iterative/{numQoI}D-fewepochs.png')\n", + "plt.savefig(f\"iterative/{numQoI}D-fewepochs.png\")\n", "plt.show()" ] }, @@ -141,25 +148,29 @@ "outputs": [], "source": [ "np.random.seed(21)\n", - "plt.figure(figsize=(10,10))\n", + "plt.figure(figsize=(10, 10))\n", "idx_choice = [np.random.randint(numQoI), np.random.randint(numQoI)]\n", - "mud_chain = funs.iterate(A, b, y, initial_mean, initial_cov, data_cov=None, num_epochs=5, idx=idx_choice)\n", - "plot.plot_contours(A, ref_param, subset=idx_choice, color='k')\n", - "plot.plotChain(mud_chain,ref_param)\n", + "mud_chain = funs.iterate(\n", + " A, b, y, initial_mean, initial_cov, data_cov=None, num_epochs=5, idx=idx_choice\n", + ")\n", + "plot.plot_contours(A, ref_param, subset=idx_choice, color=\"k\")\n", + "plot.plotChain(mud_chain, ref_param)\n", "\n", "applystyle(plt)\n", - "plt.savefig(f'iterative/{numQoI}D-fewepochs-pair.png')\n", + "plt.savefig(f\"iterative/{numQoI}D-fewepochs-pair.png\")\n", "plt.show()\n", "\n", "np.random.seed(19)\n", "idx_choice = [np.random.randint(numQoI), np.random.randint(numQoI)]\n", - "plt.figure(figsize=(10,10))\n", - "mud_chain = funs.iterate(A, b, y, initial_mean, initial_cov, data_cov=None, num_epochs=5, idx=idx_choice)\n", - "plot.plot_contours(A, ref_param, color='b', subset=idx_choice)\n", + "plt.figure(figsize=(10, 10))\n", + "mud_chain = funs.iterate(\n", + " A, b, y, initial_mean, initial_cov, data_cov=None, num_epochs=5, idx=idx_choice\n", + ")\n", + "plot.plot_contours(A, ref_param, color=\"b\", subset=idx_choice)\n", "plot.plotChain(mud_chain, ref_param)\n", "\n", "applystyle(plt)\n", - "plt.savefig(f'iterative/{numQoI}D-fewepochs-pair-alt.png')\n", + "plt.savefig(f\"iterative/{numQoI}D-fewepochs-pair-alt.png\")\n", "plt.show()" ] }, @@ -176,31 +187,37 @@ "metadata": {}, "outputs": [], "source": [ - "plt.figure(figsize=(10,10))\n", - "idx_choice = [3,8]\n", - "mud_chain = funs.iterate(A, b, y, initial_mean, initial_cov, data_cov=None, num_epochs=5, idx=idx_choice)\n", - "plot.plot_contours(A, ref_param, color='k', subset=idx_choice)\n", + "plt.figure(figsize=(10, 10))\n", + "idx_choice = [3, 8]\n", + "mud_chain = funs.iterate(\n", + " A, b, y, initial_mean, initial_cov, data_cov=None, num_epochs=5, idx=idx_choice\n", + ")\n", + "plot.plot_contours(A, ref_param, color=\"k\", subset=idx_choice)\n", "plot.plotChain(mud_chain, ref_param)\n", "\n", - "idx_choice = [1,6]\n", - "mud_chain = funs.iterate(A, b, y, initial_mean, initial_cov, data_cov=None, num_epochs=5, idx=idx_choice)\n", - "plot.plot_contours(A, ref_param, color='b', subset=idx_choice)\n", - "plot.plotChain(mud_chain, ref_param, 'b')\n", + "idx_choice = [1, 6]\n", + "mud_chain = funs.iterate(\n", + " A, b, y, initial_mean, initial_cov, data_cov=None, num_epochs=5, idx=idx_choice\n", + ")\n", + "plot.plot_contours(A, ref_param, color=\"b\", subset=idx_choice)\n", + "plot.plotChain(mud_chain, ref_param, \"b\")\n", "\n", "applystyle(plt)\n", - "plt.savefig(f'iterative/{numQoI}D-firstepoch-pair-smart.png')\n", + "plt.savefig(f\"iterative/{numQoI}D-firstepoch-pair-smart.png\")\n", "plt.show()\n", "\n", "np.random.seed(12)\n", - "plt.figure(figsize=(10,10))\n", + "plt.figure(figsize=(10, 10))\n", "idx_choice = np.arange(numQoI)\n", "np.random.shuffle(idx_choice)\n", - "mud_chain = funs.iterate(A, b, y, initial_mean, initial_cov, data_cov=None, num_epochs=1, idx=idx_choice)\n", - "plot.plot_contours(A[np.array(idx_choice),:], ref_param)\n", + "mud_chain = funs.iterate(\n", + " A, b, y, initial_mean, initial_cov, data_cov=None, num_epochs=1, idx=idx_choice\n", + ")\n", + "plot.plot_contours(A[np.array(idx_choice), :], ref_param)\n", "plot.plotChain(mud_chain, ref_param)\n", "\n", "applystyle(plt)\n", - "plt.savefig(f'iterative/{numQoI}D-firstepoch-rand.png')\n", + "plt.savefig(f\"iterative/{numQoI}D-firstepoch-rand.png\")\n", "plt.show()" ] }, @@ -223,13 +240,22 @@ "source": [ "num_trials = 100\n", "model_eval_budget = 100\n", - "num_epochs = model_eval_budget//numQoI\n", + "num_epochs = model_eval_budget // numQoI\n", "# ordered\n", "idx_choice = np.arange(numQoI)\n", "mud_chains_from_ordered = []\n", "for trial in range(num_trials):\n", - " initial_mean = np.random.rand(2,1)\n", - " mud_chain = funs.iterate(A, b, y, initial_mean, initial_cov, data_cov=None, num_epochs=num_epochs, idx=idx_choice)\n", + " initial_mean = np.random.rand(2, 1)\n", + " mud_chain = funs.iterate(\n", + " A,\n", + " b,\n", + " y,\n", + " initial_mean,\n", + " initial_cov,\n", + " data_cov=None,\n", + " num_epochs=num_epochs,\n", + " idx=idx_choice,\n", + " )\n", " mud_chains_from_ordered.append(mud_chain)\n", "\n", "# some random ordering which will be used over multiple epochs\n", @@ -237,29 +263,42 @@ "idx_choice = np.arange(numQoI)\n", "mud_chains_from_shuffles = []\n", "for trial in range(num_trials):\n", - " initial_mean = np.random.rand(2,1)\n", + " initial_mean = np.random.rand(2, 1)\n", " np.random.shuffle(idx_choice)\n", - " mud_chain = funs.iterate(A, b, y, initial_mean, initial_cov, data_cov=None, num_epochs=num_epochs, idx=idx_choice)\n", + " mud_chain = funs.iterate(\n", + " A,\n", + " b,\n", + " y,\n", + " initial_mean,\n", + " initial_cov,\n", + " data_cov=None,\n", + " num_epochs=num_epochs,\n", + " idx=idx_choice,\n", + " )\n", " mud_chains_from_shuffles.append(mud_chain)\n", "\n", "# uses all the same data same number of times, but only one epoch due to\n", "# the indices being repeated in a list before being shuffled\n", - "idx_choice = list(np.arange(numQoI))*num_epochs\n", + "idx_choice = list(np.arange(numQoI)) * num_epochs\n", "mud_chains_from_batch = []\n", "for trial in range(num_trials):\n", - " initial_mean = np.random.rand(2,1)\n", + " initial_mean = np.random.rand(2, 1)\n", " np.random.shuffle(idx_choice)\n", - " mud_chain = funs.iterate(A, b, y, initial_mean, initial_cov, data_cov=None, num_epochs=1, idx=idx_choice)\n", + " mud_chain = funs.iterate(\n", + " A, b, y, initial_mean, initial_cov, data_cov=None, num_epochs=1, idx=idx_choice\n", + " )\n", " mud_chains_from_batch.append(mud_chain)\n", "\n", "# no imposed limitation on using data equal number of times by the end\n", - "idx_choice = np.random.choice(np.arange(numQoI), size=num_epochs*numQoI)\n", + "idx_choice = np.random.choice(np.arange(numQoI), size=num_epochs * numQoI)\n", "mud_chains_from_randoms = []\n", "for trial in range(num_trials):\n", - " initial_mean = np.random.rand(2,1)\n", + " initial_mean = np.random.rand(2, 1)\n", " np.random.shuffle(idx_choice)\n", - " mud_chain = funs.iterate(A, b, y, initial_mean, initial_cov, data_cov=None, num_epochs=1, idx=idx_choice)\n", - " mud_chains_from_randoms.append(mud_chain)\n" + " mud_chain = funs.iterate(\n", + " A, b, y, initial_mean, initial_cov, data_cov=None, num_epochs=1, idx=idx_choice\n", + " )\n", + " mud_chains_from_randoms.append(mud_chain)" ] }, { @@ -268,43 +307,54 @@ "metadata": {}, "outputs": [], "source": [ - "plt.figure(figsize=(20,10))\n", + "plt.figure(figsize=(20, 10))\n", "error = []\n", "for mud_chain in mud_chains_from_ordered:\n", " approx_error = [np.linalg.norm(m - ref_param) for m in mud_chain]\n", " error.append(approx_error)\n", - " plt.plot(approx_error, 'r', alpha=0.1)\n", - "plt.plot(np.mean(np.array(error),axis=0), 'r', lw=5, label='Ordered QoI $(10\\\\times 10D)$)')\n", + " plt.plot(approx_error, \"r\", alpha=0.1)\n", + "plt.plot(\n", + " np.mean(np.array(error), axis=0), \"r\", lw=5, label=\"Ordered QoI $(10\\\\times 10D)$)\"\n", + ")\n", "\n", "error = []\n", "for mud_chain in mud_chains_from_shuffles:\n", " approx_error = [np.linalg.norm(m - ref_param) for m in mud_chain]\n", " error.append(approx_error)\n", - " plt.plot(approx_error, 'k', alpha=0.1)\n", - "plt.plot(np.mean(np.array(error),axis=0), 'k', lw=5, label='Shuffled QoI $(10\\\\times 10D)$')\n", + " plt.plot(approx_error, \"k\", alpha=0.1)\n", + "plt.plot(\n", + " np.mean(np.array(error), axis=0), \"k\", lw=5, label=\"Shuffled QoI $(10\\\\times 10D)$\"\n", + ")\n", "\n", "error = []\n", "for mud_chain in mud_chains_from_batch:\n", " approx_error = [np.linalg.norm(m - ref_param) for m in mud_chain]\n", " error.append(approx_error)\n", - " plt.plot(approx_error, 'b', alpha=0.1)\n", - "plt.plot(np.mean(np.array(error),axis=0), 'b', lw=5, label='Batch QoI $(1\\\\times 100D)$')\n", + " plt.plot(approx_error, \"b\", alpha=0.1)\n", + "plt.plot(\n", + " np.mean(np.array(error), axis=0), \"b\", lw=5, label=\"Batch QoI $(1\\\\times 100D)$\"\n", + ")\n", "\n", "error = []\n", "for mud_chain in mud_chains_from_randoms:\n", " approx_error = [np.linalg.norm(m - ref_param) for m in mud_chain]\n", " error.append(approx_error)\n", "# plt.plot(approx_error, 'orange', alpha=0.1)\n", - "plt.plot(np.mean(np.array(error),axis=0), 'orange', lw=5, label='Random QoI $(1\\\\times 100D)$')\n", + "plt.plot(\n", + " np.mean(np.array(error), axis=0),\n", + " \"orange\",\n", + " lw=5,\n", + " label=\"Random QoI $(1\\\\times 100D)$\",\n", + ")\n", "\n", "\n", - "plt.legend(loc='lower left')\n", - "plt.yscale('log')\n", - "plt.ylim(1E-16, 9E-1)\n", + "plt.legend(loc=\"lower left\")\n", + "plt.yscale(\"log\")\n", + "plt.ylim(1e-16, 9e-1)\n", "# plt.xscale('log')\n", "plt.ylabel(\"$||\\lambda - \\lambda^\\dagger||$\", fontsize=36)\n", "plt.xlabel(\"Iteration Step\", fontsize=32)\n", - "plt.savefig(f'iterative/{numQoI}D-convergence-comparison.png')\n", + "plt.savefig(f\"iterative/{numQoI}D-convergence-comparison.png\")\n", "plt.show()" ] }, diff --git a/notebooks/MUD-Discretized.ipynb b/notebooks/MUD-Discretized.ipynb index 02fb176..3394900 100644 --- a/notebooks/MUD-Discretized.ipynb +++ b/notebooks/MUD-Discretized.ipynb @@ -63,15 +63,16 @@ "from mud_examples.plotting import plot_experiment_measurements\n", "from mud_examples.helpers import maybe_fit_log_linear_regression\n", "import matplotlib\n", + "\n", "if not presentation:\n", - " matplotlib.rcParams['mathtext.fontset'] = 'stix'\n", - " matplotlib.rcParams['font.family'] = 'STIXGeneral'\n", - " fdir = 'pde_2D'\n", + " matplotlib.rcParams[\"mathtext.fontset\"] = \"stix\"\n", + " matplotlib.rcParams[\"font.family\"] = \"STIXGeneral\"\n", + " fdir = \"pde_2D\"\n", "else:\n", - "# fdir = '../presentation/figures/pde-highd'\n", - " fdir = 'pde_2D'\n", - "matplotlib.rcParams['font.size'] = 24\n", - "matplotlib.backend = 'Agg'" + " # fdir = '../presentation/figures/pde-highd'\n", + " fdir = \"pde_2D\"\n", + "matplotlib.rcParams[\"font.size\"] = 24\n", + "matplotlib.backend = \"Agg\"" ] }, { @@ -85,7 +86,11 @@ }, "outputs": [], "source": [ - "from poisson import poissonModel, poisson_sensor_model, eval_boundary_piecewise as pcwExpr\n", + "from poisson import (\n", + " poissonModel,\n", + " poisson_sensor_model,\n", + " eval_boundary_piecewise as pcwExpr,\n", + ")\n", "from poisson import eval_boundary, gamma_boundary_condition\n", "import pickle" ] @@ -102,6 +107,7 @@ "outputs": [], "source": [ "from mud_examples.models import generate_spatial_measurements as generate_sensors_pde\n", + "\n", "# from mud_examples.datasets import load_poisson" ] }, @@ -167,17 +173,17 @@ }, "outputs": [], "source": [ - "prefix = f'pde_2D/pde'\n", - "num_measure = 1000 # number of measurement (sensor) locations\n", - "fsize = 36\n", - "num_trials = 20 # realizations of synthetic data for numerical runs\n", - "tolerance = 0.1 # precision of measurement equipment\n", - "sigma = std_from_equipment(tolerance=tolerance, probability=0.99)\n", + "prefix = f\"pde_2D/pde\"\n", + "num_measure = 1000 # number of measurement (sensor) locations\n", + "fsize = 36\n", + "num_trials = 20 # realizations of synthetic data for numerical runs\n", + "tolerance = 0.1 # precision of measurement equipment\n", + "sigma = std_from_equipment(tolerance=tolerance, probability=0.99)\n", "np.random.seed(21)\n", - "lam_true = 3.0\n", - "input_dim = 2\n", + "lam_true = 3.0\n", + "input_dim = 2\n", "num_samples = 100\n", - "ftype = 'png'" + "ftype = \"png\"" ] }, { @@ -191,7 +197,7 @@ }, "outputs": [], "source": [ - "load = False\n" + "load = False" ] }, { @@ -228,7 +234,7 @@ "metadata": {}, "outputs": [], "source": [ - "fname = 'pde_2D/pde_ref.pkl'" + "fname = \"pde_2D/pde_ref.pkl\"" ] }, { @@ -244,7 +250,7 @@ "metadata": {}, "outputs": [], "source": [ - "p.plot_without_fenics(fname, 100, mode='hor', num_qoi=2)" + "p.plot_without_fenics(fname, 100, mode=\"hor\", num_qoi=2)" ] }, { @@ -272,16 +278,16 @@ "outputs": [], "source": [ "def load_poisson_from_disk(fname):\n", - " ref = pickle.load(open(fname, 'rb'))\n", - " lam = ref['lam']\n", + " ref = pickle.load(open(fname, \"rb\"))\n", + " lam = ref[\"lam\"]\n", " input_dim = lam.shape[1]\n", - " domain = np.array([[-4,0]]*input_dim)\n", - " qoi = ref['qoi']\n", - " qoi_true = ref['data']\n", - " lam_ref = ref['truth']\n", - " u = ref['plot_u']\n", - " g = ref['plot_g']\n", - " sensors = ref['sensors']\n", + " domain = np.array([[-4, 0]] * input_dim)\n", + " qoi = ref[\"qoi\"]\n", + " qoi_true = ref[\"data\"]\n", + " lam_ref = ref[\"truth\"]\n", + " u = ref[\"plot_u\"]\n", + " g = ref[\"plot_g\"]\n", + " sensors = ref[\"sensors\"]\n", " return domain, lam, qoi, qoi_true, lam_ref, u, g, sensors" ] }, @@ -352,7 +358,7 @@ "metadata": {}, "outputs": [], "source": [ - "qoi_2d_hor?" + "?qoi_2d_hor" ] }, { @@ -379,7 +385,7 @@ "metadata": {}, "outputs": [], "source": [ - "qoi_1d??" + "??qoi_1d" ] }, { @@ -388,7 +394,7 @@ "metadata": {}, "outputs": [], "source": [ - "qoi_1d(20,0.05)" + "qoi_1d(20, 0.05)" ] }, { @@ -434,7 +440,7 @@ "metadata": {}, "outputs": [], "source": [ - "plt.tricontourf(lam[:,0], lam[:,1], m._up)" + "plt.tricontourf(lam[:, 0], lam[:, 1], m._up)" ] }, { @@ -455,12 +461,13 @@ "from scipy.stats import gaussian_kde as gkde\n", "from scipy.stats import distributions as dist\n", "\n", + "\n", "def ratio_dci_sing(qoi):\n", - "# if qoi.ndim == 2:\n", - "# qoi = qoi.ravel()\n", + " # if qoi.ndim == 2:\n", + " # qoi = qoi.ravel()\n", " kde = gkde(qoi.T)\n", - " ratio_eval = dist.norm.pdf(qoi)/kde.pdf(qoi.T).ravel()\n", - " return ratio_eval\n" + " ratio_eval = dist.norm.pdf(qoi) / kde.pdf(qoi.T).ravel()\n", + " return ratio_eval" ] }, { @@ -470,6 +477,7 @@ "outputs": [], "source": [ "import mud.base as mb\n", + "\n", "importlib.reload(mb)" ] }, @@ -480,6 +488,7 @@ "outputs": [], "source": [ "import mud.funs as mf\n", + "\n", "importlib.reload(mf)" ] }, @@ -529,21 +538,21 @@ "q = wme(qoi[:, 0:num_obs], data, sd)\n", "r = ratio_dci_sing(q)\n", "print(r.shape)\n", - "plt.tricontourf(lam[:,0], lam[:,1], r)\n", + "plt.tricontourf(lam[:, 0], lam[:, 1], r)\n", "plt.show()\n", "\n", "d = mb.DensityProblem(lam, q, domain)\n", "d.estimate()\n", - "plt.tricontourf(lam[:,0], lam[:,1], d._ob / d._pr)\n", + "plt.tricontourf(lam[:, 0], lam[:, 1], d._ob / d._pr)\n", "plt.show()\n", "\n", "# up = d._up\n", "up = np.multiply(d._in, np.divide(d._ob, d._pr))\n", - "plt.tricontourf(lam[:,0], lam[:,1], up)\n", + "plt.tricontourf(lam[:, 0], lam[:, 1], up)\n", "plt.show()\n", "\n", "up = d._in\n", - "plt.tricontourf(lam[:,0], lam[:,1], up)\n", + "plt.tricontourf(lam[:, 0], lam[:, 1], up)\n", "plt.show()" ] }, @@ -562,7 +571,7 @@ "metadata": {}, "outputs": [], "source": [ - "lam[np.argmax(r),:]" + "lam[np.argmax(r), :]" ] }, { @@ -651,7 +660,7 @@ "metadata": {}, "outputs": [], "source": [ - "d = p.make_mud_wrapper(domain, lam, qoi, qoi_true)(100,0.05)" + "d = p.make_mud_wrapper(domain, lam, qoi, qoi_true)(100, 0.05)" ] }, { @@ -682,11 +691,9 @@ "\n", "\n", "# SCALAR\n", - "experiments_sing, solutions_sing = experiment_measurements(num_measurements=measurements,\n", - " sd=sigma,\n", - " num_trials=num_trials,\n", - " seed=21,\n", - " fun=wrap)\n", + "experiments_sing, solutions_sing = experiment_measurements(\n", + " num_measurements=measurements, sd=sigma, num_trials=num_trials, seed=21, fun=wrap\n", + ")\n", "# # VECTOR\n", "# def mud_wrapper(num_obs, sd):\n", "# qois = split_qoi_by_indices(qoi_indices, qoi_ref, qoi,\n", @@ -751,7 +758,7 @@ "metadata": {}, "outputs": [], "source": [ - "plot_samples(lam, qoi, qoi_ref, sols = solutions_sing, num=100, save=False)" + "plot_samples(lam, qoi, qoi_ref, sols=solutions_sing, num=100, save=False)" ] }, { @@ -761,34 +768,37 @@ "outputs": [], "source": [ "for num_plot_sensors in [20, 100]:\n", - " plt.figure(figsize=(10,10))\n", + " plt.figure(figsize=(10, 10))\n", "\n", - "# plt.subplot(111)\n", - " plt.title('MUD Estimates for $Q_{1D}$,' + f' N={num_plot_sensors}', fontsize=1.25*fsize)\n", + " # plt.subplot(111)\n", + " plt.title(\n", + " \"MUD Estimates for $Q_{1D}$,\" + f\" N={num_plot_sensors}\", fontsize=1.25 * fsize\n", + " )\n", " plt.xlabel(\"$x_2$\", fontsize=fsize)\n", " plt.ylabel(\"$g(x, \\lambda)$\", fontsize=fsize)\n", "\n", - " for _lam in solutions_sing[num_plot_sensors]: # trials\n", - " plt.plot([0]+intervals+[1], [0]+list(_lam)+[0], lw=1, c='purple', alpha=0.2)\n", - "\n", + " for _lam in solutions_sing[num_plot_sensors]: # trials\n", + " plt.plot(\n", + " [0] + intervals + [1], [0] + list(_lam) + [0], lw=1, c=\"purple\", alpha=0.2\n", + " )\n", "\n", " plt.ylim(-4, 0)\n", " plt.xlim(0, 1)\n", - " plt.legend(loc='lower left')\n", - "# plt.subplot(122)\n", - "# plt.title('Q-Q Plot', fontsize=1.25*fsize)\n", - "# for i in solutions_sing[num_plot_sensors]: # trials\n", - "# q = qoi[i,:]\n", - "# plt.scatter(q[:num_plot_sensors], qoi_ref[:num_plot_sensors] + noise[:num_plot_sensors],\n", - "# c='b', s=100, alpha=1.0/num_trials)\n", - "# plt.plot(_a,_a, c='k', lw=3)\n", - "# plt.ylabel('Collected Data', fontsize=fsize)\n", - "# plt.xlabel('Predicted Data', fontsize=fsize)\n", - "# plt.ylim(-0.5, 0.15)\n", - "# plt.xlim(-0.5, 0.15)\n", - "# plt.title(f'Solution {_r}, Index {_s}')\n", + " plt.legend(loc=\"lower left\")\n", + " # plt.subplot(122)\n", + " # plt.title('Q-Q Plot', fontsize=1.25*fsize)\n", + " # for i in solutions_sing[num_plot_sensors]: # trials\n", + " # q = qoi[i,:]\n", + " # plt.scatter(q[:num_plot_sensors], qoi_ref[:num_plot_sensors] + noise[:num_plot_sensors],\n", + " # c='b', s=100, alpha=1.0/num_trials)\n", + " # plt.plot(_a,_a, c='k', lw=3)\n", + " # plt.ylabel('Collected Data', fontsize=fsize)\n", + " # plt.xlabel('Predicted Data', fontsize=fsize)\n", + " # plt.ylim(-0.5, 0.15)\n", + " # plt.xlim(-0.5, 0.15)\n", + " # plt.title(f'Solution {_r}, Index {_s}')\n", " _fname = f\"{prefix}_mud_{input_dim}-1_N{num_plot_sensors}.{ftype}\"\n", - " plt.savefig(_fname, bbox_inches='tight')\n", + " plt.savefig(_fname, bbox_inches=\"tight\")\n", " plt.show()" ] }, @@ -866,20 +876,20 @@ "source": [ "# %%time\n", "if load:\n", - " fname = f'./{prefix}_summary_{input_dim}.pkl'\n", - " results = pickle.load(open(fname, 'rb'))\n", - " solutions_sing, solutions_mult = results['sols']\n", - " measurements = results['meas']\n", - " noise, tolerance = results['noise']\n", - " sigma = results['stdv']\n", - " lam, qoi = results['sets']\n", - " lam_ref, qoi_ref = results['true']\n", - " sensors = results['sens']\n", + " fname = f\"./{prefix}_summary_{input_dim}.pkl\"\n", + " results = pickle.load(open(fname, \"rb\"))\n", + " solutions_sing, solutions_mult = results[\"sols\"]\n", + " measurements = results[\"meas\"]\n", + " noise, tolerance = results[\"noise\"]\n", + " sigma = results[\"stdv\"]\n", + " lam, qoi = results[\"sets\"]\n", + " lam_ref, qoi_ref = results[\"true\"]\n", + " sensors = results[\"sens\"]\n", " model_list = None\n", - " \n", + "\n", " pde = pdeProblem(lam, qoi, lam_ref, qoi_ref, sensors)\n", "else:\n", - " model_list = pickle.load(open(f'res{input_dim}u.pkl', 'rb'))\n", + " model_list = pickle.load(open(f\"res{input_dim}u.pkl\", \"rb\"))\n", " sensors = generate_sensors_pde(num_measure)\n", " lam, qoi = load_poisson(sensors, model_list[0:num_samples], nx=36, ny=36)\n", " qoi_ref = poisson_sensor_model(sensors, gamma=lam_true, nx=36, ny=36)\n", @@ -899,7 +909,7 @@ }, "outputs": [], "source": [ - "noise = sigma*np.random.randn(num_measure)" + "noise = sigma * np.random.randn(num_measure)" ] }, { @@ -925,7 +935,7 @@ "metadata": {}, "outputs": [], "source": [ - "svals = [p(xi,yi) for xi,yi in sensors]" + "svals = [p(xi, yi) for xi, yi in sensors]" ] }, { @@ -935,9 +945,9 @@ "outputs": [], "source": [ "z = qoi_ref\n", - "x, y = sensors[:,0], sensors[:,1]\n", + "x, y = sensors[:, 0], sensors[:, 1]\n", "plt.tricontourf(x, y, z, levels=20, vmin=-0.5, vmax=0)\n", - "plt.scatter(x, y, c='r', s=0.5)\n", + "plt.scatter(x, y, c=\"r\", s=0.5)\n", "plt.show()" ] }, @@ -1008,7 +1018,7 @@ "source": [ "fin.plot(pn, vmin=-0.5, vmax=0)\n", "plt.title(f\"Response Surface\\n$\\\\sigma$ = {sigma:1.3E} ($\\\\tau$ = {tolerance:1.1E})\")\n", - "plt.scatter(sensors[0:100,0], sensors[0:100,1], s=100, c='k')\n", + "plt.scatter(sensors[0:100, 0], sensors[0:100, 1], s=100, c=\"k\")\n", "plt.show()" ] }, @@ -1025,7 +1035,7 @@ "source": [ "w = fin.Expression(pcwExpr(u, input_dim, d=0), degree=2)\n", "u_plot = fin.Expression(pcwExpr(u, 1000, d=0), degree=2)\n", - "domain = np.array([[-4,0]*input_dim]).reshape(-1,2)" + "domain = np.array([[-4, 0] * input_dim]).reshape(-1, 2)" ] }, { @@ -1063,37 +1073,65 @@ }, "outputs": [], "source": [ - "plt.figure(figsize=(20,10))\n", + "plt.figure(figsize=(20, 10))\n", "plt.subplot(121)\n", - "fin.plot(w, mesh=mesh, lw=5, c='k')\n", + "fin.plot(w, mesh=mesh, lw=5, c=\"k\")\n", "gt = list(lam[closest_fit_index_in, :])\n", - "plt.plot([0]+intervals+[1], [0]+gt+[0], lw=5, c='purple', alpha=1, label=f'Closest in Input: {closest_fit_index_in}')\n", + "plt.plot(\n", + " [0] + intervals + [1],\n", + " [0] + gt + [0],\n", + " lw=5,\n", + " c=\"purple\",\n", + " alpha=1,\n", + " label=f\"Closest in Input: {closest_fit_index_in}\",\n", + ")\n", "\n", "projected_line = list(lam[closest_fit_index_out, :])\n", - "plt.plot([0]+intervals+[1], [0]+projected_line+[0], lw=5, c='green', alpha=1, label=f'Closest in Output: {closest_fit_index_out}')\n", - "\n", - "plt.legend(fontsize=fsize*0.75)\n", - "plt.title(f'Parameter Space', fontsize=fsize*1.25)\n", + "plt.plot(\n", + " [0] + intervals + [1],\n", + " [0] + projected_line + [0],\n", + " lw=5,\n", + " c=\"green\",\n", + " alpha=1,\n", + " label=f\"Closest in Output: {closest_fit_index_out}\",\n", + ")\n", + "\n", + "plt.legend(fontsize=fsize * 0.75)\n", + "plt.title(f\"Parameter Space\", fontsize=fsize * 1.25)\n", "plt.ylim(-4, 0)\n", "plt.xlim(0, 1)\n", - "plt.ylabel('$u(x, \\lambda)$', fontsize=fsize)\n", - "plt.xlabel('$x_1$', fontsize=fsize)\n", + "plt.ylabel(\"$u(x, \\lambda)$\", fontsize=fsize)\n", + "plt.xlabel(\"$x_1$\", fontsize=fsize)\n", "\n", "plt.subplot(122)\n", "_plot_num = max(measurements)\n", - "q = qoi[closest_fit_index_in,:]\n", - "plt.scatter(q[:_plot_num], qoi_ref[:_plot_num] + noise[:_plot_num], c='purple', s=100, alpha=0.4, label=f'Closest in Input: {closest_fit_index_in}')\n", - "q = qoi[closest_fit_index_out,:]\n", - "plt.scatter(q[:_plot_num], qoi_ref[:_plot_num] + noise[:_plot_num], c='green', s=100, alpha=0.4, label=f'Closest in Output: {closest_fit_index_out}')\n", - "\n", - "_a = np.linspace(min(qoi_ref),max(qoi_ref), 2)\n", - "plt.plot(_a,_a, c='k', lw=3)\n", + "q = qoi[closest_fit_index_in, :]\n", + "plt.scatter(\n", + " q[:_plot_num],\n", + " qoi_ref[:_plot_num] + noise[:_plot_num],\n", + " c=\"purple\",\n", + " s=100,\n", + " alpha=0.4,\n", + " label=f\"Closest in Input: {closest_fit_index_in}\",\n", + ")\n", + "q = qoi[closest_fit_index_out, :]\n", + "plt.scatter(\n", + " q[:_plot_num],\n", + " qoi_ref[:_plot_num] + noise[:_plot_num],\n", + " c=\"green\",\n", + " s=100,\n", + " alpha=0.4,\n", + " label=f\"Closest in Output: {closest_fit_index_out}\",\n", + ")\n", + "\n", + "_a = np.linspace(min(qoi_ref), max(qoi_ref), 2)\n", + "plt.plot(_a, _a, c=\"k\", lw=3)\n", "plt.xlim(-0.5, 0.2)\n", "plt.ylim(-0.5, 0.2)\n", "# plt.legend(fontsize=fsize)\n", - "plt.xlabel('Predicted Data', fontsize=fsize)\n", - "plt.ylabel('Collected Data', fontsize=fsize)\n", - "plt.title(\"Q-Q Plot\", fontsize=fsize*1.25)\n", + "plt.xlabel(\"Predicted Data\", fontsize=fsize)\n", + "plt.ylabel(\"Collected Data\", fontsize=fsize)\n", + "plt.title(\"Q-Q Plot\", fontsize=fsize * 1.25)\n", "\n", "_fname = f\"{prefix}_proj_{input_dim}D.{ftype}\"\n", "# plt.savefig(_fname, bbox_inches='tight')\n", @@ -1117,7 +1155,7 @@ "# plot_qoi = [20, 100, 500, 1000][::-1]\n", "plot_qoi = measurements[::-2]\n", "\n", - "colors = ['xkcd:red', 'xkcd:black', 'xkcd:orange', 'xkcd:blue', 'xkcd:green']" + "colors = [\"xkcd:red\", \"xkcd:black\", \"xkcd:orange\", \"xkcd:blue\", \"xkcd:green\"]" ] }, { @@ -1126,7 +1164,9 @@ "metadata": {}, "outputs": [], "source": [ - "_intervals = np.array(intervals[1:]) + ( np.array(intervals[:-1]) - np.array(intervals[1:]) ) / 2" + "_intervals = (\n", + " np.array(intervals[1:]) + (np.array(intervals[:-1]) - np.array(intervals[1:])) / 2\n", + ")" ] }, { @@ -1135,7 +1175,7 @@ "metadata": {}, "outputs": [], "source": [ - "qoi_indices = band_qoi(sensors, num_qoi, axis=1)\n", + "qoi_indices = band_qoi(sensors, num_qoi, axis=1)\n", "qoi_indices_bad = band_qoi(sensors, num_qoi, axis=0)" ] }, @@ -1145,8 +1185,9 @@ "metadata": {}, "outputs": [], "source": [ - "fdir = '/'.join(prefix.split('/')[:-1])\n", + "fdir = \"/\".join(prefix.split(\"/\")[:-1])\n", "from mud_examples.helpers import check_dir\n", + "\n", "check_dir(fdir)" ] }, @@ -1157,36 +1198,42 @@ "outputs": [], "source": [ "# horizontal plot\n", - "plt.figure(figsize=(10,10))\n", + "plt.figure(figsize=(10, 10))\n", "fin.plot(pn, vmin=-0.5, vmax=0)\n", - "plt.title(f\"Simulated Measurement Surface\\n$\\\\sigma$ = {sigma:1.3E} ($\\\\tau$ = {tolerance:1.1E})\")\n", + "plt.title(\n", + " f\"Simulated Measurement Surface\\n$\\\\sigma$ = {sigma:1.3E} ($\\\\tau$ = {tolerance:1.1E})\"\n", + ")\n", "for i in range(0, num_qoi):\n", - " if i < num_qoi - 1: plt.axhline(_intervals[i], lw=3, c='k')\n", - " _q = qoi_indices[i][qoi_indices[i] < 100 ]\n", - " plt.scatter(sensors[_q,0], sensors[_q,1], s=100, color=colors[i%2])\n", - "plt.scatter([0]*input_dim, intervals, s=200, marker='^', c='w')\n", - "plt.xlim(0,1)\n", - "plt.ylim(0,1)\n", + " if i < num_qoi - 1:\n", + " plt.axhline(_intervals[i], lw=3, c=\"k\")\n", + " _q = qoi_indices[i][qoi_indices[i] < 100]\n", + " plt.scatter(sensors[_q, 0], sensors[_q, 1], s=100, color=colors[i % 2])\n", + "plt.scatter([0] * input_dim, intervals, s=200, marker=\"^\", c=\"w\")\n", + "plt.xlim(0, 1)\n", + "plt.ylim(0, 1)\n", "plt.xticks([])\n", "plt.yticks([])\n", "plt.xlabel(\"$x_1$\", fontsize=fsize)\n", "plt.ylabel(\"$x_2$\", fontsize=fsize)\n", "\n", "_fname = f\"{prefix}_sensors_{input_dim}D.{ftype}\"\n", - "plt.savefig(_fname, bbox_inches='tight')\n", + "plt.savefig(_fname, bbox_inches=\"tight\")\n", "plt.show()\n", "\n", "# vertical plot\n", - "plt.figure(figsize=(10,10))\n", + "plt.figure(figsize=(10, 10))\n", "fin.plot(pn, vmin=-0.5, vmax=0)\n", - "plt.title(f\"Simulated Measurement Surface\\n$\\\\sigma$ = {sigma:1.3E} ($\\\\tau$ = {tolerance:1.1E})\")\n", + "plt.title(\n", + " f\"Simulated Measurement Surface\\n$\\\\sigma$ = {sigma:1.3E} ($\\\\tau$ = {tolerance:1.1E})\"\n", + ")\n", "for i in range(0, num_qoi):\n", - " if i < num_qoi - 1: plt.axvline(_intervals[i], lw=3, c='k')\n", - " _q = qoi_indices_bad[i][qoi_indices_bad[i] < 100 ]\n", - " plt.scatter(sensors[_q,0], sensors[_q,1], s=100, color=colors[i%2])\n", - "plt.scatter([0]*input_dim,intervals, s=200, marker='^', c='w')\n", - "plt.xlim(0,1)\n", - "plt.ylim(0,1)\n", + " if i < num_qoi - 1:\n", + " plt.axvline(_intervals[i], lw=3, c=\"k\")\n", + " _q = qoi_indices_bad[i][qoi_indices_bad[i] < 100]\n", + " plt.scatter(sensors[_q, 0], sensors[_q, 1], s=100, color=colors[i % 2])\n", + "plt.scatter([0] * input_dim, intervals, s=200, marker=\"^\", c=\"w\")\n", + "plt.xlim(0, 1)\n", + "plt.ylim(0, 1)\n", "plt.xticks([])\n", "plt.yticks([])\n", "plt.xlabel(\"$x_1$\", fontsize=fsize)\n", @@ -1211,57 +1258,83 @@ "outputs": [], "source": [ "if input_dim == 2:\n", - " plt.figure(figsize=(20,10))\n", + " plt.figure(figsize=(20, 10))\n", " plt.subplot(121)\n", - " colors = ['xkcd:red', 'xkcd:black', 'xkcd:orange', 'xkcd:blue', 'xkcd:green'][::-1]\n", - "# plot_qoi = [20, 100, 500, 1000][::-1]\n", + " colors = [\"xkcd:red\", \"xkcd:black\", \"xkcd:orange\", \"xkcd:blue\", \"xkcd:green\"][::-1]\n", + " # plot_qoi = [20, 100, 500, 1000][::-1]\n", " for idx, _first in enumerate(plot_qoi):\n", - " qois = split_qoi_by_indices(qoi_indices, qoi_ref, qoi, noise, sigma, max_index=_first)\n", - " plt.scatter(qois[0], qois[1], label=f'First {_first}', s=20, c=colors[idx], alpha=1)\n", + " qois = split_qoi_by_indices(\n", + " qoi_indices, qoi_ref, qoi, noise, sigma, max_index=_first\n", + " )\n", + " plt.scatter(\n", + " qois[0], qois[1], label=f\"First {_first}\", s=20, c=colors[idx], alpha=1\n", + " )\n", " plt.legend()\n", - " plt.title(\"Horizontal Band QoI\", fontsize=1.25*fsize)\n", + " plt.title(\"Horizontal Band QoI\", fontsize=1.25 * fsize)\n", " plt.xlabel(\"$q_1$\", fontsize=fsize)\n", " plt.ylabel(\"$q_2$\", fontsize=fsize)\n", - " \n", + "\n", " plt.subplot(122)\n", " for idx, _first in enumerate(plot_qoi):\n", - " qois = split_qoi_by_indices(qoi_indices_bad, qoi_ref, qoi, noise, sigma, max_index=_first)\n", - " plt.scatter(qois[0], qois[1], label=f'First {_first}', s=20, c=colors[idx], alpha=1)\n", + " qois = split_qoi_by_indices(\n", + " qoi_indices_bad, qoi_ref, qoi, noise, sigma, max_index=_first\n", + " )\n", + " plt.scatter(\n", + " qois[0], qois[1], label=f\"First {_first}\", s=20, c=colors[idx], alpha=1\n", + " )\n", " plt.legend()\n", - " plt.title(\"Vertical Band QoI\", fontsize=1.25*fsize)\n", + " plt.title(\"Vertical Band QoI\", fontsize=1.25 * fsize)\n", " plt.xlabel(\"$q^*_1$\", fontsize=fsize)\n", " plt.ylabel(\"$q^*_2$\", fontsize=fsize)\n", "\n", "else:\n", - " plt.figure(figsize=(20,20))\n", - "# lim = 7.5/tolerance\n", - " lim = 3/tolerance\n", - " fig, axs = plt.subplots(input_dim, input_dim, figsize=(20,20))\n", + " plt.figure(figsize=(20, 20))\n", + " # lim = 7.5/tolerance\n", + " lim = 3 / tolerance\n", + " fig, axs = plt.subplots(input_dim, input_dim, figsize=(20, 20))\n", " for _i in range(input_dim):\n", " for _j in range(_i, input_dim):\n", " if _i == _j:\n", " ax = axs[_i][_i]\n", "\n", - " ax.annotate(f\"$q{_i+1}$\", (-lim/10,0), fontsize=fsize)\n", - " # ax.set_ylabel(f\"$q{_i+1}$\")\n", - " # ax.set_xlabel(f\"$q{_i+1}$\")\n", + " ax.annotate(f\"$q{_i+1}$\", (-lim / 10, 0), fontsize=fsize)\n", + " # ax.set_ylabel(f\"$q{_i+1}$\")\n", + " # ax.set_xlabel(f\"$q{_i+1}$\")\n", " ax.set_xlim(-lim, lim)\n", " ax.set_ylim(-lim, lim)\n", - " # ax.set_xticks([])\n", - " # ax.set_yticks([])\n", + " # ax.set_xticks([])\n", + " # ax.set_yticks([])\n", " else:\n", " for idx, _first in enumerate(plot_qoi):\n", " ax = axs[_i][_j]\n", - " qois = split_qoi_by_indices(qoi_indices, qoi_ref, qoi, noise, sigma, max_index=_first)\n", - " ax.scatter(qois[_i], qois[_j], label=f'First {_first}', s=20, c=colors[idx], alpha=1)\n", + " qois = split_qoi_by_indices(\n", + " qoi_indices, qoi_ref, qoi, noise, sigma, max_index=_first\n", + " )\n", + " ax.scatter(\n", + " qois[_i],\n", + " qois[_j],\n", + " label=f\"First {_first}\",\n", + " s=20,\n", + " c=colors[idx],\n", + " alpha=1,\n", + " )\n", " ax.set_xlim(-lim, lim)\n", " ax.set_ylim(-lim, lim)\n", " ax.set_xticks([])\n", " ax.set_yticks([])\n", "\n", " ax = axs[_j][_i]\n", - " qois = split_qoi_by_indices(qoi_indices_bad, qoi_ref, qoi, noise, sigma, max_index=_first)\n", - " ax.scatter(qois[_i], qois[_j], label=f'First {_first}', s=20, c=colors[idx], alpha=1)\n", + " qois = split_qoi_by_indices(\n", + " qoi_indices_bad, qoi_ref, qoi, noise, sigma, max_index=_first\n", + " )\n", + " ax.scatter(\n", + " qois[_i],\n", + " qois[_j],\n", + " label=f\"First {_first}\",\n", + " s=20,\n", + " c=colors[idx],\n", + " alpha=1,\n", + " )\n", " ax.set_xlim(-lim, lim)\n", " ax.set_ylim(-lim, lim)\n", " ax.set_xticks([])\n", @@ -1269,7 +1342,7 @@ "\n", "_fname = f\"{prefix}_geom_{input_dim}D.{ftype}\"\n", "# plt.savefig(_fname, bbox_inches='tight')\n", - "plt.show() " + "plt.show()" ] }, { @@ -1303,7 +1376,7 @@ "outputs": [], "source": [ "%%time\n", - "X = qoi[:,0:num_samps_ex_sol]\n", + "X = qoi[:, 0:num_samps_ex_sol]\n", "data = qoi_ref[0:num_samps_ex_sol] + noise[0:num_samps_ex_sol]\n", "newqoi = wme(X, data, sigma)\n", "r_sing = ratio_dci_sing(newqoi)" @@ -1316,7 +1389,7 @@ "outputs": [], "source": [ "mud_idx = np.argmax(r_sing)\n", - "mud_fun = lam[mud_idx,:]\n", + "mud_fun = lam[mud_idx, :]\n", "print(mud_idx)" ] }, @@ -1352,7 +1425,9 @@ "outputs": [], "source": [ "%%time\n", - "qois = split_qoi_by_indices(qoi_indices, qoi_ref, qoi, noise, sigma, max_index=num_samps_ex_sol)\n", + "qois = split_qoi_by_indices(\n", + " qoi_indices, qoi_ref, qoi, noise, sigma, max_index=num_samps_ex_sol\n", + ")\n", "r_mult = ratio_dci_mult(qois)" ] }, @@ -1363,7 +1438,7 @@ "outputs": [], "source": [ "mud_idx_mult = np.argmax(r_mult)\n", - "mud_fun_mult = lam[mud_idx_mult,:]" + "mud_fun_mult = lam[mud_idx_mult, :]" ] }, { @@ -1388,34 +1463,35 @@ "metadata": {}, "outputs": [], "source": [ - "plt.figure(figsize=(30,10))\n", - "colors = ['xkcd:red', 'xkcd:black', 'xkcd:orange', 'xkcd:blue', 'xkcd:green']\n", + "plt.figure(figsize=(30, 10))\n", + "colors = [\"xkcd:red\", \"xkcd:black\", \"xkcd:orange\", \"xkcd:blue\", \"xkcd:green\"]\n", "plt.subplot(131)\n", - "fin.plot(poissonModel(mud_fun, nx=36,ny=36), vmin=-0.5, vmax=0)\n", - "plt.title('MUD (Scalar-Valued)', fontsize=1.25*fsize)\n", + "fin.plot(poissonModel(mud_fun, nx=36, ny=36), vmin=-0.5, vmax=0)\n", + "plt.title(\"MUD (Scalar-Valued)\", fontsize=1.25 * fsize)\n", "\n", "plt.subplot(132)\n", "fin.plot(pn, vmin=-0.5, vmax=0)\n", - "plt.title('(Noisy) Response Surface', fontsize=1.25*fsize)\n", + "plt.title(\"(Noisy) Response Surface\", fontsize=1.25 * fsize)\n", "for i in range(0, num_qoi):\n", - " if i < num_qoi - 1: plt.axhline(_intervals[i], lw=3, c='k')\n", - " _q = qoi_indices[i][qoi_indices[i] < 100 ]\n", - " plt.scatter(sensors[_q,0], sensors[_q,1], s=100, color=colors[i%2])\n", - "plt.scatter([0]*input_dim, intervals, s=400, marker='^', c='w')\n", - "plt.xlim(0,1)\n", - "plt.ylim(0,1)\n", + " if i < num_qoi - 1:\n", + " plt.axhline(_intervals[i], lw=3, c=\"k\")\n", + " _q = qoi_indices[i][qoi_indices[i] < 100]\n", + " plt.scatter(sensors[_q, 0], sensors[_q, 1], s=100, color=colors[i % 2])\n", + "plt.scatter([0] * input_dim, intervals, s=400, marker=\"^\", c=\"w\")\n", + "plt.xlim(0, 1)\n", + "plt.ylim(0, 1)\n", "plt.xticks([])\n", "plt.yticks([])\n", "plt.xlabel(\"$x_1$\", fontsize=fsize)\n", "plt.ylabel(\"$x_2$\", fontsize=fsize)\n", "\n", "plt.subplot(133)\n", - "fin.plot(poissonModel(mud_fun_mult, nx=36,ny=36), vmin=-0.5, vmax=0)\n", - "plt.title('MUD (Vector-Valued)', fontsize=1.25*fsize)\n", + "fin.plot(poissonModel(mud_fun_mult, nx=36, ny=36), vmin=-0.5, vmax=0)\n", + "plt.title(\"MUD (Vector-Valued)\", fontsize=1.25 * fsize)\n", "\n", "_fname = f\"{prefix}_surf_exmud_{input_dim}D_m{num_samps_ex_sol}.{ftype}\"\n", "# plt.savefig(_fname, bbox_inches='tight')\n", - "plt.show()\n" + "plt.show()" ] }, { @@ -1424,31 +1500,41 @@ "metadata": {}, "outputs": [], "source": [ - "plt.figure(figsize=(10,10))\n", + "plt.figure(figsize=(10, 10))\n", "\n", "plt.subplot(111)\n", - "fin.plot(u_plot, mesh=mesh, lw=5, c='k')\n", - "fin.plot(w, mesh=mesh, lw=5, c='k', ls='--', alpha=0.5, label='Interpolant')\n", + "fin.plot(u_plot, mesh=mesh, lw=5, c=\"k\")\n", + "fin.plot(w, mesh=mesh, lw=5, c=\"k\", ls=\"--\", alpha=0.5, label=\"Interpolant\")\n", "\n", "# fin.plot(w, mesh=mesh, lw=10, c='r', ls='-', alpha=0.5)\n", - "plt.scatter(intervals, lam_ref, marker='^', s=200, c='purple', zorder=10)\n", - "plt.title(f'Ex. MUD Solution, m={num_samps_ex_sol}', fontsize=1.25*fsize)\n", - "\n", - "plt.plot(np.linspace(0,1,input_dim+2),\n", - " [0] + list(mud_fun) + [0],\n", - " alpha=0.5, c='r', lw=10, label=f'Scalar MUD, Sample {mud_idx}')\n", - "\n", - "plt.plot(np.linspace(0,1,input_dim+2),\n", - " [0] + list(mud_fun_mult) + [0],\n", - " alpha=0.5, c='b', lw=10, label=f'Vector MUD, Sample {mud_idx_mult}')\n", - "\n", - "plt.axvline(2/7, alpha=0.4, ls=':')\n", - "plt.axhline(-lam_true, alpha=0.4, ls=':')\n", - "plt.ylim(-4,0)\n", - "plt.xlim(0,1)\n", - "plt.ylabel('$u(x, \\lambda)$', fontsize=fsize)\n", - "plt.xlabel('$x_1$', fontsize=fsize)\n", - "plt.legend(fontsize=fsize*0.6, loc='lower right')\n", + "plt.scatter(intervals, lam_ref, marker=\"^\", s=200, c=\"purple\", zorder=10)\n", + "plt.title(f\"Ex. MUD Solution, m={num_samps_ex_sol}\", fontsize=1.25 * fsize)\n", + "\n", + "plt.plot(\n", + " np.linspace(0, 1, input_dim + 2),\n", + " [0] + list(mud_fun) + [0],\n", + " alpha=0.5,\n", + " c=\"r\",\n", + " lw=10,\n", + " label=f\"Scalar MUD, Sample {mud_idx}\",\n", + ")\n", + "\n", + "plt.plot(\n", + " np.linspace(0, 1, input_dim + 2),\n", + " [0] + list(mud_fun_mult) + [0],\n", + " alpha=0.5,\n", + " c=\"b\",\n", + " lw=10,\n", + " label=f\"Vector MUD, Sample {mud_idx_mult}\",\n", + ")\n", + "\n", + "plt.axvline(2 / 7, alpha=0.4, ls=\":\")\n", + "plt.axhline(-lam_true, alpha=0.4, ls=\":\")\n", + "plt.ylim(-4, 0)\n", + "plt.xlim(0, 1)\n", + "plt.ylabel(\"$u(x, \\lambda)$\", fontsize=fsize)\n", + "plt.xlabel(\"$x_1$\", fontsize=fsize)\n", + "plt.legend(fontsize=fsize * 0.6, loc=\"lower right\")\n", "\n", "# plt.subplot(122)\n", "# _data = qoi_ref[:num_samps_ex_sol] + noise[:num_samps_ex_sol]\n", @@ -1485,63 +1571,111 @@ "metadata": {}, "outputs": [], "source": [ - "labels = ['Scalar QoI', 'Vector QoI']\n", + "labels = [\"Scalar QoI\", \"Vector QoI\"]\n", "plot_top = 1000\n", - "thresh = lam.shape[0]**-1\n", + "thresh = lam.shape[0] ** -1\n", "# thresh = 1E-16\n", - "colors = ['xkcd:red', 'xkcd:black', 'xkcd:orange', 'xkcd:blue', 'xkcd:green']\n", + "colors = [\"xkcd:red\", \"xkcd:black\", \"xkcd:orange\", \"xkcd:blue\", \"xkcd:green\"]\n", "\n", "if input_dim == 2:\n", " for _i in range(input_dim):\n", - " for _j in range(_i+1, input_dim):\n", + " for _j in range(_i + 1, input_dim):\n", " for idx, ratio_eval in enumerate([r_sing, r_mult]):\n", - " _m = np.where(ratio_eval/max(ratio_eval) > thresh)[0]\n", - "# plt.scatter(lam[_m[101:5000], _i], lam[_m[101:5000], _j], c='orange', marker='^', alpha=0.2)\n", - "# plt.scatter(lam[_m[plot_top+1:1000], _i], lam[_m[plot_top+1:1000], _j], c='orange', marker='^', alpha=0.2)\n", - " plt.scatter(lam[_m[:plot_top], _i], lam[_m[:plot_top], _j], c=colors[idx], label= labels[idx] + ' (Total %d)'%len(_m), s=50)\n", + " _m = np.where(ratio_eval / max(ratio_eval) > thresh)[0]\n", + " # plt.scatter(lam[_m[101:5000], _i], lam[_m[101:5000], _j], c='orange', marker='^', alpha=0.2)\n", + " # plt.scatter(lam[_m[plot_top+1:1000], _i], lam[_m[plot_top+1:1000], _j], c='orange', marker='^', alpha=0.2)\n", + " plt.scatter(\n", + " lam[_m[:plot_top], _i],\n", + " lam[_m[:plot_top], _j],\n", + " c=colors[idx],\n", + " label=labels[idx] + \" (Total %d)\" % len(_m),\n", + " s=50,\n", + " )\n", " plt.xlabel(f\"$\\lambda_{_i+1}$\", fontsize=fsize)\n", " plt.ylabel(f\"$\\lambda_{_j+1}$\", fontsize=fsize)\n", - " plt.ylim(-4,0)\n", - " plt.xlim(-4,0)\n", - " plt.scatter(lam_ref[_i], lam_ref[_j], c='k', s=500, alpha=0.5, label='Interpolant', zorder=-10)\n", - " plt.scatter(lam[closest_fit_index_out, _i], lam[closest_fit_index_out, _j], c='g', s=500, alpha=0.8, label='Projection', zorder=15)\n", + " plt.ylim(-4, 0)\n", + " plt.xlim(-4, 0)\n", + " plt.scatter(\n", + " lam_ref[_i],\n", + " lam_ref[_j],\n", + " c=\"k\",\n", + " s=500,\n", + " alpha=0.5,\n", + " label=\"Interpolant\",\n", + " zorder=-10,\n", + " )\n", + " plt.scatter(\n", + " lam[closest_fit_index_out, _i],\n", + " lam[closest_fit_index_out, _j],\n", + " c=\"g\",\n", + " s=500,\n", + " alpha=0.8,\n", + " label=\"Projection\",\n", + " zorder=15,\n", + " )\n", " plt.legend()\n", - " plt.title(f\"Samples (m = {num_samps_ex_sol}) with\\nRelative Ratio > {thresh:1.1E}\", fontsize=fsize)\n", + " plt.title(\n", + " f\"Samples (m = {num_samps_ex_sol}) with\\nRelative Ratio > {thresh:1.1E}\",\n", + " fontsize=fsize,\n", + " )\n", "\n", "else:\n", - "\n", - " fig, axs = plt.subplots(input_dim, input_dim, figsize=(20,20))\n", + " fig, axs = plt.subplots(input_dim, input_dim, figsize=(20, 20))\n", " for _i in range(input_dim):\n", " for _j in range(_i, input_dim):\n", " if _i != _j:\n", " for idx, ratio_eval in enumerate([r_sing, r_mult]):\n", " ax = axs[_j][_i] if not idx else axs[_i][_j]\n", - " _m = np.where(ratio_eval/max(ratio_eval) > thresh)[0]\n", - " ax.scatter(lam[_m[:plot_top], _i], lam[_m[:plot_top], _j], c=colors[2+idx], label= str(len(_m)) + ' ' + labels[idx], s=10)\n", - "\n", - "\n", - " # ax.set_xlabel(f\"$\\lambda_{_i+1}$\")\n", - " # ax.set_ylabel(f\"$\\lambda_{_j+1}$\")\n", - " ax.set_ylim(-4,0)\n", - " ax.set_xlim(-4,0)\n", + " _m = np.where(ratio_eval / max(ratio_eval) > thresh)[0]\n", + " ax.scatter(\n", + " lam[_m[:plot_top], _i],\n", + " lam[_m[:plot_top], _j],\n", + " c=colors[2 + idx],\n", + " label=str(len(_m)) + \" \" + labels[idx],\n", + " s=10,\n", + " )\n", + "\n", + " # ax.set_xlabel(f\"$\\lambda_{_i+1}$\")\n", + " # ax.set_ylabel(f\"$\\lambda_{_j+1}$\")\n", + " ax.set_ylim(-4, 0)\n", + " ax.set_xlim(-4, 0)\n", " ax.set_xticks([])\n", " ax.set_yticks([])\n", - " ax.scatter(lam_ref[_i], lam_ref[_j], c='k', s=250, alpha=1, label='Interpolant', zorder=-10)\n", - " ax.scatter(lam[closest_fit_index_out, _i], lam[closest_fit_index_out, _j], c='g', s=400, alpha=0.5, label='Projection', zorder=-10)\n", - " # ax.legend()\n", + " ax.scatter(\n", + " lam_ref[_i],\n", + " lam_ref[_j],\n", + " c=\"k\",\n", + " s=250,\n", + " alpha=1,\n", + " label=\"Interpolant\",\n", + " zorder=-10,\n", + " )\n", + " ax.scatter(\n", + " lam[closest_fit_index_out, _i],\n", + " lam[closest_fit_index_out, _j],\n", + " c=\"g\",\n", + " s=400,\n", + " alpha=0.5,\n", + " label=\"Projection\",\n", + " zorder=-10,\n", + " )\n", + " # ax.legend()\n", " else:\n", " ax = axs[_i][_i]\n", - " ax.annotate(f\"$\\lambda_{_i+1}$\", (-0.6,0.5), fontsize=fsize)\n", - "# ax.set_xlabel(f\"$\\lambda_{_i+1}$\", fontsize=fsize)\n", - "# ax.set_ylabel(f\"$\\lambda_{_i+1}$\", fontsize=fsize)\n", + " ax.annotate(f\"$\\lambda_{_i+1}$\", (-0.6, 0.5), fontsize=fsize)\n", + " # ax.set_xlabel(f\"$\\lambda_{_i+1}$\", fontsize=fsize)\n", + " # ax.set_ylabel(f\"$\\lambda_{_i+1}$\", fontsize=fsize)\n", "\n", " ax.set_xticks([])\n", " ax.set_yticks([])\n", - " ax.set_xlim(-1,0)\n", + " ax.set_xlim(-1, 0)\n", "\n", "\n", - "_fname = f\"{prefix}_update_scatter_{input_dim}D_t{thresh:1.1E}\".replace('.', '-') + f\".{ftype}\"\n", - "plt.savefig(_fname, bbox_inches='tight')\n", + "_fname = (\n", + " f\"{prefix}_update_scatter_{input_dim}D_t{thresh:1.1E}\".replace(\".\", \"-\")\n", + " + f\".{ftype}\"\n", + ")\n", + "plt.savefig(_fname, bbox_inches=\"tight\")\n", "plt.show()" ] }, @@ -1565,28 +1699,43 @@ "if not load:\n", " # SCALAR\n", " def mud_wrapper(num_obs, sd):\n", - " newqoi = wme(X=qoi[:,0:num_obs], data=qoi_ref[0:num_obs] + np.random.randn(num_obs)*sd, sd=sd)\n", + " newqoi = wme(\n", + " X=qoi[:, 0:num_obs],\n", + " data=qoi_ref[0:num_obs] + np.random.randn(num_obs) * sd,\n", + " sd=sd,\n", + " )\n", " r_sing = ratio_dci_sing(newqoi)\n", " return r_sing\n", "\n", - " experiments_sing, solutions_sing = experiment_measurements_index(num_measurements=measurements,\n", - " sd=sigma,\n", - " num_trials=num_trials,\n", - " seed=21,\n", - " fun=mud_wrapper)\n", + " experiments_sing, solutions_sing = experiment_measurements_index(\n", + " num_measurements=measurements,\n", + " sd=sigma,\n", + " num_trials=num_trials,\n", + " seed=21,\n", + " fun=mud_wrapper,\n", + " )\n", + "\n", " # VECTOR\n", " def mud_wrapper(num_obs, sd):\n", - " qois = split_qoi_by_indices(qoi_indices, qoi_ref, qoi,\n", - " noise=np.random.randn(num_obs)*sd, sigma=sd, max_index=num_obs)\n", + " qois = split_qoi_by_indices(\n", + " qoi_indices,\n", + " qoi_ref,\n", + " qoi,\n", + " noise=np.random.randn(num_obs) * sd,\n", + " sigma=sd,\n", + " max_index=num_obs,\n", + " )\n", " r_mult = ratio_dci_mult(qois)\n", " return r_mult\n", - " \n", - " experiments_mult, solutions_mult = experiment_measurements_index(num_measurements=measurements,\n", - " sd=sigma,\n", - " num_trials=num_trials,\n", - " seed=21,\n", - " fun=mud_wrapper)\n", - " \n", + "\n", + " experiments_mult, solutions_mult = experiment_measurements_index(\n", + " num_measurements=measurements,\n", + " sd=sigma,\n", + " num_trials=num_trials,\n", + " seed=21,\n", + " fun=mud_wrapper,\n", + " )\n", + "\n", " del mud_wrapper" ] }, @@ -1627,9 +1776,9 @@ " np.random.seed(21)\n", " num_draws = 3\n", " for idx in range(num_draws):\n", - " i = np.random.randint(0,len(model_list))\n", - " mudU = fin.Function(V, model_list[i][i]['u'])\n", - " plt.subplot(int(f'{num_draws}{3}{1+3*idx}'))\n", + " i = np.random.randint(0, len(model_list))\n", + " mudU = fin.Function(V, model_list[i][i][\"u\"])\n", + " plt.subplot(int(f\"{num_draws}{3}{1+3*idx}\"))\n", " fin.plot(mudU, vmin=-0.5, vmax=0)\n", " plt.xticks([])\n", " plt.yticks([])\n", @@ -1638,22 +1787,22 @@ "\n", " num_plot_sensors = max(measurements)\n", " for idx in range(num_draws):\n", - " # _r = np.random.randint(0, num_trials)\n", + " # _r = np.random.randint(0, num_trials)\n", " _r = idx\n", " i = solutions[num_plot_sensors][_r]\n", - " mudU = fin.Function(V, model_list[i][i]['u'])\n", - " plt.subplot(int(f'{num_draws}{3}{2+3*idx}'))\n", + " mudU = fin.Function(V, model_list[i][i][\"u\"])\n", + " plt.subplot(int(f\"{num_draws}{3}{2+3*idx}\"))\n", " fin.plot(mudU, vmin=-0.5, vmax=0)\n", " plt.xticks([])\n", " plt.yticks([])\n", " plt.title(f\"MUD#{idx}: {i}\")\n", "\n", - " q = qoi[i,:]\n", - " plt.subplot(int(f'{num_draws}{3}{3+3*idx}'))\n", - " plt.scatter(qoi_ref + noise, q, c='b', s=50, alpha=0.05)\n", - " plt.plot(_a,_a, c='k', lw=3)\n", - " # plt.xlabel('True QoI (Noiseless)')\n", - " # plt.ylabel('Predicted Signal')\n", + " q = qoi[i, :]\n", + " plt.subplot(int(f\"{num_draws}{3}{3+3*idx}\"))\n", + " plt.scatter(qoi_ref + noise, q, c=\"b\", s=50, alpha=0.05)\n", + " plt.plot(_a, _a, c=\"k\", lw=3)\n", + " # plt.xlabel('True QoI (Noiseless)')\n", + " # plt.ylabel('Predicted Signal')\n", " plt.xticks([])\n", " plt.yticks([])\n", " plt.xlim(-0.5, 0.2)\n", @@ -1682,33 +1831,40 @@ "outputs": [], "source": [ "gt = list(lam[closest_fit_index_out, :])\n", - "plt.figure(figsize=(10,10))\n", + "plt.figure(figsize=(10, 10))\n", "\n", "# fin.plot(u_plot, mesh=mesh, lw=5, c='k', label=\"$g$\")\n", - "plt.plot([0]+intervals+[1], [0]+gt+[0], lw=5, c='green', alpha=0.6, ls='--', label='$\\hat{g}$', zorder=5)\n", + "plt.plot(\n", + " [0] + intervals + [1],\n", + " [0] + gt + [0],\n", + " lw=5,\n", + " c=\"green\",\n", + " alpha=0.6,\n", + " ls=\"--\",\n", + " label=\"$\\hat{g}$\",\n", + " zorder=5,\n", + ")\n", "\n", "# plt.scatter(intervals, lam_ref, marker='^', s=200, c='purple', zorder=10)\n", "# plt.plot([0]+intervals+[1], [0]+projected_line+[0], lw=5, c='green', alpha=1, label=\"$Proj_{g}$\")\n", "\n", "\n", "for i in range(100):\n", - " _lam= lam[i,:]\n", - " plt.plot([0]+intervals+[1], [0]+list(_lam)+[0], lw=1, c='purple', alpha=0.2)\n", - "plt.title('Samples from Initial Density', fontsize=1.25*fsize)\n", + " _lam = lam[i, :]\n", + " plt.plot([0] + intervals + [1], [0] + list(_lam) + [0], lw=1, c=\"purple\", alpha=0.2)\n", + "plt.title(\"Samples from Initial Density\", fontsize=1.25 * fsize)\n", "plt.xlabel(\"$x_2$\", fontsize=fsize)\n", "plt.ylabel(\"$g(x, \\lambda)$\", fontsize=fsize)\n", "\n", "\n", - "\n", - "\n", "# plt.axvline(2/7, alpha=0.4, ls=':')\n", "# plt.axhline(-lam_true, alpha=0.4, ls=':')\n", - "plt.ylim(-4,0)\n", - "plt.xlim(0,1)\n", + "plt.ylim(-4, 0)\n", + "plt.xlim(0, 1)\n", "plt.legend()\n", "\n", "_fname = f\"{prefix}_in_{input_dim}D.{ftype}\"\n", - "plt.savefig(_fname, bbox_inches='tight')\n", + "plt.savefig(_fname, bbox_inches=\"tight\")\n", "plt.show()" ] }, @@ -1719,35 +1875,47 @@ "outputs": [], "source": [ "for num_plot_sensors in [20, 100]:\n", - " plt.figure(figsize=(10,10))\n", - " \n", + " plt.figure(figsize=(10, 10))\n", + "\n", " plt.subplot(111)\n", - " plt.title('MUD Estimates for $Q_{1D}$,' + f' N={num_plot_sensors}', fontsize=1.25*fsize)\n", + " plt.title(\n", + " \"MUD Estimates for $Q_{1D}$,\" + f\" N={num_plot_sensors}\", fontsize=1.25 * fsize\n", + " )\n", " plt.xlabel(\"$x_2$\", fontsize=fsize)\n", " plt.ylabel(\"$g(x, \\lambda)$\", fontsize=fsize)\n", - " for i in solutions_sing[num_plot_sensors]: # trials\n", - " _lam = lam[i,:]\n", - " plt.plot([0]+intervals+[1], [0]+list(_lam)+[0], lw=1, c='purple', alpha=0.2)\n", - "# fin.plot(u_plot, mesh=mesh, lw=5, c='k', label=\"$g$\")\n", - " plt.plot([0]+intervals+[1], [0]+gt+[0], lw=5, c='green', alpha=0.6, ls='--', label='$\\hat{g}$')\n", + " for i in solutions_sing[num_plot_sensors]: # trials\n", + " _lam = lam[i, :]\n", + " plt.plot(\n", + " [0] + intervals + [1], [0] + list(_lam) + [0], lw=1, c=\"purple\", alpha=0.2\n", + " )\n", + " # fin.plot(u_plot, mesh=mesh, lw=5, c='k', label=\"$g$\")\n", + " plt.plot(\n", + " [0] + intervals + [1],\n", + " [0] + gt + [0],\n", + " lw=5,\n", + " c=\"green\",\n", + " alpha=0.6,\n", + " ls=\"--\",\n", + " label=\"$\\hat{g}$\",\n", + " )\n", "\n", " plt.ylim(-4, 0)\n", " plt.xlim(0, 1)\n", - " plt.legend(loc='lower left')\n", - "# plt.subplot(122)\n", - "# plt.title('Q-Q Plot', fontsize=1.25*fsize)\n", - "# for i in solutions_sing[num_plot_sensors]: # trials\n", - "# q = qoi[i,:]\n", - "# plt.scatter(q[:num_plot_sensors], qoi_ref[:num_plot_sensors] + noise[:num_plot_sensors],\n", - "# c='b', s=100, alpha=1.0/num_trials)\n", - "# plt.plot(_a,_a, c='k', lw=3)\n", - "# plt.ylabel('Collected Data', fontsize=fsize)\n", - "# plt.xlabel('Predicted Data', fontsize=fsize)\n", - "# plt.ylim(-0.5, 0.15)\n", - "# plt.xlim(-0.5, 0.15)\n", - "# plt.title(f'Solution {_r}, Index {_s}')\n", + " plt.legend(loc=\"lower left\")\n", + " # plt.subplot(122)\n", + " # plt.title('Q-Q Plot', fontsize=1.25*fsize)\n", + " # for i in solutions_sing[num_plot_sensors]: # trials\n", + " # q = qoi[i,:]\n", + " # plt.scatter(q[:num_plot_sensors], qoi_ref[:num_plot_sensors] + noise[:num_plot_sensors],\n", + " # c='b', s=100, alpha=1.0/num_trials)\n", + " # plt.plot(_a,_a, c='k', lw=3)\n", + " # plt.ylabel('Collected Data', fontsize=fsize)\n", + " # plt.xlabel('Predicted Data', fontsize=fsize)\n", + " # plt.ylim(-0.5, 0.15)\n", + " # plt.xlim(-0.5, 0.15)\n", + " # plt.title(f'Solution {_r}, Index {_s}')\n", " _fname = f\"{prefix}_mud_{input_dim}-1_N{num_plot_sensors}.{ftype}\"\n", - " plt.savefig(_fname, bbox_inches='tight')\n", + " plt.savefig(_fname, bbox_inches=\"tight\")\n", " plt.show()" ] }, @@ -1758,36 +1926,48 @@ "outputs": [], "source": [ "for num_plot_sensors in [20, 100]:\n", - " plt.figure(figsize=(10,10))\n", - " \n", + " plt.figure(figsize=(10, 10))\n", + "\n", " plt.subplot(111)\n", - " plt.title('MUD Estimates for $Q_{2D}$,' + f' N={num_plot_sensors}', fontsize=1.25*fsize)\n", + " plt.title(\n", + " \"MUD Estimates for $Q_{2D}$,\" + f\" N={num_plot_sensors}\", fontsize=1.25 * fsize\n", + " )\n", " plt.xlabel(\"$x_1$\", fontsize=fsize)\n", " plt.ylabel(\"$g(x, \\lambda)$\", fontsize=fsize)\n", - "# plt.plot([0]+intervals+[1], [0]+ [w(i) for i in intervals] +[0], lw=5, c='k',label=\"Interpolant\")\n", - " for i in solutions_mult[num_plot_sensors]: # trials\n", - " _lam = lam[i,:]\n", - " plt.plot([0]+intervals+[1], [0]+list(_lam)+[0], lw=1, c='purple', alpha=0.2)\n", - "# fin.plot(u_plot, mesh=mesh, lw=5, c='k', label=\"$g$\")\n", + " # plt.plot([0]+intervals+[1], [0]+ [w(i) for i in intervals] +[0], lw=5, c='k',label=\"Interpolant\")\n", + " for i in solutions_mult[num_plot_sensors]: # trials\n", + " _lam = lam[i, :]\n", + " plt.plot(\n", + " [0] + intervals + [1], [0] + list(_lam) + [0], lw=1, c=\"purple\", alpha=0.2\n", + " )\n", + " # fin.plot(u_plot, mesh=mesh, lw=5, c='k', label=\"$g$\")\n", " gt = list(lam[closest_fit_index_out, :])\n", - " plt.plot([0]+intervals+[1], [0]+gt+[0], lw=5, c='green', alpha=0.6, ls='--', label='$\\hat{g}$')\n", + " plt.plot(\n", + " [0] + intervals + [1],\n", + " [0] + gt + [0],\n", + " lw=5,\n", + " c=\"green\",\n", + " alpha=0.6,\n", + " ls=\"--\",\n", + " label=\"$\\hat{g}$\",\n", + " )\n", " plt.ylim(-4, 0)\n", " plt.xlim(0, 1)\n", - " plt.legend(loc='lower left')\n", - "# plt.subplot(122)\n", - "# plt.title('Q-Q Plot', fontsize=1.25*fsize)\n", - "# for i in solutions_mult[num_plot_sensors]: # trials\n", - "# q = qoi[i,:]\n", - "# plt.scatter(q[:num_plot_sensors], qoi_ref[:num_plot_sensors] + noise[:num_plot_sensors],\n", - "# c='b', s=100, alpha=1.0/num_trials)\n", - "# plt.plot(_a,_a, c='k', lw=3)\n", - "# plt.ylabel('Collected Data', fontsize=fsize)\n", - "# plt.xlabel('Predicted Data', fontsize=fsize)\n", - "# plt.ylim(-0.5, 0.15)\n", - "# plt.xlim(-0.5, 0.15)\n", - "# plt.title(f'Solution {_r}, Index {_s}')\n", + " plt.legend(loc=\"lower left\")\n", + " # plt.subplot(122)\n", + " # plt.title('Q-Q Plot', fontsize=1.25*fsize)\n", + " # for i in solutions_mult[num_plot_sensors]: # trials\n", + " # q = qoi[i,:]\n", + " # plt.scatter(q[:num_plot_sensors], qoi_ref[:num_plot_sensors] + noise[:num_plot_sensors],\n", + " # c='b', s=100, alpha=1.0/num_trials)\n", + " # plt.plot(_a,_a, c='k', lw=3)\n", + " # plt.ylabel('Collected Data', fontsize=fsize)\n", + " # plt.xlabel('Predicted Data', fontsize=fsize)\n", + " # plt.ylim(-0.5, 0.15)\n", + " # plt.xlim(-0.5, 0.15)\n", + " # plt.title(f'Solution {_r}, Index {_s}')\n", " _fname = f\"{prefix}_mud_{input_dim}-{num_qoi}_N{num_plot_sensors}.{ftype}\"\n", - " plt.savefig(_fname, bbox_inches='tight')\n", + " plt.savefig(_fname, bbox_inches=\"tight\")\n", " plt.show()" ] }, @@ -1799,8 +1979,19 @@ "source": [ "if not load:\n", " print(\"Saving\")\n", - " fname = f'{prefix}_summary_{input_dim}.pkl'\n", - " pickle.dump({'sets': (lam, qoi), 'sols': (solutions_sing, solutions_mult), 'meas': measurements, 'noise': (noise, tolerance), 'stdv': sigma, 'true': (lam_ref, qoi_ref), 'sens': sensors }, open(fname, 'wb'))" + " fname = f\"{prefix}_summary_{input_dim}.pkl\"\n", + " pickle.dump(\n", + " {\n", + " \"sets\": (lam, qoi),\n", + " \"sols\": (solutions_sing, solutions_mult),\n", + " \"meas\": measurements,\n", + " \"noise\": (noise, tolerance),\n", + " \"stdv\": sigma,\n", + " \"true\": (lam_ref, qoi_ref),\n", + " \"sens\": sensors,\n", + " },\n", + " open(fname, \"wb\"),\n", + " )" ] }, { @@ -1819,16 +2010,25 @@ "source": [ "# VECTOR\n", "def mud_wrapper(num_obs, sd):\n", - " qois = split_qoi_by_indices(qoi_indices_bad, qoi_ref, qoi,\n", - " noise=np.random.randn(num_obs)*sd, sigma=sd, max_index=num_obs)\n", + " qois = split_qoi_by_indices(\n", + " qoi_indices_bad,\n", + " qoi_ref,\n", + " qoi,\n", + " noise=np.random.randn(num_obs) * sd,\n", + " sigma=sd,\n", + " max_index=num_obs,\n", + " )\n", " r_mult = ratio_dci_mult(qois)\n", " return r_mult\n", "\n", - "experiments_mult_bad, solutions_mult_bad = experiment_measurements_index(num_measurements=measurements,\n", - " sd=sigma,\n", - " num_trials=num_trials,\n", - " seed=21,\n", - " fun=mud_wrapper)" + "\n", + "experiments_mult_bad, solutions_mult_bad = experiment_measurements_index(\n", + " num_measurements=measurements,\n", + " sd=sigma,\n", + " num_trials=num_trials,\n", + " seed=21,\n", + " fun=mud_wrapper,\n", + ")" ] }, { @@ -1838,35 +2038,54 @@ "outputs": [], "source": [ "for num_plot_sensors in [20, 100]:\n", - " plt.figure(figsize=(10,10))\n", - " \n", + " plt.figure(figsize=(10, 10))\n", + "\n", " plt.subplot(111)\n", - " plt.title('MUD Estimates for $Q_{2D}^\\prime$,' + f' S={num_plot_sensors}', fontsize=1.25*fsize)\n", + " plt.title(\n", + " \"MUD Estimates for $Q_{2D}^\\prime$,\" + f\" S={num_plot_sensors}\",\n", + " fontsize=1.25 * fsize,\n", + " )\n", " plt.xlabel(\"$x_1$\", fontsize=fsize)\n", " plt.ylabel(\"$g(x, \\lambda)$\", fontsize=fsize)\n", - " plt.plot([0]+intervals+[1], [0]+ [w(i) for i in intervals] +[0], lw=5, c='k',label=\"Interpolant\")\n", - " for i in solutions_mult_bad[num_plot_sensors]: # trials\n", - " gt = lam[i,:]\n", - " plt.plot([0]+intervals+[1], [0]+list(gt)+[0], lw=1, c='purple', alpha=0.2)\n", + " plt.plot(\n", + " [0] + intervals + [1],\n", + " [0] + [w(i) for i in intervals] + [0],\n", + " lw=5,\n", + " c=\"k\",\n", + " label=\"Interpolant\",\n", + " )\n", + " for i in solutions_mult_bad[num_plot_sensors]: # trials\n", + " gt = lam[i, :]\n", + " plt.plot(\n", + " [0] + intervals + [1], [0] + list(gt) + [0], lw=1, c=\"purple\", alpha=0.2\n", + " )\n", " gt = list(lam[closest_fit_index_out, :])\n", - " plt.plot([0]+intervals+[1], [0]+gt+[0], lw=5, c='green', alpha=0.6, ls='--', label=f'Closest in Output: {closest_fit_index_out}')\n", + " plt.plot(\n", + " [0] + intervals + [1],\n", + " [0] + gt + [0],\n", + " lw=5,\n", + " c=\"green\",\n", + " alpha=0.6,\n", + " ls=\"--\",\n", + " label=f\"Closest in Output: {closest_fit_index_out}\",\n", + " )\n", " plt.ylim(-4, 0)\n", " plt.xlim(0, 1)\n", - " plt.legend(loc='lower left')\n", - "# plt.subplot(122)\n", - "# plt.title('Q-Q Plot', fontsize=1.25*fsize)\n", - "# for i in solutions_mult[num_plot_sensors]: # trials\n", - "# q = qoi[i,:]\n", - "# plt.scatter(q[:num_plot_sensors], qoi_ref[:num_plot_sensors] + noise[:num_plot_sensors],\n", - "# c='b', s=100, alpha=1.0/num_trials)\n", - "# plt.plot(_a,_a, c='k', lw=3)\n", - "# plt.ylabel('Collected Data', fontsize=fsize)\n", - "# plt.xlabel('Predicted Data', fontsize=fsize)\n", - "# plt.ylim(-0.5, 0.15)\n", - "# plt.xlim(-0.5, 0.15)\n", - "# plt.title(f'Solution {_r}, Index {_s}')\n", + " plt.legend(loc=\"lower left\")\n", + " # plt.subplot(122)\n", + " # plt.title('Q-Q Plot', fontsize=1.25*fsize)\n", + " # for i in solutions_mult[num_plot_sensors]: # trials\n", + " # q = qoi[i,:]\n", + " # plt.scatter(q[:num_plot_sensors], qoi_ref[:num_plot_sensors] + noise[:num_plot_sensors],\n", + " # c='b', s=100, alpha=1.0/num_trials)\n", + " # plt.plot(_a,_a, c='k', lw=3)\n", + " # plt.ylabel('Collected Data', fontsize=fsize)\n", + " # plt.xlabel('Predicted Data', fontsize=fsize)\n", + " # plt.ylim(-0.5, 0.15)\n", + " # plt.xlim(-0.5, 0.15)\n", + " # plt.title(f'Solution {_r}, Index {_s}')\n", " _fname = f\"{prefix}_mud-alt_{input_dim}-{num_qoi}_N{num_plot_sensors}.{ftype}\"\n", - "# plt.savefig(_fname, bbox_inches='tight')\n", + " # plt.savefig(_fname, bbox_inches='tight')\n", " plt.show()" ] }, diff --git a/notebooks/MUD_SVD.ipynb b/notebooks/MUD_SVD.ipynb index 101a4d5..0f8a6ba 100644 --- a/notebooks/MUD_SVD.ipynb +++ b/notebooks/MUD_SVD.ipynb @@ -59,7 +59,9 @@ "outputs": [], "source": [ "num_samples, num_max_qoi = P.qoi.shape\n", - "print(f\"There are {num_max_qoi} sensors available and {num_samples} samples of {input_dim}-D parameter space.\")" + "print(\n", + " f\"There are {num_max_qoi} sensors available and {num_samples} samples of {input_dim}-D parameter space.\"\n", + ")" ] }, { @@ -106,7 +108,7 @@ "metadata": {}, "outputs": [], "source": [ - "mud_problem = mud_generator(500, 0.01) # simulates noisy measurements" + "mud_problem = mud_generator(500, 0.01) # simulates noisy measurements" ] }, { @@ -145,9 +147,7 @@ "outputs": [], "source": [ "# TODO : turn this snippet into a method\n", - "closest_fit_index_out = np.argmin(\n", - " np.linalg.norm(P.qoi - np.array(P.qoi_ref), axis=1)\n", - ")\n", + "closest_fit_index_out = np.argmin(np.linalg.norm(P.qoi - np.array(P.qoi_ref), axis=1))\n", "g_projected = P.lam[closest_fit_index_out, :].ravel()\n", "lam_true = g_projected\n", "print(lam_true)" @@ -226,12 +226,12 @@ "outputs": [], "source": [ "top_singular_values = 5\n", - "inds = np.arange(1, top_singular_values+1)\n", + "inds = np.arange(1, top_singular_values + 1)\n", "plt.plot(inds, singular_values[0:top_singular_values])\n", "plt.xticks(inds)\n", "plt.xlabel(\"Index\", fontsize=12)\n", "plt.ylabel(\"$\\sigma_i$\", fontsize=12)\n", - "plt.yscale('log')\n", + "plt.yscale(\"log\")\n", "plt.title(f\"Top {top_singular_values} Singular Values of Scaled Residuals\", fontsize=16)\n", "plt.show()" ] @@ -255,13 +255,27 @@ "outputs": [], "source": [ "top_singular_vectors = 2\n", - "fig, ax = plt.subplots(1, top_singular_vectors, figsize=(1+5*top_singular_vectors, 5), sharey=True)\n", + "fig, ax = plt.subplots(\n", + " 1, top_singular_vectors, figsize=(1 + 5 * top_singular_vectors, 5), sharey=True\n", + ")\n", "for i in range(top_singular_vectors):\n", - " ax[i].tricontour(P.sensors[0:num_obs, 0], P.sensors[0:num_obs, 1], singular_vectors[i, :], levels=20, alpha=1)\n", - " ax[i].scatter(P.sensors[0:num_obs, 0], P.sensors[0:num_obs, 1], s=100, c=singular_vectors[i, :])\n", + " ax[i].tricontour(\n", + " P.sensors[0:num_obs, 0],\n", + " P.sensors[0:num_obs, 1],\n", + " singular_vectors[i, :],\n", + " levels=20,\n", + " alpha=1,\n", + " )\n", + " ax[i].scatter(\n", + " P.sensors[0:num_obs, 0],\n", + " P.sensors[0:num_obs, 1],\n", + " s=100,\n", + " c=singular_vectors[i, :],\n", + " )\n", " ax[i].set_title(f\"$\\sigma_{i+1}$ = {singular_values[i]:1.2e}\")\n", " ax[i].set_xlabel(\"$x_1$\", fontsize=12)\n", - " if i == 0: ax[i].set_ylabel(\"$x_2$\", fontsize=12)\n", + " if i == 0:\n", + " ax[i].set_ylabel(\"$x_2$\", fontsize=12)\n", "fig.suptitle(\"Singular Vector Components as Weights for Sensors\", fontsize=16)\n", "fig.show()" ] @@ -273,10 +287,10 @@ "metadata": {}, "outputs": [], "source": [ - "fig, ax = plt.subplots(1, 1, figsize=(5,4))\n", + "fig, ax = plt.subplots(1, 1, figsize=(5, 4))\n", "# plt.plot((singular_vectors[0, :] + singular_vectors[1, :]))\n", - "ax.plot(singular_vectors[0, :], c='xkcd:orange', label=\"$v_1$\", lw=3)\n", - "ax.plot(singular_vectors[1, :], c='xkcd:light blue', label=\"$v_2$\", lw=3, zorder=0)\n", + "ax.plot(singular_vectors[0, :], c=\"xkcd:orange\", label=\"$v_1$\", lw=3)\n", + "ax.plot(singular_vectors[1, :], c=\"xkcd:light blue\", label=\"$v_2$\", lw=3, zorder=0)\n", "ax.legend(fontsize=12)\n", "ax.set_xlabel(\"Sensor Index\", fontsize=12)\n", "ax.set_ylabel(\"Singular Vector Component Value\", fontsize=12)\n", @@ -304,7 +318,9 @@ "source": [ "num_qoi = 2\n", "new_qoi_map = singular_vectors[0:num_qoi, :]\n", - "new_qoi = residuals @ new_qoi_map.T # ok shape, wrong solution, identical to using U@sigma" + "new_qoi = (\n", + " residuals @ new_qoi_map.T\n", + ") # ok shape, wrong solution, identical to using U@sigma" ] }, { @@ -325,14 +341,17 @@ "outputs": [], "source": [ "top_singular_vectors = 2\n", - "fig, ax = plt.subplots(1, top_singular_vectors, figsize=(1+5*top_singular_vectors, 5), sharey=True)\n", + "fig, ax = plt.subplots(\n", + " 1, top_singular_vectors, figsize=(1 + 5 * top_singular_vectors, 5), sharey=True\n", + ")\n", "for i in range(top_singular_vectors):\n", - " ax[i].scatter(P.lam[:,0], P.lam[:,1], s=50, c=new_qoi[:,i])\n", - " ax[i].tricontour(P.lam[:,0], P.lam[:,1], new_qoi[:,i], levels=20)\n", - "# ax[i].set_title(f\"$\\sigma_{i}$ = {singular_values[i]:1.2e}\")\n", + " ax[i].scatter(P.lam[:, 0], P.lam[:, 1], s=50, c=new_qoi[:, i])\n", + " ax[i].tricontour(P.lam[:, 0], P.lam[:, 1], new_qoi[:, i], levels=20)\n", + " # ax[i].set_title(f\"$\\sigma_{i}$ = {singular_values[i]:1.2e}\")\n", " ax[i].set_title(f\"$q_{i+1}$\", fontsize=12)\n", " ax[i].set_xlabel(\"$\\lambda_1$\", fontsize=12)\n", - " if i == 0: ax[i].set_ylabel(\"$\\lambda_2$\", fontsize=12)\n", + " if i == 0:\n", + " ax[i].set_ylabel(\"$\\lambda_2$\", fontsize=12)\n", "# fig.supylabel(\"$\\lambda_2$\")\n", "fig.suptitle(\"QoI Component Surface Plot\", fontsize=16)\n", "fig.show()" @@ -357,8 +376,20 @@ "outputs": [], "source": [ "P.plot()\n", - "plt.plot(np.linspace(0, 1, input_dim + 2), [0] + list(d.mud_point()) + [0], lw=5, c='blue', label='svd')\n", - "plt.plot(np.linspace(0, 1, input_dim + 2), [0] + list(mud_problem.mud_point()) + [0], lw=5, c='purple', label='split')\n", + "plt.plot(\n", + " np.linspace(0, 1, input_dim + 2),\n", + " [0] + list(d.mud_point()) + [0],\n", + " lw=5,\n", + " c=\"blue\",\n", + " label=\"svd\",\n", + ")\n", + "plt.plot(\n", + " np.linspace(0, 1, input_dim + 2),\n", + " [0] + list(mud_problem.mud_point()) + [0],\n", + " lw=5,\n", + " c=\"purple\",\n", + " label=\"split\",\n", + ")\n", "plt.legend()\n", "plt.show()" ] diff --git a/scripts/process_data.py b/scripts/process_data.py index 6c11986..5dfc17b 100644 --- a/scripts/process_data.py +++ b/scripts/process_data.py @@ -3,12 +3,12 @@ num_sensors = 500 -dist = 'u' +dist = "u" for input_dim in [2, 5, 11]: fname = mf(input_dim=input_dim, sample_dist=dist, num_measure=num_sensors) print(fname) -dist = 'n' +dist = "n" input_dim = 2 for tol in [0.9999, 0.99, 0.95]: fname = mf(input_dim=input_dim, sample_dist=dist, tol=tol, num_measure=num_sensors) diff --git a/setup.cfg b/setup.cfg index bee05ab..e6bf59e 100644 --- a/setup.cfg +++ b/setup.cfg @@ -30,7 +30,7 @@ setup_requires = setuptools_scm[toml]>=5 # Add here dependencies of your project (semicolon/line-separated), e.g. install_requires = importlib-metadata; python_version<"3.10" - mud~=0.0.24 + mud>=0.0.24,!=0.1.0 numpy scipy matplotlib diff --git a/src/mud_examples/linear/lin.py b/src/mud_examples/linear/lin.py index eda0acd..0c0e00a 100644 --- a/src/mud_examples/linear/lin.py +++ b/src/mud_examples/linear/lin.py @@ -927,7 +927,7 @@ def contour_example( inputs[:, 1].reshape(N, N), (zd).reshape(N, N), 25, - cmap=cm.viridis, + cmap="viridis", alpha=0.5, vmin=0, vmax=4, @@ -957,7 +957,7 @@ def contour_example( inputs[:, 1].reshape(N, N), z.reshape(N, N), 50, - cmap=cm.viridis, + cmap="viridis", alpha=1.0, ) elif alpha + omega > 0: @@ -966,7 +966,7 @@ def contour_example( inputs[:, 1].reshape(N, N), (alpha * zi - omega * zp).reshape(N, N), 100, - cmap=cm.viridis, + cmap="viridis", alpha=0.25, ) plt.axis("equal") @@ -1124,7 +1124,7 @@ def contour_example( new_line[0, mn:mx], new_line[1, mn:mx], lw=1, - label="projection line", + label="Projection Line", c="k", ) elif show_full: diff --git a/src/mud_examples/monomial.py b/src/mud_examples/monomial.py index c0b3969..124bc8b 100644 --- a/src/mud_examples/monomial.py +++ b/src/mud_examples/monomial.py @@ -85,7 +85,7 @@ def main(args): np.random.seed( 123456 ) # Just for reproducibility, you can comment out if you want. - data = norm.rvs(loc=mu, scale=sigma ** 2, size=num_data) + data = norm.rvs(loc=mu, scale=sigma**2, size=num_data) # We will estimate the observed distribution using a parametric estimate to keep # the assumptions involved as similar as possible between the BIP and the SIP @@ -189,7 +189,7 @@ def QoI(lam, p): """ Defines a QoI mapping function as monomials to some power p """ - q = lam ** p + q = lam**p return q diff --git a/src/mud_examples/plotting.py b/src/mud_examples/plotting.py index 60126e8..7f49883 100644 --- a/src/mud_examples/plotting.py +++ b/src/mud_examples/plotting.py @@ -276,10 +276,10 @@ def plot_scalar_poisson_summary( def plotChain(mud_chain, ref_param, color="k", s=100): num_steps = len(mud_chain) - current_point = mud_chain[0].reshape(-1,1) + current_point = mud_chain[0].reshape(-1, 1) plt.scatter(current_point[0], current_point[1], c="b", s=s) for i in range(0, num_steps): - next_point = mud_chain[i].reshape(-1,1) + next_point = mud_chain[i].reshape(-1, 1) points = np.hstack([current_point, next_point]) plt.plot(points[0, :], points[1, :], c=color) current_point = next_point