Add tutorial on univariate interpolation

Author: jdjakem

docs/source/conf.py

@@ -70,6 +70,7 @@
     # ExpDesign
     # 'plot_bayesian_oed.py', # ignore until paper published
     # Surrogates
+    'plot_univariate_interpolation.py',
     'plot_tensor_product_interpolation.py',
     'plot_sparse_grids.py',
     'plot_adaptive_leja_interpolation.py',
@@ -164,7 +165,7 @@ def __call__(self, filename):
 #     add_markdown_cell(work_notebook, first_cell)#jdj
 #     add_code_cell(work_notebook,"%matplotlib inline")#jdj
 # then add user defs like so
-sphinx_gallery_conf['first_notebook_cell'] = r"Add latex macros$$\newcommand{\V}[1]{{\boldsymbol{#1}}}\newcommand{mean}[1]{{\mathbb{E}\left[#1\right]}}\newcommand{var}[1]{{\mathbb{V}\left[#1\right]}}\newcommand{covar}[2]{\mathbb{C}\text{ov}\left[#1,#2\right]}\newcommand{corr}[2]{\mathbb{C}\text{or}\left[#1,#2\right]}\newcommand{argmin}{\mathrm{argmin}}\def\rv{z}\def\reals{\mathbb{R}}\def\rvset{{\mathcal{Z}}}\def\pdf{\rho}\def\rvdom{\Gamma}\def\coloneqq{\colon=}\newcommand{norm}{\lVert #1 \rVert}\def\argmax{\operatorname{argmax}}\def\ai{\alpha}\def\bi{\beta}\newcommand{\dx}[1]{\;\mathrm{d}#1}\newcommand{mat}[1]{{\boldsymbol{\mathrm{#1}}}}$$"
+sphinx_gallery_conf['first_notebook_cell'] = r"Add latex macros$$\newcommand{\V}[1]{{\boldsymbol{#1}}}\newcommand{mean}[1]{{\mathbb{E}\left[#1\right]}}\newcommand{var}[1]{{\mathbb{V}\left[#1\right]}}\newcommand{covar}[2]{\mathbb{C}\text{ov}\left[#1,#2\right]}\newcommand{corr}[2]{\mathbb{C}\text{or}\left[#1,#2\right]}\newcommand{argmin}{\mathrm{argmin}}\def\rv{z}\def\reals{\mathbb{R}}\def\rvset{{\mathcal{Z}}}\def\pdf{\rho}\def\rvdom{\Gamma}\def\coloneqq{\colon=}\newcommand{norm}{\lVert #1 \rVert}\def\argmax{\operatorname{argmax}}\def\ai{\alpha}\def\bi{\beta}\newcommand{\dx}[1]{\;\text{d}#1}\newcommand{\mat}[1]{{\boldsymbol{\mathrm{#1}}}}$$"
 
 # if change conf make sure to remove source/auto_examples, using make clean
 # Note sphinx can use align with single line, e.g. a=1 & & b=1 if \\ is added to the end of the line, i.e  a=1 & & b=1\\
