Skip to content

Commit

Permalink
Remove HSGP as it became a part of PyMC (#136)
Browse files Browse the repository at this point in the history
* remove HSGP as it became a part of PyMC

* isort
  • Loading branch information
ferrine authored Apr 7, 2023
1 parent f554df1 commit b730449
Show file tree
Hide file tree
Showing 2 changed files with 1 addition and 107 deletions.
6 changes: 1 addition & 5 deletions pymc_experimental/gp/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -13,8 +13,4 @@
# limitations under the License.


from pymc_experimental.gp.latent_approx import (
HSGP,
KarhunenLoeveExpansion,
ProjectedProcess,
)
from pymc_experimental.gp.latent_approx import KarhunenLoeveExpansion, ProjectedProcess
102 changes: 0 additions & 102 deletions pymc_experimental/gp/latent_approx.py
Original file line number Diff line number Diff line change
Expand Up @@ -75,108 +75,6 @@ def conditional(self, name, Xnew, jitter=1e-6, **kwargs):
return pm.MvNormal(name, mu=mu, chol=chol)


class HSGP(pm.gp.Latent):
## inputs: M, c

def __init__(
self, n_basis, c=3 / 2, *, mean_func=pm.gp.mean.Zero(), cov_func=pm.gp.cov.Constant(0.0)
):
## TODO: specify either c or L
self.M = n_basis
self.c = c
super().__init__(mean_func=mean_func, cov_func=cov_func)

def _validate_cov_func(self, cov_func):
## TODO: actually validate it. Right now this fails unless cov func is exactly
# in the form eta**2 * pm.gp.cov.Matern12(...) and will error otherwise.
cov, scaling_factor = cov_func.factor_list
return scaling_factor, cov.ls, cov.spectral_density

def prior(self, name, X, **kwargs):
f, Phi, L, spd, beta, Xmu, Xsd = self._build_prior(name, X, **kwargs)
self.X, self.f = X, f
self.Phi, self.L, self.spd, self.beta = Phi, L, spd, beta
self.Xmu, self.Xsd = Xmu, Xsd
return f

def _generate_basis(self, X, L):
indices = pt.arange(1, self.M + 1)
m1 = (np.pi / (2.0 * L)) * pt.tile(L + X, self.M)
m2 = pt.diag(indices)
Phi = pt.sin(m1 @ m2) / pt.sqrt(L)
omega = (np.pi * indices) / (2.0 * L)
return Phi, omega

def _build_prior(self, name, X, **kwargs):
n_obs = np.shape(X)[0]

# standardize input scale
X = pt.as_tensor_variable(X)
Xmu = pt.mean(X, axis=0)
Xsd = pt.std(X, axis=0)
Xz = (X - Xmu) / Xsd

# define L using Xz and c
La = pt.abs(pt.min(Xz)) # .eval()?
Lb = pt.max(Xz)
L = self.c * pt.max([La, Lb])

# make basis and omega, spectral density
Phi, omega = self._generate_basis(Xz, L)
scale, ls, spectral_density = self._validate_cov_func(self.cov_func)
spd = scale * spectral_density(omega, ls / Xsd).flatten()

beta = pm.Normal(f"{name}_coeffs_", size=self.M)
f = pm.Deterministic(name, self.mean_func(X) + pt.dot(Phi * pt.sqrt(spd), beta))
return f, Phi, L, spd, beta, Xmu, Xsd

def _build_conditional(self, Xnew, Xmu, Xsd, L, beta):
Xnewz = (Xnew - Xmu) / Xsd
Phi, omega = self._generate_basis(Xnewz, L)
scale, ls, spectral_density = self._validate_cov_func(self.cov_func)
spd = scale * spectral_density(omega, ls / Xsd).flatten()
return self.mean_func(Xnew) + pt.dot(Phi * pt.sqrt(spd), beta)

def conditional(self, name, Xnew):
# warn about extrapolation
fnew = self._build_conditional(Xnew, self.Xmu, self.Xsd, self.L, self.beta)
return pm.Deterministic(name, fnew)


class ExpQuad(pm.gp.cov.ExpQuad):
@staticmethod
def spectral_density(omega, ls):
# univariate spectral denisty, implement multi
return pt.sqrt(2 * np.pi) * ls * pt.exp(-0.5 * ls**2 * omega**2)


class Matern52(pm.gp.cov.Matern52):
@staticmethod
def spectral_density(omega, ls):
# univariate spectral denisty, implement multi
# https://arxiv.org/pdf/1611.06740.pdf
lam = pt.sqrt(5) * (1.0 / ls)
return (16.0 / 3.0) * lam**5 * (1.0 / (lam**2 + omega**2) ** 3)


class Matern32(pm.gp.cov.Matern32):
@staticmethod
def spectral_density(omega, ls):
# univariate spectral denisty, implement multi
# https://arxiv.org/pdf/1611.06740.pdf
lam = np.sqrt(3.0) * (1.0 / ls)
return 4.0 * lam**3 * (1.0 / pt.square(lam**2 + omega**2))


class Matern12(pm.gp.cov.Matern12):
@staticmethod
def spectral_density(omega, ls):
# univariate spectral denisty, implement multi
# https://arxiv.org/pdf/1611.06740.pdf
lam = 1.0 / ls
return 2.0 * lam * (1.0 / (lam**2 + omega**2))


class KarhunenLoeveExpansion(pm.gp.Latent):
def __init__(
self,
Expand Down

0 comments on commit b730449

Please sign in to comment.