import os
import sys
import re
import random

import tklib.tkimport as imp

np    = imp.import_lib("numpy", stop_by_error = False)
scipy = imp.import_lib("scipy", stop_by_error = False)
mpl   = imp.import_lib("matplotlib", stop_by_error = False)
pd    = imp.import_lib("pandas", stop_by_error = False)
imp.messages(stop_by_error = True)
#sns = imp.import_lib("seaborn", stop_by_error = False)
imp.messages(stop_by_error = True)

from numpy import sin, cos, tan, pi, exp, log, sqrt
from scipy.stats import norm
from matplotlib import pyplot as plt
import matplotlib.widgets as wg

from sklearn.model_selection import train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.utils.validation import column_or_1d
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score

from sklearn.linear_model import LassoCV, RidgeCV, ElasticNetCV
from sklearn.model_selection import GridSearchCV

from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import Pipeline

from sklearn.linear_model import LinearRegression, Ridge, Lasso, ElasticNet, SGDRegressor
from sklearn.linear_model import PassiveAggressiveRegressor, ARDRegression, RidgeCV
from sklearn.neural_network import MLPRegressor

from sklearn.gaussian_process.kernels import WhiteKernel, RBF, DotProduct, ExpSineSquared
from sklearn.gaussian_process.kernels import Matern, ConstantKernel as C
from sklearn.gaussian_process import GaussianProcessRegressor


from sklearn.model_selection import StratifiedKFold
from sklearn.model_selection import cross_val_score
from sklearn.linear_model import TheilSenRegressor, RANSACRegressor, HuberRegressor
from sklearn.svm import SVR, LinearSVR
from sklearn.neighbors import KNeighborsRegressor
from sklearn.tree import DecisionTreeRegressor
from sklearn.ensemble import RandomForestRegressor, AdaBoostRegressor, ExtraTreesRegressor, HistGradientBoostingRegressor
from sklearn.ensemble import BaggingRegressor, GradientBoostingRegressor, VotingRegressor, StackingRegressor
from sklearn.cross_decomposition import PLSRegression


from tklib.tkobject import tkObject
from tklib.tkutils import pint, pfloat, pintfloat
from tklib.tkinifile import tkIniFile
from tklib.tkparams import tkParams
from tklib.tkgraphic.tkplotevent import tkPlotEvent, RangeSelector


def extract_labels(labels, regexp_list):
    label_ret = []
    for l in labels:
        for rg in regexp_list:
            if re.search(rg, l):
                if l not in label_ret:
                    label_ret.append(l)

    return label_ret

def extract_unique(_list):
    list_ret = []
    for l in _list:
        if l not in list_ret:
            list_ret.append(l)
    return list_ret

class tkMLData(tkParams):
    def __init__(self, **kwargs):
        super().__init__(**kwargs)

        self._df = None

    @property
    def ncolumns(self):
        return len(self._df.columns)

    @property
    def ndata(self):
        return len(self._df)

    @property
    def labels(self):
        return self._df.columns

    @property
    def columns(self):
        return self._df.columns

    @property
    def df(self):
        return self._df

    @df.setter
    def df(self, dataframe):
        self._df = dataframe
#        self.labels   = dataframe.columns
#        self.columns  = dataframe.columns
#        self.ncolumns = len(dataframe.columns)
#        self.ndata    = len(dataframe)
    
    @df.deleter
    def df(self):
        del self._df
        del self.labels
        del self.columns
        del self.ncolumns
        del self.ndata
        
class tkMLR(tkMLData):
    def __init__(self, **kwargs):
        super().__init__(**kwargs)
        
        self._df = None

    def get_objective_labels(self, labels = [], keys = [], id = [], regexp_list = [], apply_default = True):
        labels = list(self.df.columns)

        objective_labels = []
        objective_labels += [label for label, _id in zip(keys, id) if _id == 0]
        if apply_default:
            objective_labels += extract_labels(labels, regexp_list = [r"^o:", r"^t:", r"^min:", r"^max:"])
        objective_labels += extract_labels(labels, regexp_list = regexp_list)
        objective_labels = extract_unique(objective_labels)

        return objective_labels

    def get_exclude_labels(self, labels = [], keys = [], id = [], regexp_list = [], apply_default = True):
        labels = list(self.df.columns)

        exclude_labels = []
        exclude_labels += [label for label, _id in zip(keys, id) if _id == 0]
        if apply_default:
            exclude_labels += extract_labels(labels, regexp_list = [r"^\-", r"^'\-"])
        exclude_labels += extract_labels(labels, regexp_list = regexp_list)
        exclude_labels = extract_unique(exclude_labels)
        
        return exclude_labels

    def get_descriptor_labels(self, labels = [], keys = [], id = [], regexp_list = [], remove_keys = []):
        labels = list(self.df.columns)

        descriptor_labels = []
        for l in labels:
            if l not in remove_keys:
                descriptor_labels.append(l)
        descriptor_labels += [label for label, _id in zip(keys, id) if _id == 0]
        descriptor_labels += extract_labels(labels, regexp_list = regexp_list)
        descriptor_labels = extract_unique(descriptor_labels)

        return descriptor_labels

    def read_data(self, path, header = 0, index_col = None):
        _, ext = os.path.splitext(path)
        if ext.lower() == '.csv':
            _df = pd.read_csv(path, header = header, index_col = index_col)
        else:
            _df = pd.read_excel(path, header = header, index_col = index_col)

        self.df = _df
