Multi-Fidelity GPR

Delta multi-fidelity GPR for simulation-plus-experiment and other low-/high-fidelity materials workflows.

matgpr.multifidelity_gpr

MultiFidelityObservationData dataclass

Validated observation table for ordered multi-fidelity datasets.

Attributes:

Name	Type	Description
X, y		Feature matrix and scalar target values for every observed row.
fidelity_index	ndarray	Integer fidelity level for each row. The order is defined by fidelity_names.
fidelity_names	tuple[str, ...]	Ordered fidelity labels, usually from cheapest/lowest to target/highest fidelity.
target_fidelity	str	Fidelity label used as the default decision target.
sample_id	ndarray | None	Optional material/sample identifiers aligned with rows of X.
noise_variance	ndarray | None	Optional known per-observation noise variances aligned with rows of X.
feature_names	tuple[str, ...] | None	Optional descriptor names aligned with columns of X.

Source code in matgpr/multifidelity_gpr.py

@dataclass(frozen=True)
class MultiFidelityObservationData:
    """Validated observation table for ordered multi-fidelity datasets.

    Attributes
    ----------
    X, y
        Feature matrix and scalar target values for every observed row.
    fidelity_index
        Integer fidelity level for each row. The order is defined by
        ``fidelity_names``.
    fidelity_names
        Ordered fidelity labels, usually from cheapest/lowest to target/highest
        fidelity.
    target_fidelity
        Fidelity label used as the default decision target.
    sample_id
        Optional material/sample identifiers aligned with rows of ``X``.
    noise_variance
        Optional known per-observation noise variances aligned with rows of
        ``X``.
    feature_names
        Optional descriptor names aligned with columns of ``X``.
    """

    X: np.ndarray
    y: np.ndarray
    fidelity_index: np.ndarray
    fidelity_names: tuple[str, ...]
    target_fidelity: str
    sample_id: np.ndarray | None = None
    noise_variance: np.ndarray | None = None
    feature_names: tuple[str, ...] | None = None

    @property
    def n_samples(self) -> int:
        """Number of observed fidelity rows."""
        return int(self.X.shape[0])

    @property
    def n_features(self) -> int:
        """Number of feature columns."""
        return int(self.X.shape[1])

    @property
    def n_fidelities(self) -> int:
        """Number of ordered fidelity levels."""
        return len(self.fidelity_names)

    @property
    def target_fidelity_index(self) -> int:
        """Integer index of the target fidelity."""
        return self.fidelity_names.index(self.target_fidelity)

    @property
    def fidelity_labels(self) -> np.ndarray:
        """Row-aligned string fidelity labels."""
        return np.asarray(
            [self.fidelity_names[index] for index in self.fidelity_index],
            dtype=object,
        )

    @property
    def fidelity_observation_counts(self) -> dict[str, int]:
        """Observation counts keyed by fidelity label."""
        return {
            name: int(np.sum(self.fidelity_index == index))
            for index, name in enumerate(self.fidelity_names)
        }

    @property
    def target_rows(self) -> np.ndarray:
        """Row indices belonging to the target fidelity."""
        return self.rows_for_fidelity(self.target_fidelity)

    def rows_for_fidelity(self, fidelity: str | int) -> np.ndarray:
        """Return row indices for a fidelity label or integer index."""
        fidelity_index = self._resolve_fidelity_index(fidelity)
        return np.flatnonzero(self.fidelity_index == fidelity_index)

    def _resolve_fidelity_index(self, fidelity: str | int) -> int:
        if isinstance(fidelity, (int, np.integer)):
            index = int(fidelity)
            if 0 <= index < self.n_fidelities:
                return index
            raise ValueError(
                f"fidelity index must be between 0 and {self.n_fidelities - 1}; got {index}"
            )

        label = str(fidelity)
        if label not in self.fidelity_names:
            raise ValueError(f"Unknown fidelity {label!r}; expected one of {self.fidelity_names}")
        return self.fidelity_names.index(label)

n_samples property

Number of observed fidelity rows.

n_features property

Number of feature columns.

n_fidelities property

Number of ordered fidelity levels.

target_fidelity_index property

Integer index of the target fidelity.

fidelity_labels property

Row-aligned string fidelity labels.

fidelity_observation_counts property

Observation counts keyed by fidelity label.

target_rows property

Row indices belonging to the target fidelity.

rows_for_fidelity(fidelity)

Return row indices for a fidelity label or integer index.

Source code in matgpr/multifidelity_gpr.py

def rows_for_fidelity(self, fidelity: str | int) -> np.ndarray:
    """Return row indices for a fidelity label or integer index."""
    fidelity_index = self._resolve_fidelity_index(fidelity)
    return np.flatnonzero(self.fidelity_index == fidelity_index)

MultiFidelityGPRPrediction dataclass

Prediction from a delta multi-fidelity Gaussian-process model.

Attributes:

Name	Type	Description
mean, std, lower, upper		High-fidelity predictive mean, standard deviation, and optional confidence interval in original high-fidelity target units.
low_fidelity_mean, low_fidelity_std		Low-fidelity values or low-fidelity surrogate predictions used in the autoregressive correction.
correction_mean, correction_std		GP correction term delta(x) and its uncertainty.
rho, intercept		Fitted linear mapping from low fidelity to high fidelity.

Source code in matgpr/multifidelity_gpr.py

@dataclass(frozen=True)
class MultiFidelityGPRPrediction:
    """Prediction from a delta multi-fidelity Gaussian-process model.

    Attributes
    ----------
    mean, std, lower, upper
        High-fidelity predictive mean, standard deviation, and optional
        confidence interval in original high-fidelity target units.
    low_fidelity_mean, low_fidelity_std
        Low-fidelity values or low-fidelity surrogate predictions used in the
        autoregressive correction.
    correction_mean, correction_std
        GP correction term ``delta(x)`` and its uncertainty.
    rho, intercept
        Fitted linear mapping from low fidelity to high fidelity.
    """

    mean: np.ndarray
    std: np.ndarray | None = None
    lower: np.ndarray | None = None
    upper: np.ndarray | None = None
    low_fidelity_mean: np.ndarray | None = None
    low_fidelity_std: np.ndarray | None = None
    correction_mean: np.ndarray | None = None
    correction_std: np.ndarray | None = None
    rho: float = 1.0
    intercept: float = 0.0

