diff --git a/gp/src/lib.rs b/gp/src/lib.rs index c536a069..4368146e 100644 --- a/gp/src/lib.rs +++ b/gp/src/lib.rs @@ -88,10 +88,15 @@ mod algorithm; pub mod correlation_models; mod errors; pub mod mean_models; +mod sgp_algorithm; + mod parameters; +mod sgp_parameters; mod utils; pub use algorithm::*; pub use errors::*; pub use parameters::*; +pub use sgp_algorithm::*; +pub use sgp_parameters::*; pub use utils::NormalizedMatrix; diff --git a/gp/src/parameters.rs b/gp/src/parameters.rs index b1068bef..64e6578d 100644 --- a/gp/src/parameters.rs +++ b/gp/src/parameters.rs @@ -7,15 +7,15 @@ use linfa::{Float, ParamGuard}; #[derive(Clone, Debug, PartialEq, Eq)] pub struct GpValidParams, Corr: CorrelationModel> { /// Parameter of the autocorrelation model - theta: Option>, + pub(crate) theta: Option>, /// Regression model representing the mean(x) - mean: Mean, + pub(crate) mean: Mean, /// Correlation model representing the spatial correlation between errors at e(x) and e(x') - corr: Corr, + pub(crate) corr: Corr, /// Optionally apply dimension reduction (KPLS) or not - kpls_dim: Option, + pub(crate) kpls_dim: Option, /// Parameter to improve numerical stability - nugget: F, + pub(crate) nugget: F, } impl Default for GpValidParams { diff --git a/gp/src/sgp_algorithm.rs b/gp/src/sgp_algorithm.rs new file mode 100644 index 00000000..dfd8dcfb --- /dev/null +++ b/gp/src/sgp_algorithm.rs @@ -0,0 +1,1175 @@ +use crate::errors::{GpError, Result}; +use crate::mean_models::*; +use crate::sgp_parameters::{SgpParams, SgpValidParams, SparseMethod}; +use crate::utils::pairwise_differences; +use crate::{correlation_models::*, Inducings}; +use egobox_doe::{Lhs, SamplingMethod}; +use linfa::prelude::{Dataset, DatasetBase, Fit, Float, PredictInplace}; +use linfa_linalg::{cholesky::*, triangular::*}; +use linfa_pls::PlsRegression; +use ndarray::{arr1, s, Array, Array1, Array2, ArrayBase, Axis, Data, Ix2, Zip}; +use ndarray_einsum_beta::*; +use ndarray_stats::QuantileExt; + +use rand_xoshiro::rand_core::SeedableRng; +use rand_xoshiro::Xoshiro256Plus; +#[cfg(feature = "serializable")] +use serde::{Deserialize, Serialize}; +use std::fmt; + +// use ndarray_rand::rand_distr::Normal; +// use ndarray_rand::RandomExt; + +// const LOG10_20: f64 = 1.301_029_995_663_981_3; //f64::log10(20.); +const N_START: usize = 10; // number of optimization restart (aka multistart) + +/// Woodbury data computed during training and used for prediction +#[cfg_attr( + feature = "serializable", + derive(Serialize, Deserialize), + serde(bound(serialize = "F: Serialize", deserialize = "F: Deserialize<'de>")) +)] +#[derive(Debug)] +pub struct WoodburyData { + vec: Array2, + inv: Array2, +} +impl Clone for WoodburyData { + fn clone(&self) -> Self { + WoodburyData { + vec: self.vec.to_owned(), + inv: self.inv.clone(), + } + } +} + +/// Structure for trained Gaussian Process model +#[derive(Debug)] +#[cfg_attr( + feature = "serializable", + derive(Serialize, Deserialize), + serde(bound(serialize = "F: Serialize", deserialize = "F: Deserialize<'de>")) +)] +pub struct SparseGaussianProcess, Corr: CorrelationModel> { + /// Regression model + #[cfg_attr( + feature = "serializable", + serde(bound(serialize = "Mean: Serialize", deserialize = "Mean: Deserialize<'de>")) + )] + mean: Mean, + /// Correlation kernel + #[cfg_attr( + feature = "serializable", + serde(bound(serialize = "Corr: Serialize", deserialize = "Corr: Deserialize<'de>")) + )] + corr: Corr, + /// Sparse method used + method: SparseMethod, + /// Parameter of the autocorrelation model + theta: Array1, + /// Gaussian process variance + sigma2: F, + /// Gaussian noise variance + noise: F, + /// Weights in case of KPLS dimension reduction coming from PLS regression (orig_dim, kpls_dim) + w_star: Array2, + /// Training inputs + xtrain: Array2, + /// Training outputs + ytrain: Array2, + /// Training outputs + inducings: Array2, + /// Data used for prediction + w_data: WoodburyData, +} + +/// Kriging as GP special case when using constant mean and squared exponential correlation +pub type SparseKriging = SgpParams; + +impl SparseKriging { + pub fn params() -> SgpParams { + SgpParams::new(ConstantMean(), SquaredExponentialCorr()) + } +} + +impl, Corr: CorrelationModel> Clone + for SparseGaussianProcess +{ + fn clone(&self) -> Self { + Self { + mean: self.mean, + corr: self.corr, + method: self.method.clone(), + theta: self.theta.to_owned(), + sigma2: self.sigma2, + noise: self.noise, + w_star: self.w_star.to_owned(), + xtrain: self.xtrain.clone(), + ytrain: self.xtrain.clone(), + inducings: self.inducings.clone(), + w_data: self.w_data.clone(), + } + } +} + +impl, Corr: CorrelationModel> fmt::Display + for SparseGaussianProcess +{ + fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { + write!(f, "SGP({}, {})", self.mean, self.corr) + } +} + +impl, Corr: CorrelationModel> + SparseGaussianProcess +{ + /// Gp parameters contructor + pub fn params, NewCorr: CorrelationModel>( + mean: NewMean, + corr: NewCorr, + ) -> SgpParams { + SgpParams::new(mean, corr) + } + + fn compute_k( + &self, + a: &ArrayBase, Ix2>, + b: &ArrayBase, Ix2>, + w_star: &Array2, + theta: &Array1, + sigma2: F, + ) -> Array2 { + // Get pairwise componentwise L1-distances to the input training set + let dx = pairwise_differences(a, b); + // Compute the correlation function + let r = self.corr.value(&dx, theta, w_star); + r.into_shape((a.nrows(), b.nrows())) + .unwrap() + .mapv(|v| v * sigma2) + } + + /// Predict output values at n given `x` points of nx components specified as a (n, nx) matrix. + /// Returns n scalar output values as (n, 1) column vector. + pub fn predict_values(&self, x: &ArrayBase, Ix2>) -> Result> { + let kx = self.compute_k(x, &self.inducings, &self.w_star, &self.theta, self.sigma2); + let mu = kx.dot(&self.w_data.vec); + Ok(mu) + } + + /// Predict variance values at n given `x` points of nx components specified as a (n, nx) matrix. + /// Returns n variance values as (n, 1) column vector. + pub fn predict_variances(&self, x: &ArrayBase, Ix2>) -> Result> { + let kx = self.compute_k(x, &self.inducings, &self.w_star, &self.theta, self.sigma2); + let kxx = Array::from_elem(x.nrows(), self.sigma2); + // let var = (kxx - np.sum(np.dot(self.woodbury_data["inv"].T, Kx) * Kx, 0))[:, None] + println!("{:?}", kx.shape()); + let tmp = self.w_data.inv.t().clone().dot(&kx) * &kx; + println!("{:?}", tmp.shape()); + let var = kxx - (self.w_data.inv.t().clone().dot(&kx) * &kx).sum_axis(Axis(0)); + let var = var.mapv(|v| { + if v < F::cast(1e-15) { + F::cast(1e-15) + self.noise + } else { + v + self.noise + } + }); + Ok(var.insert_axis(Axis(1))) + } + + /// Retrieve number of PLS components 1 <= n <= x dimension + pub fn kpls_dim(&self) -> Option { + if self.w_star.ncols() < self.xtrain.ncols() { + Some(self.w_star.ncols()) + } else { + None + } + } + + /// Retrieve input dimension before kpls dimension reduction if any + pub fn input_dim(&self) -> usize { + self.ytrain.ncols() + } + + /// Retrieve output dimension + pub fn output_dim(&self) -> usize { + self.ytrain.ncols() + } +} + +impl PredictInplace, Array2> + for SparseGaussianProcess +where + F: Float, + D: Data, + Mean: RegressionModel, + Corr: CorrelationModel, +{ + fn predict_inplace(&self, x: &ArrayBase, y: &mut Array2) { + assert_eq!( + x.nrows(), + y.nrows(), + "The number of data points must match the number of output targets." + ); + + let values = self.predict_values(x).expect("GP Prediction"); + *y = values; + } + + fn default_target(&self, x: &ArrayBase) -> Array2 { + Array2::zeros((x.nrows(), self.output_dim())) + } +} + +/// Gausssian Process adaptator to implement `linfa::Predict` trait for variance prediction. +struct GpVariancePredictor<'a, F, Mean, Corr>(&'a SparseGaussianProcess) +where + F: Float, + Mean: RegressionModel, + Corr: CorrelationModel; + +impl<'a, F, D, Mean, Corr> PredictInplace, Array2> + for GpVariancePredictor<'a, F, Mean, Corr> +where + F: Float, + D: Data, + Mean: RegressionModel, + Corr: CorrelationModel, +{ + fn predict_inplace(&self, x: &ArrayBase, y: &mut Array2) { + assert_eq!( + x.nrows(), + y.nrows(), + "The number of data points must match the number of output targets." + ); + + let values = self.0.predict_variances(x).expect("GP Prediction"); + *y = values; + } + + fn default_target(&self, x: &ArrayBase) -> Array2 { + Array2::zeros((x.nrows(), self.0.output_dim())) + } +} + +impl, Corr: CorrelationModel, D: Data> + Fit, ArrayBase, GpError> for SgpValidParams +{ + type Object = SparseGaussianProcess; + + /// Fit GP parameters using maximum likelihood + fn fit( + &self, + dataset: &DatasetBase, ArrayBase>, + ) -> Result { + let x = dataset.records(); + let y = dataset.targets(); + if let Some(d) = self.kpls_dim() { + if *d > x.ncols() { + return Err(GpError::InvalidValue(format!( + "Dimension reduction {} should be smaller than actual \ + training input dimensions {}", + d, + x.ncols() + ))); + }; + } + + let xtrain = x.to_owned(); + let ytrain = y.to_owned(); + + let mut w_star = Array2::eye(x.ncols()); + if let Some(n_components) = self.kpls_dim() { + let ds = Dataset::new(x.to_owned(), y.to_owned()); + w_star = PlsRegression::params(*n_components).fit(&ds).map_or_else( + |e| match e { + linfa_pls::PlsError::PowerMethodConstantResidualError() => { + Ok(Array2::zeros((x.ncols(), *n_components))) + } + err => Err(err), + }, + |v| Ok(v.rotations().0.to_owned()), + )?; + }; + let theta0: Array1<_> = self + .initial_theta() + .clone() + .map_or(Array1::from_elem(w_star.ncols(), F::cast(1e-2)), |v| { + Array::from_vec(v) + }); + let sigma2_0 = F::cast(1e-2); + let noise_0 = F::cast(1e-2); + + let n = theta0.len() + 2; + let mut params_0 = Array1::zeros(n); + params_0.slice_mut(s![..n - 2]).assign(&theta0); + params_0[n - 2] = sigma2_0; + params_0[n - 1] = noise_0; + + let base: f64 = 10.; + let objfn = |x: &[f64], _gradient: Option<&mut [f64]>, _params: &mut ()| -> f64 { + for v in x.iter() { + // check theta as optimizer may give nan values + if v.is_nan() { + // shortcut return worst value wrt to rlf minimization + return f64::INFINITY; + } + } + let input = Array1::from_shape_vec( + (x.len(),), + x.iter().map(|v| F::cast(base.powf(*v))).collect(), + ) + .unwrap(); + + let theta = input.slice(s![..input.len() - 2]); + let sigma2 = input[input.len() - 2]; + let noise = input[input.len() - 1]; + + let theta = theta.mapv(F::cast); + match self.reduced_likelihood( + &theta, + sigma2, + noise, + &w_star, + &xtrain, + &ytrain, + self.nugget(), + ) { + Ok(r) => unsafe { -(*(&r.0 as *const F as *const f64)) }, + Err(_) => f64::INFINITY, + } + }; + + // Multistart: user theta0 + 1e-5, 1e-4, 1e-3, 1e-2, 0.1, 1., 10. + let mut params = Array2::zeros((N_START + 1, params_0.len())); + params.row_mut(0).assign(¶ms_0.mapv(|v| F::log10(v))); + let mut xlimits: Array2 = Array2::zeros((params_0.len(), 2)); + for mut row in xlimits.rows_mut() { + row.assign(&arr1(&[F::cast(-6), F::cast(2)])); + } + // Use a seed here for reproducibility. Do we need to make it truly random + // Probably no, as it is just to get init values spread over + // [1e-6, 20] for multistart thanks to LHS method. + let seeds = Lhs::new(&xlimits) + .kind(egobox_doe::LhsKind::Maximin) + .with_rng(Xoshiro256Plus::seed_from_u64(42)) + .sample(N_START); + Zip::from(params.slice_mut(s![1.., ..]).rows_mut()) + .and(seeds.rows()) + .par_for_each(|mut theta, row| theta.assign(&row)); + + let opt_params = params.map_axis(Axis(1), |p| optimize_theta(objfn, &p.to_owned())); + let opt_index = opt_params.map(|(_, opt_f)| opt_f).argmin().unwrap(); + let opt_params = &(opt_params[opt_index]).0.mapv(F::cast); + // println!("opt_theta={}", opt_theta); + let opt_theta = opt_params.slice(s![..n - 2]).to_owned(); + let opt_sigma2 = opt_params[n - 2]; + let opt_noise = opt_params[n - 1]; + let (_, w_data) = self.reduced_likelihood( + &opt_theta, + opt_sigma2, + opt_noise, + &w_star, + &xtrain, + &ytrain, + self.nugget(), + )?; + let default_z = Array::from_elem((10, xtrain.ncols()), F::one()); // TODO + let z = match self.inducings() { + Inducings::Randomized(_n) => &default_z, + Inducings::Located(z) => z, + }; + println!("{} {} {}", opt_theta, opt_sigma2, opt_noise); + Ok(SparseGaussianProcess { + mean: *self.mean(), + corr: *self.corr(), + method: self.method().clone(), + theta: opt_theta, + sigma2: opt_sigma2, + noise: opt_noise, + w_data, + w_star, + xtrain: xtrain.to_owned(), + ytrain: ytrain.to_owned(), + inducings: z.clone(), + }) + } +} + +impl, Corr: CorrelationModel> SgpValidParams { + /// Constructor of valid params from values + #[doc(hidden)] + pub fn from( + mean: Mean, + corr: Corr, + method: SparseMethod, + theta: Array1, + sigma2: F, + noise: F, + w_star: Array2, + xtrain: Array2, + ytrain: Array2, + inducings: Array2, + w_data: WoodburyData, + ) -> Result> { + // TODO: add some consistency checks + Ok(SparseGaussianProcess { + mean, + corr, + method, + theta, + sigma2, + noise, + w_star, + xtrain, + ytrain, + inducings, + w_data, + }) + } + + /// Compute reduced likelihood function + /// nugget: factor to improve numerical stability + pub fn reduced_likelihood( + &self, + theta: &Array1, + sigma2: F, + noise: F, + w_star: &Array2, + xtrain: &Array2, + ytrain: &Array2, + nugget: F, + ) -> Result<(F, WoodburyData)> { + let (likelihood, w_data) = match self.method() { + SparseMethod::Fitc => self.fitc(theta, sigma2, noise, w_star, xtrain, ytrain, nugget), + SparseMethod::Vfe => self.vfe(theta, sigma2, noise, nugget), + }; + + Ok((likelihood, w_data)) + } + + /// Compute covariance matrix between a and b matrices + fn compute_k( + &self, + a: &ArrayBase, Ix2>, + b: &ArrayBase, Ix2>, + w_star: &Array2, + theta: &Array1, + sigma2: F, + ) -> Array2 { + // Get pairwise componentwise L1-distances to the input training set + let dx = pairwise_differences(a, b); + // Compute the correlation function + let r = self.corr().value(&dx, &theta, w_star); + r.into_shape((a.nrows(), b.nrows())) + .unwrap() + .mapv(|v| v * sigma2) + } + + /// FITC method + pub fn fitc( + &self, + theta: &Array1, + sigma2: F, + noise: F, + w_star: &Array2, + xtrain: &Array2, + ytrain: &Array2, + nugget: F, + ) -> (F, WoodburyData) { + let default_z = Array::from_elem((10, xtrain.ncols()), F::one()); // TODO + let z = match self.inducings() { + Inducings::Randomized(_n) => &default_z, + Inducings::Located(z) => z, + }; + let nz = z.nrows(); + //println!("nz = {}", nz); + // Compute: diag(Knn), Kmm and Kmn + let knn = Array1::from_elem(xtrain.nrows(), sigma2); + let kmm = self.compute_k(z, z, w_star, theta, sigma2) + Array::eye(nz) * nugget; + let kmn = self.compute_k(z, xtrain, w_star, theta, sigma2); + + // Compute (lower) Cholesky decomposition: Kmm = U U^T + let u = kmm.cholesky().unwrap(); + //U = linalg.cholesky(Kmm, lower=True) + + // Compute (upper) Cholesky decomposition: Qnn = V^T V + let ui = u + .solve_triangular(&Array::eye(u.nrows()), UPLO::Lower) + .unwrap(); + // Ui = linalg.inv(U) + let v = ui.dot(&kmn); + //V = Ui @ Kmn + + // Assumption on the gaussian noise on training outputs + let eta2 = noise; + + // Compute diagonal correction: nu = Knn_diag - Qnn_diag + \eta^2 + let nu = knn; + //println!("knn={:?}", nu.shape()); + //println!("kmm={:?}", kmm.shape()); + //println!("v={:?}", v.shape()); + let nu = nu - (v.to_owned() * &v).sum_axis(Axis(0)); + let nu = nu + eta2; + //let nu = knn - np.sum(np.square(V), 0) + eta2; + // Compute beta, the effective noise precision + let beta = nu.mapv(|v| F::one() / v); + //beta = 1.0 / nu + + // Compute (lower) Cholesky decomposition: A = I + V diag(beta) V^T = L L^T + let a = Array::eye(nz) + &(v.to_owned() * beta.to_owned().insert_axis(Axis(0))).dot(&v.t()); + // A = np.eye(self.nz) + V * beta @ V.T + + let l = a.cholesky().unwrap(); + // L = linalg.cholesky(A, lower=True) + let li = l + .solve_triangular(&Array::eye(l.nrows()), UPLO::Lower) + .unwrap(); + // Li = linalg.inv(L) + + // Compute a and b + //println!("ytrain={:?}", ytrain.shape()); + //println!("beta={:?}", beta.shape()); + let a = einsum("ij,i->ij", &[ytrain, &beta]) + .unwrap() + .into_dimensionality::() + .unwrap(); + // a = np.einsum("ij,i->ij", Y, beta) # avoid reshape for mat-vec multiplication + let tmp = li.dot(&v); + let b = tmp.dot(&a); + //println!("b={:?}", b.shape()); + // b = Li @ V @ a + + // Compute marginal log-likelihood + // constant term ignored in reduced likelihood + //let term0 = self.ytrain.nrows() * F::cast(2. * std::f64::consts::PI); + let term1 = nu.mapv(|v| v.ln()).sum(); + let term2 = F::cast(2.) * l.diag().mapv(|v| v.ln()).sum(); + let term3 = (a.t().to_owned()).dot(ytrain)[[0, 0]]; + //let term4 = einsum("ij,ij->", &[&b, &b]).unwrap(); + let term4 = -(b.to_owned() * &b).sum(); + let likelihood = -F::cast(0.5) * (term1 + term2 + term3 + term4); + // likelihood = -0.5 * ( + // # num_data * np.log(2.0 * np.pi) # constant term ignored in reduced likelihood + // +np.sum(np.log(nu)) + // + 2.0 * np.sum(np.log(np.diag(L))) + // + a.T @ Y + // - np.einsum("ij,ij->", b, b) + // ) + + // Store Woodbury vectors for prediction step + // LiUi = Li @ Ui + // LiUiT = LiUi.T + // woodbury_vec = LiUiT @ b + // woodbury_inv = Ui.T @ Ui - LiUiT @ LiUi + let li_ui = li.dot(&ui); + let li_ui_t = li_ui.t(); + + //println!("li_ui={:?}", li_ui.shape()); + //println!("li_ui_t={:?}", li_ui_t); + //println!("b={:?}", b); + + let w_data = WoodburyData { + vec: li_ui_t.dot(&b), + inv: (ui.t()).dot(&ui) - li_ui_t.dot(&li_ui), + }; + + //println!("wdata = {:?}", w_data); + (likelihood, w_data) + } + + fn vfe(&self, _theta: &Array1, _variance: F, _noise: F, _nugget: F) -> (F, WoodburyData) { + todo!() + } +} + +/// Optimize gp hyper parameter theta given an initial guess `theta0` +#[cfg(not(feature = "nlopt"))] +fn optimize_theta(objfn: ObjF, theta0: &Array1) -> (Array1, f64) +where + ObjF: Fn(&[f64], Option<&mut [f64]>, &mut ()) -> f64, + F: Float, +{ + use cobyla::{minimize, Func, StopTols}; + + let base: f64 = 10.; + let cons: Vec<&dyn Func<()>> = vec![]; + let theta_init = theta0 + .map(|v| unsafe { *(v as *const F as *const f64) }) + .into_raw_vec(); + + let initial_step = 0.5; + let ftol_rel = 1e-4; + let maxeval = 15 * theta0.len(); + let bounds = vec![(-6., 2.); theta0.len()]; + + match minimize( + |x, u| objfn(x, None, u), + &theta_init, + &bounds, + &cons, + (), + maxeval, + cobyla::RhoBeg::All(initial_step), + Some(StopTols { + ftol_rel, + ..StopTols::default() + }), + ) { + Ok((_, x_opt, fval)) => { + let thetas_opt = arr1(&x_opt).mapv(|v| base.powf(v)); + let fval = if f64::is_nan(fval) { + f64::INFINITY + } else { + fval + }; + (thetas_opt, fval) + } + Err((status, x_opt, _)) => { + println!("ERROR Cobyla optimizer in GP status={:?}", status); + (arr1(&x_opt).mapv(|v| base.powf(v)), f64::INFINITY) + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + + use linfa::ParamGuard; + use ndarray::{array, Array}; + use ndarray_npy::{read_npy, write_npy}; + use ndarray_rand::rand::seq::SliceRandom; + use ndarray_rand::rand::SeedableRng; + use ndarray_rand::rand_distr::{Normal, Uniform}; + use ndarray_rand::RandomExt; + use rand_xoshiro::Xoshiro256Plus; + + const PI: f64 = std::f64::consts::PI; + + fn f_obj(x: &ArrayBase, Ix2>) -> Array2 { + x.mapv(|v| (3. * PI * v).sin() + 0.3 * (9. * PI * v).cos() + 0.5 * (7. * PI * v).sin()) + } + + #[test] + fn test_sgp() { + let mut rng = Xoshiro256Plus::seed_from_u64(42); + // Generate training data + let nt = 200; + // Variance of the gaussian noise on our trainingg data + let eta2: f64 = 0.01; + let normal = Normal::new(0., eta2.sqrt()).unwrap(); + let gaussian_noise = Array::::random_using((nt, 1), normal, &mut rng); + let xt = 2. * Array::::random_using((nt, 1), Uniform::new(0., 1.), &mut rng) - 1.; + let yt = f_obj(&xt) + gaussian_noise; + + let test_dir = "target/tests"; + std::fs::create_dir_all(test_dir).ok(); + + let file_path = format!("{}/smt_xt.npy", test_dir); + let xt: Array2 = read_npy(file_path).unwrap(); + let file_path = format!("{}/smt_yt.npy", test_dir); + let yt: Array2 = read_npy(file_path).unwrap(); + + let xplot = Array::linspace(-1., 1., 100).insert_axis(Axis(1)); + + let n_inducing = 30; + let mut indices = (0..nt).collect::>(); + indices.shuffle(&mut rng); + let mut z = Array2::zeros((n_inducing, xt.ncols())); + let idx = indices[..n_inducing].to_vec(); + Zip::from(z.rows_mut()) + .and(&Array1::from_vec(idx)) + .for_each(|mut zi, i| zi.assign(&xt.row(*i))); + + let file_path = format!("{}/smt_z.npy", test_dir); + let z: Array2 = read_npy(file_path).unwrap(); + println!("NZ = {}", z.nrows()); + + let sgp = SparseGaussianProcess::::params( + ConstantMean::default(), + SquaredExponentialCorr::default(), + ) + .inducings(z.clone()) + .initial_theta(Some(vec![0.1])) + .fit(&Dataset::new(xt.clone(), yt.clone())) + .expect("GP fit error"); + + let sgp_vals = sgp.predict_values(&xplot).unwrap(); + + let sgp_vars = sgp.predict_variances(&xplot).unwrap(); + + let file_path = format!("{}/sgp_xt.npy", test_dir); + write_npy(file_path, &xt).expect("xt saved"); + let file_path = format!("{}/sgp_yt.npy", test_dir); + write_npy(file_path, &yt).expect("yt saved"); + let file_path = format!("{}/sgp_z.npy", test_dir); + write_npy(file_path, &z).expect("z saved"); + + let file_path = format!("{}/sgp_x.npy", test_dir); + write_npy(file_path, &xplot).expect("x saved"); + + let file_path = format!("{}/sgp_vals.npy", test_dir); + write_npy(file_path, &sgp_vals).expect("sgp vals saved"); + let file_path = format!("{}/sgp_vars.npy", test_dir); + write_npy(file_path, &sgp_vars).expect("sgp vars saved"); + } + + #[test] + fn test_fitc() { + let xt = array![ + [0.92387276], + [-0.41570495], + [-0.51834244], + [-0.79941212], + [-0.96714074], + [0.85905863], + [0.33983309], + [0.57030582], + [-0.43653979], + [0.17282033], + [-0.87208947], + [-0.02874481], + [0.95499028], + [0.75301049], + [-0.3236821], + [0.92314031], + [-0.53659675], + [0.89863764], + [0.88275541], + [0.59840517], + [0.26089587], + [0.74857593], + [-0.41395943], + [0.69788711], + [0.23575338], + [-0.97352628], + [-0.30553296], + [-0.70371828], + [0.96365878], + [-0.04325939], + [-0.00521727], + [0.27894503], + [-0.26283079], + [-0.72619946], + [0.64423547], + [-0.62030418], + [0.02263797], + [-0.55136594], + [-0.80431103], + [0.72438303], + [0.94583898], + [0.92166932], + [0.813111], + [0.54809467], + [-0.3337097], + [-0.83779722], + [-0.18551766], + [-0.53553172], + [-0.73502473], + [-0.89314564], + [0.45118873], + [-0.97714508], + [0.5411615], + [-0.70610671], + [-0.84095583], + [-0.82079393], + [0.34409561], + [-0.50926558], + [-0.15892107], + [0.11473758], + [0.72110235], + [0.45408853], + [-0.45934419], + [-0.7370344], + [-0.88925136], + [-0.39680273], + [-0.4757637], + [-0.08771887], + [0.36656267], + [0.39125089], + [-0.43296231], + [-0.24014609], + [-0.63769808], + [0.57709102], + [-0.88630385], + [0.39399448], + [0.55739079], + [0.55481512], + [-0.48115487], + [-0.25237372], + [0.17519927], + [-0.4543562], + [-0.2582944], + [-0.60589144], + [-0.08028823], + [-0.9107754], + [0.59959177], + [-0.84608711], + [0.0376703], + [-0.3863798], + [0.1550859], + [0.91886668], + [0.29114049], + [-0.92927513], + [-0.13919512], + [0.0200337], + [0.07235499], + [0.36278502], + [-0.4448078], + [-0.74227887], + [-0.21464865], + [0.91281145], + [-0.62573822], + [0.80796791], + [0.0876119], + [-0.08617716], + [0.76408282], + [-0.08279208], + [0.44833527], + [-0.20194936], + [0.80808879], + [0.38005004], + [0.39924411], + [-0.3445592], + [0.51355729], + [0.27212211], + [-0.51995945], + [-0.67892236], + [0.59278295], + [0.91833321], + [-0.08372235], + [0.18196833], + [0.71544529], + [-0.08555309], + [0.90374895], + [0.15150232], + [0.64153424], + [0.81768744], + [0.63104764], + [-0.68117107], + [0.25779688], + [-0.20313148], + [-0.8745741], + [-0.1519355], + [-0.48263187], + [0.69807662], + [-0.93339075], + [0.91796544], + [-0.2892623], + [-0.28658622], + [-0.96734299], + [-0.62953535], + [-0.197481], + [0.85858283], + [-0.80077014], + [0.89060307], + [0.73897706], + [-0.09167521], + [-0.34659824], + [-0.53451174], + [0.22892941], + [-0.93385082], + [-0.96878787], + [-0.14240856], + [-0.86385185], + [-0.49611802], + [-0.55767817], + [-0.49361761], + [-0.73788954], + [-0.97592755], + [-0.76903141], + [0.23696052], + [0.94851243], + [0.98069], + [-0.18189181], + [-0.67409115], + [0.27752351], + [-0.01938931], + [0.97881955], + [-0.86939159], + [0.56646888], + [-0.42320301], + [-0.51716276], + [0.32500914], + [-0.50787363], + [0.33171824], + [0.03461703], + [-0.15182202], + [0.10937562], + [-0.42589696], + [0.41314941], + [-0.17028626], + [-0.27890888], + [0.65731383], + [0.84993382], + [-0.90798538], + [-0.53474601], + [-0.30296126], + [0.62993296], + [0.97098286], + [0.93794341], + [0.80989669], + [-0.40688747], + [0.98402249], + [-0.50115992], + [-0.78818769], + [0.90190522], + [-0.53315949], + [0.37953653], + [-0.88328728] + ]; + let yt = array![ + [1.49604938e+00], + [8.07376194e-01], + [1.39413512e+00], + [-4.93356947e-01], + [-6.28255039e-01], + [1.09355377e+00], + [2.03049694e-01], + [-1.09015072e+00], + [1.19579331e+00], + [7.85179169e-01], + [-8.14238069e-01], + [-2.11363794e-01], + [8.17530692e-01], + [1.33482589e-01], + [-7.05981901e-01], + [1.36245744e+00], + [1.17921756e+00], + [1.48018011e+00], + [1.48724123e+00], + [-5.11232534e-01], + [2.53873172e-01], + [2.17184128e-01], + [8.12652757e-01], + [5.82742457e-01], + [8.55379663e-01], + [-8.87098517e-01], + [-6.77529315e-01], + [-3.26126672e-01], + [6.92440375e-01], + [-5.54453333e-01], + [2.05843685e-01], + [4.44202508e-01], + [-3.41246114e-01], + [-6.34343524e-01], + [4.96800247e-01], + [7.61669265e-02], + [8.14149554e-01], + [9.20941252e-01], + [-7.62253363e-01], + [3.58868632e-01], + [8.36275513e-01], + [1.20535170e+00], + [2.36963332e-01], + [-1.24274435e+00], + [-7.82468114e-01], + [-7.98801480e-01], + [-5.53868793e-01], + [1.12397933e+00], + [-6.48021291e-01], + [-9.24508060e-01], + [-9.29583831e-01], + [-6.55417766e-01], + [-1.56103604e+00], + [-4.38709243e-01], + [-7.61831419e-01], + [-6.95887015e-01], + [9.87355308e-02], + [1.43839544e+00], + [-9.52964322e-01], + [8.37765720e-01], + [3.58410881e-01], + [-9.20954343e-01], + [1.43306148e+00], + [-6.70439602e-01], + [-8.71022078e-01], + [2.67949598e-01], + [1.43151602e+00], + [-1.39444012e+00], + [-8.65527160e-02], + [-1.28119627e-01], + [1.21239970e+00], + [-7.31427195e-02], + [9.16685909e-02], + [-1.05492474e+00], + [-8.37628001e-01], + [-2.21800752e-01], + [-1.39583732e+00], + [-1.40595026e+00], + [1.56317991e+00], + [-1.53443204e-01], + [6.25542276e-01], + [1.55583567e+00], + [-1.62573924e-01], + [5.51829523e-04], + [-1.22134415e+00], + [-7.73497897e-01], + [-2.78758355e-01], + [-7.84554800e-01], + [7.54301251e-01], + [1.63398796e-01], + [7.24398852e-01], + [1.50172983e+00], + [3.56941530e-01], + [-8.95805245e-01], + [-1.18154204e+00], + [7.24720645e-01], + [9.94038969e-01], + [1.98859988e-01], + [1.35516697e+00], + [-4.86054345e-01], + [8.19676635e-02], + [1.30177136e+00], + [-9.52225114e-02], + [4.30007289e-01], + [8.50212505e-01], + [-1.23385671e+00], + [3.08275807e-02], + [-1.47156589e+00], + [-6.03830756e-01], + [-6.34039326e-02], + [5.21492866e-01], + [2.80907213e-02], + [-2.80967099e-01], + [-4.69381286e-01], + [-1.60879187e+00], + [5.26206431e-01], + [1.36922961e+00], + [-1.99397948e-01], + [-5.02335739e-01], + [1.47741916e+00], + [-1.36834256e+00], + [6.26506988e-01], + [5.17911964e-01], + [-1.29019806e+00], + [1.41953538e+00], + [7.55630184e-01], + [4.49004982e-01], + [6.64235015e-01], + [3.81324072e-01], + [-1.53179711e-01], + [4.48859419e-01], + [-1.45218771e-01], + [-9.04086969e-01], + [-1.00940699e+00], + [1.52849017e+00], + [7.23115442e-01], + [-9.34641068e-01], + [1.36823110e+00], + [-4.98603556e-01], + [-6.19304014e-01], + [-9.62020510e-01], + [5.74046486e-02], + [-2.45671541e-01], + [1.24759100e+00], + [-4.80924667e-01], + [1.58753127e+00], + [1.43471070e-01], + [-1.35583870e+00], + [-7.72660570e-01], + [1.01574139e+00], + [6.46210065e-01], + [-8.20611577e-01], + [-8.71859416e-01], + [-1.25156140e+00], + [-8.14348346e-01], + [1.46402410e+00], + [8.17953850e-01], + [1.43912533e+00], + [-6.17272004e-01], + [-7.54370054e-01], + [-6.95865216e-01], + [8.16619046e-01], + [9.79639090e-01], + [1.39359009e-01], + [-6.08669223e-01], + [-7.86833410e-02], + [3.14633053e-01], + [-2.86960334e-01], + [2.93868680e-01], + [-7.88614939e-01], + [-1.05828298e+00], + [9.69780217e-01], + [1.39750383e+00], + [1.01809826e-01], + [1.32011089e+00], + [2.07555561e-01], + [7.52491896e-01], + [-1.08349928e+00], + [8.48313905e-01], + [9.97911173e-01], + [-3.62317505e-01], + [-8.22624958e-01], + [-4.90143818e-01], + [4.54140638e-01], + [1.10706979e+00], + [-1.11504798e+00], + [9.47867772e-01], + [-6.58281827e-01], + [2.18874692e-01], + [5.17658963e-01], + [9.67056863e-01], + [3.86897104e-01], + [5.51768589e-01], + [-6.46088554e-02], + [1.54226825e+00], + [-7.15830378e-01], + [1.57188675e+00], + [1.16492529e+00], + [1.56678108e-01], + [-7.33200252e-01] + ]; + let z = array![ + [-0.18551766], + [0.90190522], + [-0.97352628], + [-0.2892623], + [-0.82079393], + [0.59959177], + [0.23575338], + [0.81768744], + [0.51355729], + [-0.15892107], + [-0.43653979], + [0.39125089], + [-0.15182202], + [0.55481512], + [-0.42320301], + [-0.55767817], + [0.91886668], + [0.85858283], + [0.27894503], + [-0.4448078], + [0.29114049], + [-0.34659824], + [0.64423547], + [-0.01938931], + [0.07235499], + [-0.02874481], + [-0.49361761], + [0.91796544], + [-0.93339075], + [-0.62030418] + ]; + // th[0.01]; 0.7575623927583945 0.01 2.220446049250313e-13 + //[[-6330.83001996]] + let sgp = SparseGaussianProcess::::params( + ConstantMean::default(), + SquaredExponentialCorr::default(), + ) + .inducings(z.clone()) + .initial_theta(Some(vec![0.01])); + + let theta = array![0.01]; + let sigma2 = 0.7575623927583945; + let noise = 0.01; + let w_star = array![[1.]]; + let nugget = 1000.0 * f64::EPSILON; + let (likelihood, _) = sgp + .check() + .unwrap() + .fitc(&theta, sigma2, noise, &w_star, &xt, &yt, nugget); + println!("{}", likelihood); + } +} diff --git a/gp/src/sgp_parameters.rs b/gp/src/sgp_parameters.rs new file mode 100644 index 00000000..473e0144 --- /dev/null +++ b/gp/src/sgp_parameters.rs @@ -0,0 +1,200 @@ +use crate::correlation_models::{CorrelationModel, SquaredExponentialCorr}; +use crate::errors::{GpError, Result}; +use crate::mean_models::{ConstantMean, RegressionModel}; +use crate::parameters::GpValidParams; +use linfa::{Float, ParamGuard}; +use ndarray::Array2; +#[cfg(feature = "serializable")] +use serde::{Deserialize, Serialize}; + +#[derive(Clone, Debug, PartialEq, Eq)] +pub enum NoiseVarianceConfig { + Constant(F), + Estimated { + initial_guess: F, + lower_bound: F, + upper_bound: F, + }, +} +impl Default for NoiseVarianceConfig { + fn default() -> NoiseVarianceConfig { + Self::Constant(F::one()) + } +} + +#[derive(Clone, Debug, PartialEq, Eq)] +#[cfg_attr(feature = "serializable", derive(Serialize, Deserialize))] +pub enum Inducings { + Randomized(usize), + Located(Array2), +} +impl Default for Inducings { + fn default() -> Inducings { + Self::Randomized(10) + } +} + +#[derive(Clone, Debug, PartialEq, Eq, Default)] +#[cfg_attr(feature = "serializable", derive(Serialize, Deserialize))] +pub enum SparseMethod { + #[default] + Fitc, + Vfe, +} + +/// A set of validated SGP parameters. +#[derive(Clone, Debug, PartialEq, Eq)] +pub struct SgpValidParams, Corr: CorrelationModel> { + /// gp + gp_params: GpValidParams, + /// Gaussian homeoscedastic noise variance + noise: NoiseVarianceConfig, + /// Inducing points + z: Inducings, + /// Method + method: SparseMethod, +} + +impl Default for SgpValidParams { + fn default() -> SgpValidParams { + SgpValidParams { + gp_params: GpValidParams::default(), + noise: NoiseVarianceConfig::default(), + z: Inducings::default(), + method: SparseMethod::default(), + } + } +} + +impl, Corr: CorrelationModel> SgpValidParams { + /// Get starting theta value for optimization + pub fn initial_theta(&self) -> &Option> { + &self.gp_params.theta + } + + /// Get mean model + pub fn mean(&self) -> &Mean { + &self.gp_params.mean + } + + /// Get correlation corr k(x, x') + pub fn corr(&self) -> &Corr { + &self.gp_params.corr + } + + /// Get number of components used by PLS + pub fn kpls_dim(&self) -> &Option { + &self.gp_params.kpls_dim + } + + /// Get number of components used by PLS + pub fn nugget(&self) -> F { + self.gp_params.nugget + } + + /// Get used sparse method + pub fn method(&self) -> &SparseMethod { + &self.method + } + + /// Get inducing points + pub fn inducings(&self) -> &Inducings { + &self.z + } +} + +#[derive(Clone, Debug)] +/// The set of hyperparameters that can be specified for the execution of +/// the [GP algorithm](struct.GaussianProcess.html). +pub struct SgpParams, Corr: CorrelationModel>( + SgpValidParams, +); + +impl, Corr: CorrelationModel> SgpParams { + /// A constructor for GP parameters given mean and correlation models + pub fn new(mean: Mean, corr: Corr) -> SgpParams { + Self(SgpValidParams { + gp_params: GpValidParams { + theta: None, + mean, + corr, + kpls_dim: None, + nugget: F::cast(1000.0) * F::epsilon(), + }, + noise: NoiseVarianceConfig::default(), + z: Inducings::default(), + method: SparseMethod::default(), + }) + } + + /// Set initial value for theta hyper parameter. + /// + /// During training process, the internal optimization is started from `initial_theta`. + pub fn initial_theta(mut self, theta: Option>) -> Self { + self.0.gp_params.theta = theta; + self + } + + /// Set mean model. + pub fn mean(mut self, mean: Mean) -> Self { + self.0.gp_params.mean = mean; + self + } + + /// Set correlation model. + pub fn corr(mut self, corr: Corr) -> Self { + self.0.gp_params.corr = corr; + self + } + + /// Set number of PLS components. + /// Should be 0 < n < pb size (i.e. x dimension) + pub fn kpls_dim(mut self, kpls_dim: Option) -> Self { + self.0.gp_params.kpls_dim = kpls_dim; + self + } + + /// Set nugget. + /// + /// Nugget is used to improve numerical stability + pub fn nugget(mut self, nugget: F) -> Self { + self.0.gp_params.nugget = nugget; + self + } + + pub fn inducings(mut self, z: Array2) -> Self { + self.0.z = Inducings::Located(z); + self + } +} + +impl, Corr: CorrelationModel> ParamGuard + for SgpParams +{ + type Checked = SgpValidParams; + type Error = GpError; + + fn check_ref(&self) -> Result<&Self::Checked> { + if let Some(d) = self.0.gp_params.kpls_dim { + if d == 0 { + return Err(GpError::InvalidValue("`kpls_dim` canot be 0!".to_string())); + } + if let Some(theta) = self.0.initial_theta() { + if theta.len() > 1 && d > theta.len() { + return Err(GpError::InvalidValue(format!( + "Dimension reduction ({}) should be smaller than expected + training input size infered from given initial theta length ({})", + d, + theta.len() + ))); + }; + } + } + Ok(&self.0) + } + + fn check(self) -> Result { + self.check_ref()?; + Ok(self.0) + } +}