#        self.ncolumns = len(_df.columns)
#        self.ndata    = len(_df)
#        self.labels = self.df.columns

        return self.df

    def extract_data(self, df = None, delete_na = True, limit_min_list = [], limit_max_list = []):
        if df is None:
            df = self.df

        _df = self.df.select_dtypes(include=[np.number])
        
        for cond in limit_max_list:
            key, max_val = cond
            _df = _df[_df[key] <= max_val]

        for cond in limit_min_list:
            key, min_val = cond
            _df = _df[_df[key] >= min_val]

        return _df

    def convert_df(self, df, scale_dict = {}):
        self.scale_dict = scale_dict
        _df = pd.DataFrame()
        keys = list(scale_dict.keys())
        for i, label in enumerate(df.columns):
            if label not in keys:
#                print(f"Error: Can not find descriptor [{label}] in fit.varname, ", fit.varname)
                _df[label] = df[label].copy()
            else:
                scale = scale_dict[label]
#                print(f"  {label}: {scale}")
                if scale == 'log':
                    _df[label] = log(df[label])
                else:
                    _df[label] = df[label].copy()

        return _df

    def standardize(self, df):
        if df is None: return None

        self.scaler = StandardScaler()
        if isinstance(df, np.ndarray):
            ret = self.scaler.fit_transform(df)
            self.mean_standardize   = self.scaler.mean_ 
            self.stddev_standardize = self.scaler.scale_
            return ret
        else:
            ret = self.scaler.fit_transform(df)
            self.mean_standardize   = self.scaler.mean_ 
            self.stddev_standardize = self.scaler.scale_
            return pd.DataFrame(ret, columns = df.columns)

    def inverse_standardize(self, df_std):
        if df_std is None: return None

        _df = self.scaler.inverse_transform(df_std)
        return _df

    def mean_std(self, df):
        return np.mean(df, 0), np.std(df, 0)

    def split_data(self, descriptor_df, objective_df, test_size):
        d_train, d_test, o_train, o_test =train_test_split(descriptor_df, objective_df, test_size = test_size)

        return d_train, d_test, o_train, o_test

    def extract_order(self, method, pos):
        if len(method) <= pos:
            return None

        order = pint(method[pos:], None)
        if order is None:
            print(f"\nError in tkmlr.extract_order(): Invalid order=[{method[pos:]}] given by [{method}]")
#            return None
            exit()

        return order

    def get_model(self, descriptor_df, objective_df, method = 'mlp', 
            use_cv = False, cv = 5,
            n_restarts_optimizer = 50,
            alpha = 0.1, l1_ratio = 0.5, 
            hidden_layer_sizes = (5, 5), max_iter_mlp = 1000, random_state = 1, 
            max_num_probes = 1, interval = 50, num_rand_basis = 400, seed = None, score_mode = 'EI', 
            print_level = 1):

        if type(cv) is str and cv == '':
            cv = None

        self.method = method
        self.model = None
