7607107083

Committed 22 Jan 2024 06:01AM UTC coverage: 75.216% (+2.5%) from 72.738%

Build # 7607107083

Build Type

Pull #71

github

Committed by

web-flow

Commit Message

Merge 7d7cd3b4f into a134ca7cc

Pull Request Pull Request #71: housekeeping

Run Details

664 of 944 branches covered (0.0%)

Branch coverage included in aggregate %.

160 of 164 new or added lines in 6 files covered. (97.56%)

5 existing lines in 3 files now uncovered.

1597 of 2062 relevant lines covered (77.45%)

3.09 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

84.0

/src/mud/examples/poisson.py

"""
Poisson Problem Example

The functions here implement the poisson problem example found in [ref] section
6.1.
"""
import logging
import random
from typing import List, Optional, Union

import matplotlib.pyplot as plt  # type: ignore
import numpy as np
import pandas as pd  # type: ignore
from scipy.interpolate import interp1d  # type: ignore
from scipy.optimize import minimize  # type: ignore

from mud.base import SpatioTemporalProblem
from mud.plot import mud_plot_params, plot_vert_line, save_figure

_logger = logging.getLogger(__name__)

# Matplotlib plotting options
plt.rcParams.update(mud_plot_params)


def load_poisson_prob(
    prob_data_file,
    std_dev=0.05,
    seed=21,
):
    """
    Load Poisson Problem Dataset


    Parameters
    ----------
    prob_data_file : pathlib.Path, str
        Path to pickled data file storing data from solving poisson problem.
    std_dev : float, default = 0.05
        Standard deviation of noise to add to true values stored in dataset
        to create fake "measurements".
    seed : int, default = 21
        For reproducible results, leave seed as is or set to a value. Otherwise,
        set to None.

    Returns
    -------
    data, poisson_prob : Tuple(dict, mud.base.SpatioTemporalProblem)
        Returns a pair of the raw dictionary loaded from the pickled dataset
        and the associated SpatioTemporalProblem class object associated with it.
    """

    poisson_prob = SpatioTemporalProblem()
    data = poisson_prob.load(
        str(prob_data_file),
        lam="lam",
        data="data",
        true_vals="true_vals",
        measurements=None,
        std_dev=0.05,
        sample_dist="uniform",
        domain="domain",
        lam_ref=None,
        sensors="sensors",
        times=np.array([0]),
    )

    poisson_prob.measurements_from_reference(std_dev=0.05, seed=seed)

    return data, poisson_prob


def plot_solution_spline(lam, aff=1000, plot_true=True, ax=None, **kwargs):
    """
    Poisson Problem Plot Boundary Solution Spline

    Parameters
    ----------

    Returns
    ----------
    """

    if ax is None:
        fig = plt.figure(figsize=(5, 5))
        ax = fig.add_subplot(1, 1, 1)

    if plot_true:
        x = np.linspace(0, 1, aff)
        a = 197.65032
        g = a * np.power(x, 2) * np.power((x - 1), 5)
        ax.plot(x, g, lw=5, c="k", label="$g$")

    intervals = list(np.linspace(0, 1, 2 + len(lam))[1:-1])
    ax.plot([0] + intervals + [1], [0] + list(lam) + [0], **kwargs)

    ax.set_xlabel("$x_2$")
    ax.set_ylabel(r"$g(x_2,\lambda)$")
    ax.set_xlim([0, 1])
    ax.set_ylim([-4, 0])

    return ax


def spline_objective_function_2d(lam, aff=10000):
    """
    Spline Objective Function

    Parameters
    ----------
    lam1: float
        $lambda_1$ parameter value - location of first knot at x=1/3.
    lam2: float
        $lambda_2$ parameter value - location of second knot not at x=2/3.

    Returns
    -------

    """
    x = np.linspace(0, 1, aff)
    a = 197.65032
    g = a * np.power(x, 2) * np.power((x - 1), 5)
    spline = interp1d([0, 1 / 3, 2 / 3, 1], [0, lam[0], lam[1], 0])
    vals = spline(x)

    return np.linalg.norm(g - vals)


def plot_qoi_ex(mud_prob, num_components: int, **kwargs):
    fig = plt.figure(figsize=(10, 5))
    for i in range(num_components):
        ax = fig.add_subplot(1, 2, i + 1)
        mud_prob.plot_params_2d(ax=ax, y=i, contours=True, colorbar=True)
        if i == 1:
            ax.set_ylabel("")
    if kwargs.get("save_path"):
        save_figure("learned_qoi", **kwargs)
    return ax


