From c06d12ba2c85a086003b55577a75a3570a67d1f2 Mon Sep 17 00:00:00 2001 From: Marcus Date: Wed, 27 Nov 2024 10:48:41 -0800 Subject: [PATCH] x_out optional even for multi-task and bug fix with entropy related posterior quantities --- docs/source/examples/MultiTaskTest.ipynb | 34 +++-- .../examples/NonEuclideanInputSpaces.ipynb | 5 +- docs/source/examples/SingleTaskTest.ipynb | 8 +- docs/source/examples/gp2ScaleTest.ipynb | 3 +- examples/MultiTaskTest.ipynb | 21 ++- examples/NonEuclideanInputSpaces.ipynb | 5 +- examples/SingleTaskTest.ipynb | 4 +- examples/gp2ScaleTest.ipynb | 3 +- fvgp/fvgp.py | 2 - fvgp/gp.py | 31 +--- fvgp/gp_posterior.py | 136 +++++++----------- tests/test_fvgp.py | 40 +++++- 12 files changed, 146 insertions(+), 146 deletions(-) diff --git a/docs/source/examples/MultiTaskTest.ipynb b/docs/source/examples/MultiTaskTest.ipynb index f919224..a9e0358 100644 --- a/docs/source/examples/MultiTaskTest.ipynb +++ b/docs/source/examples/MultiTaskTest.ipynb @@ -16,8 +16,9 @@ "metadata": {}, "outputs": [], "source": [ - "##first install the newest version of fvgp\n", - "#!pip install fvgp~=4.5.3" + "##first, install the newest version of fvgp\n", + "#!pip install fvgp~=4.5.6\n", + "#!pip install torch" ] }, { @@ -133,9 +134,10 @@ "outputs": [], "source": [ "#mean and standard deviation\n", - "mean = my_gp2.posterior_mean(x_pred=x_pred1d.reshape(50,1), x_out=np.array([0,1]))[\"f(x)\"]\n", + "mean = my_gp2.posterior_mean(x_pred=x_pred1d.reshape(50,1))[\"f(x)\"]\n", "std = np.sqrt(my_gp2.posterior_covariance(x_pred=x_pred1d.reshape(50,1), x_out=np.array([0,1]))[\"v(x)\"])\n", "\n", + "\n", "plt.plot(x_pred1d.reshape(50,1),mean[:,0], label = \"mean task 1\")\n", "plt.plot(x_pred1d.reshape(50,1),mean[:,1], label = \"mean task 2\")\n", "plt.scatter(x_data,y_data1) \n", @@ -411,9 +413,9 @@ "\n", "\n", "bounds = np.zeros((n.number_of_hps+2,2))\n", - "bounds[0] = np.array([0.001,10.])\n", - "bounds[1] = np.array([0.001,10.])\n", - "bounds[2:] = np.array([-1,1])\n", + "bounds[0] = np.array([0.1,10.])\n", + "bounds[1] = np.array([0.1,10.])\n", + "bounds[2:] = np.array([0.01,1])\n", "my_gp2.train(hyperparameter_bounds=bounds,max_iter = 2)" ] }, @@ -508,13 +510,29 @@ "fig.show()\n", "\n" ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "811a3337-a377-44bd-85b7-9a699c299cd0", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "069a28fb-fecf-42e9-bc27-16fcccec3080", + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": "fvgp_dev", "language": "python", - "name": "python3" + "name": "fvgp_dev" }, "language_info": { "codemirror_mode": { diff --git a/docs/source/examples/NonEuclideanInputSpaces.ipynb b/docs/source/examples/NonEuclideanInputSpaces.ipynb index dca3739..547d27b 100644 --- a/docs/source/examples/NonEuclideanInputSpaces.ipynb +++ b/docs/source/examples/NonEuclideanInputSpaces.ipynb @@ -19,7 +19,7 @@ "outputs": [], "source": [ "#install the newest version of fvgp\n", - "#!pip install fvgp~=4.5.3" + "#!pip install fvgp~=4.5.6" ] }, { @@ -214,7 +214,8 @@ "outputs": [], "source": [ "x_pred = [\"dwed\",\"dwe\"]\n", - "my_gp2.posterior_mean(x_pred, x_out = np.array([0,1,2,3]))" + "my_gp2.posterior_mean(x_pred, x_out = np.array([0,1,2,3]))\n", + "my_gp2.posterior_mean(x_pred)" ] }, { diff --git a/docs/source/examples/SingleTaskTest.ipynb b/docs/source/examples/SingleTaskTest.ipynb index a4cdfc1..5606110 100644 --- a/docs/source/examples/SingleTaskTest.ipynb +++ b/docs/source/examples/SingleTaskTest.ipynb @@ -13,7 +13,7 @@ "id": "a1b1c026", "metadata": {}, "source": [ - "This is the new test for fvgp version 4.3.0 and later." + "This is the new test for fvgp version 4.5.5 and later." ] }, { @@ -23,7 +23,7 @@ "metadata": {}, "outputs": [], "source": [ - "#!pip install fvgp~=4.5.3" + "#!pip install fvgp~=4.5.5" ] }, { @@ -170,7 +170,7 @@ " [0.001,0.1], #noise\n", " [0.01,1.] #mean\n", " ])\n", - "my_gp1.update_gp_data(x_data, y_data, noise_variances_new=np.ones(y_data.shape) * 0.01)\n", + "my_gp1.update_gp_data(x_data, y_data, noise_variances_new=np.ones(y_data.shape) * 0.01) #this is just for testing, not needed\n", "print(\"Standard Training\")\n", "my_gp1.train(hyperparameter_bounds=hps_bounds)\n", "print(\"Global Training\")\n", @@ -316,7 +316,6 @@ "my_gp1.gp_entropy(x_test)\n", "my_gp1.gp_entropy_grad(x_test, 0)\n", "my_gp1.gp_kl_div(x_test, np.ones((len(x_test))), np.identity((len(x_test))))\n", - "my_gp1.gp_kl_div_grad(x_test, np.ones((len(x_test))), np.identity((len(x_test))), 0)\n", "my_gp1.gp_relative_information_entropy(x_test)\n", "my_gp1.gp_relative_information_entropy_set(x_test)\n", "my_gp1.posterior_covariance(x_test)\n", @@ -324,7 +323,6 @@ "my_gp1.posterior_mean(x_test)\n", "my_gp1.posterior_mean_grad(x_test)\n", "my_gp1.posterior_probability(x_test, np.ones((len(x_test))), np.identity((len(x_test))))\n", - "my_gp1.posterior_probability_grad(x_test, np.ones((len(x_test))), np.identity((len(x_test))),0)\n", "\n" ] }, diff --git a/docs/source/examples/gp2ScaleTest.ipynb b/docs/source/examples/gp2ScaleTest.ipynb index a767a8b..45c4805 100644 --- a/docs/source/examples/gp2ScaleTest.ipynb +++ b/docs/source/examples/gp2ScaleTest.ipynb @@ -19,7 +19,8 @@ "outputs": [], "source": [ "##first install the newest version of fvgp\n", - "#!pip install fvgp~=4.5.3" + "#!pip install fvgp~=4.5.6\n", + "#!pip install imate" ] }, { diff --git a/examples/MultiTaskTest.ipynb b/examples/MultiTaskTest.ipynb index 199df1a..a9e0358 100644 --- a/examples/MultiTaskTest.ipynb +++ b/examples/MultiTaskTest.ipynb @@ -17,7 +17,8 @@ "outputs": [], "source": [ "##first, install the newest version of fvgp\n", - "#!pip install fvgp~=4.5.5" + "#!pip install fvgp~=4.5.6\n", + "#!pip install torch" ] }, { @@ -412,9 +413,9 @@ "\n", "\n", "bounds = np.zeros((n.number_of_hps+2,2))\n", - "bounds[0] = np.array([0.001,10.])\n", - "bounds[1] = np.array([0.001,10.])\n", - "bounds[2:] = np.array([-1,1])\n", + "bounds[0] = np.array([0.1,10.])\n", + "bounds[1] = np.array([0.1,10.])\n", + "bounds[2:] = np.array([0.01,1])\n", "my_gp2.train(hyperparameter_bounds=bounds,max_iter = 2)" ] }, @@ -517,13 +518,21 @@ "metadata": {}, "outputs": [], "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "069a28fb-fecf-42e9-bc27-16fcccec3080", + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": "fvgp_dev", "language": "python", - "name": "python3" + "name": "fvgp_dev" }, "language_info": { "codemirror_mode": { diff --git a/examples/NonEuclideanInputSpaces.ipynb b/examples/NonEuclideanInputSpaces.ipynb index dca3739..547d27b 100644 --- a/examples/NonEuclideanInputSpaces.ipynb +++ b/examples/NonEuclideanInputSpaces.ipynb @@ -19,7 +19,7 @@ "outputs": [], "source": [ "#install the newest version of fvgp\n", - "#!pip install fvgp~=4.5.3" + "#!pip install fvgp~=4.5.6" ] }, { @@ -214,7 +214,8 @@ "outputs": [], "source": [ "x_pred = [\"dwed\",\"dwe\"]\n", - "my_gp2.posterior_mean(x_pred, x_out = np.array([0,1,2,3]))" + "my_gp2.posterior_mean(x_pred, x_out = np.array([0,1,2,3]))\n", + "my_gp2.posterior_mean(x_pred)" ] }, { diff --git a/examples/SingleTaskTest.ipynb b/examples/SingleTaskTest.ipynb index b3daf2d..5606110 100644 --- a/examples/SingleTaskTest.ipynb +++ b/examples/SingleTaskTest.ipynb @@ -170,7 +170,7 @@ " [0.001,0.1], #noise\n", " [0.01,1.] #mean\n", " ])\n", - "my_gp1.update_gp_data(x_data, y_data, noise_variances_new=np.ones(y_data.shape) * 0.01)\n", + "my_gp1.update_gp_data(x_data, y_data, noise_variances_new=np.ones(y_data.shape) * 0.01) #this is just for testing, not needed\n", "print(\"Standard Training\")\n", "my_gp1.train(hyperparameter_bounds=hps_bounds)\n", "print(\"Global Training\")\n", @@ -316,7 +316,6 @@ "my_gp1.gp_entropy(x_test)\n", "my_gp1.gp_entropy_grad(x_test, 0)\n", "my_gp1.gp_kl_div(x_test, np.ones((len(x_test))), np.identity((len(x_test))))\n", - "my_gp1.gp_kl_div_grad(x_test, np.ones((len(x_test))), np.identity((len(x_test))), 0)\n", "my_gp1.gp_relative_information_entropy(x_test)\n", "my_gp1.gp_relative_information_entropy_set(x_test)\n", "my_gp1.posterior_covariance(x_test)\n", @@ -324,7 +323,6 @@ "my_gp1.posterior_mean(x_test)\n", "my_gp1.posterior_mean_grad(x_test)\n", "my_gp1.posterior_probability(x_test, np.ones((len(x_test))), np.identity((len(x_test))))\n", - "my_gp1.posterior_probability_grad(x_test, np.ones((len(x_test))), np.identity((len(x_test))),0)\n", "\n" ] }, diff --git a/examples/gp2ScaleTest.ipynb b/examples/gp2ScaleTest.ipynb index a767a8b..45c4805 100644 --- a/examples/gp2ScaleTest.ipynb +++ b/examples/gp2ScaleTest.ipynb @@ -19,7 +19,8 @@ "outputs": [], "source": [ "##first install the newest version of fvgp\n", - "#!pip install fvgp~=4.5.3" + "#!pip install fvgp~=4.5.6\n", + "#!pip install imate" ] }, { diff --git a/fvgp/fvgp.py b/fvgp/fvgp.py index 2666c07..2f12bc0 100755 --- a/fvgp/fvgp.py +++ b/fvgp/fvgp.py @@ -266,8 +266,6 @@ class provides all the methods described for the GP (:py:class:`fvgp.GP`) class. :py:meth:`fvgp.GP.posterior_probability` - :py:meth:`fvgp.GP.posterior_probability_grad` - Validation Methods: :py:meth:`fvgp.GP.crps` diff --git a/fvgp/gp.py b/fvgp/gp.py index ad66b9b..bbd8939 100755 --- a/fvgp/gp.py +++ b/fvgp/gp.py @@ -1232,7 +1232,7 @@ def posterior_probability(self, x_pred, comp_mean, comp_cov, x_out=None): comp_mean: np.ndarray A vector of mean values, same length as x_pred. comp_cov: np.nparray - Covariance matrix, in R^{len(x_pred) times len(x_pred)} + Covariance matrix, in R^{len(x_pred) x len(x_pred)} x_out : np.ndarray, optional Output coordinates in case of multi-task GP use; a numpy array of size (N), where N is the number evaluation points in the output direction. @@ -1245,35 +1245,6 @@ def posterior_probability(self, x_pred, comp_mean, comp_cov, x_out=None): """ return self.posterior.posterior_probability(x_pred, comp_mean, comp_cov, x_out=x_out) - def posterior_probability_grad(self, x_pred, comp_mean, comp_cov, direction, x_out=None): - """ - Function to compute the gradient of the probability of a probabilistic quantity of interest, - given the GP posterior at a given point. - - Parameters - ---------- - x_pred : np.ndarray or list - A numpy array of shape (V x D), interpreted as an array of input point positions, or a list for - GPs on non-Euclidean input spaces. - comp_mean: np.ndarray - A vector of mean values, same length as x_pred. - comp_cov: np.nparray - Covariance matrix, in R^{len(x_pred) times len(x_pred)} - direction : int - The direction to compute the gradient in. - x_out : np.ndarray, optional - Output coordinates in case of multi-task GP use; a numpy array of size (N), - where N is the number evaluation points in the output direction. - Usually this is np.ndarray([0,1,2,...]). - - Return - ------ - Solution : dict - The gradient of the probability of a probabilistic quantity of interest, - given the GP posterior at a given point. - """ - return self.posterior.posterior_probability_grad(x_pred, comp_mean, comp_cov, direction, x_out=x_out) - #################################################################################### #################################################################################### #######################VALIDATION################################################### diff --git a/fvgp/gp_posterior.py b/fvgp/gp_posterior.py index 69ac463..6783f71 100755 --- a/fvgp/gp_posterior.py +++ b/fvgp/gp_posterior.py @@ -34,15 +34,18 @@ def posterior_mean(self, x_pred, hyperparameters=None, x_out=None): self._perform_input_checks(x_pred, x_out) x_orig = x_pred.copy() - if x_out is not None: x_pred = self.cartesian_product(x_pred, x_out) + if isinstance(x_out, np.ndarray): x_pred = self.cartesian_product(x_pred, x_out) k = self.kernel(x_data, x_pred, hyperparameters) A = k.T @ KVinvY posterior_mean = self.mean_function(x_pred, hyperparameters) + A - if x_out is not None: posterior_mean = posterior_mean.reshape(len(x_orig), len(x_out), order='F') + if isinstance(x_out, np.ndarray): posterior_mean_re = posterior_mean.reshape(len(x_orig), len(x_out), order='F') + else: posterior_mean_re = posterior_mean return {"x": x_orig, - "f(x)": posterior_mean} + "f(x)": posterior_mean_re, + "f(x)_flat": posterior_mean, + "x_pred": x_pred} def posterior_mean_grad(self, x_pred, hyperparameters=None, x_out=None, direction=None): x_data, y_data, KVinvY = \ @@ -59,7 +62,7 @@ def posterior_mean_grad(self, x_pred, hyperparameters=None, x_out=None, directio if x_out is None: x_out = self.x_out self._perform_input_checks(x_pred, x_out) x_orig = x_pred.copy() - if x_out is not None: x_pred = self.cartesian_product(x_pred, x_out) + if isinstance(x_out, np.ndarray): x_pred = self.cartesian_product(x_pred, x_out) f = self.mean_function(x_pred, hyperparameters) eps = 1e-6 @@ -69,7 +72,8 @@ def posterior_mean_grad(self, x_pred, hyperparameters=None, x_out=None, directio mean_der = (self.mean_function(x1, hyperparameters) - f) / eps k_g = self.d_kernel_dx(x_pred, x_data, direction, hyperparameters) posterior_mean_grad = mean_der + (k_g @ KVinvY) - if x_out is not None: posterior_mean_grad = posterior_mean_grad.reshape(len(x_orig), len(x_out), order='F') + if isinstance(x_out, np.ndarray): + posterior_mean_grad = posterior_mean_grad.reshape(len(x_orig), len(x_out), order='F') else: posterior_mean_grad = np.zeros((len(x_pred), x_orig.shape[1])) for direction in range(len(x_orig[0])): @@ -79,7 +83,7 @@ def posterior_mean_grad(self, x_pred, hyperparameters=None, x_out=None, directio k_g = self.d_kernel_dx(x_pred, x_data, direction, hyperparameters) posterior_mean_grad[:, direction] = mean_der + (k_g @ KVinvY) direction = "ALL" - if x_out is not None: + if isinstance(x_out, np.ndarray): posterior_mean_grad = posterior_mean_grad.reshape(len(x_orig), len(x_orig[0]), len(x_out), order='F') return {"x": x_orig, @@ -92,7 +96,7 @@ def posterior_covariance(self, x_pred, x_out=None, variance_only=False, add_nois if x_out is None: x_out = self.x_out self._perform_input_checks(x_pred, x_out) x_orig = x_pred.copy() - if x_out is not None: x_pred = self.cartesian_product(x_pred, x_out) + if isinstance(x_out, np.ndarray): x_pred = self.cartesian_product(x_pred, x_out) k = self.kernel(x_data, x_pred, self.prior_obj.hyperparameters) kk = self.kernel(x_pred, x_pred, self.prior_obj.hyperparameters) @@ -121,20 +125,27 @@ def posterior_covariance(self, x_pred, x_out=None, variance_only=False, add_nois if add_noise: v, S = self.add_noise(x_pred, v, S) - if x_out is not None: - v = v.reshape(len(x_orig), len(x_out), order='F') - if S is not None: S = S.reshape(len(x_orig), len(x_orig), len(x_out), len(x_out), order='F') + if isinstance(x_out, np.ndarray): + v_re = v.reshape(len(x_orig), len(x_out), order='F') + if S is not None: S_re = S.reshape(len(x_orig), len(x_orig), len(x_out), len(x_out), order='F') + else: S_re = None + else: + v_re = v + S_re = S return {"x": x_orig, - "v(x)": v, - "S": S} + "x_pred": x_pred, + "v(x)": v_re, + "S": S_re, + "S_flat": S, + "v_flat": v} def posterior_covariance_grad(self, x_pred, x_out=None, direction=None): x_data = self.data_obj.x_data.copy() if x_out is None: x_out = self.x_out self._perform_input_checks(x_pred, x_out) x_orig = x_pred.copy() - if x_out is not None: x_pred = self.cartesian_product(x_pred, x_out) + if isinstance(x_out, np.ndarray): x_pred = self.cartesian_product(x_pred, x_out) k = self.kernel(x_data, x_pred, self.prior_obj.hyperparameters) k_covariance_prod = self.marginal_density_obj.KVlinalg.solve(k) @@ -148,7 +159,7 @@ def posterior_covariance_grad(self, x_pred, x_out=None, direction=None): self.kernel(x2, x2, self.prior_obj.hyperparameters)) / eps dSdx = kk_g - (2.0 * k_g.T @ k_covariance_prod) a = np.diag(dSdx) - if x_out is not None: + if isinstance(x_out, np.ndarray): a = a.reshape(len(x_orig), len(x_out), order='F') dSdx = dSdx.reshape(len(x_orig), len(x_orig), len(x_out), len(x_out), order='F') return {"x": x_orig, @@ -166,7 +177,8 @@ def posterior_covariance_grad(self, x_pred, x_out=None, direction=None): self.kernel(x2, x2, self.prior_obj.hyperparameters)) / eps grad_v[:, direction] = np.diag(kk_g - (2.0 * k_g.T @ k_covariance_prod)) - if x_out is not None: grad_v = grad_v.reshape(len(x_orig), len(x_orig[0]), len(x_out), order='F') + if isinstance(x_out, np.ndarray): + grad_v = grad_v.reshape(len(x_orig), len(x_orig[0]), len(x_out), order='F') return {"x": x_orig, "dv/dx": grad_v} @@ -178,7 +190,7 @@ def joint_gp_prior(self, x_pred, x_out=None): self.prior_obj.m.copy()) if x_out is None: x_out = self.x_out self._perform_input_checks(x_pred, x_out) - if x_out is not None: x_pred = self.cartesian_product(x_pred, x_out) + if isinstance(x_out, np.ndarray): x_pred = self.cartesian_product(x_pred, x_out) k = self.kernel(x_data, x_pred, self.prior_obj.hyperparameters) kk = self.kernel(x_pred, x_pred, self.prior_obj.hyperparameters) @@ -200,7 +212,7 @@ def joint_gp_prior_grad(self, x_pred, direction, x_out=None): self.prior_obj.m.copy()) if x_out is None: x_out = self.x_out self._perform_input_checks(x_pred, x_out) - if x_out is not None: x_pred = self.cartesian_product(x_pred, x_out) + if isinstance(x_out, np.ndarray): x_pred = self.cartesian_product(x_pred, x_out) k_g = self.d_kernel_dx(x_pred, x_data, direction, self.prior_obj.hyperparameters).T x1 = np.array(x_pred) @@ -248,11 +260,8 @@ def gp_entropy(self, x_pred, x_out=None): ------ Entropy : float """ - if x_out is None: x_out = self.x_out - self._perform_input_checks(x_pred, x_out) - if x_out is not None: x_pred = self.cartesian_product(x_pred, x_out) - priors = self.joint_gp_prior(x_pred, x_out=None) + priors = self.joint_gp_prior(x_pred, x_out=x_out) S = priors["S"] dim = len(S[0]) ldet = calculate_logdet(S) @@ -260,17 +269,14 @@ def gp_entropy(self, x_pred, x_out=None): ########################################################################### def gp_entropy_grad(self, x_pred, direction, x_out=None): - if x_out is None: x_out = self.x_out - self._perform_input_checks(x_pred, x_out) - if x_out is not None: x_pred = self.cartesian_product(x_pred, x_out) - - priors1 = self.joint_gp_prior(x_pred, x_out=None) - priors2 = self.joint_gp_prior_grad(x_pred, direction, x_out=None) + priors1 = self.joint_gp_prior(x_pred, x_out=x_out) + priors2 = self.joint_gp_prior_grad(x_pred, direction, x_out=x_out) S1 = priors1["S"] S2 = priors2["dS/dx"] return 0.5 * np.trace(calculate_inv(S1) @ S2) ########################################################################### + @staticmethod def kl_div_grad(self, mu1, dmu1dx, mu2, S1, dS1dx, S2): x1 = solve(S2, dS1dx) mu = np.subtract(mu2, mu1) @@ -299,12 +305,11 @@ def kl_div(self, mu1, mu2, S1, S2): ########################################################################### def gp_kl_div(self, x_pred, comp_mean, comp_cov, x_out=None): if x_out is None: x_out = self.x_out - self._perform_input_checks(x_pred, x_out) - if x_out is not None: x_pred = self.cartesian_product(x_pred, x_out) - res = self.posterior_mean(x_pred, x_out=None) - gp_mean = res["f(x)"] - gp_cov = self.posterior_covariance(x_pred, x_out=None)["S"] + np.identity(len(x_pred)) * 1e-9 + res = self.posterior_mean(x_pred, x_out=x_out) + gp_mean = res["f(x)_flat"] + gp_cov = self.posterior_covariance(x_pred, x_out=x_out)["S_flat"] + gp_cov = gp_cov + + np.identity(len(gp_cov)) * 1e-9 comp_cov = comp_cov + np.identity(len(comp_cov)) * 1e-9 return {"x": x_pred, "gp posterior mean": gp_mean, @@ -313,26 +318,6 @@ def gp_kl_div(self, x_pred, comp_mean, comp_cov, x_out=None): "given covariance": comp_cov, "kl-div": self.kl_div(gp_mean, comp_mean, gp_cov, comp_cov)} - ########################################################################### - def gp_kl_div_grad(self, x_pred, comp_mean, comp_cov, direction, x_out=None): - if x_out is None: x_out = self.x_out - self._perform_input_checks(x_pred, x_out) - if x_out is not None: x_pred = self.cartesian_product(x_pred, x_out) - - gp_mean = self.posterior_mean(x_pred, x_out=None)["f(x)"] - gp_mean_grad = self.posterior_mean_grad(x_pred, direction=direction, x_out=None)["df/dx"] - gp_cov = self.posterior_covariance(x_pred, x_out=None)["S"] + np.identity(len(x_pred)) * 1e-9 - gp_cov_grad = self.posterior_covariance_grad(x_pred, direction=direction, x_out=None)["dS/dx"] - comp_cov = comp_cov + np.identity(len(comp_cov)) * 1e-9 - return {"x": x_pred, - "gp posterior mean": gp_mean, - "gp posterior mean grad": gp_mean_grad, - "gp posterior covariance": gp_cov, - "gp posterior covariance grad": gp_cov_grad, - "given mean": comp_mean, - "given covariance": comp_cov, - "kl-div grad": self.kl_div_grad(gp_mean, gp_mean_grad, comp_mean, gp_cov, gp_cov_grad, comp_cov)} - ########################################################################### def mutual_information(self, joint, m1, m2): return self.entropy(m1) + self.entropy(m2) - self.entropy(joint) @@ -343,7 +328,7 @@ def gp_mutual_information(self, x_pred, x_out=None, add_noise=False): if x_out is None: x_out = self.x_out self._perform_input_checks(x_pred, x_out) x_orig = x_pred.copy() - if x_out is not None: x_pred = self.cartesian_product(x_pred, x_out) + if isinstance(x_out, np.ndarray): x_pred = self.cartesian_product(x_pred, x_out) k = self.kernel(x_data, x_pred, self.prior_obj.hyperparameters) kk = self.kernel(x_pred, x_pred, self.prior_obj.hyperparameters) + (np.identity(len(x_pred)) * 1e-9) @@ -359,7 +344,7 @@ def gp_total_correlation(self, x_pred, x_out=None, add_noise=False): if x_out is None: x_out = self.x_out self._perform_input_checks(x_pred, x_out) x_orig = x_pred.copy() - if x_out is not None: x_pred = self.cartesian_product(x_pred, x_out) + if isinstance(x_out, np.ndarray): x_pred = self.cartesian_product(x_pred, x_out) k = self.kernel(x_data, x_pred, self.prior_obj.hyperparameters) kk = self.kernel(x_pred, x_pred, self.prior_obj.hyperparameters) + (np.identity(len(x_pred)) * 1e-9) @@ -377,27 +362,27 @@ def gp_relative_information_entropy(self, x_pred, x_out=None, add_noise=False): if x_out is None: x_out = self.x_out self._perform_input_checks(x_pred, x_out) x_orig = x_pred.copy() - if x_out is not None: x_pred = self.cartesian_product(x_pred, x_out) - kk = self.kernel(x_pred, x_pred, self.prior_obj.hyperparameters) + (np.identity(len(x_pred)) * 1e-9) - post_cov = self.posterior_covariance(x_pred, x_out=None)["S"] + (np.identity(len(x_pred)) * 1e-9) - if add_noise: v, post_cov = self.add_noise(x_pred, np.diag(post_cov), post_cov) - + if isinstance(x_out, np.ndarray): x_pred_aux = self.cartesian_product(x_pred, x_out) + else: x_pred_aux = x_pred + kk = self.kernel(x_pred_aux, x_pred_aux, self.prior_obj.hyperparameters) + (np.identity(len(x_pred_aux)) * 1e-9) + post_cov = self.posterior_covariance(x_pred, x_out=x_out, add_noise=add_noise)["S_flat"] + post_cov = post_cov + (np.identity(len(post_cov)) * 1e-9) hyperparameters = self.prior_obj.hyperparameters - post_mean = self.posterior_mean(x_pred, x_out=None)["f(x)"] - prio_mean = self.mean_function(x_pred, hyperparameters) + post_mean = self.posterior_mean(x_pred, x_out=x_out)["f(x)_flat"] + prio_mean = self.mean_function(x_pred_aux, hyperparameters) return {"x": x_orig, "RIE": self.kl_div(prio_mean, post_mean, kk, post_cov)} ########################################################################### def gp_relative_information_entropy_set(self, x_pred, x_out=None, add_noise=False): if x_out is None: x_out = self.x_out - self._perform_input_checks(x_pred, x_out) + #self._perform_input_checks(x_pred, x_out) x_orig = x_pred.copy() - if x_out is not None: x_pred = self.cartesian_product(x_pred, x_out) + #if isinstance(x_out, np.ndarray): x_pred = self.cartesian_product(x_pred, x_out) RIE = np.zeros((len(x_pred))) for i in range(len(x_pred)): RIE[i] = self.gp_relative_information_entropy( - x_pred[i].reshape(1, len(x_pred[i])), x_out=None, add_noise=add_noise)["RIE"] + x_pred[i].reshape(1, len(x_pred[i])), x_out=x_out, add_noise=add_noise)["RIE"] return {"x": x_orig, "RIE": RIE} @@ -406,11 +391,10 @@ def gp_relative_information_entropy_set(self, x_pred, x_out=None, add_noise=Fals def posterior_probability(self, x_pred, comp_mean, comp_cov, x_out=None): if x_out is None: x_out = self.x_out self._perform_input_checks(x_pred, x_out) - if x_out is not None: x_pred = self.cartesian_product(x_pred, x_out) + #if isinstance(x_out, np.ndarray): x_pred = self.cartesian_product(x_pred, x_out) - res = self.posterior_mean(x_pred, x_out=None) - gp_mean = res["f(x)"] - gp_cov = self.posterior_covariance(x_pred, x_out=None)["S"] + (np.identity(len(x_pred)) * 1e-9) + gp_mean = self.posterior_mean(x_pred, x_out=x_out)["f(x)_flat"] + gp_cov = self.posterior_covariance(x_pred, x_out=x_out, add_noise=True)["S_flat"] gp_cov_inv = calculate_inv(gp_cov) comp_cov_inv = calculate_inv(comp_cov) cov = calculate_inv(gp_cov_inv + comp_cov_inv) @@ -429,20 +413,6 @@ def posterior_probability(self, x_pred, comp_mean, comp_cov, x_out=None): np.exp(ln_p) } - def posterior_probability_grad(self, x_pred, comp_mean, comp_cov, direction, x_out=None): - if x_out is None: x_out = self.x_out - self._perform_input_checks(x_pred, x_out) - if x_out is not None: x_pred = self.cartesian_product(x_pred, x_out) - - x1 = np.array(x_pred) - x2 = np.array(x_pred) - x1[:, direction] = x1[:, direction] + 1e-6 - x2[:, direction] = x2[:, direction] - 1e-6 - - probability_grad = (self.posterior_probability(x1, comp_mean, comp_cov, x_out=None)["probability"] - - self.posterior_probability(x2, comp_mean, comp_cov, x_out=None)["probability"]) / 2e-6 - return {"probability grad": probability_grad} - def add_noise(self, x_pred, v, S): if callable(self.likelihood_obj.noise_function): noise = self.likelihood_obj.noise_function(x_pred, self.prior_obj.hyperparameters) @@ -462,7 +432,7 @@ def add_noise(self, x_pred, v, S): ########################################################################### @staticmethod - def _int_gauss(self, S): + def _int_gauss(S): return ((2.0 * np.pi) ** (len(S) / 2.0)) * np.sqrt(np.linalg.det(S)) def _perform_input_checks(self, x_pred, x_out): diff --git a/tests/test_fvgp.py b/tests/test_fvgp.py index fda5ba0..778f749 100755 --- a/tests/test_fvgp.py +++ b/tests/test_fvgp.py @@ -167,9 +167,7 @@ def noiseC(x,hps): B = A.T @ A res = my_gp1.entropy(B) res = my_gp1.gp_kl_div(np.random.rand(10,len(x_data[0])), np.random.rand(10), B) - res = my_gp1.gp_kl_div_grad(np.random.rand(10,len(x_data[0])), np.random.rand(10), B,0) res = my_gp1.posterior_probability(np.random.rand(10,len(x_data[0])), np.random.rand(10), B) - res = my_gp1.posterior_probability_grad(np.random.rand(10,len(x_data[0])), np.random.rand(10), B, direction = 0) res = squared_exponential_kernel(1.,1.) res = squared_exponential_kernel_robust(1.,1.) @@ -244,7 +242,43 @@ def mkernel(x1,x2,hps): my_fvgp.posterior_mean_grad(np.random.rand(10,5))["df/dx"] my_fvgp.posterior_covariance(np.random.rand(10,5), x_out = np.array([0,1]))["v(x)"] my_fvgp.posterior_covariance(np.random.rand(10,5))["v(x)"] - + my_fvgp.posterior_covariance_grad(np.random.rand(10,5)) + my_fvgp.posterior_covariance_grad(np.random.rand(10,5), x_out = np.array([0,1])) + + my_fvgp.joint_gp_prior(np.random.rand(10,5)) + my_fvgp.joint_gp_prior(np.random.rand(10,5), x_out = np.array([0,1])) + + my_fvgp.joint_gp_prior_grad(np.random.rand(10,5), 0) + my_fvgp.joint_gp_prior_grad(np.random.rand(10,5), 0, x_out = np.array([0,1])) + + my_fvgp.gp_entropy(np.random.rand(10,5)) + my_fvgp.gp_entropy_grad(np.random.rand(10,5), 0) + my_fvgp.gp_entropy(np.random.rand(10,5), x_out = np.array([0,1])) + my_fvgp.gp_entropy_grad(np.random.rand(10,5),0, x_out = np.array([0,1])) + + A = np.random.rand(20,20) + B = A.T @ A + + + my_fvgp.gp_kl_div(np.random.rand(10,5), np.random.rand(20), B) + my_fvgp.gp_kl_div(np.random.rand(10,5), np.random.rand(20), B ,x_out = np.array([0,1])) + + my_fvgp.gp_mutual_information(np.random.rand(10,5)) + my_fvgp.gp_mutual_information(np.random.rand(10,5), x_out = np.array([0,1])) + + + my_fvgp.gp_total_correlation(np.random.rand(10,5)) + my_fvgp.gp_total_correlation(np.random.rand(10,5), x_out = np.array([0,1])) + + + my_fvgp.gp_relative_information_entropy(np.random.rand(10,5)) + my_fvgp.gp_relative_information_entropy(np.random.rand(10,5), x_out = np.array([0,1])) + + my_fvgp.gp_relative_information_entropy_set(np.random.rand(10,5)) + my_fvgp.gp_relative_information_entropy_set(np.random.rand(10,5), x_out = np.array([0,1])) + + my_fvgp.posterior_probability(np.random.rand(10,5), np.random.rand(20), B) + my_fvgp.posterior_probability(np.random.rand(10,5), np.random.rand(20), B, x_out = np.array([0,1])) my_fvgp = fvGP(np.random.rand(3,5), np.random.rand(3,2), noise_variances = None, init_hyperparameters = np.array([1, 1]), gp_kernel_function=mkernel)