# CVをした結果を受け取るためのdict変数のkey
        self.named_steps_key = None
        self.mean = None
        self.std  = None
        self.variance = None
        self.score_mode = score_mode

        if print_level:
            print()
            print(f"Fitting by [{method}]")

        if method == 'linear':
            model = LinearRegression()

        elif len(method) >= 5 and method[:4] == 'poly':
            order = pint(method[4:], None)
            if order is None:
                print(f"\nError in tkmlr.get_model(): Invalid order=[method[4:]] given by [{method}]")
            model = Pipeline([('poly',   PolynomialFeatures(degree = order)),
                              ('linear', LinearRegression())])

        elif len(method) >= 5 and method[:5] == 'ridge':
            if print_level:
                print("Ridge Regression:")

            order = self.extract_order(method, 5)
            if order is None:
                if not use_cv:
                    model = Ridge(alpha = alpha)
                else:
                    model = RidgeCV()
                    self.named_steps_key = 'ridgecv'
            else:
                if not use_cv:
                    model = make_pipeline(
                        PolynomialFeatures(degree = order),
                        Ridge(alpha=0.1)
                        )
                else:
                    model = make_pipeline(
                        PolynomialFeatures(degree = order),
                        RidgeCV()
                        )
                    self.named_steps_key = 'ridgecv'
            
        elif len(method) >= 5 and method[:5] == 'lasso':
            if print_level:
                print("Lasso Regression:")

            order = self.extract_order(method, 5)
            if order is None:
                if not use_cv:
                    model = Lasso(alpha = alpha, max_iter = max_iter, tol = tol)
                else:
                    model = LassoCV(max_iter = max_iter, tol = tol)
                    self.named_steps_key = 'lassocv'
            else:
                if not use_cv:
                    model = make_pipeline(
                        PolynomialFeatures(degree = order),
                        Lasso(alpha = alpha, max_iter = max_iter, tol = tol)
                        )
                else:
                    model = make_pipeline(
                        PolynomialFeatures(degree = order),
                        LassoCV(max_iter = max_iter, tol = tol)
                        )
                    self.named_steps_key = 'lassocv'

        elif len(method) >= 5 and method[:7] == 'elastic':
            model = ElasticNet(alpha = alpha, l1_ratio = l1_ratio, max_iter = max_iter, tol = tol)
            order = self.extract_order(method, 7)
            if order is None:
                if not use_cv:
                    model = ElasticNet(alpha = alpha, l1_ratio = l1_ratio, max_iter = max_iter, tol = tol)
                else:
                    model = ElasticNetCV(cv = cv, max_iter = max_iter, tol = tol)
                    self.named_steps_key = 'lassocv'
            else:
                if not use_cv:
                    model = make_pipeline(
                        PolynomialFeatures(degree = order),
                        ElasticNet(alpha = alpha, l1_ratio = l1_ratio, max_iter = max_iter, tol = tol)
                        )
                else:
                    model = make_pipeline(
                        PolynomialFeatures(degree = order),
                        ElasticNetCV(cv = cv, max_iter = max_iter, tol = tol)
                        )
                    self.named_steps_key = 'elasticnetcv'

        elif method == 'mlp':
            if print_level:
                print("Multi-layer prceptron regression:")

            if not use_cv:
                if type(hidden_layer_sizes) is str:
                    hidden_layer_sizes = hidden_layer_sizes.split(',')
                model = MLPRegressor(hidden_layer_sizes = hidden_layer_sizes, max_iter = max_iter_mlp, random_state = 1, alpha = 3) 
            else:
                param_grid = {
                    'hidden_layer_sizes': [(5, 5), (10, 10), (5, 10, 5)],
                    'alpha': [0.0001, 0.001, 0.01, 0.1, 1.0],
                    'max_iter': [max_iter]
                    }
                grid_search = GridSearchCV(MLPRegressor(), param_grid, cv=5)
                grid_search.fit(descriptor_df, objective_df)

                best_params = grid_search.best_params_
                best_model = grid_search.best_estimator_

                if print_level:
                    print("Best Hyperparameters:", best_params)
                    print("Best Model:", best_model)

                model = MLPRegressor(**best_params)

        elif method == 'physbo':
            print("")
            print("Bayes search by PHYSBO:")
            self.physbo = imp.import_lib("physbo", stop_by_error = True)

            descriptor_copy = descriptor_df.iloc[0:3].copy()
            descriptor_df = pd.concat([descriptor_df, descriptor_copy], ignore_index = True, axis = 0)

            descriptors_np     = descriptor_df.to_numpy()
#            objective_train_np = objective_df.index.to_numpy()
            objective_train_np = objective_df.to_numpy().T[0]
            idx_train          = np.array(range(objective_train_np.shape[0]))
#            idx_train          = np.array(objective_df.index)

            self.policy = self.physbo.search.discrete.policy(
                    test_X = descriptors_np, initial_data = (idx_train, objective_train_np))

            model = None

        elif len(method) >= 8 and method[:7] == 'GPy:gpr':
            GPy = imp.import_lib("GPy", stop_by_error = True)
            order = pint(method[7:], None)
            if order is None:
                print(f"\nError in tkmlr.fit(): Invalid order=[method[3:]] given by [{method}]")

            if print_level:
                print(f"GPy:Gaussian process regression with white and order {order} polynomial kernels :")

            ndim = len(descriptor_df.columns)
            custom_kernel = GPy.kern.Bias(input_dim = ndim)
            for norder in range(1, order):
                custom_kernel += GPy.kern.Poly(input_dim = ndim, order = order)

            x = descriptor_df.values
            y = objective_df.values
#            print("x=", x.shape)
#            print("y=", y.shape)
            model = GPy.models.GPRegression(x, y, kernel = custom_kernel)

        elif len(method) >= 4 and method[:3] == 'gpr':
            order = pint(method[3:], None)
            if order is None:
                print(f"\nError in tkmlr.fit(): Invalid order=[method[3:]] given by [{method}]")

            if print_level:
                print(f"Gaussian process regression with white and order {order} polynomial kernels :")

            custom_kernel = C() #+ WhiteKernel() + RBF()
            for norder in range(1, order + 1):
                custom_kernel += DotProduct(1.0)**norder

            model = GaussianProcessRegressor(kernel = custom_kernel, alpha = alpha)

        elif method == 'gpr' or method == 'gp':
            if print_level:
                print("Gaussian process regression:")