def plot_response_ex(
    poisson_prob,
    closest,
    raw_data,
    sensors_mask,
    group_idxs: List[List[int]] = [[0, 5], [5, 50], [50, -1]],
    markers: List[str] = ["*", "+", "."],
    colors: List[str] = ["red", "white", "k"],
    **kwargs,
):
    fig = plt.figure(figsize=(10, 5))
    ax = fig.add_subplot(1, 2, 1)

    # Plot response surface from solving Eq. 30.
    # u stores the mesh and values as solved by fenics
    mesh, vals = raw_data["u"]
    tcf = ax.tricontourf(mesh[:, 0], mesh[:, 1], vals, levels=20, vmin=-0.5, vmax=0)
    fig.colorbar(tcf)

    # Plot points used for each ordering
    for idx, oi in enumerate(group_idxs):
        start_idx, end_idx = group_idxs[idx][0], group_idxs[idx][1]
        mask = sensors_mask[start_idx:end_idx]
        poisson_prob.sensor_scatter_plot(
            ax=ax,
            mask=mask,
            color=colors[idx],
            marker=markers[idx],
        )

    # Label and format figure
    _ = plt.xlim(0, 1)
    _ = plt.ylim(0, 1)
    _ = plt.title("Response Surface")
    _ = plt.xlabel("$x_1$")
    _ = plt.ylabel("$x_2$")

    ax = fig.add_subplot(1, 2, 2)
    # Plot closest solution in sample set to the reference solution
    plot_solution_spline(
        closest,
        ax=ax,
        lw=5,
        c="green",
        ls="--",
        label=r"$\hat{g}(\lambda^\dagger)$",
        zorder=50,
    )
    # Plot first 100 lambda initial
    for i, lam in enumerate(poisson_prob.lam[0:50]):
        plot_solution_spline(
            lam,
            plot_true=False,
            ax=ax,
            lw=1,
            c="purple",
            alpha=0.1,
        )

    plot_vert_line(ax, 1 / 3.0, color="b", linestyle="--", alpha=0.6)
    plot_vert_line(ax, 2 / 3.0, color="b", linestyle="--", alpha=0.6)

    ax.set_title("Boundary Condition")
    ax.set_ylabel("")
    ax.legend(
        ["$g(x_2)$", r"$\hat{g}(x_2,\lambda^\dagger)$", r"$\hat{g}(x_2,\lambda_i)$"]
    )

    fig.tight_layout()
    if kwargs.get("save_path"):
        save_figure("response_surface", **kwargs)
    return ax