@@ -265,11 +266,10 @@ def __call__(self, filename):
 \def\rvdom{\Gamma}
 \def\coloneqq{\colon=}
 \newcommand{\norm}[1]{\lVert #1 \rVert}
-%\def\argmax{\operatorname{argmax}}
 \def\ai{\alpha}
 \def\bi{\beta}
-\newcommand{\dx}[1]{\;\mathrm{d}#1}
-\newcommand{mat}[1]{{\boldsymbol{\mathrm{#1}}}}
+\newcommand{\dx}[1]{\;\text{d}#1}
+\newcommand{\mat}[1]{{\boldsymbol{\mathrm{#1}}}}
 ''',
 }
 
@@ -327,7 +327,7 @@ def __call__(self, filename):
           "corr": [r'\mathbb{C}\text{or}\left[#1,#2\right]', 2],
           "ai": r'\alpha',
           "bi": r'\beta',
-          "dx": [r'\;\mathrm{d}#1', 1],
+          "dx": [r'\;\text{d}#1', 1],
           "mat": [r'{\boldsymbol{\mathrm{#1}}}', 1],
     },
   },

pyapprox/surrogates/interp/barycentric_interpolation.py

@@ -513,54 +513,3 @@ def tensor_product_barycentric_lagrange_interpolation(
     if not return_all:
         return interp_vals
     return interp_vals, grid_samples
-
-
-def plot_tensor_product_lagrange_basis_2d(
-        univariate_quad_rule, npts, ii, jj, ax=None):
-    abscissa = univariate_quad_rule(npts)[0]
-    abscissa_1d = [abscissa, abscissa]
-    barycentric_weights_1d = [compute_barycentric_weights_1d(abscissa_1d[0]),
-                              compute_barycentric_weights_1d(abscissa_1d[1])]
-    training_samples = cartesian_product(abscissa_1d, 1)
-    fn_vals = np.zeros((training_samples.shape[1], 1))
-    idx = jj*abscissa_1d[1].shape[0]+ii
-    fn_vals[idx] = 1.
-
-    def f(samples): return multivariate_barycentric_lagrange_interpolation(
-        samples, abscissa_1d, barycentric_weights_1d, fn_vals,
-        np.array([0, 1]))
-
-    plot_limits = [-1, 1, -1, 1]
-    num_pts_1d = 101
-    X, Y, Z = get_meshgrid_function_data(f, plot_limits, num_pts_1d)
-    if ax is None:
-        ax = create_3d_axis()
-    cmap = mpl.cm.coolwarm
-
-    plot_surface(X, Y, Z, ax, axis_labels=None, limit_state=None,
-                 alpha=0.3, cmap=mpl.cm.coolwarm, zorder=3, plot_axes=False)
-    num_contour_levels = 30
-    offset = -(Z.max()-Z.min())/2
-    cmap = mpl.cm.gray
-    ax.contourf(
-        X, Y, Z, zdir='z', offset=offset,
-        levels=np.linspace(Z.min(), Z.max(), num_contour_levels),
-        cmap=cmap, zorder=-1)
-    ax.plot(training_samples[0, :], training_samples[1, :],
-            offset*np.ones(training_samples.shape[1]), 'o',
-            zorder=100, color='b')
-
-    x = np.linspace(-1, 1, 100)
-    y = training_samples[1, idx]*np.ones((x.shape[0]))
-    z = f(np.vstack((x[np.newaxis, :], y[np.newaxis, :])))[:, 0]
-    ax.plot(x, Y.max()*np.ones((x.shape[0])), z, '-r')
-    ax.plot(abscissa_1d[0], Y.max()*np.ones(
-        (abscissa_1d[0].shape[0])), np.zeros(abscissa_1d[0].shape[0]), 'or')
-
-    y = np.linspace(-1, 1, 100)
-    x = training_samples[0, idx]*np.ones((y.shape[0]))
-    z = f(np.vstack((x[np.newaxis, :], y[np.newaxis, :])))[:, 0]
-    ax.plot(X.min()*np.ones((x.shape[0])), y, z, '-r')
-    ax.plot(X.min()*np.ones(
-        (abscissa_1d[1].shape[0])), abscissa_1d[1],
-        np.zeros(abscissa_1d[1].shape[0]), 'or')

pyapprox/surrogates/interp/tensorprod.py

@@ -504,15 +504,11 @@ def get_univariate_interpolation_basis(basis_type):
     return basis_dict[basis_type]()
 
 
-class TensorProductInterpolant():
+class TensorProductBasis():
     def __init__(self, bases_1d):
         self._bases_1d = bases_1d
 
-        self._nodes_1d = None
-        self._nnodes_1d = None
-        self._values = None
-
-    def basis(self, nodes_1d, samples):
+    def __call__(self, nodes_1d, samples):
         # assumes each array in nodes_1d is in ascending order
         nvars = len(nodes_1d)
         nnodes_1d = np.array([n.shape[0] for n in nodes_1d])
@@ -535,11 +531,65 @@ def basis(self, nodes_1d, samples):
         indices = cartesian_product([np.arange(N) for N in nnodes_1d_active])
         nindices = indices.shape[1]
         temp2 = temp1[indices.ravel()+np.repeat(
-            np.arange(nactive_vars)*basis_vals_1d.shape[1], nindices), :].reshape(
-                nactive_vars, nindices, nsamples)
+            np.arange(nactive_vars)*basis_vals_1d.shape[1],
+            nindices), :].reshape(nactive_vars, nindices, nsamples)
         basis_vals = np.prod(temp2, axis=0).T
         return basis_vals
 
+    def _single_basis_fun(self, nodes_1d, idx, xx):
+        return self(nodes_1d, xx)[:, idx]
+
+    def plot_single_basis(self, ax, nodes_1d, ii, jj, nodes=None,
+                          plot_limits=[-1, 1, -1, 1],
+                          num_pts_1d=101, surface_cmap="coolwarm",
+                          contour_cmap="gray"):
+        from pyapprox.util.visualization import (
+            get_meshgrid_function_data, plot_surface)
+        idx = jj*nodes_1d[0].shape[0]+ii
+        single_basis_fun = partial(self._single_basis_fun, nodes_1d, idx)
+        X, Y, Z = get_meshgrid_function_data(
+            single_basis_fun, plot_limits, num_pts_1d)
+        if surface_cmap is not None:
+            plot_surface(X, Y, Z, ax, axis_labels=None, limit_state=None,
+                         alpha=0.3, cmap="coolwarm", zorder=3, plot_axes=False)
+        if contour_cmap is not None:
+            num_contour_levels = 30
+            offset = -(Z.max()-Z.min())/2
+            ax.contourf(
+                X, Y, Z, zdir='z', offset=offset,
+                levels=np.linspace(Z.min(), Z.max(), num_contour_levels),
+                cmap="gray", zorder=-1)
+
+        if nodes is None:
+            return
+        ax.plot(nodes[0, :], nodes[1, :],
+                offset*np.ones(nodes.shape[1]), 'o',
+                zorder=100, color='b')
+
+        x = np.linspace(-1, 1, 100)
+        y = nodes[1, idx]*np.ones((x.shape[0]))
+        z = single_basis_fun(np.vstack((x[None, :], y[None, :])))
+        ax.plot(x, Y.max()*np.ones((x.shape[0])), z, '-r')
+        ax.plot(nodes_1d[0], Y.max()*np.ones(
+            (nodes_1d[0].shape[0])), np.zeros(nodes_1d[0].shape[0]), 'or')
+
+        y = np.linspace(-1, 1, 100)
+        x = nodes[0, idx]*np.ones((y.shape[0]))
+        z = single_basis_fun(np.vstack((x[None, :], y[None, :])))
+        ax.plot(X.min()*np.ones((x.shape[0])), y, z, '-r')
+        ax.plot(X.min()*np.ones(
+            (nodes_1d[1].shape[0])), nodes_1d[1],
+                np.zeros(nodes_1d[1].shape[0]), 'or')
+
+
+class TensorProductInterpolant():
+    def __init__(self, bases_1d):
+        self.basis = TensorProductBasis(bases_1d)
+
+        self._nodes_1d = None
+        self._nnodes_1d = None
+        self._values = None
+
     def tensor_product_grid(self, nodes_1d):
         return cartesian_product(nodes_1d)
 

tutorials/multi_fidelity/plot_ensemble_selection.py

@@ -13,7 +13,7 @@
 from pyapprox.interface.wrappers import WorkTrackingModel, TimerModel
 from pyapprox.util.visualization import mathrm_label
 
-np.random.seed()
+np.random.seed(1)
 
 #%%
 #Configure the benchmark
@@ -41,8 +41,6 @@
 # Add wraper to compute the time each model takes to run
 funs = [WorkTrackingModel(
     TimerModel(fun), base_model=fun) for fun in benchmark.funs]
-for fun in funs:
-    print(fun)
 
 #%%
 #Run the pilot study
@@ -66,17 +64,22 @@
 ax = plt.subplots(1, 1, figsize=(8, 6))[1]
 multifidelity.plot_model_costs(model_costs, ax=ax)
 
+ax = plt.subplots(1, 1, figsize=(16, 12))[1]
+_ = multifidelity.plot_correlation_matrix(
+    multifidelity.get_correlation_from_covariance(stat._cov.numpy()), ax=ax,
+    format_string=None)
+
 #%%
-#Find the best estimator
-#-----------------------
-# The following code finds the best estimator by iterating over all possible subsets of models. Note this code can easily be modified to also iterate over a list of estimator types.
+#Find the best subset of of low-fidelity models
+#----------------------------------------------
+# The following code finds the subset of models that minimizes the variance of a single estimator by iterating over all possible subsets of models and computing the optimal sample allocation and assocated estimator variance. Here we restrict our attention to model subsets that contain at most 4 models.
 
 # Some MFMC estimators will fail because the models
 # do not satisfy its hierarchical condition so set
 # allow_failures=True
 best_est = multifidelity.get_estimator(
     "mfmc", stat, model_costs, allow_failures=True,
-    max_nmodels=4, save_candidate_estimators=True)
+    max_nmodels=5, save_candidate_estimators=True)
 
 target_cost = 1e2
 best_est.allocate_samples(target_cost)
@@ -110,3 +113,48 @@
     best_ests, est_labels, ax)
 ax.set_xlabel(mathrm_label("Low fidelity models"))
 plt.show()
+
+#%%
+#Find the best estimator
+#-----------------------
+#We can also find the best estimator from a list of estimator types while still determining the best model subset. This code chooses the best estimator from two possible parameterized ACV estimator classes. Specifically it chooses from all possible generalized recursive difference (GRD) estimators and genearlized multifidelity estimators that use at most 3 models and have a maximum tree depth of 3.
+best_est = multifidelity.get_estimator(
+    ["grd", "gmf"], stat, model_costs, allow_failures=True,
+    max_nmodels=3, tree_depth=3,
+    save_candidate_estimators=True)
+
+target_cost = 1e2
+best_est.allocate_samples(target_cost)
+print("Predicted variance",
+      best_est._covariance_from_npartition_samples(
+          best_est._rounded_npartition_samples))
+
+# Sort candidate estimators into lists with the same estimator type
+from collections import defaultdict
+est_dict = defaultdict(list)
+for result in best_est._candidate_estimators:
+    if result[0] is None:
+        # skip failures
+        continue
+    est_dict[result[0].__class__.__name__].append(result)
+
+
+est_name_list = list(est_dict.keys())
+est_name_list.sort()
+best_est_indices = [
+    np.argmin([result[0]._optimized_criteria for result in est_dict[name]])
+    for name in est_name_list]
+best_ests = [est_dict[est_name_list[ii]][best_est_indices[ii]][0]
+             for ii in range(len(est_name_list))]
+est_labels = [
+    "{0}({1}, {2})".format(
+        est_name_list[ii],
+        est_dict[est_name_list[ii]][best_est_indices[ii]][1],
+        est_dict[est_name_list[ii]][best_est_indices[ii]][0]._recursion_index)
+    for ii in range(len(est_name_list))]
+
+ax = plt.subplots(1, 1, figsize=(8, 6))[1]
+_ = multifidelity.plot_estimator_variance_reductions(
+    best_ests, est_labels, ax)
+ax.set_xlabel(mathrm_label("Estimator types"))
+plt.show()

tutorials/multi_fidelity/plot_multioutput_acv.py

@@ -17,8 +17,14 @@
 nmodels = 3
 
 cov = benchmark.covariance
-W = benchmark.fun.covariance_of_centered_values_kronker_product()
-B = benchmark.fun.covariance_of_mean_and_variance_estimators()
+
+labels = ([r"$(f_{0})_{%d}$" % (ii+1) for ii in range(benchmark.nqoi)] +
+          [r"$(f_{2})_{%d}$" % (ii+1) for ii in range(benchmark.nqoi)] +
+          [r"$(f_{2})_{%d}$" % (ii+1) for ii in range(benchmark.nqoi)])
+ax = plt.subplots(1, 1, figsize=(8, 6))[1]
+_ = mf.plot_correlation_matrix(
+    mf.get_correlation_from_covariance(cov), ax=ax, model_names=labels,
+    label_fontsize=20)
 
 target_cost = 10
 stat = mf.multioutput_stats["mean"](benchmark.nqoi)

tutorials/multi_fidelity/plot_multioutput_monte_carlo.py

@@ -51,7 +51,7 @@
 #Now construct the estimator. The following requires the costs of the model
 #the covariance and some other quantities W and B. These four quantities are not
 #needed to compute the value of the estimator from a sample set, however they are needed to compute the MSE of the estimator and the number of samples of that model for a given target cost. We load these quantities but ignore there meaning for the moment.
-costs = [1]
+costs = np.array([1])
 nqoi = 3
 cov = benchmark.covariance[:3, :3]
 W = benchmark.fun.covariance_of_centered_values_kronker_product()[:9, :9]

tutorials/multi_fidelity/plot_pilot_studies.py

@@ -183,7 +183,7 @@ def compute_mse(build_acv, funs, variable, target_cost, npilot_samples,
 oracle_mse = oracle_est._optimized_covariance[0, 0].item()
 ax = plt.subplots(1, 1, figsize=(8, 6))[1]
 ax.axhline(y=oracle_mse, ls='--', color='k', label=mathrm_label("Oracle MSE"))
-ax.axhline(y=mc_mse, ls=':', color='r', label=mathrm_label("MC MSE"))
+ax.axhline(y=mc_mse, ls=':', color='r', label=mathrm_label("Oracle MC MSE"))
 ax.plot(npilot_samples_list, mse_list, '-o', label=mathrm_label("Pilot MSE"))
 ax.set_xlabel(mathrm_label("Number of pilot samples"))
 _ = ax.legend()
@@ -212,7 +212,7 @@ def compute_mse(build_acv, funs, variable, target_cost, npilot_samples,
 
 ax = plt.subplots(1, 1, figsize=(8, 6))[1]
 ax.axhline(y=oracle_mse, ls='--', color='k', label=mathrm_label("Oracle MSE"))
-ax.axhline(y=mc_mse, ls=':', color='r', label=mathrm_label("MC MSE"))
+ax.axhline(y=mc_mse, ls=':', color='r', label=mathrm_label("Oracle MC MSE"))
 ax.plot(npilot_samples_list, mse_list, '-o', label=mathrm_label("Pilot MSE"))
 ax.set_xlabel(mathrm_label("Number of pilot samples"))
 _ = ax.legend()