#           kernel = Matern(length_scale = 2.0, nu = 2.5)
#            kernel = C(1.0, (1.0e-3, 1.0e3)) * Matern(length_scale = 1.0, length_scale_bounds = (1.0e-2, 1.0e1), nu = 2.5)
#           kernel = RBF(length_scale=2.0)
            kernel = C(1.0, (1.0e-3, 1.0e3)) * RBF(1.0, (1.0e-4, 1.0e2))
            model = GaussianProcessRegressor(kernel = kernel, n_restarts_optimizer = n_restarts_optimizer, alpha = alpha)

            '''
            if not use_cv:
                model = GaussianProcessRegressor()
            else:
                param_grid = {
                    'alpha': [0.1, 0.25, 0.5, 0.75, 1.0],
                    'kernel': [PolynomialKernel(degree=2), PolynomialKernel(degree=3), PolynomialKernel(degree=4)]
                }
                grid_search = GridSearchCV(GaussianProcessRegressor(), param_grid, cv=5)
                grid_search.fit(descriptor_df, column_or_1d(objective_df))

                best_params = grid_search.best_params_
                best_model = grid_search.best_estimator_

                if print_level:
                    print("Best Hyperparameters:", best_params)
                    print("Best Model:", best_model)

                model = GaussianProcessRegressor(**best_params)
            '''
        else:
            print(f"\nError in tkmlr.get_model(): Invalide method [{method}]")
            exit()

        self.method = method
        self.model = model
        
        return self.model

    def fit(self, descriptor_df, objective_df, method = 'mlp', model = None, 
            use_cv = False, cv = 5, 
            n_restarts_optimizer = 50,
            alpha = 0.1, l1_ratio = 0.5, 
            hidden_layer_sizes = (5, 5), max_iter_mlp = 1000, random_state = 1, 
            max_num_probes = 1, interval = 50, num_rand_basis = 400, seed = None, score_mode = 'EI', 
            max_iter = 1000, tol = 1.0e-4, print_level = 1):

        if method is None: method = self.method

        if model is None:
            self.model = self.get_model(descriptor_df, objective_df, method = method, 
                                use_cv = use_cv, cv = cv, 
                                n_restarts_optimizer = n_restarts_optimizer,
                                alpha = alpha, l1_ratio = l1_ratio, 
                                hidden_layer_sizes = hidden_layer_sizes, max_iter_mlp = max_iter_mlp, random_state = random_state, 
                                max_num_probes = max_num_probes, interval = interval, 
                                num_rand_basis = num_rand_basis, seed = seed, score_mode = score_mode, 
                                print_level = print_level)
            model = self.model
        elif model is None:
            model = self.model

        if method == 'physbo':
            descriptors_np = descriptor_df.to_numpy()

            if type(seed) is int: self.policy.set_seed(seed)

# bayes_searchは、simulator = Noneでは actions が返り、simulatorに関数を渡すと Hisotry object が返ってくる
            print("")
            print("Start Bayes search:")
            actions = self.policy.bayes_search(max_num_probes = max_num_probes, simulator = None, 
                        score = score_mode, interval = interval, num_rand_basis = num_rand_basis)

            print("")
            print("Bayse search:")
            print("show_search_results")
            self.physbo.search.utility.show_search_results(self.policy.history, 10)

# Hisotry objectの取得
            self.res = self.policy.export_history()
            self.best_fx, self.best_action = self.res.export_all_sequence_best_fx()
            self.bayes_x = self.res.chosen_actions
            self.x_bayes = descriptors_np[self.bayes_x]
            self.y_bayes = self.res.fx

# 獲得関数
            self.score = self.policy.get_score(mode = score_mode, xs = descriptors_np)

# 回帰。事後分布の平均値、分散
            self.mean      = self.policy.get_post_fmean(descriptors_np)
            self.variance  = self.policy.get_post_fcov(descriptors_np)
            self.std       = np.sqrt(self.variance)
            self.mean_m_sigma = self.mean - self.std
            self.mean_p_sigma = self.mean + self.std
#            print("objective_df=", objective_df)
#            print("self.mean=", self.mean)

        elif method == 'mlp':
            model.fit(descriptor_df, column_or_1d(objective_df))

        elif len(method) >= 8 and method[:7] == 'GPy:gpr':
            model.optimize(messages = True)
            '''
#            X_test = np.linspace(0, 10, 100)[:, np.newaxis]
#           y_pred, var = model.predict(X_test)
            y_pred, var = model.predict(x)

            plt.figure(figsize=(8, 6))
            
            plt.scatter(y, y_pred, label = 'train')
#            plt.scatter(x, y, c='r', label='Data')
#            plt.plot(X_test, y_pred, 'b-', label='Prediction')
#            plt.fill_between(X_test.ravel(), y_pred[:, 0] - np.sqrt(var[:, 0]), y_pred[:, 0] + np.sqrt(var[:, 0]), color='blue', alpha=0.2)
            plt.title('Gaussian Process Regression with Polynomial and Constant Kernel (GPy)')
            _legend = plt.legend()
            _legend.set_draggable(True)
            plt.pause(1.0e-5)
            input(">>")
            '''
        else:
            model.fit(descriptor_df, objective_df)

        if self.named_steps_key is not None:
            optimul_vals = model.named_steps[self.named_steps_key]
            print("  Cross validation resutls:")
#            print("  Optimized hyper parameters:", optimul_vals)
            for key in ["alpha_", "l1_ratio_", "best_score_", "cv_values_", "intercept_", 
                        "alphas_", "mse_path_", 
                        "cv_alphas_", "cv_alphas_std_"]:
                if hasattr(optimul_vals, key):
                    print(f"    Optimized {key}:", getattr(optimul_vals, key, '-'))


        self.model = model

        return self.model
    
    def predict(self, df, columns, method = None, score_mode = None):
        if method is None:
            method = self.method
        if score_mode is None:
            score_mode = self.score_mode    