# TODO: Document group_idxs, markers, colors
def run_2d_poisson_sol(
    data_file: str,
    sigma: float = 0.05,
    seed: Optional[int] = None,
    plot_fig: Union[List[str], str] = "all",
    save_path: Optional[str] = None,
    dpi: int = 500,
    close_fig: bool = False,
    order: str = "random",
    group_idxs: List[List[int]] = [[0, 5], [5, 50], [50, -1]],
    markers: List[str] = ["*", "+", "."],
    colors: List[str] = ["red", "white", "k"],
    param1_kwargs: dict = {},
    param2_kwargs: dict = {},
):
    """
    Run 2D Poisson problem

    Parameters
    ----------
    data_file : str
        Path to pickled data file storing data from solving poisson problem.
    sigma: float, default=0.05
        N(0, sigma) error added to true solution surface to produce
        "measurements".
    seed: int, default = None
        To fix results for reproducibility. Set to None to randomize results.
    plot_fig : str, default='all'
        Which figures to produce. Possible options are response, qoi, or all to
        produce all figures.
    save_path: str, optional
        If provided, path to save the resulting figure to.
    dpi: int, default=500
        Resolution of images saved
    close_fig: bool, default=False
        Set to True to close figure and only save it.

    Returns
    -------
    posson_prob, mud_prob, axes: Tuple[mud.base.SpatioTemporalProblem,
                                       mud.base.DensityProblem,
                                       List[matplotlib.pyplot.Axes]]
        Tuple of SpatioTemporalProblem data structure storing poisson problem
        as configured by inputs, and DensityProblem object used to solve the
        inverse problem, along with the list of axes objects where desired
        plots were plotted on.
    """
    res = minimize(spline_objective_function_2d, x0=[-3, -1])
    closest = res["x"]

    raw_data, poisson_prob = load_poisson_prob(data_file, std_dev=sigma, seed=seed)
    if order == "random":
        idx_o = list(
            random.sample(range(poisson_prob.n_sensors), poisson_prob.n_sensors)
        )
    elif order == "sorted":
        idx_o = list(
            np.lexsort((poisson_prob.sensors[:, 1], poisson_prob.sensors[:, 0]))
        )
    else:
        idx_o = list(np.arange(0, poisson_prob.n_sensors, 1))
    num_components = 2
    mud_prob = poisson_prob.mud_problem(
        method="pca", num_components=num_components, sensors_mask=idx_o
    )
    _ = mud_prob.estimate()

    plot_fig = [plot_fig] if not isinstance(plot_fig, list) else plot_fig

    axes = []
    if "response" in plot_fig or "all" in plot_fig:
        ax = plot_response_ex(
            poisson_prob,
            closest=closest,
            raw_data=raw_data,
            group_idxs=group_idxs,
            markers=markers,
            sensors_mask=idx_o,
            colors=colors,
            save_path=save_path,
            close_fig=close_fig,
            dpi=dpi,
            bbox_inches="tight",
        )
        axes.append(ax)
    if "qoi" in plot_fig or "all" in plot_fig:
        ax = plot_qoi_ex(
            mud_prob,
            num_components=num_components,
            save_path=save_path,
            close_fig=close_fig,
            dpi=dpi,
            bbox_inches="tight",
        )
        axes.append(ax)
    if "densities" in plot_fig or "all" in plot_fig:
        ax1 = mud_prob.plot_param_space(param_idx=0, **param1_kwargs)
        if save_path is not None:
            save_figure(
                "lam1", save_path, close_fig=close_fig, dpi=dpi, bbox_inches="tight"
            )

        ax2 = mud_prob.plot_param_space(param_idx=1, **param2_kwargs)
        if save_path is not None:
            save_figure(
                "lam2", save_path, close_fig=close_fig, dpi=dpi, bbox_inches="tight"
            )
        axes.append([ax1, ax2])

    return (poisson_prob, mud_prob, axes)