CoKrigingGPRPrediction dataclass

Prediction from a two-level autoregressive co-kriging model.

Source code in matgpr/multifidelity_gpr.py

@dataclass(frozen=True)
class CoKrigingGPRPrediction:
    """Prediction from a two-level autoregressive co-kriging model."""

    mean: np.ndarray
    std: np.ndarray | None = None
    lower: np.ndarray | None = None
    upper: np.ndarray | None = None
    low_fidelity_mean: np.ndarray | None = None
    low_fidelity_std: np.ndarray | None = None
    scaled_low_fidelity_mean: np.ndarray | None = None
    scaled_low_fidelity_std: np.ndarray | None = None
    discrepancy_mean: np.ndarray | None = None
    discrepancy_std: np.ndarray | None = None
    fidelity: str = ""
    rho: float = 1.0

CoKrigingGPRResult dataclass

Fitted two-level autoregressive co-kriging GPR result.

Source code in matgpr/multifidelity_gpr.py

@dataclass(frozen=True)
class CoKrigingGPRResult:
    """Fitted two-level autoregressive co-kriging GPR result."""

    model: "ExactTwoLevelCoKrigingGPRModel"
    likelihood: gpytorch.likelihoods.GaussianLikelihood
    observation_data: MultiFidelityObservationData
    loss_history: list[float]
    target_mean: float
    target_std: float
    standardize_y: bool
    low_fidelity: str
    target_fidelity: str
    low_fidelity_index: int
    target_fidelity_index: int
    fidelity_observation_counts: dict[str, int]
    low_fidelity_kernel: str
    discrepancy_kernel: str
    ard: bool
    noise_mode: str
    noise_variance: float
    rho: float
    device: str
    dtype: torch.dtype

    def predict(
        self,
        X,
        *,
        target_fidelity: str | int | None = None,
        return_std: bool = True,
        confidence_level: float | None = None,
        include_observation_noise: bool = True,
    ) -> CoKrigingGPRPrediction:
        """Predict one modeled fidelity for new samples."""
        return _predict_cokriging_gpr(
            self,
            X,
            target_fidelity=target_fidelity,
            return_std=return_std,
            confidence_level=confidence_level,
            include_observation_noise=include_observation_noise,
        )

predict(X, *, target_fidelity=None, return_std=True, confidence_level=None, include_observation_noise=True)

Predict one modeled fidelity for new samples.

Source code in matgpr/multifidelity_gpr.py

def predict(
    self,
    X,
    *,
    target_fidelity: str | int | None = None,
    return_std: bool = True,
    confidence_level: float | None = None,
    include_observation_noise: bool = True,
) -> CoKrigingGPRPrediction:
    """Predict one modeled fidelity for new samples."""
    return _predict_cokriging_gpr(
        self,
        X,
        target_fidelity=target_fidelity,
        return_std=return_std,
        confidence_level=confidence_level,
        include_observation_noise=include_observation_noise,
    )

ExactTwoLevelCoKrigingGPRModel

Bases: ExactGP

Exact two-level autoregressive co-kriging model.

The covariance follows

f_high(x) = rho * f_low(x) + delta(x),

with independent low-fidelity and discrepancy kernels. The first implementation supports exactly two fidelity labels and a shared learned Gaussian likelihood noise term.

Source code in matgpr/multifidelity_gpr.py

class ExactTwoLevelCoKrigingGPRModel(gpytorch.models.ExactGP):
    """Exact two-level autoregressive co-kriging model.

    The covariance follows

    ``f_high(x) = rho * f_low(x) + delta(x)``,

    with independent low-fidelity and discrepancy kernels. The first
    implementation supports exactly two fidelity labels and a shared learned
    Gaussian likelihood noise term.
    """

    def __init__(
        self,
        train_x: torch.Tensor,
        train_fidelity_indices: torch.Tensor,
        train_y: torch.Tensor,
        likelihood: gpytorch.likelihoods.GaussianLikelihood,
        *,
        low_fidelity_index: int,
        target_fidelity_index: int,
        low_fidelity_kernel: str = "matern",
        discrepancy_kernel: str = "matern",
        ard_num_dims: int | None = None,
        initial_rho: float = 1.0,
    ):
        super().__init__((train_x, train_fidelity_indices), train_y, likelihood)
        self.low_fidelity_index = int(low_fidelity_index)
        self.target_fidelity_index = int(target_fidelity_index)
        self.register_parameter(
            name="rho",
            parameter=torch.nn.Parameter(
                torch.as_tensor(float(initial_rho), dtype=train_x.dtype, device=train_x.device)
            ),
        )
        self.register_parameter(
            name="low_mean",
            parameter=torch.nn.Parameter(
                torch.zeros((), dtype=train_x.dtype, device=train_x.device)
            ),
        )
        self.register_parameter(
            name="discrepancy_mean",
            parameter=torch.nn.Parameter(
                torch.zeros((), dtype=train_x.dtype, device=train_x.device)
            ),
        )
        self.low_covar_module = gpytorch.kernels.ScaleKernel(
            _make_gpytorch_base_kernel(low_fidelity_kernel, ard_num_dims=ard_num_dims)
        )
        self.discrepancy_covar_module = gpytorch.kernels.ScaleKernel(
            _make_gpytorch_base_kernel(discrepancy_kernel, ard_num_dims=ard_num_dims)
        )
        # Buffers so target standardization survives ``state_dict()``.
        self.register_buffer(
            "target_mean", torch.zeros((), dtype=train_x.dtype, device=train_x.device)
        )
        self.register_buffer(
            "target_std", torch.ones((), dtype=train_x.dtype, device=train_x.device)
        )

    def forward(
        self,
        x: torch.Tensor,
        fidelity_indices: torch.Tensor,
    ) -> gpytorch.distributions.MultivariateNormal:
        fidelity_indices = fidelity_indices.long().reshape(-1)
        if x.shape[0] != fidelity_indices.shape[0]:
            raise ValueError("x and fidelity_indices must contain the same number of rows")

        mean_x = self._mean_for_fidelity(fidelity_indices, reference=x)
        covar_x = self._covariance_for_fidelity(x, fidelity_indices)
        return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)

    def _mean_for_fidelity(
        self,
        fidelity_indices: torch.Tensor,
        *,
        reference: torch.Tensor,
    ) -> torch.Tensor:
        rho = self.rho.to(dtype=reference.dtype, device=reference.device)
        low_mean = self.low_mean.to(dtype=reference.dtype, device=reference.device)
        discrepancy_mean = self.discrepancy_mean.to(dtype=reference.dtype, device=reference.device)
        target_mean = rho * low_mean + discrepancy_mean
        return torch.where(
            fidelity_indices == self.target_fidelity_index,
            target_mean.expand(fidelity_indices.shape[0]),
            low_mean.expand(fidelity_indices.shape[0]),
        )

    def _covariance_for_fidelity(
        self,
        x: torch.Tensor,
        fidelity_indices: torch.Tensor,
    ):
        rho = self.rho.to(dtype=x.dtype, device=x.device)
        is_target = fidelity_indices == self.target_fidelity_index
        low_weights = torch.where(is_target, rho.expand_as(x[:, 0]), torch.ones_like(x[:, 0]))
        discrepancy_weights = torch.where(
            is_target,
            torch.ones_like(x[:, 0]),
            torch.zeros_like(x[:, 0]),
        )
        low_weight_matrix = low_weights.unsqueeze(-1) * low_weights.unsqueeze(-2)
        discrepancy_weight_matrix = discrepancy_weights.unsqueeze(
            -1
        ) * discrepancy_weights.unsqueeze(-2)

        low_covar = self.low_covar_module(x).mul(low_weight_matrix)
        discrepancy_covar = self.discrepancy_covar_module(x).mul(discrepancy_weight_matrix)
        return low_covar + discrepancy_covar