#        print("method=", method)
        
        if len(method) >= 8 and method[:7] == 'GPy:gpr':
            y_pred, var = self.model.predict(df.values)
            self.variation = var
            self.std = sqrt(var)
            return pd.DataFrame(y_pred, columns = ['predicted'])

        elif len(method) >= 3 and method[:3] == 'gpr':
            predict_np, std = self.model.predict(df, return_std = True)
            df = pd.DataFrame(predict_np, columns = columns)
            self.std = std
            self.variation = pow(std, 2)
            self.mean         = predict_np
            self.mean_m_sigma = self.mean - self.std
            self.mean_p_sigma = self.mean + self.std
            return df

        elif method == 'physbo':
# 回帰。事後分布の平均値、分散
            X_std = df.to_numpy()
            self.score    = self.policy.get_score(mode = score_mode, xs = X_std)
#            print("x=", X_std, self.score)
            self.variance = self.policy.get_post_fcov(X_std)
            self.std      = np.sqrt(self.variance)
            self.mean     = self.policy.get_post_fmean(X_std)
            self.mean_m_sigma = self.mean - self.std
            self.mean_p_sigma = self.mean + self.std
            
            return pd.DataFrame(self.mean, columns = columns)

        else:
            predict_np = self.model.predict(df)
            df = pd.DataFrame(predict_np, columns = columns)
            self.variation = None
            self.std = None
            return df

    def cv_available(self, method):
        if "gp" in method or "lasso" in method or "ridge" in method or "elastic" in method or cfg.method in ["mlp"]:
            return True 
        else:
            return False

    def fit_cv_alpha(self, d_train, o_train, d_test, o_test,
                method = None,
                surrogate_model = "", surrogate_method = '', surrogate_scores = "mean|max", nmaxdata = 20,
                nmaxiter_surrogate = 10, tol_surrogate = 1.0e-5,
                n_restarts_optimizer = 50, alpha = 1.0e-4, l1_ratio = 0.5,
                hidden_layer_sizes = [3, 3, 3], max_iter_mlp = 1000,
                max_num_probes = 10, interval = 50, num_rand_basis = 400, score_mode = "EI", seed = None,
                log_alpha_min = 10, log_alpha_max = 4,
                 max_iter = 20, tol = 1.0e-5, print_level = 1):
        use_cv = self.cv_available(method)
        if use_cv:
            irange = range(log_alpha_min, log_alpha_max + 1)
        else:
            irange = [0]

        for log_alpha in irange:
            alpha = pow(10.0, log_alpha)
            self.fit(
                 descriptor_df = d_train,
                 objective_df  = o_train,
                 method = method, model = None,
                 n_restarts_optimizer = n_restarts_optimizer,
                 hidden_layer_sizes = hidden_layer_sizes, max_iter_mlp = max_iter_mlp,
                 alpha = alpha, l1_ratio = l1_ratio, 
                 max_num_probes = max_num_probes, interval = interval,
                 num_rand_basis = num_rand_basis, score_mode = score_mode, seed = seed,
                 max_iter = max_iter, tol = tol, print_level = print_level)

            predict_train = self.predict(d_train, ["predicted (train)"])
            predict_test  = self.predict(d_test,  ["predicted(test)"])
            mae_train, mse_train, r2_train = self.error(o_train, predict_train)
            mae_test,  mse_test,  r2_test  = self.error(o_test,  predict_test)

            if use_cv:
                print(f"  alpha={alpha:12.6g}")
                print(f"  train: MAE={mae_train:10.4g} MSE={mse_train:10.4g} R2={r2_train:10.4g}")
                print(f"  test : MAE={mae_test:10.4g} MSE={mse_test:10.4g} R2={r2_test:10.4g}")

        return predict_train, predict_test

    def error(self, df1, df2):
        mae = mean_absolute_error(df1, df2)
        mse = mean_squared_error(df1, df2)
        r2  = r2_score(df1, df2)

        return mae, mse, r2

#f_minは既知の最小値
    def EI(UCB, mean, std, f_min):
        Z = (mean - f_min) / std
        ei = std * (Z * norm.cdf(Z) + norm.pdf(Z))
        return ei

    def PI(mean, std, f_min):
        Z = (mean - f_min) / std
        pi = norm.cdf(Z)
        return pi