def run_2d_poisson_trials(
    data_file: pd.DataFrame,
    N_vals: List[int] = [5, 50, 500],
    order: str = "random",
    ylim1: List[float] = [-0.1, 3.5],
    ylim2: List[float] = [-0.1, 2.5],
    xlim1: List[float] = [-4.5, 0.5],
    xlim2: List[float] = [-4.5, 0.5],
    annotate_location_1: List[float] = [-2.6, 1.2, 0.8],
    annotate_location_2: List[float] = [-3.5, 0.83, 0.53],
    sigma: float = 0.05,
    seed: Optional[int] = None,
    save_path: Optional[str] = None,
    dpi: int = 500,
    close_fig: bool = False,
):
    """
    Run Poisson Problem 2D Solution Trials

    Run trials for solving Poisson boundary estimation problem using a 2D linear
    spline with increasing number of sensors over the solution space.

    Parameters
    ----------
    data_file : str
        Path to pickled data file storing data from solving poisson problem.
    N_vals: List[int], default=[5, 50, 500]
        Number of sensors to use in each trial. Each N value will produce two
        figures, one for the updated parameter distribution and MUD solution for
        each parameter value.
    ylim1: List[float], default=[-0.1, 3.5]
        y limits to use in plot for first parameter space. Set to something
        consistent so that plots for all N values will be clearly visible.
    ylim2: List[float], default=[-0.1, 2.5]
        y limits to use in plot for second parameter space. Set to something
        consistent so that plots for all N values will be clearly visible.
    xlim1: List[float], default=[-4.5, 0.5]
        x limits to use in plot for second parameter space. Set to something
        consistent so that plots for all N values will be clearly visible.
    xlim2: List[float], default=[-4.5, 0.5]
        x limits to use in plot for second parameter space. Set to something
        consistent so that plots for all N values will be clearly visible.
    sigma: float, default=1e-1
        N(0, sigma) error added to produce "measurements" from linear operator.
    seed: int, default = None
        To fix results for reproducibility. Set to None to randomize results.
    save_path: str, optional
        If provided, path to save the resulting figure to.
    dpi: int, default=500
        Resolution of images saved
    close_fig: bool, default=False
        Set to True to close figure and only save it.

    Returns
    ----------
    poisson_prob, probs, axes : Tuple[mud.base.SpatioTemporalProblem,
                                     List[mud.base.DensityProblem],
                                     List[Tuple[matplotlib.pyplot.Axes,
                                                matplotlib.pyplot.Axes]]]
        Tuple of poisson problem from loaded data, list of mud problems solved
        for each N value, and the pair of axes objects that each mud problem
        was plotted to.
    """
    if seed is not None:
        random.seed(seed)

    res = minimize(spline_objective_function_2d, x0=[-3, -1])
    closest = res["x"]
    raw_data, poisson_prob = load_poisson_prob(data_file, std_dev=sigma, seed=seed)
    x_range = np.array([xlim1, xlim2])
    axes = []
    probs = []
    if order == "random":
        idx_o = list(
            random.sample(range(poisson_prob.n_sensors), poisson_prob.n_sensors)
        )
    elif order == "sorted":
        idx_o = list(
            np.lexsort((poisson_prob.sensors[:, 1], poisson_prob.sensors[:, 0]))
        )
    else:
        idx_o = list(np.arange(0, poisson_prob.n_sensors, 1))
    for N in N_vals:
        mud_prob = poisson_prob.mud_problem(
            method="pca", num_components=2, sensors_mask=idx_o[0:N]
        )
        _ = mud_prob.estimate()

        fig = plt.figure(figsize=(9, 4))

        ax1 = fig.add_subplot(1, 2, 1)
        mud_prob.plot_param_space(
            ax=ax1, x_range=x_range, param_idx=0, mud_opts=None, true_opts=None
        )
        ax1.set_ylim(*ylim1)
        mud_prob.plot_param_space(ax=ax1, true_val=closest, in_opts=None, up_opts=None)
        ax1.set_xlabel(r"$\lambda_1$")
        # annotate_location_1 = [-2.8, 1.2, 0.8]
        if annotate_location_1 is not None:
            x = annotate_location_1[0]
            y1 = annotate_location_1[1]
            y2 = annotate_location_1[2]
            ax1.text(x, y1, f"$N = {N}$", fontsize=18)
            e_r = mud_prob.expected_ratio()
            ax1.text(x, y2, rf"$\mathbb{{E}}(r) = {e_r:0.4}$", fontsize=18)
        ax1.legend()

        ax2 = fig.add_subplot(1, 2, 2)
        mud_prob.plot_param_space(
            ax=ax2,
            x_range=x_range,
            param_idx=1,
            mud_opts=None,
            true_opts=None,
        )
        ax2.set_ylim(*ylim2)
        mud_prob.plot_param_space(
            ax=ax2, param_idx=1, true_val=closest, in_opts=None, up_opts=None
        )
        ax2.set_xlabel(r"$\lambda_2$")
        # annotate_location_2 = [-3.5, 0.83, 0.53]
        if annotate_location_2 is not None:
            x = annotate_location_2[0]
            y1 = annotate_location_2[1]
            y2 = annotate_location_2[2]
            ax2.text(x, y1, f"$N = {N}$", fontsize=18)
            e_r = mud_prob.expected_ratio()
            ax2.text(x, y2, rf"$\mathbb{{E}}(r) = {e_r:0.4}$", fontsize=18)
        ax2.legend()
        if save_path is not None:
            save_figure(
                f"solution_n{N}",
                save_path,
                close_fig=close_fig,
                dpi=dpi,
                bbox_inches="tight",
            )

        axes.append([ax1, ax2])
        probs.append(mud_prob)

    return (poisson_prob, probs, axes)