DeltaMultiFidelityGPRResult dataclass

Fitted two-stage delta multi-fidelity GPR result.

Source code in matgpr/multifidelity_gpr.py

@dataclass(frozen=True)
class DeltaMultiFidelityGPRResult:
    """Fitted two-stage delta multi-fidelity GPR result."""

    correction_model: GPyTorchGPRResult
    low_fidelity_model: GPyTorchGPRResult | None
    rho: float
    intercept: float
    high_fidelity_target_mean: float
    high_fidelity_target_std: float
    correction_target: np.ndarray
    low_fidelity_at_high: np.ndarray
    fit_intercept: bool
    include_low_fidelity_uncertainty: bool
    device: str
    dtype: torch.dtype

    def predict(
        self,
        X,
        *,
        low_fidelity=None,
        return_std: bool = True,
        confidence_level: float | None = None,
        include_observation_noise: bool = True,
        include_low_fidelity_uncertainty: bool | None = None,
    ) -> MultiFidelityGPRPrediction:
        """Predict high-fidelity values for new samples."""
        return _predict_delta_multifidelity_gpr(
            self,
            X,
            low_fidelity=low_fidelity,
            return_std=return_std,
            confidence_level=confidence_level,
            include_observation_noise=include_observation_noise,
            include_low_fidelity_uncertainty=include_low_fidelity_uncertainty,
        )

predict(X, *, low_fidelity=None, return_std=True, confidence_level=None, include_observation_noise=True, include_low_fidelity_uncertainty=None)

Predict high-fidelity values for new samples.

Source code in matgpr/multifidelity_gpr.py

def predict(
    self,
    X,
    *,
    low_fidelity=None,
    return_std: bool = True,
    confidence_level: float | None = None,
    include_observation_noise: bool = True,
    include_low_fidelity_uncertainty: bool | None = None,
) -> MultiFidelityGPRPrediction:
    """Predict high-fidelity values for new samples."""
    return _predict_delta_multifidelity_gpr(
        self,
        X,
        low_fidelity=low_fidelity,
        return_std=return_std,
        confidence_level=confidence_level,
        include_observation_noise=include_observation_noise,
        include_low_fidelity_uncertainty=include_low_fidelity_uncertainty,
    )

MultiFidelityGPRRegressor

Bases: RegressorMixin, BaseEstimator

Scikit-learn-style delta multi-fidelity GPR estimator.

The estimator models scarce high-fidelity measurements by combining low-fidelity values with a GPR correction. Pass low_fidelity to :meth:fit and :meth:predict when low-fidelity values are externally available. If X_low and y_low are passed to :meth:fit, the estimator fits an internal low-fidelity surrogate and can predict without explicit low-fidelity values.

Source code in matgpr/multifidelity_gpr.py

class MultiFidelityGPRRegressor(RegressorMixin, BaseEstimator):
    """Scikit-learn-style delta multi-fidelity GPR estimator.

    The estimator models scarce high-fidelity measurements by combining
    low-fidelity values with a GPR correction. Pass ``low_fidelity`` to
    :meth:`fit` and :meth:`predict` when low-fidelity values are externally
    available. If ``X_low`` and ``y_low`` are passed to :meth:`fit`, the
    estimator fits an internal low-fidelity surrogate and can predict without
    explicit low-fidelity values.
    """

    def __init__(
        self,
        *,
        fit_intercept: bool = True,
        correction_kernel: str = "matern",
        low_fidelity_kernel: str = "matern",
        ard: bool = True,
        lr: float = 0.01,
        training_iter: int = 1000,
        low_fidelity_training_iter: int | None = None,
        correction_initial_noise: float | None = 0.1,
        low_fidelity_initial_noise: float | None = 0.1,
        standardize_y: bool = True,
        include_low_fidelity_uncertainty: bool = True,
        device: str = "cpu",
        dtype: str | torch.dtype = "float64",
        verbose: bool = False,
        log_every: int = 100,
        include_observation_noise: bool = True,
        random_state: int | None = None,
    ):
        self.fit_intercept = fit_intercept
        self.correction_kernel = correction_kernel
        self.low_fidelity_kernel = low_fidelity_kernel
        self.ard = ard
        self.lr = lr
        self.training_iter = training_iter
        self.low_fidelity_training_iter = low_fidelity_training_iter
        self.correction_initial_noise = correction_initial_noise
        self.low_fidelity_initial_noise = low_fidelity_initial_noise
        self.standardize_y = standardize_y
        self.include_low_fidelity_uncertainty = include_low_fidelity_uncertainty
        self.device = device
        self.dtype = dtype
        self.verbose = verbose
        self.log_every = log_every
        self.include_observation_noise = include_observation_noise
        self.random_state = random_state

    def fit(self, X, y, *, low_fidelity=None, X_low=None, y_low=None):
        """Fit the multi-fidelity GPR model."""
        _seed_torch(self.random_state)
        X_checked = check_array(X, ensure_2d=True, dtype="numeric", ensure_all_finite=True)
        y_checked = column_or_1d(
            check_array(
                y,
                ensure_2d=False,
                dtype="numeric",
                ensure_all_finite=True,
                input_name="y",
            ),
            warn=True,
        )
        check_consistent_length(X_checked, y_checked)
        self.n_features_in_ = X_checked.shape[1]

        self.result_ = fit_delta_multifidelity_gpr(
            X_checked,
            y_checked,
            low_fidelity_high=low_fidelity,
            X_low=X_low,
            y_low=y_low,
            fit_intercept=self.fit_intercept,
            correction_kernel=self.correction_kernel,
            low_fidelity_kernel=self.low_fidelity_kernel,
            ard=self.ard,
            lr=self.lr,
            training_iter=self.training_iter,
            low_fidelity_training_iter=self.low_fidelity_training_iter,
            correction_initial_noise=self.correction_initial_noise,
            low_fidelity_initial_noise=self.low_fidelity_initial_noise,
            standardize_y=self.standardize_y,
            include_low_fidelity_uncertainty=self.include_low_fidelity_uncertainty,
            device=self.device,
            dtype=_resolve_torch_dtype(self.dtype),
            verbose=self.verbose,
            log_every=self.log_every,
        )
        self.correction_model_ = self.result_.correction_model
        self.low_fidelity_model_ = self.result_.low_fidelity_model
        self.rho_ = self.result_.rho
        self.intercept_ = self.result_.intercept
        self.correction_target_ = self.result_.correction_target.copy()
        self.low_fidelity_at_high_ = self.result_.low_fidelity_at_high.copy()
        return self

    def predict(
        self,
        X,
        *,
        low_fidelity=None,
        return_std: bool = False,
        include_observation_noise: bool | None = None,
        include_low_fidelity_uncertainty: bool | None = None,
    ):
        """Predict high-fidelity values."""
        prediction = self.predict_distribution(
            X,
            low_fidelity=low_fidelity,
            return_std=return_std,
            include_observation_noise=include_observation_noise,
            include_low_fidelity_uncertainty=include_low_fidelity_uncertainty,
        )
        if return_std:
            return prediction.mean, prediction.std
        return prediction.mean

    def predict_distribution(
        self,
        X,
        *,
        low_fidelity=None,
        return_std: bool = True,
        confidence_level: float | None = None,
        include_observation_noise: bool | None = None,
        include_low_fidelity_uncertainty: bool | None = None,
    ) -> MultiFidelityGPRPrediction:
        """Return high-fidelity predictive mean, uncertainty, and components."""
        check_is_fitted(self, "result_")
        X_checked = check_array(X, ensure_2d=True, dtype="numeric", ensure_all_finite=True)
        if X_checked.shape[1] != self.n_features_in_:
            raise ValueError(
                f"X has {X_checked.shape[1]} features, but this estimator was fitted with "
                f"{self.n_features_in_} features"
            )
        if include_observation_noise is None:
            include_observation_noise = self.include_observation_noise

        return self.result_.predict(
            X_checked,
            low_fidelity=low_fidelity,
            return_std=return_std,
            confidence_level=confidence_level,
            include_observation_noise=include_observation_noise,
            include_low_fidelity_uncertainty=include_low_fidelity_uncertainty,
        )

    def score(self, X, y, sample_weight=None, *, low_fidelity=None) -> float:
        """Return R2 for high-fidelity predictions."""
        return r2_score(
            y,
            self.predict(X, low_fidelity=low_fidelity),
            sample_weight=sample_weight,
        )

fit(X, y, *, low_fidelity=None, X_low=None, y_low=None)

Fit the multi-fidelity GPR model.

Source code in matgpr/multifidelity_gpr.py

def fit(self, X, y, *, low_fidelity=None, X_low=None, y_low=None):
    """Fit the multi-fidelity GPR model."""
    _seed_torch(self.random_state)
    X_checked = check_array(X, ensure_2d=True, dtype="numeric", ensure_all_finite=True)
    y_checked = column_or_1d(
        check_array(
            y,
            ensure_2d=False,
            dtype="numeric",
            ensure_all_finite=True,
            input_name="y",
        ),
        warn=True,
    )
    check_consistent_length(X_checked, y_checked)
    self.n_features_in_ = X_checked.shape[1]

    self.result_ = fit_delta_multifidelity_gpr(
        X_checked,
        y_checked,
        low_fidelity_high=low_fidelity,
        X_low=X_low,
        y_low=y_low,
        fit_intercept=self.fit_intercept,
        correction_kernel=self.correction_kernel,
        low_fidelity_kernel=self.low_fidelity_kernel,
        ard=self.ard,
        lr=self.lr,
        training_iter=self.training_iter,
        low_fidelity_training_iter=self.low_fidelity_training_iter,
        correction_initial_noise=self.correction_initial_noise,
        low_fidelity_initial_noise=self.low_fidelity_initial_noise,
        standardize_y=self.standardize_y,
        include_low_fidelity_uncertainty=self.include_low_fidelity_uncertainty,
        device=self.device,
        dtype=_resolve_torch_dtype(self.dtype),
        verbose=self.verbose,
        log_every=self.log_every,
    )
    self.correction_model_ = self.result_.correction_model
    self.low_fidelity_model_ = self.result_.low_fidelity_model
    self.rho_ = self.result_.rho
    self.intercept_ = self.result_.intercept
    self.correction_target_ = self.result_.correction_target.copy()
    self.low_fidelity_at_high_ = self.result_.low_fidelity_at_high.copy()
    return self

predict(X, *, low_fidelity=None, return_std=False, include_observation_noise=None, include_low_fidelity_uncertainty=None)

Predict high-fidelity values.

Source code in matgpr/multifidelity_gpr.py

def predict(
    self,
    X,
    *,
    low_fidelity=None,
    return_std: bool = False,
    include_observation_noise: bool | None = None,
    include_low_fidelity_uncertainty: bool | None = None,
):
    """Predict high-fidelity values."""
    prediction = self.predict_distribution(
        X,
        low_fidelity=low_fidelity,
        return_std=return_std,
        include_observation_noise=include_observation_noise,
        include_low_fidelity_uncertainty=include_low_fidelity_uncertainty,
    )
    if return_std:
        return prediction.mean, prediction.std
    return prediction.mean

predict_distribution(X, *, low_fidelity=None, return_std=True, confidence_level=None, include_observation_noise=None, include_low_fidelity_uncertainty=None)

Return high-fidelity predictive mean, uncertainty, and components.

Source code in matgpr/multifidelity_gpr.py

def predict_distribution(
    self,
    X,
    *,
    low_fidelity=None,
    return_std: bool = True,
    confidence_level: float | None = None,
    include_observation_noise: bool | None = None,
    include_low_fidelity_uncertainty: bool | None = None,
) -> MultiFidelityGPRPrediction:
    """Return high-fidelity predictive mean, uncertainty, and components."""
    check_is_fitted(self, "result_")
    X_checked = check_array(X, ensure_2d=True, dtype="numeric", ensure_all_finite=True)
    if X_checked.shape[1] != self.n_features_in_:
        raise ValueError(
            f"X has {X_checked.shape[1]} features, but this estimator was fitted with "
            f"{self.n_features_in_} features"
        )
    if include_observation_noise is None:
        include_observation_noise = self.include_observation_noise

    return self.result_.predict(
        X_checked,
        low_fidelity=low_fidelity,
        return_std=return_std,
        confidence_level=confidence_level,
        include_observation_noise=include_observation_noise,
        include_low_fidelity_uncertainty=include_low_fidelity_uncertainty,
    )

score(X, y, sample_weight=None, *, low_fidelity=None)

Return R2 for high-fidelity predictions.

Source code in matgpr/multifidelity_gpr.py

def score(self, X, y, sample_weight=None, *, low_fidelity=None) -> float:
    """Return R2 for high-fidelity predictions."""
    return r2_score(
        y,
        self.predict(X, low_fidelity=low_fidelity),
        sample_weight=sample_weight,
    )

CoKrigingGPRRegressor

Bases: RegressorMixin, BaseEstimator

Scikit-learn-style two-level autoregressive co-kriging estimator.

Source code in matgpr/multifidelity_gpr.py

class CoKrigingGPRRegressor(RegressorMixin, BaseEstimator):
    """Scikit-learn-style two-level autoregressive co-kriging estimator."""

    def __init__(
        self,
        *,
        fidelity_order: Sequence[str] | None = None,
        low_fidelity: str | None = None,
        target_fidelity: str | None = None,
        low_fidelity_kernel: str = "matern",
        discrepancy_kernel: str = "matern",
        ard: bool = True,
        initial_rho: float | None = None,
        lr: float = 0.01,
        training_iter: int = 1000,
        initial_noise: float | None = 0.1,
        standardize_y: bool = True,
        noise_mode: str = "shared",
        min_observations_per_fidelity: int = 2,
        device: str = "cpu",
        dtype: str | torch.dtype = "float64",
        verbose: bool = False,
        log_every: int = 100,
        include_observation_noise: bool = True,
        random_state: int | None = None,
    ):
        self.fidelity_order = fidelity_order
        self.low_fidelity = low_fidelity
        self.target_fidelity = target_fidelity
        self.low_fidelity_kernel = low_fidelity_kernel
        self.discrepancy_kernel = discrepancy_kernel
        self.ard = ard
        self.initial_rho = initial_rho
        self.lr = lr
        self.training_iter = training_iter
        self.initial_noise = initial_noise
        self.standardize_y = standardize_y
        self.noise_mode = noise_mode
        self.min_observations_per_fidelity = min_observations_per_fidelity
        self.device = device
        self.dtype = dtype
        self.verbose = verbose
        self.log_every = log_every
        self.include_observation_noise = include_observation_noise
        self.random_state = random_state

    def fit(self, X, y, *, fidelity, sample_id=None):
        """Fit a two-level co-kriging GPR model."""
        _seed_torch(self.random_state)
        X_checked = check_array(X, ensure_2d=True, dtype="numeric", ensure_all_finite=True)
        y_checked = column_or_1d(
            check_array(
                y,
                ensure_2d=False,
                dtype="numeric",
                ensure_all_finite=True,
                input_name="y",
            ),
            warn=True,
        )
        check_consistent_length(X_checked, y_checked)
        self.n_features_in_ = X_checked.shape[1]

        observation_data = prepare_multifidelity_observations(
            X_checked,
            y_checked,
            fidelity,
            fidelity_order=self.fidelity_order,
            target_fidelity=self.target_fidelity,
            sample_id=sample_id,
            min_observations_per_fidelity=self.min_observations_per_fidelity,
        )
        self.result_ = fit_cokriging_gpr(
            observation_data,
            low_fidelity=self.low_fidelity,
            target_fidelity=self.target_fidelity,
            low_fidelity_kernel=self.low_fidelity_kernel,
            discrepancy_kernel=self.discrepancy_kernel,
            ard=self.ard,
            initial_rho=self.initial_rho,
            lr=self.lr,
            training_iter=self.training_iter,
            initial_noise=self.initial_noise,
            standardize_y=self.standardize_y,
            noise_mode=self.noise_mode,
            min_observations_per_fidelity=self.min_observations_per_fidelity,
            device=self.device,
            dtype=_resolve_torch_dtype(self.dtype),
            verbose=self.verbose,
            log_every=self.log_every,
        )
        self.model_ = self.result_.model
        self.likelihood_ = self.result_.likelihood
        self.observation_data_ = self.result_.observation_data
        self.fidelity_names_ = self.result_.observation_data.fidelity_names
        self.low_fidelity_ = self.result_.low_fidelity
        self.target_fidelity_ = self.result_.target_fidelity
        self.fidelity_observation_counts_ = self.result_.fidelity_observation_counts
        self.loss_history_ = list(self.result_.loss_history)
        self.rho_ = self.result_.rho
        self.noise_variance_ = self.result_.noise_variance
        return self

    def predict(
        self,
        X,
        *,
        target_fidelity: str | int | None = None,
        return_std: bool = False,
        include_observation_noise: bool | None = None,
    ):
        """Predict one modeled fidelity."""
        prediction = self.predict_distribution(
            X,
            target_fidelity=target_fidelity,
            return_std=return_std,
            include_observation_noise=include_observation_noise,
        )
        if return_std:
            return prediction.mean, prediction.std
        return prediction.mean

    def predict_distribution(
        self,
        X,
        *,
        target_fidelity: str | int | None = None,
        return_std: bool = True,
        confidence_level: float | None = None,
        include_observation_noise: bool | None = None,
    ) -> CoKrigingGPRPrediction:
        """Return predictive mean, uncertainty, and interval for one fidelity."""
        check_is_fitted(self, "result_")
        X_checked = check_array(X, ensure_2d=True, dtype="numeric", ensure_all_finite=True)
        if X_checked.shape[1] != self.n_features_in_:
            raise ValueError(
                f"X has {X_checked.shape[1]} features, but this estimator was fitted with "
                f"{self.n_features_in_} features"
            )
        if include_observation_noise is None:
            include_observation_noise = self.include_observation_noise
        return self.result_.predict(
            X_checked,
            target_fidelity=target_fidelity,
            return_std=return_std,
            confidence_level=confidence_level,
            include_observation_noise=include_observation_noise,
        )

    def score(self, X, y, sample_weight=None, *, target_fidelity: str | int | None = None) -> float:
        """Return R2 for predictions at one fidelity."""
        return r2_score(
            y,
            self.predict(X, target_fidelity=target_fidelity),
            sample_weight=sample_weight,
        )