# typically kappa in [0, 2]. larger value would explore more challenging region
    def UCB(self, mean, std, kappa):
        ucb = mean + kappa * std
        return ucb

    def TS(self, means, covs):
        samples = [np.random.multivariate_normal(mean, cov) for mean, cov in zip(means, covs)]
        return np.argmax(samples)

    def probability_of_improvement_for_multiple_points(self, means, stds, f_min):
        z_values = (means - f_min) / stds
        pis = norm.cdf(z_values)
        return np.argmax(pis)
    
    def max_value_entropy_search(self, means, stds):
        entropies = np.log(stds)
        return np.argmax(entropies)

    def expected_improvement_per_second(means, stds, f_min, times):
        z_values = (means - f_min) / stds
        eis = stds * (z_values * norm.cdf(z_values) + norm.pdf(z_values))
        eis_per_second = eis / times
        return np.argmax(eis_per_second)
    
    def knowledge_gradient(self, means, stds):
        gradients = means + stds * log(stds)
        return np.argmax(gradients)
    
    def aquisition(self, model = None, score_mode = 'EI'):
        if model is None:
            model = self.model

        if len(method) >= 8 and method[:7] == 'GPy:gpr':
            GPyOpt = imp.import_lib("GPyOpt", stop_by_error = True)
            
            if score_mode == 'EI':
                acquisition = GPyOpt.acquisitions.AcquisitionEI(self.model)
            elif score_mode == 'PI':
                acquisition = GPyOpt.acquisitions.AcquisitionPI(self.model)
            elif score_mode == 'UCB':
                acquisition = GPyOpt.acquisitions.AcquisitionUCB(self.model)
            else:
                print(f"\nError in tkmlr.aquisition(): Invalid aquisition function [{acquisition_func}] for [{model}]")
                return None
            
            next_point = GPyOpt.methods.BayesianOptimization(model, acquisition).suggest_next_locations()
            return next_point

    '''
    def approximate_gradient(gpr, x, epsilon=1e-5):
#    ガウス過程モデルから勾配情報を近似する関数
    # xの各次元ごとに微小な変化を加える
        x_plus = x + epsilon
        x_minus = x - epsilon

    # 各次元ごとの微小な変化における予測値を取得
        y_plus = gpr.predict(x_plus.reshape(1, -1))
        y_minus = gpr.predict(x_minus.reshape(1, -1))

    # 勾配を計算
        gradient = (y_plus - y_minus) / (2 * epsilon)

        return gradient
    
    def stein_score(gpr, x):
#    Stein Scoreを計算する関数
# 勾配情報をサンプリング
        gradients = approximate_gradient(gpr, x)

# Stein Scoreを計算
        self.score = np.mean(np.dot(gradients, gradients.T))

        return self.score

# ダミーデータの生成
    X = np.random.uniform(0, 10, (20, 1))
    y = np.sin(X)

# ガウス過程モデルの構築
    kernel = RBF(length_scale=1.0)
    gpr = GaussianProcessRegressor(kernel=kernel)
    gpr.fit(X, y)

# Stein Scoreを計算
    x_next = np.array([5.0])  # 次の評価ポイントの候補
    self.score = stein_score(gpr, x_next)
    print("Stein Score for next evaluation point:", self.score)    
    '''

    def inf_format_func(self, plotevent, sel, iline, idata, label, inf, x = 0, y = 0, opt = None):
        if inf is None:
            return

        x_label = inf["x_label"]
        y_label = inf["y_label"]
        input_list    = inf["x_list"]
        predict_list  = inf["y_list"]
        if isinstance(input_list, pd.DataFrame):
            input_list = input_list.to_numpy().T[0]
        if isinstance(predict_list, pd.DataFrame):
            predict_list = predict_list.to_numpy().T[0]

        descriptor_df = inf.get("descriptor_df", None)
        objective_df  = inf.get("objective_df", None)

        scale_keys = self.scale_dict.keys()

        print(f"{label}: line #{iline}  data#{idata}:")
#        print("columns=", descriptor_columns)
#        print("scale=", self.scale_dict)

        def invert_scale(v, scale = ''):
            if scale == 'log':
                v_conv = exp(v)
                str_scale = 'log(x)'
            else:
                v_conv = v
                str_scale = 'x'

            return v_conv, str_scale

        print(f"  {x_label}={input_list[idata]}  {y_label}={predict_list[idata]:g}")
#        print("descriptor_df=", descriptor_df)
        if descriptor_df is None:
            return
        else:
            ferror = False
            descriptor_df_org_np = self.inverse_standardize(descriptor_df)
            descriptor_columns = descriptor_df.columns
            for icol, l in enumerate(descriptor_columns):
                try:
                    val_std   = descriptor_df[l][idata]
                    val_nostd = descriptor_df_org_np[idata][icol]
                    scale = ''
                    if l in scale_keys:
                        scale = self.scale_dict[l] 
                        if scale == 'log':
                            val_original = exp(val_nostd)

                    val_original, str_scale = invert_scale(val_nostd, scale)
                    print(f"  {l:10}: x={val_original:12.4g}  {str_scale:6}:{val_nostd:12.4g}  standardized:{val_std:12.4g}")
                except:
                    ferror = True
#                print(f"  Error in icol={icol} label={l} for idata#{idata}")

            if ferror:
                print(f"  ** Dataframe label error(s) occured for idata#{idata}")
                _seriese = descriptor_df.iloc[idata, :]
#            print(_seriese)
                for icol, l in enumerate(descriptor_columns):
                    val_std      = _seriese[l]
                    val_nostd    = descriptor_df_org_np[idata][icol]
                    scale = ''
                    if l in scale_keys:
                        scale = self.scale_dict[l] 
                        if scale == 'log':
                            val_original = exp(val_nostd)

                    val_original, str_scale = invert_scale(val_nostd, scale)
                    print(f"  {l:10}: x={val_original:12.4g}  {str_scale:6}:{val_nostd:12.4g}  standardized:{val_std:12.4g}")

    def annotation_format_func(self, plotevent, sel, iline, idata, label, inf, x = 0, y = 0, opt = None):
        x_label = inf["x_label"]
        y_label = inf["y_label"]
        input_list    = inf["x_list"]
        predict_list  = inf["y_list"]
        if isinstance(input_list, pd.DataFrame):
            input_list = input_list.to_numpy().T[0]
        if isinstance(predict_list, pd.DataFrame):
            predict_list = predict_list.to_numpy().T[0]

