# coding: utf-8
# pylint: disable = invalid-name, C0111
import gpboost as gpb
import numpy as np
from sklearn.metrics import mean_squared_error
import matplotlib.pyplot as plt
plt.style.use('ggplot')

#--------------------Cross validation for tree-boosting without GP or random effects----------------
print('Simulating data...')
# Simulate and create your dataset
def f1d(x):
    """Non-linear function for simulation"""
    return (1.7 * (1 / (1 + np.exp(-(x - 0.5) * 20)) + 0.75 * x))
x = np.linspace(0, 1, 200, endpoint=True)
plt.plot(x, f1d(x), linewidth=2, color="r")
plt.title("Mean function")
plt.show()
def sim_data(n):
    """Function that simulates data. Two covariates of which only one has an effect"""
    X = np.random.rand(n, 2)
    # mean function plus noise
    y = f1d(X[:, 0]) + np.random.normal(scale=0.1, size=n)
    return ([X, y])

# Simulate data
n = 1000
data = sim_data(2 * n)
# create dataset for gpb.train
data_train = gpb.Dataset(data[0][0:n, :], data[1][0:n])

# specify your configurations as a dict
params = {
    'objective': 'regression_l2',
    'metric': {'l2', 'l1'},
    'learning_rate': 0.1,
    'max_depth': 6,
    'min_data_in_leaf': 5,
    'verbose': 0
}

print('Starting cross-validation...')
# do cross-validation
cvbst = gpb.cv(params=params, train_set=data_train,
               num_boost_round=100, early_stopping_rounds=5,
               nfold=2, verbose_eval=True, show_stdv=False, seed=1)
print("Best number of iterations: " + str(np.argmin(cvbst['l2-mean'])))


# --------------------Combine tree-boosting and grouped random effects model----------------
print('Simulating data...')
# Simulate data
def f1d(x):
    """Non-linear function for simulation"""
    return (1.7 * (1 / (1 + np.exp(-(x - 0.5) * 20)) + 0.75 * x))
x = np.linspace(0, 1, 200, endpoint=True)
plt.figure("Mean function")
plt.plot(x, f1d(x), linewidth=2, color="r")
plt.title("Mean function")
plt.show()
n = 1000  # number of samples
np.random.seed(1)
X = np.random.rand(n, 2)
F = f1d(X[:, 0])
# Simulate grouped random effects
m = 25  # number of categories / levels for grouping variable
group = np.arange(n)  # grouping variable
for i in range(m):
    group[int(i * n / m):int((i + 1) * n / m)] = i
# incidence matrix relating grouped random effects to samples
Z1 = np.zeros((n, m))
for i in range(m):
    Z1[np.where(group == i), i] = 1
sigma2_1 = 1 ** 2  # random effect variance
sigma2 = 0.1 ** 2  # error variance
b1 = np.sqrt(sigma2_1) * np.random.normal(size=m)  # simulate random effects
eps = Z1.dot(b1)
xi = np.sqrt(sigma2) * np.random.normal(size=n)  # simulate error term
y = F + eps + xi  # observed data

# define GPModel
gp_model = gpb.GPModel(group_data=group)
gp_model.set_optim_params(params={"optimizer_cov": "fisher_scoring"})
# create dataset for gpb.train
data_train = gpb.Dataset(X, y)
# specify your configurations as a dict
params = {
    'objective': 'regression_l2',
    'learning_rate': 0.05,
    'max_depth': 6,
    'min_data_in_leaf': 5,
    'verbose': 0
}

print('Starting cross-validation...')
# do cross-validation
cvbst = gpb.cv(params=params, train_set=data_train,
               gp_model=gp_model, use_gp_model_for_validation=False,
               num_boost_round=100, early_stopping_rounds=5,
               nfold=2, verbose_eval=True, show_stdv=False, seed=1)
print("Best number of iterations: " + str(np.argmin(cvbst['l2-mean'])))

# Include random effect predictions for validation (observe the lower test error)
gp_model = gpb.GPModel(group_data=group)
print("Running cross validation for GPBoost model and use_gp_model_for_validation = TRUE")
cvbst = gpb.cv(params=params, train_set=data_train,
               gp_model=gp_model, use_gp_model_for_validation=True,
               num_boost_round=100, early_stopping_rounds=5,
               nfold=2, verbose_eval=True, show_stdv=Falsem, seed=1)
print("Best number of iterations: " + str(np.argmin(cvbst['l2-mean'])))
cvbst.best_iteration