fit(X, y, *, fidelity, sample_id=None)

Fit a two-level co-kriging GPR model.

Source code in matgpr/multifidelity_gpr.py

def fit(self, X, y, *, fidelity, sample_id=None):
    """Fit a two-level co-kriging GPR model."""
    _seed_torch(self.random_state)
    X_checked = check_array(X, ensure_2d=True, dtype="numeric", ensure_all_finite=True)
    y_checked = column_or_1d(
        check_array(
            y,
            ensure_2d=False,
            dtype="numeric",
            ensure_all_finite=True,
            input_name="y",
        ),
        warn=True,
    )
    check_consistent_length(X_checked, y_checked)
    self.n_features_in_ = X_checked.shape[1]

    observation_data = prepare_multifidelity_observations(
        X_checked,
        y_checked,
        fidelity,
        fidelity_order=self.fidelity_order,
        target_fidelity=self.target_fidelity,
        sample_id=sample_id,
        min_observations_per_fidelity=self.min_observations_per_fidelity,
    )
    self.result_ = fit_cokriging_gpr(
        observation_data,
        low_fidelity=self.low_fidelity,
        target_fidelity=self.target_fidelity,
        low_fidelity_kernel=self.low_fidelity_kernel,
        discrepancy_kernel=self.discrepancy_kernel,
        ard=self.ard,
        initial_rho=self.initial_rho,
        lr=self.lr,
        training_iter=self.training_iter,
        initial_noise=self.initial_noise,
        standardize_y=self.standardize_y,
        noise_mode=self.noise_mode,
        min_observations_per_fidelity=self.min_observations_per_fidelity,
        device=self.device,
        dtype=_resolve_torch_dtype(self.dtype),
        verbose=self.verbose,
        log_every=self.log_every,
    )
    self.model_ = self.result_.model
    self.likelihood_ = self.result_.likelihood
    self.observation_data_ = self.result_.observation_data
    self.fidelity_names_ = self.result_.observation_data.fidelity_names
    self.low_fidelity_ = self.result_.low_fidelity
    self.target_fidelity_ = self.result_.target_fidelity
    self.fidelity_observation_counts_ = self.result_.fidelity_observation_counts
    self.loss_history_ = list(self.result_.loss_history)
    self.rho_ = self.result_.rho
    self.noise_variance_ = self.result_.noise_variance
    return self

predict(X, *, target_fidelity=None, return_std=False, include_observation_noise=None)

Predict one modeled fidelity.

Source code in matgpr/multifidelity_gpr.py

def predict(
    self,
    X,
    *,
    target_fidelity: str | int | None = None,
    return_std: bool = False,
    include_observation_noise: bool | None = None,
):
    """Predict one modeled fidelity."""
    prediction = self.predict_distribution(
        X,
        target_fidelity=target_fidelity,
        return_std=return_std,
        include_observation_noise=include_observation_noise,
    )
    if return_std:
        return prediction.mean, prediction.std
    return prediction.mean

predict_distribution(X, *, target_fidelity=None, return_std=True, confidence_level=None, include_observation_noise=None)

Return predictive mean, uncertainty, and interval for one fidelity.

Source code in matgpr/multifidelity_gpr.py

def predict_distribution(
    self,
    X,
    *,
    target_fidelity: str | int | None = None,
    return_std: bool = True,
    confidence_level: float | None = None,
    include_observation_noise: bool | None = None,
) -> CoKrigingGPRPrediction:
    """Return predictive mean, uncertainty, and interval for one fidelity."""
    check_is_fitted(self, "result_")
    X_checked = check_array(X, ensure_2d=True, dtype="numeric", ensure_all_finite=True)
    if X_checked.shape[1] != self.n_features_in_:
        raise ValueError(
            f"X has {X_checked.shape[1]} features, but this estimator was fitted with "
            f"{self.n_features_in_} features"
        )
    if include_observation_noise is None:
        include_observation_noise = self.include_observation_noise
    return self.result_.predict(
        X_checked,
        target_fidelity=target_fidelity,
        return_std=return_std,
        confidence_level=confidence_level,
        include_observation_noise=include_observation_noise,
    )

score(X, y, sample_weight=None, *, target_fidelity=None)

Return R2 for predictions at one fidelity.

Source code in matgpr/multifidelity_gpr.py

def score(self, X, y, sample_weight=None, *, target_fidelity: str | int | None = None) -> float:
    """Return R2 for predictions at one fidelity."""
    return r2_score(
        y,
        self.predict(X, target_fidelity=target_fidelity),
        sample_weight=sample_weight,
    )

prepare_multifidelity_observations(X, y, fidelity, *, fidelity_order=None, target_fidelity=None, sample_id=None, noise_variance=None, feature_names=None, min_observations_per_fidelity=1)

Validate and encode row-wise observations from ordered fidelities.

fidelity_order should normally be supplied from low/cheap to high/target fidelity. If omitted, the order is inferred from first appearance in fidelity and the last level becomes the default target.

Source code in matgpr/multifidelity_gpr.py

def prepare_multifidelity_observations(
    X,
    y,
    fidelity,
    *,
    fidelity_order: Sequence[str] | None = None,
    target_fidelity: str | None = None,
    sample_id=None,
    noise_variance=None,
    feature_names: Sequence[str] | None = None,
    min_observations_per_fidelity: int = 1,
) -> MultiFidelityObservationData:
    """Validate and encode row-wise observations from ordered fidelities.

    ``fidelity_order`` should normally be supplied from low/cheap to
    high/target fidelity. If omitted, the order is inferred from first
    appearance in ``fidelity`` and the last level becomes the default target.
    """
    inferred_feature_names = _infer_feature_names(X)
    X_array = _to_2d_numpy(X, "X")
    y_array = _to_1d_numpy(y, "y")
    if X_array.shape[0] != y_array.shape[0]:
        raise ValueError("X and y must contain the same number of samples")

    labels = _to_1d_fidelity_labels(fidelity, "fidelity")
    if X_array.shape[0] != labels.shape[0]:
        raise ValueError("X and fidelity must contain the same number of samples")

    min_count = _validate_min_observations_per_fidelity(min_observations_per_fidelity)
    fidelity_names = _resolve_fidelity_names(labels, fidelity_order, min_count)
    target_fidelity_resolved = _resolve_target_fidelity(target_fidelity, fidelity_names)
    index_by_name = {name: index for index, name in enumerate(fidelity_names)}
    fidelity_index = np.asarray([index_by_name[label] for label in labels], dtype=int)

    sample_id_array = _resolve_sample_id(sample_id, X_array.shape[0])
    noise_variance_array = _resolve_noise_variance(noise_variance, X_array.shape[0])
    feature_names_resolved = _resolve_feature_names(
        feature_names if feature_names is not None else inferred_feature_names,
        X_array.shape[1],
    )

    return MultiFidelityObservationData(
        X=X_array.copy(),
        y=y_array.astype(float, copy=True),
        fidelity_index=fidelity_index,
        fidelity_names=fidelity_names,
        target_fidelity=target_fidelity_resolved,
        sample_id=sample_id_array,
        noise_variance=noise_variance_array,
        feature_names=feature_names_resolved,
    )

fit_cokriging_gpr(observation_data, *, low_fidelity=None, target_fidelity=None, low_fidelity_kernel='matern', discrepancy_kernel='matern', ard=True, initial_rho=None, lr=0.01, training_iter=1000, initial_noise=0.1, standardize_y=True, noise_mode='shared', min_observations_per_fidelity=2, device='cpu', dtype=torch.float64, verbose=False, log_every=100)

Fit a two-level autoregressive co-kriging GPR model.

This first implementation supports exactly two fidelity levels with one shared learned Gaussian noise term. Use :func:prepare_multifidelity_observations to create observation_data.

Source code in matgpr/multifidelity_gpr.py

def fit_cokriging_gpr(
    observation_data: MultiFidelityObservationData,
    *,
    low_fidelity: str | None = None,
    target_fidelity: str | None = None,
    low_fidelity_kernel: str = "matern",
    discrepancy_kernel: str = "matern",
    ard: bool = True,
    initial_rho: float | None = None,
    lr: float = 0.01,
    training_iter: int = 1000,
    initial_noise: float | None = 0.1,
    standardize_y: bool = True,
    noise_mode: str = "shared",
    min_observations_per_fidelity: int = 2,
    device: str = "cpu",
    dtype: torch.dtype | str = torch.float64,
    verbose: bool = False,
    log_every: int = 100,
) -> CoKrigingGPRResult:
    """Fit a two-level autoregressive co-kriging GPR model.

    This first implementation supports exactly two fidelity levels with one
    shared learned Gaussian noise term. Use
    :func:`prepare_multifidelity_observations` to create ``observation_data``.
    """
    if not isinstance(observation_data, MultiFidelityObservationData):
        raise ValueError("observation_data must be a MultiFidelityObservationData instance")
    dtype = _resolve_torch_dtype(dtype)
    _validate_training_options(lr=lr, training_iter=training_iter, log_every=log_every)
    noise_mode = _validate_cokriging_noise_mode(noise_mode)
    _validate_cokriging_known_noise(observation_data, noise_mode=noise_mode)

    low_name, target_name, low_index, target_index = _resolve_two_level_fidelities(
        observation_data,
        low_fidelity=low_fidelity,
        target_fidelity=target_fidelity,
        min_observations_per_fidelity=min_observations_per_fidelity,
    )
    initial_rho_value = _resolve_initial_cokriging_rho(
        observation_data,
        low_fidelity_index=low_index,
        target_fidelity_index=target_index,
        initial_rho=initial_rho,
    )

    train_x = _to_tensor(observation_data.X, device=device, dtype=dtype)
    train_fidelity_indices = torch.as_tensor(
        observation_data.fidelity_index,
        dtype=torch.long,
        device=device,
    )
    train_y = _to_tensor(observation_data.y, device=device, dtype=dtype).reshape(-1)

    if standardize_y:
        target_mean = train_y.mean()
        target_std = train_y.std(unbiased=False)
        if target_std.item() <= 0:
            raise ValueError("y has zero standard deviation")
        train_y_model = (train_y - target_mean) / target_std
    else:
        target_mean = torch.tensor(0.0, dtype=dtype, device=device)
        target_std = torch.tensor(1.0, dtype=dtype, device=device)
        train_y_model = train_y

    likelihood = gpytorch.likelihoods.GaussianLikelihood().to(device=device, dtype=dtype)
    if initial_noise is not None:
        if initial_noise <= 0:
            raise ValueError("initial_noise must be positive")
        likelihood.initialize(noise=float(initial_noise))

    model = ExactTwoLevelCoKrigingGPRModel(
        train_x,
        train_fidelity_indices,
        train_y_model,
        likelihood,
        low_fidelity_index=low_index,
        target_fidelity_index=target_index,
        low_fidelity_kernel=low_fidelity_kernel,
        discrepancy_kernel=discrepancy_kernel,
        ard_num_dims=train_x.shape[1] if ard else None,
        initial_rho=initial_rho_value,
    ).to(device=device, dtype=dtype)
    model.target_mean = target_mean.detach()
    model.target_std = target_std.detach()
    model.standardize_y = standardize_y
    model.fidelity_names = observation_data.fidelity_names
    model.low_fidelity = low_name
    model.target_fidelity = target_name

    model.train()
    likelihood.train()

    optimizer = torch.optim.Adam(model.parameters(), lr=lr)
    marginal_log_likelihood = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)
    loss_history: list[float] = []

    for iteration in range(training_iter):
        optimizer.zero_grad()
        output = model(train_x, train_fidelity_indices)
        loss = -marginal_log_likelihood(output, train_y_model)
        loss.backward()
        optimizer.step()
        loss_history.append(float(loss.detach().cpu().item()))

        if verbose and _should_log_iteration(iteration, training_iter, log_every):
            _print_cokriging_training_status(iteration, training_iter, loss, likelihood, model)

    model.training_loss_history = loss_history
    rho = float(model.rho.detach().cpu().item())
    standardized_noise = float(likelihood.noise.detach().cpu().reshape(-1)[0])
    noise_variance = standardized_noise * float(target_std.detach().cpu().item()) ** 2

    return CoKrigingGPRResult(
        model=model,
        likelihood=likelihood,
        observation_data=observation_data,
        loss_history=loss_history,
        target_mean=float(target_mean.detach().cpu().item()),
        target_std=float(target_std.detach().cpu().item()),
        standardize_y=standardize_y,
        low_fidelity=low_name,
        target_fidelity=target_name,
        low_fidelity_index=low_index,
        target_fidelity_index=target_index,
        fidelity_observation_counts=observation_data.fidelity_observation_counts,
        low_fidelity_kernel=low_fidelity_kernel,
        discrepancy_kernel=discrepancy_kernel,
        ard=ard,
        noise_mode=noise_mode,
        noise_variance=noise_variance,
        rho=rho,
        device=device,
        dtype=dtype,
    )

fit_delta_multifidelity_gpr(X_high, y_high, *, low_fidelity_high=None, X_low=None, y_low=None, fit_intercept=True, correction_kernel='matern', low_fidelity_kernel='matern', ard=True, lr=0.01, training_iter=1000, low_fidelity_training_iter=None, correction_initial_noise=0.1, low_fidelity_initial_noise=0.1, standardize_y=True, include_low_fidelity_uncertainty=True, device='cpu', dtype=torch.float64, verbose=False, log_every=100)

Fit a delta multi-fidelity GPR model.

The model uses the autoregressive correction

y_high(x) = rho * y_low(x) + intercept + delta(x),

where delta(x) is modeled by a standard exact GPR. Provide low_fidelity_high when low-fidelity values are already available at the high-fidelity training points. Alternatively, provide X_low and y_low to fit a low-fidelity GPR surrogate first.

Source code in matgpr/multifidelity_gpr.py

def fit_delta_multifidelity_gpr(
    X_high,
    y_high,
    *,
    low_fidelity_high=None,
    X_low=None,
    y_low=None,
    fit_intercept: bool = True,
    correction_kernel: str = "matern",
    low_fidelity_kernel: str = "matern",
    ard: bool = True,
    lr: float = 0.01,
    training_iter: int = 1000,
    low_fidelity_training_iter: int | None = None,
    correction_initial_noise: float | None = 0.1,
    low_fidelity_initial_noise: float | None = 0.1,
    standardize_y: bool = True,
    include_low_fidelity_uncertainty: bool = True,
    device: str = "cpu",
    dtype: torch.dtype | str = torch.float64,
    verbose: bool = False,
    log_every: int = 100,
) -> DeltaMultiFidelityGPRResult:
    """Fit a delta multi-fidelity GPR model.

    The model uses the autoregressive correction

    ``y_high(x) = rho * y_low(x) + intercept + delta(x)``,

    where ``delta(x)`` is modeled by a standard exact GPR. Provide
    ``low_fidelity_high`` when low-fidelity values are already available at the
    high-fidelity training points. Alternatively, provide ``X_low`` and
    ``y_low`` to fit a low-fidelity GPR surrogate first.
    """
    dtype = _resolve_torch_dtype(dtype)
    X_high_array = _to_2d_numpy(X_high, "X_high")
    y_high_array = _to_1d_numpy(y_high, "y_high")
    check_consistent_length(X_high_array, y_high_array)
    _validate_feature_width_match(X_high_array, X_low, "X_high", "X_low")

    low_model: GPyTorchGPRResult | None = None
    if low_fidelity_high is None:
        if X_low is None or y_low is None:
            raise ValueError(
                "Provide low_fidelity_high, or provide both X_low and y_low to fit "
                "a low-fidelity surrogate"
            )
        low_model = fit_gpytorch_gpr(
            X_low,
            y_low,
            kernel=low_fidelity_kernel,
            ard=ard,
            lr=lr,
            training_iter=low_fidelity_training_iter or training_iter,
            initial_noise=low_fidelity_initial_noise,
            standardize_y=standardize_y,
            device=device,
            dtype=dtype,
            verbose=verbose,
            log_every=log_every,
        )
        low_prediction = low_model.predict(
            X_high_array,
            return_std=True,
            include_observation_noise=False,
        )
        low_fidelity_high_array = low_prediction.mean
    else:
        if X_low is not None or y_low is not None:
            raise ValueError("Use either low_fidelity_high or X_low/y_low, not both")
        low_fidelity_high_array = _to_1d_numpy(low_fidelity_high, "low_fidelity_high")

    check_consistent_length(X_high_array, low_fidelity_high_array)

    rho, intercept = _fit_linear_fidelity_map(
        low_fidelity_high_array,
        y_high_array,
        fit_intercept=fit_intercept,
    )
    correction_target = y_high_array - (rho * low_fidelity_high_array + intercept)
    standardize_correction = standardize_y and np.std(correction_target) > 0.0
    correction_model = fit_gpytorch_gpr(
        X_high_array,
        correction_target,
        kernel=correction_kernel,
        ard=ard,
        lr=lr,
        training_iter=training_iter,
        initial_noise=correction_initial_noise,
        standardize_y=standardize_correction,
        device=device,
        dtype=dtype,
        verbose=verbose,
        log_every=log_every,
    )

    return DeltaMultiFidelityGPRResult(
        correction_model=correction_model,
        low_fidelity_model=low_model,
        rho=float(rho),
        intercept=float(intercept),
        high_fidelity_target_mean=float(y_high_array.mean()),
        high_fidelity_target_std=float(y_high_array.std(ddof=0)),
        correction_target=correction_target.copy(),
        low_fidelity_at_high=low_fidelity_high_array.copy(),
        fit_intercept=fit_intercept,
        include_low_fidelity_uncertainty=include_low_fidelity_uncertainty,
        device=device,
        dtype=dtype,
    )