#        text = f"{label}: line #{iline}  data#{idata}:\n x={x:g} y={y:g}"
        text = f"{label}: line #{iline}  data#{idata}:\n {x_label}={input_list[idata]:g} {y_label}={predict_list[idata]:g}"
        sel.annotation.set(text = text)

    def plot_accuracy(self, ax, fig, data_list = [], 
                xlabel = 'Input', ylabel = 'Prediction', 
                plot_event = None, register_events = False,
                show_legend = True, fontsize = None, legend_fontsize = None):

        self.fig = fig
        if plot_event:
            self.plot_event = plot_event
        else:
            self.plot_event = tkPlotEvent(plt)

        if register_events:
            if "stop_button" in register_events:
                self.plot_event.add_stop_button()
            if "annotation" in register_events:
                self.plot_event.prepare_annotation()
            if "move_text" in register_events:
                self.plot_event.prepare_move_text(fig)

        inf_format_func        = lambda plotevent, sel, iline, idata, label, inf, x = 0, y = 0, opt = None: \
                                    self.inf_format_func(plotevent, sel, iline, idata, label, inf, x = x, y = y, opt = opt)
        annotation_format_func = lambda plotevent, sel, iline, idata, label, inf, x = 0, y = 0, opt = None: \
                                    self.annotation_format_func(plotevent, sel, iline, idata, label, inf, x = x, y = y, opt = opt)

        ax.tick_params(labelsize = fontsize)

        xmin =  1.0e300
        xmax = -1.0e300
        ymin =  1.0e300
        ymax = -1.0e300
        for inf in data_list:
            label         = inf["label"]
            x_list        = inf["x_df"].to_numpy().T
            y_list        = inf["y_df"].to_numpy().T
            descriptor_df = inf["descriptor_df"]
            objective_df  = inf["objective_df"]
            color         = inf["color"]

            nx, ny = x_list.shape
            for i in range(nx):
                xl = x_list[i]
                yl = y_list[i]
                _label = f"{label}: object var#{i}"
                line = ax.scatter(x_list, y_list, label = _label, color = color)

                if register_events:
                    if "annotation" in register_events:
                        self.plot_event.annotation.add_line(_label, ax, ax, xl, yl, line, 
                                inf_list = {"x_label": "input", "y_label": "predict",
                                            "x_list": xl, "y_list": yl, "descriptor_df": descriptor_df, "objective_df": objective_df},
                                annotation_format = annotation_format_func, inf_format = inf_format_func)

                    if "move_text" in register_events:
                        self.plot_event.move_text.add_annotation(ax, ax, x_list = x_list, y_list = y_list, frac = None, 
                            text = label, fontsize = 16, ha = 'left', va = 'top', color = color, fc = "w", ec = "none")

            x_flat = x_list.flatten()
            _xmin = min(x_flat)
            _xmax = max(x_flat)
            y_flat = y_list.flatten()
            _ymin = min(y_flat)
            _ymax = max(y_flat)
#            print("range:", _xmin, _xmax, _ymin, _ymax)

            if type(_xmin) is list or type(_xmin) is tuple or isinstance(_xmin, np.ndarray):
                _xmin = min(_xmin)
                _xmax = max(_xmax)
            if type(_ymin) is list or type(_ymin) is tuple or isinstance(_ymin, np.ndarray):
                _ymin = min(_ymin)
                _ymin = max(_ymax)

            xmin = min([xmin, _xmin])#, _ymin])
            xmax = max([xmax, _xmax])#, _ymax])
            ymin = min([ymin, _ymin])#, _ymin])
            ymax = max([ymax, _ymax])#, _ymax])

        if type(xmin) is list or type(xmin) is tuple or isinstance(xmin, np.ndarray):
            _xmin = min(xmin)
        if type(xmax) is list or type(xmax) is tuple or isinstance(xmax, np.ndarray):
            xmax = min(xmax)
