# Path: sdcit/hsic.py
# def HSIC(K: np.ndarray, L: np.ndarray, p_val_method='bootstrap', num_boot=1000) -> float:
#     if p_val_method == 'bootstrap':
#         return HSIC_boot(K, L, num_boot)
#     elif p_val_method == 'gamma':
#         return HSIC_gamma_approx(K, L)
#     else:
#         raise ValueError('unknown p value computation method: {}'.format(p_val_method))
# 
# Path: sdcit/utils.py
# def rbf_kernel_median(data: np.ndarray, *args, without_two=False):
#     """A list of RBF kernel matrices for data sets in arguments based on median heuristic"""
#     if args is None:
#         args = []
# 
#     outs = []
#     for x in [data, *args]:
#         D_squared = euclidean_distances(x, squared=True)
#         # masking upper triangle and the diagonal.
#         mask = np.triu(np.ones(D_squared.shape), 0)
#         median_squared_distance = ma.median(ma.array(D_squared, mask=mask))
#         if without_two:
#             kx = exp(-D_squared / median_squared_distance)
#         else:
#             kx = exp(-0.5 * D_squared / median_squared_distance)
#         outs.append(kx)
# 
#     if len(outs) == 1:
#         return outs[0]
#     else:
#         return outs
# 
# def residual_kernel(K_Y: np.ndarray, K_X: np.ndarray, use_expectation=True, with_gp=True, sigma_squared=1e-3, return_learned_K_X=False):
#     """Kernel matrix of residual of Y given X based on their kernel matrices, Y=f(X)"""
#     import gpflow
#     from gpflow.kernels import White, Linear
#     from gpflow.models import GPR
# 
#     K_Y, K_X = centering(K_Y), centering(K_X)
#     T = len(K_Y)
# 
#     if with_gp:
#         eig_Ky, eiy = truncated_eigen(*eigdec(K_Y, min(100, T // 4)))
#         eig_Kx, eix = truncated_eigen(*eigdec(K_X, min(100, T // 4)))
# 
#         X = eix @ diag(sqrt(eig_Kx))  # X @ X.T is close to K_X
#         Y = eiy @ diag(sqrt(eig_Ky))
#         n_feats = X.shape[1]
# 
#         linear = Linear(n_feats, ARD=True)
#         white = White(n_feats)
#         gp_model = GPR(X, Y, linear + white)
#         gpflow.train.ScipyOptimizer().minimize(gp_model)
# 
#         K_X = linear.compute_K_symm(X)
#         sigma_squared = white.variance.value
# 
#     P = pdinv(np.eye(T) + K_X / sigma_squared)  # == I-K @ inv(K+Sigma) in Zhang et al. 2011
#     if use_expectation:  # Flaxman et al. 2016 Gaussian Processes for Independence Tests with Non-iid Data in Causal Inference.
#         RK = (K_X + P @ K_Y) @ P
#     else:  # Zhang et al. 2011. Kernel-based Conditional Independence Test and Application in Causal Discovery.
#         RK = P @ K_Y @ P
# 
#     if return_learned_K_X:
#         return RK, K_X
#     else:
#         return RK
# 
# def residualize(Y, X=None, gp_kernel=None):
#     """Residual of Y given X. Y_i - E[Y_i|X_i]"""
#     import gpflow
#     from gpflow.models import GPR
# 
#     if X is None:
#         return Y - np.mean(Y)  # nothing is residualized!
# 
#     if gp_kernel is None:
#         gp_kernel = default_gp_kernel(X)
# 
#     m = GPR(X, Y, gp_kernel)
#     gpflow.train.ScipyOptimizer().minimize(m)
# 
#     Yhat, _ = m.predict_y(X)
#     return Y - Yhat
# 
# def columnwise_normalizes(*Xs) -> typing.List[Union[None, np.ndarray]]:
#     """normalize per column for multiple data"""
#     return [columnwise_normalize(X) for X in Xs]


import numpy as np
from sdcit.hsic import HSIC
from sdcit.utils import rbf_kernel_median, residual_kernel, residualize, columnwise_normalizes

    if Kz is None:
        return HSIC(Kx, Ky, **hsic_kws)

    RX_Z = residual_kernel(Kx, Kz, use_expectation=use_expectation, with_gp=with_gp, sigma_squared=sigma_squared)
    RY_Z = residual_kernel(Ky, Kz, use_expectation=use_expectation, with_gp=with_gp, sigma_squared=sigma_squared)

    return HSIC(RX_Z, RY_Z, **hsic_kws)


def FCIT(X, Y, Z=None, kern=rbf_kernel_median, normalize=False, seed=None, hsic_kws=None):
    """Flaxman et al. Residualization-based CI Test, X_||_Y | Z

    References
    ----------
    Flaxman, S. R., Neill, D. B., & Smola, A. J. (2016). Gaussian Processes for Independence Tests with Non-iid Data in Causal Inference.
    ACM Transactions on Intelligent Systems and Technology, 7(2), 1–23.
    """
    if seed is not None:
        np.random.seed(seed)

    if hsic_kws is None:
        hsic_kws = {}

    if normalize:
        X, Y, Z = columnwise_normalizes(X, Y, Z)

    if Z is None:
        return HSIC(kern(X), kern(Y), **hsic_kws)


# Path: sdcit/hsic.py
# def HSIC(K: np.ndarray, L: np.ndarray, p_val_method='bootstrap', num_boot=1000) -> float:
#     if p_val_method == 'bootstrap':
#         return HSIC_boot(K, L, num_boot)
#     elif p_val_method == 'gamma':
#         return HSIC_gamma_approx(K, L)
#     else:
#         raise ValueError('unknown p value computation method: {}'.format(p_val_method))
# 
# Path: sdcit/utils.py
# def rbf_kernel_median(data: np.ndarray, *args, without_two=False):
#     """A list of RBF kernel matrices for data sets in arguments based on median heuristic"""
#     if args is None:
#         args = []
# 
#     outs = []
#     for x in [data, *args]:
#         D_squared = euclidean_distances(x, squared=True)
#         # masking upper triangle and the diagonal.
#         mask = np.triu(np.ones(D_squared.shape), 0)
#         median_squared_distance = ma.median(ma.array(D_squared, mask=mask))
#         if without_two:
#             kx = exp(-D_squared / median_squared_distance)
#         else:
#             kx = exp(-0.5 * D_squared / median_squared_distance)
#         outs.append(kx)
# 
#     if len(outs) == 1:
#         return outs[0]
#     else:
#         return outs
# 
# def residual_kernel(K_Y: np.ndarray, K_X: np.ndarray, use_expectation=True, with_gp=True, sigma_squared=1e-3, return_learned_K_X=False):
#     """Kernel matrix of residual of Y given X based on their kernel matrices, Y=f(X)"""
#     import gpflow
#     from gpflow.kernels import White, Linear
#     from gpflow.models import GPR
# 
#     K_Y, K_X = centering(K_Y), centering(K_X)
#     T = len(K_Y)
# 
#     if with_gp:
#         eig_Ky, eiy = truncated_eigen(*eigdec(K_Y, min(100, T // 4)))
#         eig_Kx, eix = truncated_eigen(*eigdec(K_X, min(100, T // 4)))
# 
#         X = eix @ diag(sqrt(eig_Kx))  # X @ X.T is close to K_X
#         Y = eiy @ diag(sqrt(eig_Ky))
#         n_feats = X.shape[1]
# 
#         linear = Linear(n_feats, ARD=True)
#         white = White(n_feats)
#         gp_model = GPR(X, Y, linear + white)
#         gpflow.train.ScipyOptimizer().minimize(gp_model)
# 
#         K_X = linear.compute_K_symm(X)
#         sigma_squared = white.variance.value
# 
#     P = pdinv(np.eye(T) + K_X / sigma_squared)  # == I-K @ inv(K+Sigma) in Zhang et al. 2011
#     if use_expectation:  # Flaxman et al. 2016 Gaussian Processes for Independence Tests with Non-iid Data in Causal Inference.
#         RK = (K_X + P @ K_Y) @ P
#     else:  # Zhang et al. 2011. Kernel-based Conditional Independence Test and Application in Causal Discovery.
#         RK = P @ K_Y @ P
# 
#     if return_learned_K_X:
#         return RK, K_X
#     else:
#         return RK
# 
# def residualize(Y, X=None, gp_kernel=None):
#     """Residual of Y given X. Y_i - E[Y_i|X_i]"""
#     import gpflow
#     from gpflow.models import GPR
# 
#     if X is None:
#         return Y - np.mean(Y)  # nothing is residualized!
# 
#     if gp_kernel is None:
#         gp_kernel = default_gp_kernel(X)
# 
#     m = GPR(X, Y, gp_kernel)
#     gpflow.train.ScipyOptimizer().minimize(m)
# 
#     Yhat, _ = m.predict_y(X)
#     return Y - Yhat
# 
# def columnwise_normalizes(*Xs) -> typing.List[Union[None, np.ndarray]]:
#     """normalize per column for multiple data"""
#     return [columnwise_normalize(X) for X in Xs]


# Path: sdcit/flaxman.py
import numpy as np
from sdcit.hsic import HSIC
from sdcit.utils import rbf_kernel_median, residual_kernel, residualize, columnwise_normalizes


    if Kz is None:
        return HSIC(Kx, Ky, **hsic_kws)

    RX_Z = residual_kernel(Kx, Kz, use_expectation=use_expectation, with_gp=with_gp, sigma_squared=sigma_squared)
    RY_Z = residual_kernel(Ky, Kz, use_expectation=use_expectation, with_gp=with_gp, sigma_squared=sigma_squared)

    return HSIC(RX_Z, RY_Z, **hsic_kws)


def FCIT(X, Y, Z=None, kern=rbf_kernel_median, normalize=False, seed=None, hsic_kws=None):
    """Flaxman et al. Residualization-based CI Test, X_||_Y | Z

    References
    ----------
    Flaxman, S. R., Neill, D. B., & Smola, A. J. (2016). Gaussian Processes for Independence Tests with Non-iid Data in Causal Inference.
    ACM Transactions on Intelligent Systems and Technology, 7(2), 1–23.
    """
    if seed is not None:
        np.random.seed(seed)

    if hsic_kws is None:
        hsic_kws = {}

    if normalize:
        X, Y, Z = columnwise_normalizes(X, Y, Z)

    if Z is None:
        return HSIC(kern(X), kern(Y), **hsic_kws)


    e_YZ = residualize(Y, Z)