[((7, 0, 7, 23), 'matplotlib.pyplot.style.use', 'plt.style.use', ({(7, 14, 7, 22): '"""ggplot"""'}, {}), "('ggplot')", True, 'import matplotlib.pyplot as plt\n'), ((15, 4, 15, 41), 'numpy.linspace', 'np.linspace', (), '', True, 'import numpy as np\n'), ((17, 0, 17, 26), 'matplotlib.pyplot.title', 'plt.title', ({(17, 10, 17, 25): '"""Mean function"""'}, {}), "('Mean function')", True, 'import matplotlib.pyplot as plt\n'), ((18, 0, 18, 10), 'matplotlib.pyplot.show', 'plt.show', ({}, {}), '()', True, 'import matplotlib.pyplot as plt\n'), ((30, 13, 30, 55), 'gpboost.Dataset', 'gpb.Dataset', ({(30, 25, 30, 40): 'data[0][0:n, :]', (30, 42, 30, 54): 'data[1][0:n]'}, {}), '(data[0][0:n, :], data[1][0:n])', True, 'import gpboost as gpb\n'), ((44, 8, 46, 67), 'gpboost.cv', 'gpb.cv', (), '', True, 'import gpboost as gpb\n'), ((56, 4, 56, 41), 'numpy.linspace', 'np.linspace', (), '', True, 'import numpy as np\n'), ((57, 0, 57, 27), 'matplotlib.pyplot.figure', 'plt.figure', ({(57, 11, 57, 26): '"""Mean function"""'}, {}), "('Mean function')", True, 'import matplotlib.pyplot as plt\n'), ((59, 0, 59, 26), 'matplotlib.pyplot.title', 'plt.title', ({(59, 10, 59, 25): '"""Mean function"""'}, {}), "('Mean function')", True, 'import matplotlib.pyplot as plt\n'), ((60, 0, 60, 10), 'matplotlib.pyplot.show', 'plt.show', ({}, {}), '()', True, 'import matplotlib.pyplot as plt\n'), ((62, 0, 62, 17), 'numpy.random.seed', 'np.random.seed', ({(62, 15, 62, 16): '(1)'}, {}), '(1)', True, 'import numpy as np\n'), ((63, 4, 63, 24), 'numpy.random.rand', 'np.random.rand', ({(63, 19, 63, 20): 'n', (63, 22, 63, 23): '2'}, {}), '(n, 2)', True, 'import numpy as np\n'), ((67, 8, 67, 20), 'numpy.arange', 'np.arange', ({(67, 18, 67, 19): 'n'}, {}), '(n)', True, 'import numpy as np\n'), ((71, 5, 71, 21), 'numpy.zeros', 'np.zeros', ({(71, 14, 71, 20): '(n, m)'}, {}), '((n, m))', True, 'import numpy as np\n'), ((82, 11, 82, 40), 'gpboost.GPModel', 'gpb.GPModel', (), '', True, 'import gpboost as gpb\n'), ((85, 13, 85, 30), 'gpboost.Dataset', 'gpb.Dataset', ({(85, 25, 85, 26): 'X', (85, 28, 85, 29): 'y'}, {}), '(X, y)', True, 'import gpboost as gpb\n'), ((97, 8, 100, 67), 'gpboost.cv', 'gpb.cv', (), '', True, 'import gpboost as gpb\n'), ((104, 11, 104, 40), 'gpboost.GPModel', 'gpb.GPModel', (), '', True, 'import gpboost as gpb\n'), ((106, 8, 109, 68), 'gpboost.cv', 'gpb.cv', (), '', True, 'import gpboost as gpb\n'), ((21, 8, 21, 28), 'numpy.random.rand', 'np.random.rand', ({(21, 23, 21, 24): 'n', (21, 26, 21, 27): '2'}, {}), '(n, 2)', True, 'import numpy as np\n'), ((76, 5, 76, 22), 'numpy.sqrt', 'np.sqrt', ({(76, 13, 76, 21): 'sigma2_1'}, {}), '(sigma2_1)', True, 'import numpy as np\n'), ((76, 25, 76, 49), 'numpy.random.normal', 'np.random.normal', (), '', True, 'import numpy as np\n'), ((78, 5, 78, 20), 'numpy.sqrt', 'np.sqrt', ({(78, 13, 78, 19): 'sigma2'}, {}), '(sigma2)', True, 'import numpy as np\n'), ((78, 23, 78, 47), 'numpy.random.normal', 'np.random.normal', (), '', True, 'import numpy as np\n'), ((23, 23, 23, 58), 'numpy.random.normal', 'np.random.normal', (), '', True, 'import numpy as np\n'), ((47, 42, 47, 69), 'numpy.argmin', 'np.argmin', ({(47, 52, 47, 68): "cvbst['l2-mean']"}, {}), "(cvbst['l2-mean'])", True, 'import numpy as np\n'), ((101, 42, 101, 69), 'numpy.argmin', 'np.argmin', ({(101, 52, 101, 68): "cvbst['l2-mean']"}, {}), "(cvbst['l2-mean'])", True, 'import numpy as np\n'), ((110, 42, 110, 69), 'numpy.argmin', 'np.argmin', ({(110, 52, 110, 68): "cvbst['l2-mean']"}, {}), "(cvbst['l2-mean'])", True, 'import numpy as np\n'), ((73, 7, 73, 27), 'numpy.where', 'np.where', ({(73, 16, 73, 26): '(group == i)'}, {}), '(group == i)', True, 'import numpy as np\n'), ((14, 28, 14, 51), 'numpy.exp', 'np.exp', ({(14, 35, 14, 50): '(-(x - 0.5) * 20)'}, {}), '(-(x - 0.5) * 20)', True, 'import numpy as np\n'), ((55, 28, 55, 51), 'numpy.exp', 'np.exp', ({(55, 35, 55, 50): '(-(x - 0.5) * 20)'}, {}), '(-(x - 0.5) * 20)', True, 'import numpy as np\n')]