#        print("x range:", xmin, xmax)

        ax.plot([xmin, xmax], [xmin, xmax], color = 'red', linestyle = '--', linewidth = 0.5)
        ax.set_xlabel(xlabel, fontsize = fontsize)
        ax.set_ylabel(ylabel, fontsize = fontsize)
        if show_legend:
            _legend = ax.legend(fontsize = legend_fontsize)
            _legend.set_draggable(True)

        if register_events:
            if "stop_button" in register_events:
                self.plot_event.add_stop_button()

            if "annotation" in register_events:
                self.plot_event.register_annotation_event(fig, activate = True, print_level = 0)

            if "move_text" in register_events:
                self.plot_event.register_move_text_event(fig)

        return ax

    def plot_index_fmin(self, ax, fig, data_list = [],
                objective_label = 'y', score_mode = None, 
                plot_event = None, register_events = False,
                fontsize = None, legend_fontsize = None):

        if plot_event:
            self.plot_event = plot_event
        else:
            self.plot_event = tkPlotEvent(plt)

        if register_events:
            if "stop_button" in register_events:
                self.plot_event.add_stop_button()
            if "annotation" in register_events:
                self.plot_event.prepare_annotation()
            if "move_text" in register_events:
                self.plot_event.prepare_move_text(fig)

        if score_mode is None:
            score_mode = self.score_mode

        inf_format_func = lambda plotevent, sel, iline, idata, label, inf, x = 0, y = 0, opt = None: \
                                    self.inf_format_func(plotevent, sel, iline, idata, label, inf, x = x, y = y, opt = opt)
        annotation_format_func = lambda plotevent, sel, iline, idata, label, inf, x = 0, y = 0, opt = None: \
                                    self.annotation_format_func(plotevent, sel, iline, idata, label, inf, x = x, y = y, opt = opt)

        ins = None
        offset = 0
        axb = None
        for inf in data_list:
            label         = inf.get("label",       None)
            input_o       = inf.get("input_df",    None)
            predict       = inf.get("predict_df",  None)
            color         = inf.get("color",       None)
            mean_m_sigma  = inf.get("-sigma_df",   None)
            mean_p_sigma  = inf.get("+sigma_df",   None)
            color_sigma   = inf.get("color_sigma", None)
            score         = inf.get("score_df",    None)
            color_score   = inf.get("color_score", None)
            descriptor_df = inf.get("descriptor_df", None)
            objective_df  = inf.get("objective_df",  None)
            alpha         = inf.get("alpha", 0.5)

            ntrain  = len(input_o)
            index_o = list(range(offset, offset + ntrain))

            ins1 = ax.plot(index_o, input_o, label = f'input ({label})',
                        linestyle = '',  marker = 'o', markerfacecolor = color, markeredgecolor = color, markersize = 8.0)
            ins2 = ax.plot(index_o, predict, label = f'predict ({label})', 
                        linestyle = '-', color = color, linewidth = 0.5)

            if ins is None:
                ins = ins1 + ins2
            else:
                ins += ins1 + ins2

            if mean_m_sigma is not None:
                ax.fill_between(index_o, mean_m_sigma, mean_p_sigma, color = color_sigma, alpha = alpha)

            ins6 = None
            if score is not None:
                if axb is None:
                   axb = ax.twinx()
                   axb.tick_params(labelsize = fontsize)

                ins6 = axb.plot(index_o, score, label = f'score {score_mode} ({label})', 
                            linestyle = 'dashed', linewidth = 0.5, color = color_score)
#               ins7 = ax3.plot(idx_best, mean[idx_best], label = 'best candidate', linestyle = '', marker = '*', markersize = 15.0)
                ins += ins6

                axb.set_ylabel(f"score {score_mode}", fontsize = fontsize)

            if register_events:
                if axb is None:
                    ax2 = ax
                else:
                    ax2 = axb
                
                if "annotation" in register_events:
                    self.plot_event.annotation.add_line(f"input {label}", ax, ax2, index_o, input_o, ins1[0], 
                        inf_list = {"x_label": "index", "y_label": objective_label,
                                    "x_list": index_o, "y_list": input_o, "input_df": input_o, "predict_df": predict,
                                    "descriptor_df": descriptor_df, "objective_df": objective_df},
                        annotation_format = annotation_format_func, inf_format = inf_format_func)

                    self.plot_event.annotation.add_line(f"predict {label}", ax, ax2, index_o, predict, ins2[0], 
                        inf_list = {"x_label": "index", "y_label": objective_label,
                                    "x_list": index_o, "y_list": input_o, "input_df": input_o, "predict_df": predict,
                                    "descriptor_df": descriptor_df, "objective_df": objective_df},
                        annotation_format = annotation_format_func, inf_format = inf_format_func)

                    if score is not None:
                        self.plot_event.annotation.add_line(f'score {score_mode} ({label})', ax2, ax2, index_o, score, ins6[0], 
                            inf_list = {"x_label": "index", "y_label": "score {score_mode}",
                                        "x_list": index_o, "y_list": score, "input_df": input_o, "predict_df": predict,
                                        "descriptor_df": descriptor_df, "objective_df": objective_df},
                            annotation_format = annotation_format_func, inf_format = inf_format_func)

                if  "move_text" in register_events:
                    self.plot_event.move_text.add_annotation(ax, axb, x_list = index_o, y_list = predict, 
                           x_offset = 0, frac = 0.5, 
                           text = label, fontsize = 16, ha = 'left', va = 'top', color = color, fc = "w", ec = "none")

            offset += ntrain

        _legend = ax.legend(ins, [l.get_label() for l in ins], fontsize = legend_fontsize, loc = 'best')
        _legend.set_draggable(True)
        ax.minorticks_on()
        ax.grid(which = "major", axis = "x", color = "green", alpha = 0.5, linestyle = '--', linewidth = 0.5)
        ax.grid(which = "minor", axis = "x", color = "green", alpha = 0.5, linestyle = '--', linewidth = 0.1)
        ax.set_xlabel ( "index", fontsize = fontsize, fontname = "MS Gothic")
        ax.set_ylabel (objective_label, fontsize = fontsize, fontname = "MS Gothic")

        if register_events:
            if "stop_button" in register_events:
                self.plot_event.add_stop_button()

            if "annotation" in register_events:
                self.plot_event.register_annotation_event(fig, activate = True, print_level = 0)

            if "move_text" in register_events:
                self.plot_event.register_move_text_event(fig)