1	"""
2	Poisson Problem Example
3
4	The functions here implement the poisson problem example found in [ref] section
5	6.1.
6	"""
7	import logging	4✔
8	import random	4✔
9	from typing import List, Optional, Union	4✔
10
11	import matplotlib.pyplot as plt # type: ignore	4✔
12	import numpy as np	4✔
13	import pandas as pd # type: ignore	4✔
14	from scipy.interpolate import interp1d # type: ignore	4✔
15	from scipy.optimize import minimize # type: ignore	4✔
16
17	from mud.base import SpatioTemporalProblem	4✔
18	from mud.plot import mud_plot_params, plot_vert_line, save_figure	4✔
19
20	_logger = logging.getLogger(__name__)	4✔
21
22	# Matplotlib plotting options
23	plt.rcParams.update(mud_plot_params)	4✔
24
25
26	def load_poisson_prob(	4✔
27	prob_data_file,
28	std_dev=0.05,
29	seed=21,
30	):
31	"""
32	Load Poisson Problem Dataset
33
34
35	Parameters
36	----------
37	prob_data_file : pathlib.Path, str
38	Path to pickled data file storing data from solving poisson problem.
39	std_dev : float, default = 0.05
40	Standard deviation of noise to add to true values stored in dataset
41	to create fake "measurements".
42	seed : int, default = 21
43	For reproducible results, leave seed as is or set to a value. Otherwise,
44	set to None.
45
46	Returns
47	-------
48	data, poisson_prob : Tuple(dict, mud.base.SpatioTemporalProblem)
49	Returns a pair of the raw dictionary loaded from the pickled dataset
50	and the associated SpatioTemporalProblem class object associated with it.
51	"""
52
53	poisson_prob = SpatioTemporalProblem()	4✔
54	data = poisson_prob.load(	4✔
55	str(prob_data_file),
56	lam="lam",
57	data="data",
58	true_vals="true_vals",
59	measurements=None,
60	std_dev=0.05,
61	sample_dist="uniform",
62	domain="domain",
63	lam_ref=None,
64	sensors="sensors",
65	times=np.array([0]),
66	)
67
68	poisson_prob.measurements_from_reference(std_dev=0.05, seed=seed)	4✔
69
70	return data, poisson_prob	4✔
71
72
73	def plot_solution_spline(lam, aff=1000, plot_true=True, ax=None, **kwargs):	4✔
74	"""
75	Poisson Problem Plot Boundary Solution Spline
76
77	Parameters
78	----------
79
80	Returns
81	----------
82	"""
83
84	if ax is None:	4!
85	fig = plt.figure(figsize=(5, 5))	×
86	ax = fig.add_subplot(1, 1, 1)	×
87
88	if plot_true:	4✔
89	x = np.linspace(0, 1, aff)	4✔
90	a = 197.65032	4✔
91	g = a * np.power(x, 2) * np.power((x - 1), 5)	4✔
92	ax.plot(x, g, lw=5, c="k", label="$g$")	4✔
93
94	intervals = list(np.linspace(0, 1, 2 + len(lam))[1:-1])	4✔
95	ax.plot([0] + intervals + [1], [0] + list(lam) + [0], **kwargs)	4✔
96
97	ax.set_xlabel("$x_2$")	4✔
98	ax.set_ylabel(r"$g(x_2,\lambda)$")	4✔
99	ax.set_xlim([0, 1])	4✔
100	ax.set_ylim([-4, 0])	4✔
101
102	return ax	4✔
103
104
105	def spline_objective_function_2d(lam, aff=10000):	4✔
106	"""
107	Spline Objective Function
108
109	Parameters
110	----------
111	lam1: float
112	$lambda_1$ parameter value - location of first knot at x=1/3.
113	lam2: float
114	$lambda_2$ parameter value - location of second knot not at x=2/3.
115
116	Returns
117	-------
118
119	"""
120	x = np.linspace(0, 1, aff)	4✔
121	a = 197.65032	4✔
122	g = a * np.power(x, 2) * np.power((x - 1), 5)	4✔
123	spline = interp1d([0, 1 / 3, 2 / 3, 1], [0, lam[0], lam[1], 0])	4✔
124	vals = spline(x)	4✔
125
126	return np.linalg.norm(g - vals)	4✔
127
128
129	def plot_qoi_ex(mud_prob, num_components: int, **kwargs):	4✔
130	fig = plt.figure(figsize=(10, 5))	4✔
131	for i in range(num_components):	4✔
132	ax = fig.add_subplot(1, 2, i + 1)	4✔
133	mud_prob.plot_params_2d(ax=ax, y=i, contours=True, colorbar=True)	4✔
134	if i == 1:	4✔
135	ax.set_ylabel("")	4✔
136	if kwargs.get("save_path"):	4!
137	save_figure("learned_qoi", **kwargs)	4✔
138	return ax	4✔
139
140
141	def plot_response_ex(	4✔
142	poisson_prob,
143	closest,
144	raw_data,
145	sensors_mask,
146	group_idxs: List[List[int]] = [[0, 5], [5, 50], [50, -1]],
147	markers: List[str] = ["*", "+", "."],
148	colors: List[str] = ["red", "white", "k"],
149	**kwargs,
150	):
151	fig = plt.figure(figsize=(10, 5))	4✔
152	ax = fig.add_subplot(1, 2, 1)	4✔
153
154	# Plot response surface from solving Eq. 30.
155	# u stores the mesh and values as solved by fenics
156	mesh, vals = raw_data["u"]	4✔
157	tcf = ax.tricontourf(mesh[:, 0], mesh[:, 1], vals, levels=20, vmin=-0.5, vmax=0)	4✔
158	fig.colorbar(tcf)	4✔
159
160	# Plot points used for each ordering
161	for idx, oi in enumerate(group_idxs):	4✔
162	start_idx, end_idx = group_idxs[idx][0], group_idxs[idx][1]	4✔
163	mask = sensors_mask[start_idx:end_idx]	4✔
164	poisson_prob.sensor_scatter_plot(	4✔
165	ax=ax,
166	mask=mask,
167	color=colors[idx],
168	marker=markers[idx],
169	)
170
171	# Label and format figure
172	_ = plt.xlim(0, 1)	4✔
173	_ = plt.ylim(0, 1)	4✔
174	_ = plt.title("Response Surface")	4✔
175	_ = plt.xlabel("$x_1$")	4✔
176	_ = plt.ylabel("$x_2$")	4✔
177
178	ax = fig.add_subplot(1, 2, 2)	4✔
179	# Plot closest solution in sample set to the reference solution
180	plot_solution_spline(	4✔
181	closest,
182	ax=ax,
183	lw=5,
184	c="green",
185	ls="--",
186	label=r"$\hat{g}(\lambda^\dagger)$",
187	zorder=50,
188	)
189	# Plot first 100 lambda initial
190	for i, lam in enumerate(poisson_prob.lam[0:50]):	4✔
191	plot_solution_spline(	4✔
192	lam,
193	plot_true=False,
194	ax=ax,
195	lw=1,
196	c="purple",
197	alpha=0.1,
198	)
199
200	plot_vert_line(ax, 1 / 3.0, color="b", linestyle="--", alpha=0.6)	4✔
201	plot_vert_line(ax, 2 / 3.0, color="b", linestyle="--", alpha=0.6)	4✔
202
203	ax.set_title("Boundary Condition")	4✔
204	ax.set_ylabel("")	4✔
205	ax.legend(	4✔
206	["$g(x_2)$", r"$\hat{g}(x_2,\lambda^\dagger)$", r"$\hat{g}(x_2,\lambda_i)$"]
207	)
208
209	fig.tight_layout()	4✔
210	if kwargs.get("save_path"):	4!
211	save_figure("response_surface", **kwargs)	4✔
212	return ax	4✔
213
214
215	# TODO: Document group_idxs, markers, colors
216	def run_2d_poisson_sol(	4✔
217	data_file: str,
218	sigma: float = 0.05,
219	seed: Optional[int] = None,
220	plot_fig: Union[List[str], str] = "all",
221	save_path: Optional[str] = None,
222	dpi: int = 500,
223	close_fig: bool = False,
224	order: str = "random",
225	group_idxs: List[List[int]] = [[0, 5], [5, 50], [50, -1]],
226	markers: List[str] = ["*", "+", "."],
227	colors: List[str] = ["red", "white", "k"],
228	param1_kwargs: dict = {},
229	param2_kwargs: dict = {},
230	):
231	"""
232	Run 2D Poisson problem
233
234	Parameters
235	----------
236	data_file : str
237	Path to pickled data file storing data from solving poisson problem.
238	sigma: float, default=0.05
239	N(0, sigma) error added to true solution surface to produce
240	"measurements".
241	seed: int, default = None
242	To fix results for reproducibility. Set to None to randomize results.
243	plot_fig : str, default='all'
244	Which figures to produce. Possible options are response, qoi, or all to
245	produce all figures.
246	save_path: str, optional
247	If provided, path to save the resulting figure to.
248	dpi: int, default=500
249	Resolution of images saved
250	close_fig: bool, default=False
251	Set to True to close figure and only save it.
252
253	Returns
254	-------
255	posson_prob, mud_prob, axes: Tuple[mud.base.SpatioTemporalProblem,
256	mud.base.DensityProblem,
257	List[matplotlib.pyplot.Axes]]
258	Tuple of SpatioTemporalProblem data structure storing poisson problem
259	as configured by inputs, and DensityProblem object used to solve the
260	inverse problem, along with the list of axes objects where desired
261	plots were plotted on.
262	"""
263	res = minimize(spline_objective_function_2d, x0=[-3, -1])	4✔
264	closest = res["x"]	4✔
265
266	raw_data, poisson_prob = load_poisson_prob(data_file, std_dev=sigma, seed=seed)	4✔
267	if order == "random":	4!
268	idx_o = list(	4✔
269	random.sample(range(poisson_prob.n_sensors), poisson_prob.n_sensors)
270	)
271	elif order == "sorted":	×
272	idx_o = list(	×
273	np.lexsort((poisson_prob.sensors[:, 1], poisson_prob.sensors[:, 0]))
274	)
275	else:
276	idx_o = list(np.arange(0, poisson_prob.n_sensors, 1))	×
277	num_components = 2	4✔
278	mud_prob = poisson_prob.mud_problem(	4✔
279	method="pca", num_components=num_components, sensors_mask=idx_o
280	)
281	_ = mud_prob.estimate()	4✔
282
283	plot_fig = [plot_fig] if not isinstance(plot_fig, list) else plot_fig	4✔
284
285	axes = []	4✔
286	if "response" in plot_fig or "all" in plot_fig:	4✔
287	ax = plot_response_ex(	4✔
288	poisson_prob,
289	closest=closest,
290	raw_data=raw_data,
291	group_idxs=group_idxs,
292	markers=markers,
293	sensors_mask=idx_o,
294	colors=colors,
295	save_path=save_path,
296	close_fig=close_fig,
297	dpi=dpi,
298	bbox_inches="tight",
299	)
300	axes.append(ax)	4✔
301	if "qoi" in plot_fig or "all" in plot_fig:	4✔
302	ax = plot_qoi_ex(	4✔
303	mud_prob,
304	num_components=num_components,
305	save_path=save_path,
306	close_fig=close_fig,
307	dpi=dpi,
308	bbox_inches="tight",
309	)
310	axes.append(ax)	4✔
311	if "densities" in plot_fig or "all" in plot_fig:	4!
UNCOV 312	ax1 = mud_prob.plot_param_space(param_idx=0, **param1_kwargs)	×
NEW 313	if save_path is not None:	×
NEW 314	save_figure(	×
315	"lam1", save_path, close_fig=close_fig, dpi=dpi, bbox_inches="tight"
316	)
317
UNCOV 318	ax2 = mud_prob.plot_param_space(param_idx=1, **param2_kwargs)	×
NEW 319	if save_path is not None:	×
NEW 320	save_figure(	×
321	"lam2", save_path, close_fig=close_fig, dpi=dpi, bbox_inches="tight"
322	)
UNCOV 323	axes.append([ax1, ax2])	×
324
325	return (poisson_prob, mud_prob, axes)	4✔
326
327
328	def run_2d_poisson_trials(	4✔
329	data_file: pd.DataFrame,
330	N_vals: List[int] = [5, 50, 500],
331	order: str = "random",
332	ylim1: List[float] = [-0.1, 3.5],
333	ylim2: List[float] = [-0.1, 2.5],
334	xlim1: List[float] = [-4.5, 0.5],
335	xlim2: List[float] = [-4.5, 0.5],
336	annotate_location_1: List[float] = [-2.6, 1.2, 0.8],
337	annotate_location_2: List[float] = [-3.5, 0.83, 0.53],
338	sigma: float = 0.05,
339	seed: Optional[int] = None,
340	save_path: Optional[str] = None,
341	dpi: int = 500,
342	close_fig: bool = False,
343	):
344	"""
345	Run Poisson Problem 2D Solution Trials
346
347	Run trials for solving Poisson boundary estimation problem using a 2D linear
348	spline with increasing number of sensors over the solution space.
349
350	Parameters
351	----------
352	data_file : str
353	Path to pickled data file storing data from solving poisson problem.
354	N_vals: List[int], default=[5, 50, 500]
355	Number of sensors to use in each trial. Each N value will produce two
356	figures, one for the updated parameter distribution and MUD solution for
357	each parameter value.
358	ylim1: List[float], default=[-0.1, 3.5]
359	y limits to use in plot for first parameter space. Set to something
360	consistent so that plots for all N values will be clearly visible.
361	ylim2: List[float], default=[-0.1, 2.5]
362	y limits to use in plot for second parameter space. Set to something
363	consistent so that plots for all N values will be clearly visible.
364	xlim1: List[float], default=[-4.5, 0.5]
365	x limits to use in plot for second parameter space. Set to something
366	consistent so that plots for all N values will be clearly visible.
367	xlim2: List[float], default=[-4.5, 0.5]
368	x limits to use in plot for second parameter space. Set to something
369	consistent so that plots for all N values will be clearly visible.
370	sigma: float, default=1e-1
371	N(0, sigma) error added to produce "measurements" from linear operator.
372	seed: int, default = None
373	To fix results for reproducibility. Set to None to randomize results.
374	save_path: str, optional
375	If provided, path to save the resulting figure to.
376	dpi: int, default=500
377	Resolution of images saved
378	close_fig: bool, default=False
379	Set to True to close figure and only save it.
380
381	Returns
382	----------
383	poisson_prob, probs, axes : Tuple[mud.base.SpatioTemporalProblem,
384	List[mud.base.DensityProblem],
385	List[Tuple[matplotlib.pyplot.Axes,
386	matplotlib.pyplot.Axes]]]
387	Tuple of poisson problem from loaded data, list of mud problems solved
388	for each N value, and the pair of axes objects that each mud problem
389	was plotted to.
390	"""
391	if seed is not None:	4✔
392	random.seed(seed)	4✔
393
394	res = minimize(spline_objective_function_2d, x0=[-3, -1])	4✔
395	closest = res["x"]	4✔
396	raw_data, poisson_prob = load_poisson_prob(data_file, std_dev=sigma, seed=seed)	4✔
397	x_range = np.array([xlim1, xlim2])	4✔
398	axes = []	4✔
399	probs = []	4✔
400	if order == "random":	4!
401	idx_o = list(	4✔
402	random.sample(range(poisson_prob.n_sensors), poisson_prob.n_sensors)
403	)
404	elif order == "sorted":	×
405	idx_o = list(	×
406	np.lexsort((poisson_prob.sensors[:, 1], poisson_prob.sensors[:, 0]))
407	)
408	else:
409	idx_o = list(np.arange(0, poisson_prob.n_sensors, 1))	×
410	for N in N_vals:	4✔
411	mud_prob = poisson_prob.mud_problem(	4✔
412	method="pca", num_components=2, sensors_mask=idx_o[0:N]
413	)
414	_ = mud_prob.estimate()	4✔
415
416	fig = plt.figure(figsize=(9, 4))	4✔
417
418	ax1 = fig.add_subplot(1, 2, 1)	4✔
419	mud_prob.plot_param_space(	4✔
420	ax=ax1, x_range=x_range, param_idx=0, mud_opts=None, true_opts=None
421	)
422	ax1.set_ylim(*ylim1)	4✔
423	mud_prob.plot_param_space(ax=ax1, true_val=closest, in_opts=None, up_opts=None)	4✔
424	ax1.set_xlabel(r"$\lambda_1$")	4✔
425	# annotate_location_1 = [-2.8, 1.2, 0.8]
426	if annotate_location_1 is not None:	4!
427	x = annotate_location_1[0]	4✔
428	y1 = annotate_location_1[1]	4✔
429	y2 = annotate_location_1[2]	4✔
430	ax1.text(x, y1, f"$N = {N}$", fontsize=18)	4✔
431	e_r = mud_prob.expected_ratio()	4✔
432	ax1.text(x, y2, rf"$\mathbb{{E}}(r) = {e_r:0.4}$", fontsize=18)	4✔
433	ax1.legend()	4✔
434
435	ax2 = fig.add_subplot(1, 2, 2)	4✔
436	mud_prob.plot_param_space(	4✔
437	ax=ax2,
438	x_range=x_range,
439	param_idx=1,
440	mud_opts=None,
441	true_opts=None,
442	)
443	ax2.set_ylim(*ylim2)	4✔
444	mud_prob.plot_param_space(	4✔
445	ax=ax2, param_idx=1, true_val=closest, in_opts=None, up_opts=None
446	)
447	ax2.set_xlabel(r"$\lambda_2$")	4✔
448	# annotate_location_2 = [-3.5, 0.83, 0.53]
449	if annotate_location_2 is not None:	4!
450	x = annotate_location_2[0]	4✔
451	y1 = annotate_location_2[1]	4✔
452	y2 = annotate_location_2[2]	4✔
453	ax2.text(x, y1, f"$N = {N}$", fontsize=18)	4✔
454	e_r = mud_prob.expected_ratio()	4✔
455	ax2.text(x, y2, rf"$\mathbb{{E}}(r) = {e_r:0.4}$", fontsize=18)	4✔
456	ax2.legend()	4✔
457	if save_path is not None:	4!
458	save_figure(	4✔
459	f"solution_n{N}",
460	save_path,
461	close_fig=close_fig,
462	dpi=dpi,
463	bbox_inches="tight",
464	)
465
466	axes.append([ax1, ax2])	4✔
467	probs.append(mud_prob)	4✔
468
469	return (poisson_prob, probs, axes)	4✔

mathematicalmichael / mud / 7607107083

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous