import os
import random
import tqdm
import numpy as np
from numpy import sin, cos, tan, pi, exp, log, sqrt
from scipy.optimize import minimize
import openpyxl
import pandas as pd
from matplotlib import pyplot as plt
import matplotlib.widgets as wg
from functools import lru_cache
#from sklearn.gaussian_process import GaussianProcessRegressor
#from sklearn.gaussian_process.kernels import Matern, RBF, ConstantKernel as C


from tklib.tkutils import terminate, pint, pfloat, getarg, getintarg, getfloatarg, pconv_by_type, str2val
#analyze_varstr, save_csv
from tklib.tkutils import format_strlist, sort_lists, is_exist, is_file, is_dir, val2str
from tklib.tkinifile import tkIniFile
from tklib.tkparams import tkParams
from tklib.tkvariousdata import tkVariousData
from tklib.tkgraphic.tkplotevent import tkPlotEvent

from tklib.tkvariousdata import tkVariousData
from tklib.tkfilter import tkFilter
from tklib.tkexcel import tkExcel
from tklib.tksci.tksci import pi, h, hbar,c, e, kB, me, Ea_Arrhenius, log10
from tklib.tksci.tkoptimize import tkOptimize, mlsq_general
from tklib.tksci.tkFit_object import tkFit_object
from tklib.tksci.tkFit_lib import save_data, add_history, inf_format_func, convert_scale, P
from tklib.tksci.tkoptimize_direct import particle_swarm_optimization, genetic_algorithm
from tklib.tksci.tkoptimize_direct import simulated_annealing, replica_exchange_monte_carlo
from tklib.tksci.tkmlr import tkMLData, tkMLR


'''
Fitting and interactive matplotlb helper library
'''


_global = tkParams()

#@lru_cache(maxsize = None)
def minimize_func_cachable(xk):
        print_level = _global.print_level
        fit = _global.fit
        x_list = _global.x_list
        y_list = _global.y_list
        w_list = _global.w_list
        run    = _global.run

        if fit.icall > fit.nmaxcall:
            if print_level:
                print(f"\nMessage: # of minimize_func calls [{fit.icall}] exceesds nmaxcall [{fit.nmaxcall}].")
                print(f"           Will stop optimization\n")
            fit.stop_flag = True

        if fit.stop(fit.cfg.stopfile): return 1.0e300

        if x_list is None: x_list = fit.x_list
        if y_list is None: y_list = fit.y_list
        if w_list is None: w_list = fit.w_list

        if print_level >= 2:
            print()
            print(f"tkFit_mxy.minimize_func_cachable(): iter: {fit.iter}   icall: {fit.icall}")

        nx = len(x_list[0])
        xk_all = fit.recover_parameters(xk, set_member = False)
        fmin   = fit.cal_total_fmin(xk_all, x_list, y_list, w_list, run)

        if fmin is None: return 1.0e10

        xdiff = fit.normalized_difference(xk_all, fit.pk_prev, fit.optid)

        if print_level:
#            print(f"\n  fmin={fmin:.6g}")
            if xdiff is not None:
                print(f"  Normalized diffrence of parameters: {xdiff:.6g}")

        xatol = 1.0e-4 if not fit.hasattr('xatol') or fit.xatol is None else fit.xatol
        if xdiff is not None and xdiff > 0 and xdiff < xatol:
            if print_level:
                print(f"Warning: Normalized difference of parameters [{xdiff}] is smaller than eps={xatol:.4g}. Skip recording")
#        elif abs(fmin - fit.fmin_prev) < fit.tol:
#            print(f"Warning: fmin={fmin:g} is close to the previous value {fit.fmin_prev:g} within eps={fit.tol:.4g}. Skip recording")
        elif fmin < fit.fmax_record:
            if xdiff is None: xdiff = 0
            fit.save_history(xk_all, fmin, xdiff, save_time = True, print_level = print_level)
        else:
            if print_level:
                print(f"Warning: fmin={fmin:g} is larger than fmax_record={fit.fmax_record:g}. Skip recording")
        
        fit.icall += 1
        fit.fmin = fmin
        fit.fmin_prev = fmin
        fit.xdiff = xdiff
        fit.pk_prev = xk_all.copy()

        return fmin


class tkFit_mxy(tkFit_object):
    def __init__(self, app, parameter_file = None, 
                    method = 'nelder-mead', jac = '3-points', tol = 1.0e-5, xatol = None, nmaxiter = 100, nmaxcall = 1000,
                    fplot = 1, fdisconnect_backward_data = False,
                    fhistory = 0, ffitfiles = 0, print_interval = 1, plot_interval = 1, **args):
        super().__init__(**args)
        
#        super(tkParams, self).__init__(parameter_file, app, args)
#        self.update(**args)

        self.app = app

        self.is_fitting = True

        self.cal_ylist = self._cal_ylist
        self.cal_penalty = self._cal_penalty
        self.cal_fmin = self._cal_fmin
        self.callback = self._callback
        self.minimize_func = self._minimize_func

# tkPlotEvent is initialized in tkFit_object.initial_plot()
#        self.plot_event = tkPlotEvent(plt)
        self.x_labels  = None
        self.x_list    = None
        self.y_labels  = None
        self.y_list    = None
        self.y_scale   = None
        self.w_list    = None
        self.func      = None
        self.iter      = 0
        self.iter_surrogate = None
        self.f         = None

        self.ix_plot   = ''

        self.varname  = []
        self.unit     = []
        self.pk       = []
        self.pk_scale = []
        self.dx       = []
        self.optid    = []
        self.kmin     = []
        self.kmax     = []
        self.kpenalty = []
        self.pk_corrected   = None

        self.fmin           = 1.0e300
        self.fmin_prev      = 1.0e300
        self.pk_prev        = None
        self.xk_all_prev    = None
        self.best_xk_all    = None
        self.best_fmin      = None
        self.best_xk_all    = None
        self.best_sm_func_mean = None
        self.best_sm_func_std = None
        self.last_message   = "No message"
        
        self.__method       = method
        self.jac            = jac
        self.tol            = tol
        self.xatol          = xatol
        self.nmaxiter       = nmaxiter
        self.nmaxcall       = nmaxcall
        self.icall          = 0

        self.fmax_record    = 1.0e300

        self.print_interval = print_interval
        self.plot_interval  = plot_interval
        self.fplot          = fplot
        self._fdisconnect_backward_data = fdisconnect_backward_data
        self.fhistory       = fhistory
        self.ffitfiles      = ffitfiles
        self.stop_flag      = False

        self.cached_history = {}
        selfhistoryfile = None
        self.fp_history = None
        self.wb_history = None
        self.ws_history = None

        self.iter_list = []
        self.fmin_list = []
        self.iter_list_exact = []
        self.fmin_list_exact = []
        self.xdiff_list = []
        self.std_iter_list = []
        self.std_list = []
        self.sm_func_std = None


    @property
    def fdisconnect_backward_data(self):
        return self._fdisconnect_backward_data
        
    @fdisconnect_backward_data.setter
    def fdisconnect_backward_data(self, fdisconnect_backward_data):
        self._fdisconnect_backward_data = fdisconnect_backward_data

    @property
    def method(self):
        return self.__method
        
    @method.setter
    def method(self, method):
        if method == 'simplex':
            self.__method = 'nelder-mead'
        else:
            self.__method = method

    def convert_scale(self, x, y_scale = 'log', eps = None):
        return convert_scale(x, y_scale, eps)

    def open_history(self, historyfile = None):
        if historyfile is None: historyfile = self.historyfile
        else: self.historyfile = historyfile

        if not os.path.isfile(historyfile):
            self.wb = openpyxl.Workbook()
        else:
            self.wb = openpyxl.load_workbook(historyfile)

        self.ws = self.wb.active

        return self.wb, self.ws

    def add_history(self, historyfile, varnames, xk, fmin, xdiff = None, 
                    wb = None, ws = None, save_time = False, print_level = 1):
        if wb is None:
            if self.get("wb", None):
                wb = self.wb
            else:
                if historyfile is None:
                    historyfile = self.historyfile

                if not os.path.isfile(historyfile):
                    if print_level:
                        print(f"tkFit_mxy_flex.add_history() 265 **create [{historyfile}]")
                    wb = openpyxl.Workbook()
                else:
                    if print_level:
                        print(f"tkFit_mxy_flex.add_history() 267 **append to [{historyfile}]")
                    wb = openpyxl.load_workbook(historyfile)

        if wb is None: 
            print("Error in tkFit_mxy.add_history(): None wb")
            exit()

        if ws is None:
            ws = wb.active

        return add_history(historyfile, varnames, xk, fmin, xdiff = xdiff, wb = wb, ws = ws, save_time = save_time, print_level = print_level)

    def close_history(self, historyfile = None, wb = None):
        if historyfile is None: historyfile = self.historyfile
        if wb is None: wb = self.wb

        try:
            wb.save(historyfile)
        except:
            return False

        self.fp_history = None
        self.wb_history = None
        self.ws_history = None
        
        return True

    def save_history(self, xk_all, fmin, xdiff = None, 
                    save_func = None, save_time = False, print_level = 1):
        xk_tuple = tuple(xk_all)
        if xk_tuple in self.cached_history.keys(): return

        self.cached_history[xk_tuple] = True

        if save_func is None: save_func = save_data

        if self.cfg.fhistory:
            if print_level:
                print(f"tkFit_mxy_flex.save_history(): Add fmin and parameters to [{self.historyfile}]")
            iter = self.add_history(self.historyfile, self.varname, xk_all, fmin, xdiff = xdiff, save_time = save_time, print_level = print_level)
            if iter is None:
                print(f"  Warning in tkFit_mxy_flex.save_history(): Failed to add to {self.historyfile}]")
                print(f"           Check write permission / disk space etc")
            elif print_level:
                if xdiff is None:
                    print(f"DB #{iter}: fmin={fmin:12.4g} ", end = '')
                else:
                    print(f"DB #{iter}: fmin={fmin:12.4g} xdiff={xdiff:12.4g} ", end = '')
                for v in xk_all:
                    print(f" {v:12.6g}".strip() + " ", end = '')
                if print_level: print()

        if self.cfg.ffitfiles and iter is not None:
            save_path = f'fit{iter:04}.xlsx'
            if print_level:
                print(f"Save last input and calculation data  to [{save_path}]")
            if hasattr(self, 'yini_list'):
                labels = self.x_labels + self.y_labels + ['initial', 'final']
                data_list = self.x_list + self.y_list + [self.yini_list[0], self.yc_list[0]]
#                labels = ['x', 'input', 'initial', 'final']
#                data_list = [self.x_list[0], self.y_list[0], self.yini_list[0], self.yc_list[0]]
                ret = save_func(save_path, labels, data_list)
            else:
                labels = self.x_labels + self.y_labels + ['final']
                data_list = self.x_list + self.y_list + [self.yc_list[0]]
#                labels = ['x', 'input', 'final']
#                data_list = [self.x_list[0], self.y_list[0], self.yc_list[0]]
#                for i in range(len(labels)):
#                    print(f"{i}: {len(data_list[i])}:", data_list[i])
                ret = save_func(save_path, labels, data_list)
            if not ret:
                print(f"  Warning: Failed to save to [{save_path}]")
                print(f"           Check write permission / disk space etc")

    def read_data_from_filter(self, infile, module_dir, module_name, imodule = 0, iregion = 0,
                x_indexes = [], y_indexes = [], app = None, cfg = None, is_print = False):
#        print("")
#        print(f"Load module by {module_dir}/{module_name}.py")
        filter = tkFilter(app = app, cfg = cfg, plugin_dir = module_dir, module_file = f"{module_name}.py")
        self.module_names, self.modules = filter.load(target = "read_data", is_print = True)
        self.inf_input = filter.read_data(imodule, infile, app = app, cfg = cfg, is_print = is_print)

        self.x_labels = []
        self.x_list   = []
        for idx in x_indexes:
            label, data = filter.find_data_array(idx)
            self.x_labels.append(label)
            self.x_list.append(data)

        self.y_labels = []
        self.y_list  = []
        for idx in y_indexes:
            label, data = filter.find_data_array(idx)
            self.y_labels.append(label)
            self.y_list.append(data)


    def read_data(self, infile, x_labels = [0], y_label = 1, y_labels = [1], xmins = [None], xmaxs = [None], usage = None,
                module_dir = None, module_name = None, imodule = 0, iregion = 0,
                x_indexes = [], y_indexes = [], app = None, cfg = None, is_print = False):
        if module_dir is not None and module_name is not None:
            return self.read_data_from_filter(infile = infile, 
                    module_dir = module_dir, module_name = module_name, imodule = imodule, iregion = iregion,
                    x_indexes = x_indexes, y_indexes = y_indexes, app = app, cfg = cfg, is_print = is_print)

        print("")
        print(f"Read [{infile}]")
        datafile = tkVariousData(infile)
        labels, datalist = datafile.Read_minimum_matrix(close_fp = True, force_numeric = False, usage = usage)
        
        _x_labels = []
        xins      = []
        for l in x_labels:
            x_label, xin = datafile.FindDataArray(l, flag = 'i')
#            print("x l=", l, x_label)
            _x_labels.append(x_label)
            xins.append(xin)
#        xlabel, xin = datafile.FindDataArray(xlabel, flag = 'i')

        _y_labels = []
        yins      = []
        for l in y_labels:
            y_label, yin = datafile.FindDataArray(l, flag = 'i')
            _y_labels.append(y_label)
            yins.append(yin)
        ylabel, yin = y_labels[0], yins[0]

        self.datafile  = datafile
        self.labels    = labels
        self.x_labels  = _x_labels
        self.y_label   = y_label
        self.y_labels  = _y_labels
        self.datalist  = datalist
        self.xdata_list = xins
        self.ydata_list = yins

        self.nx        = len(x_labels)
        self.ndata_all = len(xin)
#        print("ndata_all=", self.ndata_all)
        self.x_list = []
        for il in range(self.nx):
            self.x_list.append([])
        self.y = []
        self.included_index = []
        for i in range(self.ndata_all):
            is_included = True
            for il in range(len(_x_labels)):
                if xmins[il] is not None and xmins[il] > xins[il][i]:
                    is_included = False
                    continue
                if xmaxs[il] is not None and xmaxs[il] < xins[il][i]:
                    is_included = False
                    continue

                if is_included:
                    self.x_list[il].append(xins[il][i])

            if is_included:
                self.y.append(yin[i])
                self.included_index.append(i)

        self.ndata = len(self.x_list[0])
        self.index = range(self.ndata)

    def _cal_ylist(self, pk, x_list = None, run = True, print_level = 1):
        if x_list is None:
            x_list = self.x_list

        ny = len(self.y_list)

#        print("205 func: len(pk)=", len(pk))
        nvars = len(x_list)
        ndata = len(x_list[0])
        y_list = []
        for i in range(ny):
            y_list.append([])

        for i in range(ndata):
# i番目のxの値セット xs を取り出す
            xs = [x_list[j][i] for j in range(nvars)]
# xs に対するyのリストを取り出す
            ret = self.func(xs, pk)
            for i, v in enumerate(ret):
                y_list[i].append(v)

        return y_list

    def cal_difference(self, ymeas, ysim, y_scale = '', eps = None):
        if y_scale == 'log':
            if eps is not None:
                if ymeas <= 0.0:
                    ymeas = eps
                if ysim <= 0.0:
                    ysim = eps

            return log(ymeas) - log(ysim)

        return ymeas - ysim

    def correct_pk_range(self, xk_all, print_level = 1):
        _xk_all = xk_all.copy()
        kmin = self.kmin
        kmax = self.kmax
        nkp = len(xk_all)
        for i in range(nkp):
            if xk_all[i] < kmin[i]  :
                print(f"**Warning: [{self.varname[i]}={xk_all[i]}] is smaller than [{kmin[i]}].")
                _xk_all[i] = kmin[i]
            elif xk_all[i] > kmax[i]:
                print(f"**Warning: [{self.varname[i]}={xk_all[i]}] is larger than [{kmax[i]}].")
                _xk_all[i] = kmax[i]

        return _xk_all

    def _cal_penalty(self, xk_all, print_level = 1):
        cfg = self.app.get('cfg', None)
        if cfg: eps = cfg.get('eps', 1.0e-10)
        else: eps = 1.0e-10
        y_scale = self.app.cfg.get('y_scale', '')

        self.pk_given = xk_all.copy()

        kp   = self.kpenalty
        if kp is None: return 0.0

        kmin = self.kmin
        kmax = self.kmax
        nkp = len(kp)
        p_tot = 0.0
        for i in range(nkp):
            if xk_all[i] < kmin[i]  : origin_value = kmin[i]
            elif xk_all[i] > kmax[i]: origin_value = kmax[i]
            else: continue

            d = convert_scale(xk_all[i], y_scale, eps) - convert_scale(origin_value, y_scale, eps)
            p = kp[i] * d * d

            if xk_all[i] < kmin[i]  :
                print(f"**Warning: [{self.varname[i]}={xk_all[i]}] is smaller than [{kmin[i]}]. Add penalty [{p:10.3g}] to fmin")
#                print("  p=", y_scale, eps, "  conv=", convert_scale(xk_all[i], y_scale, eps), convert_scale(origin_value, y_scale, eps), "  d=", d, kp[i], p)
            elif xk_all[i] > kmax[i]:
                print(f"**Warning: [{self.varname[i]}={xk_all[i]}] is larger than [{kmax[i]}]. Add penalty [{p:10.3g}] to fmin")
#                print("  p=", y_scale, eps, "  conv=", convert_scale(xk_all[i], y_scale, eps), convert_scale(origin_value, y_scale, eps), "  d=", d, kp[i], p)

            p_tot += p
            xk_all[i] = origin_value

        self.pk_corrected = xk_all.copy()

        return p_tot

    def _cal_fmin(self, xk_all, x_list, y_list, w_list, run = True, print_level = 1):
        eps = self.app.cfg.get('eps', 1.0e-10)
        y_scale = self.app.cfg.get('y_scale', '')

        yc_list = self.cal_ylist(xk_all, x_list, run = run, print_level = print_level)
        if yc_list is None: 
            print("\nError in tkFit_mxy_flex._cal_fmin(): self.cal_ylist returned None\n")
            return None

        nx = len(x_list[0])
        nc = len(yc_list[0])
        if nc < nx: 
            print(f"\nWarning in cal_fmin(): # of data in yc_list[0]={nc} is smaller than that of x_list[0]={nx}\n")
            return None

        fmin = 0.0
        wtot = 0.0
        yerr_list = []
        for i in range(nx):
            ymeas = self.y_list[0][i]
            ysim  = yc_list[0][i]
            d = convert_scale(ymeas, y_scale, eps) - convert_scale(ysim, y_scale, eps)
    
            yerr_list.append(d)
            if w_list:
                fmin += w_list[0][i] * d * d
                wtot += w_list[0][i]
            else:
                fmin += d * d
                wtot += 1.0

        self.yerr2_sum = fmin
        fmin = sqrt(fmin / wtot)

        self.yc_list = yc_list
        self.yerr_list = yerr_list

        return fmin

    def cal_total_fmin(self, xk_all, x_list, y_list, w_list, run):
        p_tot = self.cal_penalty(xk_all)
        fmin  = self.cal_fmin(xk_all, x_list, y_list, w_list, run)
        if fmin is None: return None

        fmin  = fmin + p_tot

        return fmin

    def stop(self, stopfile):
        if self.stop_flag: return True
        if hasattr(self, "plot_event"):
            if self.plot_event.stop_button.status == 'stop': return True

        exist_stopfile = os.path.isfile(stopfile)
#        if exist_stopfile:
#            print(f"stop file [{stopfile}] found. terminating...")
        if exist_stopfile:
            print()
            print(f"Message: Found stop file [{stopfile}] for early stop.")
            print()
            self.stop_flag = True
            return True

        return False

    def normalized_difference(self, xk0, xk1, optid, defval = None):
        if xk0 is None or xk1 is None: return defval

        xdiff = 0.0
        npk = 0
        for x0, x1, id in zip(xk0, xk1, optid):
            if id == 0: continue
            
            maxabs = max([abs(x0), abs(x1)])
            if maxabs == 0.0: continue

            xdiff += pow(abs(x0 - x1) / maxabs, 2)
            npk += 1

        if npk == 0: return defval

        return sqrt(xdiff) / npk

    def _minimize_func(self, xk, x_list = None, y_list = None, w_list = None, run = True, print_level = 1):
        _global.fit = self
        _global.x_list = x_list
        _global.y_list = y_list
        _global.w_list = w_list
        _global.run = run
        _global.print_level = print_level

        return minimize_func_cachable(tuple(xk))

    def mutate(self, individual, bounds, mutation_rate = 0.01):
        for i in range(len(individual)):
            if random.random() < mutation_rate:
                individual[i] = random.uniform(bounds[i][0], bounds[i][1])

    def crossover(self, parent1, parent2):
        crossover_point = random.randint(1, len(parent1) - 1)
        child1 = np.concatenate((parent1[:crossover_point], parent2[crossover_point:]))
        child2 = np.concatenate((parent2[:crossover_point], parent1[crossover_point:]))
        return list(child1), list(child2)

    def choose_random(self, population):
        n = len(population) // 2
        i = random.randint(0, n-1)
        return population[i]

    def get_model(self, surrogate_model, n_restarts_optimizer, alpha):
        if surrogate_model == 'gp':
#           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-2, 1.0e2))
            model = GaussianProcessRegressor(kernel = kernel, n_restarts_optimizer = n_restarts_optimizer, alpha = alpha)
        else:
            print(f"\nError in tkFit_mxy_flex.get_model(): Invalid surrogate model [{surrogate_model}]")
            input("Press ENTER to terminate>>")
            exit()
            
        return model

    def surrogate_model_func(self, x, rettype = 'mean', print_level = 1):
        self.icall_model += 1

#        pred, sigma = self.model.predict([x], return_std = True)
        x_std = self.mlr.xscaler.transform(pd.DataFrame([x]))

        df = self.mlr.predict(x_std, columns = ['predicted'])
        pred_org = self.mlr.yscaler.inverse_transform(df)
        pred_org = pred_org[0]
        if type(pred_org) is not float: pred_org = pred_org[0]

        sigma = self.mlr.std
        if sigma is not None: 
            sigma = sigma[0] * self.mlr.scaler.scale_
            if type(sigma) is not float: sigma = sigma[0]

        self.sm_func_mean = pred_org
        self.sm_func_std  = sigma

        if sigma is None: return pred_org

        if rettype == 'mean'  : return pred_org
        if rettype == 'min'   : return pred_org - sigma
        if rettype == 'max'   : return pred_org + sigma
        if rettype == 'sigma' : return sigma
        if rettype == '-sigma': return -sigma

        print(f"Error in tkFit_mxy_flex.surrogate_model_func(): Invalid return type = {rettype}")

    def surrogate_optimization(self, X_sample, Y_sample, x0, pk, optid, dx, bounds,
                minimize_func = None, callback = None,
                tol = 1.0e-5, xatol = 1.0e-4, nmaxiter = 1000,
                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):

        self.mlr = tkMLR()
        cfg = self.cfg
        self.icall_model = 0

        hidden_layer_sizes = cfg.get("hidden_layer_sizes", (3, 3))
        nmax_surrogate_trial = cfg.get("nmax_surrogate_trial", 5)
        print(f"  surrogate_model   : {surrogate_model}")
        print(f"  surrogate_method  : {surrogate_method}")
        print(f"  surrogate_scores  : {surrogate_scores}")
        print(f"  nmaxiter_surrogate: {nmaxiter_surrogate}")
        print(f"  tol_surrogate     : {tol_surrogate}")
        print(f"  nmax_surrogate_trial: {nmax_surrogate_trial}")

        self.minimize_func_exact = self.minimize_func
        self.minimize_func = lambda x, run = True: self.surrogate_model_func(x, score)

        x_prev = None
        self.model = None
        i_surrogate_trial = 0
        
        ymin = min(Y_sample)
        idx_min = Y_sample.index(ymin)
        x0 = X_sample[idx_min]
#        ymin_list = [ymin]
#        xmin_list = [x0]
        for i in range(nmaxiter):
            self.iter_surrogate = i

            print()
            print("global iter:", i)

            if ":test" in cfg.method:
                print()
                print(f"Enter test mode (method={cfg.method}):")
                print(f"  surrogate_model   : {surrogate_model}")
                print(f"  surrogate_method  : {surrogate_method}")
                print(f"  surrogate_scores  : {surrogate_scores}")
                print(f"  nmaxdata          : {nmaxdata}")
                print(f"  n_restarts_optimizer: {n_restarts_optimizer}")
                print(f"  test_data_ratio   : {cfg.test_data_ratio}")
                print(f"  hidden_layer_sizes: {cfg.hidden_layer_sizes}")
                print(f"  alpha         : {alpha}")
                print(f"  tol           : {tol}")
                print(f"  xatol         : {xatol}")

            self.descriptor_df = pd.DataFrame(X_sample)
            self.objective_df  = pd.DataFrame(Y_sample)
            self.descriptor_df_std = self.mlr.standardize(self.descriptor_df)
            self.mlr.xscaler = self.mlr.scaler
            self.objective_df_std  = self.mlr.standardize(self.objective_df)
            self.mlr.yscaler = self.mlr.scaler

            if ":test" in cfg.method:
                print()
                print("  Standardize:")
                print("descriptor data (standardized)")
                print("  ", self.descriptor_df_std)
                print("objective data (standardized)")
                print("  ", self.objective_df_std)
                self.objective_df_org = self.mlr.yscaler.inverse_transform(self.objective_df_std)
                print("objective data (inverse standardized)")
                print("  ", self.objective_df_org)

            print("test_data_ratio=", cfg.test_data_ratio)
            print("hidden_layer_sizes=", cfg.hidden_layer_sizes)
            d_train, d_test, o_train, o_test = self.mlr.split_data(self.descriptor_df_std, self.objective_df_std, test_size = cfg.test_data_ratio)

            max_iter_mlp = cfg.get("max_iter_mlp", 1000)
            self.mlr.fit(
                    descriptor_df = d_train,
                    objective_df  = o_train,
                    method = surrogate_model, model = None,
                    n_restarts_optimizer = n_restarts_optimizer,
                    hidden_layer_sizes = hidden_layer_sizes, max_iter_mlp = max_iter_mlp,
                    alpha = cfg.alpha, l1_ratio = cfg.l1_ratio, 
                    use_cv = cfg.use_cv, cv = cfg.cv,
                    max_num_probes = cfg.max_num_probes, interval = cfg.interval,
                    num_rand_basis = cfg.num_rand_basis, score_mode = cfg.score_mode, seed = cfg.seed,
                    max_iter = 20, tol = cfg.tol, print_level = 1)

            if ":test" in cfg.method:
                use_cv = self.mlr.cv_available(cfg.method)

                print()
                if use_cv:
                    log_alpha_min = -10
                    log_alpha_max = 4
                    print("Cross validation:")
                else:
                    print("Calculate scores:")

                predict_train, predict_test = self.mlr.fit_cv_alpha(
                    d_train, o_train, d_test, o_test,
                    method = surrogate_model, surrogate_model = "", surrogate_method = '', 
                    surrogate_scores = "mean|max", nmaxdata = 20,
                    nmaxiter_surrogate = 10, tol_surrogate = 1.0e-5,
                    n_restarts_optimizer = n_restarts_optimizer, 
                    alpha = cfg.alpha, l1_ratio = cfg.l1_ratio,
                    hidden_layer_sizes = hidden_layer_sizes, max_iter_mlp = max_iter_mlp,
                    max_num_probes = cfg.max_num_probes, interval = cfg.interval, num_rand_basis = cfg.num_rand_basis, 
                    score_mode = cfg.score_mode, seed = cfg.seed,
                    log_alpha_min = log_alpha_min, log_alpha_max = log_alpha_max,
                    max_iter = 20, tol = cfg.tol, print_level = 1)

                plt.clf()
                plt.plot(o_train, predict_train, label = "train", linestyle = '', marker = 'o')
                plt.plot(o_test,  predict_test, label = "test", linestyle = '', marker = 'x')
                plt.pause(0.01)
                self.app.terminate("", pause = True)

            updated = False
            status = ''
            xstr = format_strlist(x0, "{:.6g}")
            print(f"  Predict by x0_std and sm: {xstr}, fsm: {self.surrogate_model_func(x0)}")
#            nmaxiter_surrogate = 1000
            options = {'maxiter': nmaxiter_surrogate, "disp": True, 'return_all': True}

            exact_list = []
            for score in surrogate_scores.split('|'):
                if surrogate_method == 'L-BFGS-B':
                    res = minimize(lambda x: self.surrogate_model_func(x, score), 
                            callback = lambda x: self.callback(x, use_so = True, print_level = 0), 
                            x0 = x0, bounds = bounds, method = surrogate_method,
                            tol = tol_surrogate, options = options)
                else:
                    res = minimize(lambda x: self.surrogate_model_func(x, score), 
                            callback = lambda x: self.callback(x, use_so = True, print_level = 0), 
                            x0 = x0, method = surrogate_method,
                            tol = tol_surrogate, options = options)

                xfin = list(res.x)
                print(f"*** surrogate optimization #{i} finished.")
                xstr = format_strlist(xfin, "{:.6g}")
                print(f"  evaluate eaxct minimize function:")
                exact = self.minimize_func_exact(xfin)
                print(f"    locally optimized by {score:4}: x = {xstr}, fsm(x) = {res.fun:.6g}, exact = {exact:.6g}")

                if xfin not in X_sample:
                    X_sample.append(xfin)
                    Y_sample.append(exact)
                    updated = True

                exact_list.append(exact)

            print()
            print("###################################")
            print(f"Surrogate itration #{i} finished:")

            if not updated:
#                status = "  Surrogate X was not updated. Terminate global iteration."
#                break
                status = "  Surrogate X was not updated. Add random data."
                xlist = self.ga_sampling(cfg.bounds_init, 
                            cfg.nparents, cfg.nmaxdata, cfg.ga_nmaxiter, cfg.mutation_rate, print_level = 0)
                X_sample.extend(xlist)
                Y_sample = self.fill_Y_sample(X_sample, Y_sample, self.minimize_func_exact)

            Y_sample, X_sample = sort_lists([Y_sample, X_sample])
            X_sample = list(X_sample)
            Y_sample = list(Y_sample)

# extract minimum nmaxdata
            if len(X_sample) > nmaxdata:
                print(f"    # of data exceeds {nmaxdata}. Extract minimum {nmaxdata} data.")
                X_sample = X_sample[:nmaxdata]
                Y_sample = Y_sample[:nmaxdata]

# minimum data
            xmin = X_sample[0]
            ymin = Y_sample[0]

            if x_prev is None:
                max_dx = None
            else:
                max_dx = -1.0e100
                for _xfin, _xprev in zip(xfin, x_prev):
                    dx = abs(_xfin - _xprev)
                    if dx > max_dx:
                        max_dx = dx
                        break

            xstr = format_strlist(xmin, "{:.6g}")
            print(f"    Surrogate minimum result: X={xstr}  Ymin(exact)={ymin:.6g}")
            xstr = format_strlist(xfin, "{:.6g}")
            if max_dx is None:
                print(f"    max_dx={max_dx}   x_prev={xstr}")
            else:
                print(f"    max_dx={max_dx:.6g}   x_prev={xstr}")

            self.iter_list_exact.append(self.iter)
            self.fmin_list_exact.append(ymin)
            if self.error_data_exact:
                self.error_data_exact[0].set_data(self.iter_list_exact, self.fmin_list_exact)

            if ymin <= tol:
                status = f"Reach the tolerance (Ymin={ymin} <= tol={tol})"
                print(status)
                break

            if max_dx is not None and max_dx < xatol:
                i_surrogate_trial += 1
                print(f"Reach the tolerance for X (max_dx={max_dx:.6g} < xatol={xatol:.4g})")
                print(f"  i_surrogate_trial={i_surrogate_trial} < nmax_surrogate_trial={nmax_surrogate_trial}: continue")
                if i_surrogate_trial > nmax_surrogate_trial:
                    status = f"Reach the tolerance for X (max_dx={max_dx:.6g} < xatol={xatol:.4g})"
                    print(status)
                    break
            else:
                i_surrogate_trial = 0

            '''
            xmin_list.append(xmin)
            ymin_list.append(ymin)
            print(f"    i_surrogate_trial={i_surrogate_trial}")
            print("xmin_list:")
            for _i, _xmin in enumerate(xmin_list):
                print(f"  {_i}: ", format_strlist(_xmin, "{:.6g}"))
            print(f"ymin_list: ", ymin_list)
            '''

#            x_prev = xfin
#            x0 = xfin.copy()
            x_prev = xmin
            x0 = xmin.copy()
            
            self.values = x0
            self.save_parameter_files(cfg.parameterfile, self.fitfile_path, cfg, print_level = 1)

        print()
        print("##########################################")
        print("Surrogate global iteration finished")
        print("status:", status)
        print("##########################################")
        
        ymin = min(Y_sample)
        idx_min = Y_sample.index(ymin)
        xmin = X_sample[idx_min]

        res.x   = xmin
        res.fun = ymin
        res.iter_surrogate = i
        res.status_surrogate = status
        res.iter_callback = self.iter
        res.icall_exact = self.icall
        res.icall_model = self.icall_model

        return res

    def random_sampling(self, bounds, print_level = 1):
        if print_level:
            print()
            print(f"Random sampling in the range:")
            for b in bounds:
                print("  ", b)

        return [random.uniform(low, high) for low, high in bounds]

    def generate_by_ga(self, population, bounds, ga_nmaxiter, nmaxdata = 20, mutation_rate = 0.01):
        population_size = len(population)

        for generation in range(ga_nmaxiter):
            while True:
                c = len(population)
                if c >= nmaxdata: break

                nmax_generation = len(population) // 2
#                if c > nmax_generation: break

                parent1 = self.choose_random(population)
                parent2 = self.choose_random(population)
                child1, child2 = self.crossover(parent1, parent2)
                self.mutate(child1, bounds, mutation_rate)
                self.mutate(child2, bounds, mutation_rate)
                if child1 not in population: population.append(child1)
                if child2 not in population: population.append(child2)
            
                c += 1

        return population

    def ga_sampling(self, bounds, nparents = 5, nmaxdata = 20, ga_nmaxiter = 10, mutation_rate = 0.01, print_level = 1):
        if print_level:
            print()
            print(f"Genetic algorism sampling in the range [{bounds}]: nparents={nparents} nmaxdata={nmaxdata} mutation_rate={mutation_rate}")

        xlist = [self.random_sampling(bounds, print_level = not i) for i in range(nparents)] #min([nparents, nmaxdata]))]
        xlist = self.generate_by_ga(xlist, bounds, ga_nmaxiter, nmaxdata, mutation_rate)
        return xlist

    def history_file_sampling(self, infile, print_level = 1):
        if print_level:
            print()
            print(f"History file sampling from [{infile}]")

        mlr = tkMLR()
        mlr.read_data(infile)
        print("Original data:")
        print(f"  Number of columns: {mlr.ncolumns}")
        print(f"  Number of data   : {mlr.ndata}")
#        print("df:", mlr.df)

        exclude_labels    = mlr.get_exclude_labels(keys = self.varname, id = self.optid, regexp_list = [r"^\-", r"^'\-"])
        objective_labels  = mlr.get_objective_labels(regexp_list = [r"^o:", r"^t:", r"^min:", r"^max:"])
        descriptor_labels = mlr.get_descriptor_labels(labels = mlr.columns, 
                            remove_keys = exclude_labels + objective_labels)
        print("labels:")
        print("  descriptors   :", descriptor_labels)
        print("  objective vars:", objective_labels)
        print("  exclude       :", exclude_labels)
        descriptors = mlr.df[descriptor_labels].to_numpy().tolist()
        objectives  = mlr.df[objective_labels].to_numpy().T.tolist()[0]
#        print("descriptors=", descriptors)
#        print("objectives=", objectives)

        return descriptors, objectives

    def fill_Y_sample(self, X_sample, Y_sample, func = None):
        if func is None: func = lambda x, print_level = 0: self.minimize_func(x, print_level = 0)

        print(f"fill Y_sample for {len(X_sample)} data")
        for i, x in enumerate(tqdm.tqdm(X_sample)):
            if i >= len(Y_sample) or Y_sample[i] is None:
                y = func(x, print_level = 0)
                if i >= len(Y_sample):
                    Y_sample.append(y)
                else:
                    Y_sample[i] = y
        return Y_sample

    def sampling(self, method = "ga", bounds = None, 
                infile = None,
                nparents = 5, nmaxdata = 20, ga_nmaxiter = 10, mutation_rate = 0.01, 
                X_sample = None, Y_sample = None, print_level = 1):
        if X_sample is None: X_sample = []
        if Y_sample is None: Y_sample = []

        if method == 'random':
            xlist = self.random_sampling(bounds, print_level = print_level)
            X_sample.extend(xlist)
            Y_sample = self.fill_Y_sample(X_sample, Y_sample)
        elif method == 'file':
            xlist, ylist = self.history_file_sampling(infile, print_level = print_level)
            Y_sample = self.fill_Y_sample(X_sample, Y_sample)
            X_sample.extend(xlist)
            Y_sample.extend(ylist)
        elif method == 'ga':
            xlist = self.ga_sampling(bounds, nparents, nmaxdata, ga_nmaxiter, mutation_rate, print_level = print_level)
            X_sample.extend(xlist)
            Y_sample = self.fill_Y_sample(X_sample, Y_sample)
        elif method in ['simplex', 'nelder-mead']:
            xlist = self.build_initial_simplex()
            X_sample.extend(xlist)
            Y_sample = self.fill_Y_sample(X_sample, Y_sample)
        else:
            print(f"Error in tkFit_mxy_flx.sampling(): Invalid sampling method [{method}]")
            self.app.terminate("", pause = True)

        print(f"Extract unique data from the obtained {len(X_sample)} data")
        xlist = []
        ylist = []
        for i, (_x, _y) in enumerate(tqdm.tqdm(zip(X_sample, Y_sample))):
            if _x in X_sample[:i]: continue
            
            xlist.append(_x)
            ylist.append(_y)

        return xlist, ylist

    def minimize(self, method = None, jac = None, tol = None, xatol = None, nmaxiter = None, initial_simplex = None, print_level = True):
        cfg = self.cfg

        if method is not None: self.method = method
        if jac is None: jac = self.get('jac', '3-point')
        if tol is None: tol = self.tol
        if xatol is None: xatol = self.xatol
        if nmaxiter is None: nmaxiter = self.nmaxiter

        if not self.hasattr('icall'): self.icall = 1
        if not self.hasattr('nmaxcall'): self.nmaxcall = -1
        if not self.hasattr('fmax_record'): self.fmax_record = 1.0e300

        print("")
        print("Optimizing parameters:")
        optpk = self.extract_parameters()
        print("   method  =", self.method)
        print("   nmaxiter=", nmaxiter)
        print("   tol     =", tol)
        print("   xatol   =", xatol)
        print("   optpk=", optpk)
        print("")
        print("Start minimization:")
        if self.method == 'tkcg':
            '''
            opt = tkOptimize(
                        callback = self.callback, ycalfunc = self.cal_ylist, func = self.minimize_func,
                        diff1func = jac,
                        method = 'cg', nmaxiter = nmaxiter, tolf = tol, tolx = xatol,
                        lsmode = 'wolfe', ls_nmaxiter = 100, ls_alpha = 1.0, ls_h = 1.0e-3, 
                        ls_xrange = 0.5, ls_dump = 0.3, 
                        ls_alphaeps = 1.0e-5, ls_feps = 1.0e-5,
                        print_level = 4, iprintinterval = 10)
            '''

            def diff1func(i, xk_all, optdata = None):
                h = 1.0e-5
                xh = xk_all.copy()
                xh[i] += h
                return (self.minimize_func(xh) - self.minimize_func(xk_all)) / h

            opt = tkOptimize()
            for i in range(len(self.pk)):
                opt.add_parameter(self.varname[i], self.pk[i], self.dx[i], self.optid[i])
            opt.set_method('cg', 'armijo', initial_simplex = None)
            opt.set_functions(fitfunc = self.cal_ylist, func = self.minimize_func, 
                        diff1func = diff1func, diff2func = None)
            opt.initialize(callback = self.callback)
            var = opt.make_optdata(tsleep = 0.0)

            ci, S2, optdata = opt.optimize()
            res = tkParams(x = ci, fun = S2, iter = optdata.iter)
        elif self.method == 'pso':   # particle swarm optimization
            nparents   = self.get('nparents', 20)
            inertia    = self.get('inertia', 0.5)
            cognitive  = self.get('cognitive', 1.5)
            social     = self.get('social', 1.5)
            res = particle_swarm_optimization(self.minimize_func, 
                                    x0 = optpk, pk = self.pk, optid = self.optid,
                                    dx = self.dx, kmin = self.kmin, kmax = self.kmax,
                                    callback = lambda xk: self.callback(xk), 
                                    nparticles = nparents,
                                    inertia = inertia, cognitive = cognitive, social = social,
                                    tol = tol, nmaxiter = nmaxiter)
        elif method == 'remc':
            nREMC_T = self.get('ntempeartures', 5)
            REMC_Tmin = self.get('REMC_Tmin', 1)
            REMC_Tmax = self.get('REMC_Tmax', nREMC_T)
            res = replica_exchange_monte_carlo(self.minimize_func, 
                                    x0 = optpk, pk = self.pk, optid = self.optid,
                                    dx = self.dx, kmin = self.kmin, kmax = self.kmax,
                                    callback = lambda xk: self.callback(xk), 
                                    ntemperatures = nREMC_T, REMC_Tmin = REMC_Tmin, REMC_Tmax = REMC_Tmax,
                                    tol = tol, nmaxiter = nmaxiter)
        elif method == 'ga':   # generic algorism
            nparents = self.get('nparents', 5)
            mutation_rate = self.get('mutation_rate', 0.01)
            res = genetic_algorithm(self.minimize_func, 
                                    x0 = optpk, pk = self.pk, optid = self.optid,
                                    dx = self.dx, kmin = self.kmin, kmax = self.kmax,
                                    callback = lambda xk: self.callback(xk), 
                                    population_size = nparents, mutation_rate = mutation_rate,
                                    tol = tol, nmaxiter = nmaxiter)
        elif 'sa:' in self.method:   # simulated annealing
            _aa = self.method.split(':')
            if len(_aa) == 1:
                method2 = 'rand'
            else:
                method2 = _aa[1]

            T_sa = self.get('T_sa', None)
            cooling_rate = self.get('cooling_rate', None)
            step_size = self.get('step_size', 1.0)

            if method2 == 'nelder-mead':
                tol2 = self.get('tol2', 0.1)
                nmaxiter2 = 1000
            else:
                tol2 = self.get('tol2', 0.1)
                nmaxiter2 = self.get('nmaxiter2', 10)
            
            res = simulated_annealing(self.minimize_func, 
                            x0 = optpk, pk = self.pk, optid = self.optid,
                            dx = self.dx, kmin = self.kmin, kmax = self.kmax,
                            T = T_sa, cooling_rate = cooling_rate, step_size = step_size,
                            method2 = method2, tol2 = tol2, nmaxiter2 = nmaxiter2,
                            callback = lambda xk: self.callback(xk), 
                            tol = tol, nmaxiter = nmaxiter)
        elif self.method == 'nelder-mead':
            if initial_simplex and len(initial_simplex) > 0:
                options = {'maxiter': nmaxiter, "disp": True, "initial_simplex": initial_simplex, 'return_all': True}
            else:
                options = {'maxiter': nmaxiter, "disp": True, 'return_all': True}
            if xatol is not None: options["xatol"] = xatol 

            res = minimize(self.minimize_func, optpk, 
                    method = self.method, 
                    callback = lambda xk: self.callback(xk),
                    tol = tol, options = options)
                    
        elif len(self.method) >= 2 and self.method[:2] == 'so':  # surrogate optimization
            print("Surrogate model")
            initialize_method    = cfg.get('initialize_method', '')
            surrogate_model      = cfg.get('surrogate_model'  , '')
            surrogate_method     = cfg.get('surrogate_method' , '')

            options = self.method.split(':')
            nopt = len(options)
            if nopt > 1: initialize_method = options[1]
            if initialize_method == '': initialize_method = 'ga'
            if nopt > 2: surrogate_model = options[2]
            if surrogate_model == '': surrogate_model = 'gpr'
            if nopt > 3: surrogate_method = options[3]
            if surrogate_method == '': surrogate_method = 'cg'

            nmaxiter_surrogate   = cfg.get('nmaxiter_surrogate', 10)
            tol_surrogate        = cfg.get('tol_surrogate', 1.0e-5)
            surrogate_scores     = cfg.get('surrogate_scores', "mean|max")
            cfg.nparents         = cfg.get('nparents', 5)
            cfg.ga_nmaxiter      = cfg.get('ga_nmaxiter', 3)
            cfg.mutation_rate    = cfg.get('mutation_rate', 0.01) 
            cfg.nmaxdata         = cfg.get('nmaxdata', 20)
            xatol                = cfg.get('xatol', 1.0e-4)
            mutation_rate        = cfg.get('mutation_rate', 0.01)
            n_restarts_optimizer = cfg.get('n_restarts_optimizer', 50)
            alpha                = cfg.get('alpha', 1.0e-4)

            print(f"  initialize_method : {initialize_method}")
            print(f"  surrogate_model   : {surrogate_model}")
            print(f"  surrogate_method  : {surrogate_method}")
            print(f"  nmaxiter_surrogate: {nmaxiter_surrogate}")
            print(f"  tol_surrogate     : {tol_surrogate}")
            print(f"  surrogate_scores  : {surrogate_scores}")
            print(f"  nparents          : {cfg.nparents}")
            print(f"  mutation_rate     : {cfg.mutation_rate}")
            print(f"  ga_nmaxiter       : {cfg.ga_nmaxiter}")
            print(f"  nmaxdata          : {cfg.nmaxdata}")
            print(f"  n_restarts_optimizer: {n_restarts_optimizer}")
            print(f"  alpha             : {alpha}")
            print(f"  tol               : {tol}")
            print(f"  xatol             : {xatol}")

            bounds = []
            for i in range(len(self.kmin)):
                if self.optid[i] == 0: continue
                bounds.append([self.kmin[i], self.kmax[i]])

            bounds_init = []
            for i in range(len(self.kmin)):
                if self.optid[i] == 0: continue
                
                min_bound = max([self.pk[i] - self.dx[i], self.kmin[i]])
                max_bound = min([self.pk[i] + self.dx[i], self.kmax[i]])
                bounds_init.append([min_bound, max_bound])
            cfg.bounds = bounds
            cfg.bounds_init = bounds_init

            print()
            print("initial candidates:", optpk)
            print("bounds:")
            for b in bounds:
                print(f"  ", b)
            print("bounds_init:")
            for b in bounds_init:
                print(f"  ", b)

            X_sample = [optpk.copy()]
            X_sample, Y_sample = self.sampling(initialize_method, bounds = bounds_init, infile = cfg.trainfile,
                                    nparents = cfg.nparents, nmaxdata = cfg.nmaxdata, ga_nmaxiter = cfg.ga_nmaxiter, 
                                    mutation_rate = cfg.mutation_rate,
                                    X_sample = X_sample, Y_sample = None)
            print("X_sample(initial):")
            for i, (_x, _y) in enumerate(zip(X_sample, Y_sample)):
                xstr = format_strlist(_x, "{:12.6g}")
                print(f"  {i:03d}: ({xstr}): {_y:12.6g}")

            Y_sample, X_sample = sort_lists([Y_sample, X_sample])
            X_sample = list(X_sample)
            Y_sample = list(Y_sample)
            if len(X_sample) > cfg.nmaxdata:
                print(f"    # of data exceeds {cfg.nmaxdata}. Extract minimum {cfg.nmaxdata} data.")
                X_sample = X_sample[:cfg.nmaxdata]
                Y_sample = Y_sample[:cfg.nmaxdata]

            print()
            print(f"Initial values extracted for {cfg.nmaxdata} data:")
            for i, (_x, _y) in enumerate(zip(X_sample, Y_sample)):
                xstr = format_strlist(_x, "{:12.6g}")
                print(f"  {i:03d}: ({xstr}): y={_y:12.6g}")

            res = self.surrogate_optimization(X_sample, Y_sample, x0 = optpk, pk = self.pk, optid = self.optid, 
                            dx = self.dx, bounds = bounds,
                            minimize_func = self.minimize_func, 
                            callback = lambda xk: self.callback(xk),
                            tol = tol, xatol = xatol, nmaxiter = nmaxiter,
                            surrogate_model = surrogate_model, surrogate_method = surrogate_method,
                            nmaxiter_surrogate = nmaxiter_surrogate, tol_surrogate = tol_surrogate,
                            surrogate_scores = surrogate_scores, nmaxdata = cfg.nmaxdata,
                            n_restarts_optimizer = n_restarts_optimizer, alpha = alpha)
        else:
            res = minimize(self.minimize_func, optpk, 
                    method = self.method, jac = jac, 
                    callback = self.callback,
                    tol = tol, options = {'maxiter': nmaxiter, "disp": True})

        self.finalpk = self.recover_parameters(res.x)
        self.iter = res.nit
        if self.get('fun', None) is not None:
            self.ffin = res.fun
        else:
            self.ffin = self.minimize_func(res.x)

        pk_ret = self.get("pk_corrected", None)
        if pk_ret is None:
            pk_ret = self.finalpk

        if self.best_fmin is not None:
            if self.best_fmin != res.fun:
                print(f"\n*** Warning: Return values of minimize() conflicts.")
                print(f"        res.fun={res.fun}   fit.best_fmin={self.best_fmin}")
                print(f"        res.x={res.x}   fit.best_xk_all={self.best_xk_all}")
                print(f"   Use fit.best values")
                res.fun = self.best_fmin
                res.x   = self.best_xk_all
                res.success = self.success

        return pk_ret, res.fun, res.success, res
#        return self.finalpk, res.fun, res.success, res

    def _callback(self, pk, run = True, savefig_path = None, use_so = False, print_level = 1):
        cfg = self.app.cfg
        if self.stop(cfg.stopfile): return False

        fplot = self.get('fplot', 1)
        fdisconnect_backward_data = self.fdisconnect_backward_data #if self.hasattr("fdisconnect_backward_data") else False

        if fplot >= 0 and hasattr(self, "plt") and hasattr(self.plt, "get_current_fig_manager"):
            w = self.plt.get_current_fig_manager().window
            if hasattr(self, "window") and w != self.window:
                print()
                self.last_message = f"Message in tkFit_mxy_flex._callback(): Graph window is closed for early stop"
                print(f"{self.last_message}.")
                print()
                self.stop_flag = True
                return False

        xk_all = self.recover_parameters(pk, set_member = False)
        n = len(xk_all)

        xdiff = self.normalized_difference(xk_all, self.xk_all_prev, self.optid, 0.0)
        xatol = 1.0e-4 if not self.hasattr('xatol') else self.xatol

        if xatol > 0.0 and xdiff > 0.0 and xdiff < xatol:
            self.last_message = f"Message in tkFit_mxy_flex._callback(): xdiff={xdiff:.6g} becomes smaller than xatol={xatol}"
            print(f"\n{self.last_message}.")
            print(f"   Terminating iteration\n")
            self.stop_flag = True
            return False

        self.yerr_list = None
        self.yc_list = None
        self.success = False
        fmin = None
#        if self.iter % self.print_interval == 0:
        if True:
            print()
            if self.iter_surrogate is not None:
                print(f"tkFit_mxy_flex._callback(): iter(so): {self.iter_surrogate}  iter: {self.iter}   icall: {self.icall}")
            else:
                print(f"tkFit_mxy_flex._callback(): iter: {self.iter}   icall: {self.icall}")

            if print_level:
                for i in range(n):
                    print(f"  {self.varname[i]:10}: {xk_all[i]:10.4g} {self.unit[i]}")

#            print("  X=", pk)
            fmin = self.minimize_func(pk)
            print(f"  tkFit_mxy_flex._callback(): f={fmin:12.6g}   xdiff={xdiff:12.6g}")

            self.iter_list.append(self.iter + 1)
            self.fmin_list.append(fmin)

            if xdiff == 0.0:
                self.xdiff_list.append(xatol)
            else:
                self.xdiff_list.append(xdiff)

            ycal_list = self.yc_list
            if self.is_fitting and not use_so and ycal_list is None:
                print("\nWarning in callback(): Failed to calculate ycal_list: continue")
                self.success = False
                return 1.0e10

            if self.sm_func_std is not None:
                self.std_iter_list.append(self.iter + 1)
                self.std_list.append(self.sm_func_std)

        if True:
#        if self.iter % self.plot_interval == 0:
#            if not use_so and self.yc_list is None:
#                fmin = self.minimize_func(pk)
#                print(f"  tkFit_mxy_flex._callback(): f={fmin:12.6g}   xdiff={xdiff:12.6g}")

            if fplot >= 0:
                x_list   = self.x_list
                x_labels = self.x_labels
                nx = len(x_list)
                if not use_so:
                    if len(self.data_axes) >= 1 and (nx == 1 or type(self.ix_plot) is int):
                        if type(self.ix_plot) is int and self.ix_plot < nx:
                            ix_plot = self.ix_plot
                        else:
                            ix_plot = 0
                        for i in range(len(self.y_list)):
                            axis = self.data_axes[i]
                            xlist, ylist = self.disconnect_backward_data(x_list[ix_plot], ycal_list[i], fdisconnect_backward_data)
                            self.fit_data_list[i][0].set_data(xlist, ylist)
                    elif len(self.data_axes) >= 1:
                        for i in range(len(self.y_list)):
                            axis = self.data_axes[i]
                            self.fit_data_list[i][0].set_data(range(len(ycal_list[i])), ycal_list[i])

                if self.fmin_list is not None:
                    self.error_data[0].set_data(self.iter_list, self.fmin_list)
                    self.error_axis.set_xlim([min(self.iter_list), max(self.iter_list)])
                    self.error_axis.set_ylim([min(self.fmin_list), max(self.fmin_list)])

                if self.xdiff_list is not None:
                    self.xdiff_data[0].set_data(self.iter_list, self.xdiff_list)
                    self.xdiff_axis.set_xlim([min(self.iter_list), max(self.iter_list)])
                    self.xdiff_axis.set_ylim([min(self.xdiff_list), max(self.xdiff_list) * 1.2])

                if self.std_axis and self.std_list is not None:
                    self.std_data[0].set_data(self.std_iter_list, self.std_list)
                    self.std_axis.set_xlim([min(self.std_iter_list), max(self.std_iter_list)])
                    self.std_axis.set_ylim([min(self.std_list), max(self.std_list) * 1.2])

            if fplot >= 1:
                for ax in self.data_axes:
                    ax.relim()
                    ax.autoscale_view()

                if savefig_path is not None and savefig_path != '':
                    self.plt.savefig(savefig_path)

                self.plt.tight_layout()
                self.plt.subplots_adjust(top = self.plot_region[0], bottom = self.plot_region[1])
                self.plt.pause(0.01)

        if self.best_fmin is None or self.best_fmin > fmin:
            self.best_fmin   = fmin
            self.best_xk_all = xk_all

        self.iter += 1
        self.xk_all_prev = xk_all.copy()
    
        self.retrieve_parameter_list(xk_all, self.varname, target = cfg, print_level = print_level)
#    self.print_variables(heading = "selfting parameters (callback):")

        self.pk = xk_all.copy()
        self.save_parameter_files(cfg.parameterfile, self.fitfile_path, cfg, print_level = print_level)

        if fmin < self.tol:
            self.stop_flag = True
            if print_level:
                self.last_message = f"Converged: fmin={fmin} reaches tol={self.tol}"
                self.success = True
                print(f"\n** {self.last_message}\n")
                return False

        return True

    def button_click(self, e):
        self.stop_flag = True
        self.plot_event.stop_button.status = 'stop'
        self.plot_event.stop_button.set_text('terminating')
        print("\ntkfit_mxy_flex.button_click(): Stop button pressed. Terminating...\n")

    def initial_plot(self, data_axes, error_axis = None, xdiff_axis = None, std_axis = None,
                        use_so = False,
                        label_input = None, yini_list = 'initial', label_ini = 'input',
                        label_fit = 'fit',
                        ix_plot = '', x_scale = '',
                        label_error = 'error', fmin = None,
                        label_error_exact = 'error(exact)', fmin_exact = None,
                        label_xdiff = 'xdiff', xdiff = None,
                        fontsize = 16, legend_fontsize = None,
                        fplot = None, fdisconnect_backward_data = None,
                        plt = None, fig = None, use_pause = 0.1, 
                        savefig_path = None,
                        button_region = [0.15, 0.95, 0.10, 0.03], plot_region = [0.92, 0.15]):
#                        button_region = [0.15, 0.95, 0.10, 0.03], plot_region = [0.90, 0.10]):
        if fplot is None: fplot = self.get("fplot", 1)
        if fplot < 0: return
        if legend_fontsize is None: legend_fontsize = fontsize
        if fdisconnect_backward_data is None: 
            fdisconnect_backward_data = self.fdisconnect_backward_data if self.hasattr("fdisconnect_backward_data") else False

        self.plt = plt
        self.fig = fig
        self.use_so = use_so
        self.ix_plot = ix_plot
        self.x_scale = x_scale
        self.use_pause     = use_pause
        self.button_region = button_region
        self.plot_region   = plot_region

        self.plot_event = tkPlotEvent(plt)
        self.plot_event.prepare_annotation()

        self.yini_list  = yini_list
        self.data_axes  = data_axes
        self.error_axis = error_axis
        self.xdiff_axis = xdiff_axis
        self.std_axis   = std_axis
        for axis in data_axes:
            axis.tick_params(labelsize = fontsize)

        nx    = len(self.x_list)
        index = range(len(self.x_list[0]))
        if len(data_axes) >= 1 and (nx == 1 or type(ix_plot) is int):
            if type(ix_plot) is int and ix_plot < nx:
                pass
            else:
                ix_plot = 0
            self.fit_data_list = []
            for i in range(len(self.y_list)):
                axis = data_axes[i]
#                self.input_data   = axis.plot(self.x_list[ix_plot], self.y_list[i], label = self.y_labels[i], linestyle = '',  marker = 'o', markersize = 5.0)
#                self.initial_data = axis.plot(self.x_list[ix_plot], yini_list[i],   label = label_ini,       linestyle = '-', linewidth = 0.5, color = 'red')
                line_input,   = axis.plot(self.x_list[ix_plot], self.y_list[i], label = self.y_labels[i], linestyle = '',  
                                          marker = 'o', markersize = 5.0, markerfacecolor = 'w', markeredgecolor = 'black')
                xlist, ylist = self.disconnect_backward_data(self.x_list[ix_plot], yini_list[i], fdisconnect_backward_data)
                line_initial, = axis.plot(xlist, ylist,   label = label_ini,       linestyle = '-', linewidth = 0.5, color = 'red')

                format = "{label}: line#{iline} data#{idata}: x={x:g} y={y:g}"
                self.plot_event.annotation.add_line(f"input y#{i}", axis, axis, self.x_list[ix_plot], self.y_list[i], line_input, 
                        inf_list = {"x": self.x_list[ix_plot], "y": self.y_list[i]},
                        annotation_format = format, inf_format = format)
                self.plot_event.annotation.add_line(f"initial y#{i}", axis, axis, self.x_list[ix_plot], yini_list[i], line_initial, 
                        inf_list = {"x": self.x_list[ix_plot], "y": self.y_list[i]},
                        annotation_format = format, inf_format = format)

                fit_data = axis.plot([], [], label = label_fit, linestyle = '-', linewidth = 0.5, color = 'blue')
                self.fit_data_list.append(fit_data)

            axis.set_xlabel(self.x_labels[ix_plot], fontsize = fontsize)
            axis.set_ylabel(self.y_labels[0], fontsize = fontsize)
            if self.x_scale == 'log':
                axis.set_xscale('log')
            if self.y_scale == 'log':
                axis.set_yscale('log')

            axis.legend(fontsize = legend_fontsize)

#            self.plot_event_input_data   = self.plot_event.add_data({"label": label_input, "plot_type": "2D", "axis": axis, 
#                                        "data": self.input_data})
#            self.plot_event_initial_data = self.plot_event.add_data({"label": label_ini,   "plot_type": "2D", "axis": axis, 
#                                        "data": self.initial_data})
#            self.plot_event_fit_data     = self.plot_event.add_data({"label": label_fit,   "plot_type": "2D", "axis": axis, 
#                                        "data": fit_data})
        elif len(data_axes) >= 1:
            self.fit_data_list = []
            for i in range(len(self.y_list)):
                axis = data_axes[i]
#                self.input_data   = axis.plot(index, self.y_list[i], label = self.y_labels[i], linestyle = '',  marker = 'o', markersize = 5.0)
#                self.initial_data = axis.plot(index, yini_list[i],   label = label_ini,   linestyle = '-', linewidth = 0.5, color = 'red')
                line_input,   = axis.plot(index, self.y_list[i], label = self.y_labels[i], linestyle = '',  
                                            marker = 'o', markersize = 5.0, markerfacecolor = 'w', markeredgecolor = 'black')
                line_initial, = axis.plot(index, yini_list[i],   label = label_ini,   linestyle = '-', linewidth = 0.5, color = 'red')

                format = "{label}: line#{iline} data#{idata}: index={index} x={x:g} y={y:g}"
                self.plot_event.annotation.add_line(f"input", axis, axis, index, self.y_list[i], line_input, 
                        inf_list = {"index": index, "x": self.x_list[i], "y": self.y_list[i], "fit": self},
                        annotation_format = format, inf_format = inf_format_func)
                self.plot_event.annotation.add_line(f"initial", axis, axis, index, yini_list[i], line_initial, 
                        inf_list = {"index": index, "x": self.x_list[i], "y": yini_list[i], "fit": self},
                        annotation_format = format,
                        inf_format = inf_format_func)

                fit_data     = axis.plot([], [], label = label_fit, linestyle = '-', linewidth = 0.5, color = 'blue')
                self.fit_data_list.append(fit_data)
       
            axis.set_xlabel('index', fontsize = fontsize)
            axis.set_ylabel(self.y_labels[0], fontsize = fontsize)
            if self.y_scale == 'log':
                axis.set_yscale('log')
            axis.legend(fontsize = legend_fontsize)
#            self.plot_event_input_data   = self.plot_event.add_data({"label": label_input, "plot_type": "2D", "axis": axis, 
#                                        "data": self.input_data,   "x_list": self.x_list, "x_labels": self.x_labels})
#            self.plot_event_initial_data = self.plot_event.add_data({"label": label_ini,   "plot_type": "2D", "axis": axis, 
#                                        "data": self.initial_data, "x_list": self.x_list, "x_labels": self.x_labels})
#            self.plot_event_fit_data     = self.plot_event.add_data({"label": label_fit,   "plot_type": "2D", "axis": axis, 
#                                        "data": fit_data,     "x_list": self.x_list, "x_labels": self.x_labels})

        if len(data_axes) >= 1:
            self.view_xlim_data_list = []
            self.view_ylim_data_list = []
            for i in range(len(self.y_list)):
                axis = data_axes[i]
                self.view_xlim_data_list.append(axis.get_xlim())
                self.view_ylim_data_list.append(axis.get_ylim())

        if error_axis:
            self.error_axis = error_axis
#            error_axis.tick_params(labelsize = fontsize)
            if fmin is None:
                fmin = self.minimize_func(self.pk)
            self.iter_list = [0]
            self.fmin_list = [fmin]
            self.iter_list_exact = [0]
            self.fmin_list_exact = [fmin]

            self.error_data = error_axis.plot(self.iter_list, self.fmin_list, label = label_error, 
                        linestyle = '-', linewidth = 0.5, color = 'black',
                        marker = 'o', markersize = 5.0, markerfacecolor = 'black', markeredgecolor = 'black')
            if use_so:
                self.error_data_exact = error_axis.plot(self.iter_list_exact, self.fmin_list_exact, label = label_error_exact, 
                        linestyle = '-', linewidth = 0.5, color = 'red',
                        marker = '*', markersize = 10.0, markerfacecolor = 'red', markeredgecolor = 'red')
            else:
                self.error_data_exact = None

            error_axis.set_xlabel('iteration', fontsize = fontsize)
            error_axis.set_ylabel('error', fontsize = fontsize)
            error_axis.set_yscale('log')
#            error_axis.legend(fontsize = legend_fontsize)

            if xdiff_axis:
                self.xdiff_axis = xdiff_axis
#                xdiff_axis.tick_params(labelsize = fontsize)
                if xdiff is None or xdiff == 0.0:
                    xatol = 1.0e-4 if not self.hasattr('xatol') else self.xatol
                    self.xdiff_list = [xatol]
                else:
                    self.xdiff_list = [xdiff]

                self.xdiff_data = xdiff_axis.plot(self.iter_list, self.xdiff_list, label = label_xdiff, 
                        linestyle = '', linewidth = 0.5, color = 'blue',
                        marker = '^', markersize = 5.0, markerfacecolor = 'blue', markeredgecolor = 'blue')
                xdiff_axis.set_ylabel('xdiff', fontsize = fontsize)
                xdiff_axis.set_yscale('log')
#                xdiff_axis.legend(fontsize = legend_fontsize)

            if self.std_axis:
                self.std_iter_list = [0]
                self.std_list = [0.0]
                self.std_data = std_axis.plot(self.std_iter_list, self.std_list, label = "std",
                        linestyle = '', linewidth = 0.5, color = 'blue',
                        marker = 'o', markersize = 5.0, markerfacecolor = 'blue', markeredgecolor = 'blue')
                std_axis.set_xlabel('iteration', fontsize = fontsize)
                std_axis.set_ylabel('std', fontsize = fontsize)
                std_axis.set_yscale('log')

            if xdiff_axis:
                h1, l1 = error_axis.get_legend_handles_labels()
                h2, l2 = xdiff_axis.get_legend_handles_labels()
                error_axis.legend(h1 + h2, l1 + l2, fontsize = legend_fontsize)
            else:
                error_axis.legend(fontsize = legend_fontsize)

#            self.plot_event_error_data = self.plot_event.add_data(
#                        {"label": label_error, "plot_type": "2D", "axis": self.error_axis, "data": self.error_data})


#            self.plot_event_error_data = self.plot_event.add_data(
#                        {"label": label_error, "plot_type": "2D", "axis": self.error_axis, "data": self.error_data})

#        self.plt.tight_layout()
        self.plt.subplots_adjust(top = plot_region[0], bottom = plot_region[1])

#        self.plot_event.add_stop_button()
        self.plot_event.add_stop_button(on_stop_clicked = self.button_click)
#        self.button_region = button_region
#        self.ax_button = plt.axes(button_region)
#        self.stop_button = wg.Button(self.ax_button, 'stop', color = '#f8e58c', hovercolor = '#38b48b')
#        self.stop_button.on_clicked(self.button_click)

        if fplot >= 1:
            for ax in self.data_axes:
                ax.relim()
                ax.autoscale_view()

        if savefig_path is not None and savefig_path != '':
            plt.savefig(savefig_path)

        if fplot >= 1:
            if use_pause is not None:
                self.plt.pause(use_pause)
            else:
                self.plt.show()
        self.window   = plt.get_current_fig_manager().window

    def finalize_plot(self, yfin_list, iter = None, fmin = None, 
                savefig_path = 'final.png', fdisconnect_backward_data = None,
                fplot = None, use_pause = 0.01, button_region = None, plot_region = None):
        if fplot is None: fplot = self.get("fplot", 1)
        if fplot < 0: return

        if button_region is None: button_region = self.button_region
        if plot_region is None: plot_region = self.plot_region
        if fdisconnect_backward_data is None: 
            fdisconnect_backward_data = self.fdisconnect_backward_data if self.hasattr("fdisconnect_backward_data") else False

        x_list   = self.__dict__.get('x_list', None)
        x_labels = self.__dict__.get('x_labels', None)
        nx = len(x_list)
        if yfin_list and (nx == 1 or type(self.ix_plot) is int):
            if type(self.ix_plot) is int and self.ix_plot < nx:
                ix_plot = self.ix_plot
            else:
                ix_plot = 0
            for i in range(len(self.y_list)):
                xlist, ylist = self.disconnect_backward_data(x_list[ix_plot], yfin_list[i], fdisconnect_backward_data)
                self.fit_data_list[i][0].set_data(xlist, ylist)

                axis = self.data_axes[ix_plot]
                format = "{label}: line#{iline} data#{idata}: x={x:g} y={y:g}"
                self.plot_event.annotation.add_line(f"final", axis, axis, x_list[ix_plot], yfin_list[i], self.fit_data_list[i][0], 
                        inf_list = {"x": x_list[ix_plot], "y": yfin_list[i], "fit": self},
                        annotation_format = format, inf_format = inf_format_func)

                self.data_axes[ix_plot].set_xlim([min(self.x_list[0]), max(self.x_list[0])])
                self.data_axes[ix_plot].set_ylim([min([self.view_ylim_data_list[i][0], min(yfin_list[i])]), 
                                 max([self.view_ylim_data_list[i][1], max(yfin_list[i])])])
        elif yfin_list:
            self.index = list(range(len(yfin_list[0])))
            for i in range(len(self.y_list)):
                self.fit_data_list[0][0].set_data(self.index, yfin_list[i])

                axis = self.data_axes[0]
                format = "{label}: line#{iline} data#{idata}: index={index} x={x:g} y={y:g}"
                self.plot_event.annotation.add_line(f"final", axis, axis, x_list[i], yfin_list[i], self.fit_data_list[i][0], 
                        inf_list = {"index": self.index, "x": self.x_list[i], "y": self.y_list[i], "fit": self},
                        annotation_format = format, inf_format = inf_format_func)

                self.data_axes[0].set_xlim([min(self.index), max(self.index)])
                self.data_axes[0].set_ylim([min([self.view_ylim_data_list[i][0], min(yfin_list[i])]), 
                                 max([self.view_ylim_data_list[i][1], max(yfin_list[i])])])

        if self.error_axis:
            format = "error: iter={iteration}  fmin={fmin:g}"
            self.plot_event.annotation.add_line(f"error", self.error_axis, self.error_axis, self.iter_list, self.fmin_list, self.error_data[0], 
                        inf_list = {"iteration": self.iter_list, "fmin": self.fmin_list},
                        annotation_format = format, inf_format = format)

        if self.y_scale == 'log':
            self.data_axes[0].set_yscale('log')

        if iter is None:
            iter = self.iter
        if fmin is None:
            fmin = self.ffin
            self.iter_list.append(iter)
            self.fmin_list.append(ffin)

        if self.get('error_data', None) is not None:
            self.error_data[0].set_data(self.iter_list, self.fmin_list)
            self.error_axis.set_xlim([min(self.iter_list), max(self.iter_list)])
            self.error_axis.set_ylim([min(self.fmin_list), max(self.fmin_list)])

        self.plot_event.stop_button.set_text('finished')
        self.plot_event.register_annotation_event(self.fig, activate = True, print_level = 0)

#        plt.tight_layout()
        self.plt.subplots_adjust(top = plot_region[0], bottom = plot_region[1])

        if savefig_path is not None and savefig_path != '':
            plt.savefig(savefig_path)

        if fplot >= 1:
            if use_pause is not None:
                self.plt.pause(use_pause)
            else:
                self.plt.show()

    def S2_with_x_ivar(self, x, ivar, variance, run = True, print_level = 0):
        pk = self.pk.copy()
        if ivar is not None:
            pk[ivar] = x
        fmin = self.cal_fmin(pk, self.x_list, self.y_list, self.w_list, run = run, print_level = print_level)
        
        return self.yerr2_sum / variance

    def search_xlim_p(self, variance, ivar, xstep, xlim_max, S2_0, S2_target, nmax_iter, S2_eps = 1.0e-5, kstep = 2.0):
        x0 = self.pk[ivar]
        x_min  = x0
        S2_min = S2_0
        while True:
            x_new = x0 + xstep
            if xlim_max <= x_new:
                x_new = xlim_max
            S2_new = self.S2_with_x_ivar(x_new, ivar, variance)
#            print("x_new=", x_new, S2_new, S2_target)
            if S2_target < S2_new or xlim_max == x_new:
                x_max = x_new
                S2_max = S2_new
                break

            x_min = x_new
            S2_min = S2_new
            xstep *= kstep
        
        if abs(S2_max - S2_target) < S2_eps:
            return x_max, S2_max
        else:
            for iter in range(nmax_iter):
                x_half = (x_min + x_max) / 2.0
                S2_half = self.S2_with_x_ivar(x_half, ivar, variance)
                if abs(S2_half - S2_target) < S2_eps:
                    break

                if (S2_min - S2_target) * (S2_half - S2_target) <= 0:
                    x_max = x_half
                    S2_max = S2_half
                else:
                    x_min = x_half
                    S2_min = S2_half
#                print("iter=", x_min, x_half, x_max)

        return x_half, S2_half

    def search_xlim_m(self, variance, ivar, xstep, xlim_min, S2_0, S2_target, nmax_iter, S2_eps = 1.0e-5, kstep = 2.0):
        x0 = self.pk[ivar]
        x_max  = x0
        S2_max = S2_0
        while True:
            x_new = x0 - xstep
            if x_new <= xlim_min:
                x_new = xlim_min
            S2_new = self.S2_with_x_ivar(x_new, ivar, variance)
            if S2_target < S2_new or x_new == xlim_min:
                x_min = x_new
                S2_min = S2_new
                break

            x_max = x_new
            S2_max = S2_new
            xstep *= kstep

        if abs(S2_max - S2_target) < S2_eps:
            return x_min, S2_min
        else:
            for iter in range(nmax_iter):
                x_half = (x_min + x_max) / 2.0
                S2_half = self.S2_with_x_ivar(x_half, ivar, variance)
                if abs(S2_half - S2_target) < S2_eps:
                    break
            
                if (S2_min - S2_target) * (S2_half - S2_target) <= 0:
                    x_max = x_half
                    S2_max = S2_half
                else:
                    x_min = x_half
                    S2_min = S2_half

        return x_half, S2_half

    def build_likelihood_func(self, xk, ivar, 
                x_list = None, y_list = None, w_list = None,
                xlim_min = None, xlim_max = None, kvariance = 6.0, ncal = 101, 
                nmax_iter = 10, S2_eps = 0.1, scale = None, run = True, print_level = 0):
        if x_list is None:
            x_list = self.x_list
        if y_list is None:
            y_list = self.y_list
        if w_list is None:
            w_list = self.w_list
        if xlim_min is None:
            xlim_min = self.kmin[ivar]
        if xlim_max is None:
            xlim_max = self.kmax[ivar]
        if scale is None:
            scale = self.pk_scale[ivar]

        fmin = self.cal_fmin(xk, x_list, y_list, w_list, run = run, print_level = print_level)
        yerr_list = self.yerr_list

        nx = len(yerr_list)
        yerr2_sum0 = sum([d * d for d in yerr_list])
        yerr_mean = np.mean(yerr_list)
        variance  = np.var(yerr_list)

        if print_level:
            print(f"Error statistics for {self.varname[ivar]}:")
            print(f"  # of data   : {nx}")
            print(f"  mean        : {yerr_mean:12.6g}")
            print(f"  variance    : {variance:12.6g}")
            print(f"  S2=sum(error^2): {yerr2_sum0:12.6g}")
#    print(f"  fmin=sqrt({yerr2_sum0:12.6g}/{nx}={fmin:12.6g}")

        x0 = xk[ivar]
        S2_0 = self.S2_with_x_ivar(x0, ivar, variance)
        p0 = exp(-S2_0 / 2.0)
        S2_target = S2_0 + kvariance
#       p0 = exp(-yerr2_sum0 / variance / 2.0)
#        print(f"  p0={p0:12.6g}")
#       print(f"  x0={x0:g}")
        if print_level:
            print(f"  xrange: {xlim_min:g} - {xlim_max:g}")
            print(f"  S2={S2_0:12.6g}")
            print(f"  p0={p0:12.6g}")
            print(f"  S2_target={S2_target:g}")

        x_list = [x0]
        p_list = [1.0]

        xstep = (xlim_max - x0) / 10.0
        kstep = 2.0
        nmaxiter_new = nmax_iter
        if abs(xlim_max * pow(kstep, nmax_iter)) > 1.0e-5:
            print(f"\n*** Warning in tkFit_mxy_flex.build_likelihood_func(): xlim_max={xlim_max} is too large for nmax_iter={nmax_iter}")
            nmaxiter_new = int(abs(log(1.0e-5) * log(abs(xlim_max)) / log(kstep)))
            kstep = nmaxiter_new / 20.0 * 20.0
            nmaxiter_new = int(abs(log(1.0e-5) * log(abs(xlim_max)) / log(kstep)))
            print(f"    nmaxiter is changed to {nmaxiter_new}, kstep is changed to {kstep}\n")
        xmax, S2_max = self.search_xlim_p(variance, ivar, xstep, xlim_max, S2_0, S2_target, nmaxiter_new, S2_eps)
        if print_level:
            print(f"    calculate in xrange: {x0:g} - {xmax:g} at S2={S2_max:g}")

        if scale == 'log':
            xstep = (log(xmax) - log(x0)) / (ncal - 1)
        else:
            xstep = (xmax - x0) / (ncal - 1)
        for i in range(1, ncal):
            if scale == 'log':
                x = x0 * exp(i * xstep)
            else:
                x = x0 + i * xstep
            _S2 = self.S2_with_x_ivar(x, ivar, variance)
            x_list.append(x)
            p_list.append(exp(-_S2 / 2.0) / p0)

        xstep = (x0 - xlim_min) / 10.0
        nmaxiter_new = nmax_iter
        kstep = 2.0
        print("xlim_min=", xlim_min)
        if abs(xlim_min * pow(kstep, nmax_iter)) > 0.01:
            print(f"\n*** Warning in tkFit_mxy_flex.build_likelihood_func(): xlim_min={xlim_min} is too large for nmax_iter={nmax_iter}")
            kstep = 10.0
            nmaxiter_new = int(abs(log(1.0e-5)) * abs(log(abs(xlim_min)) / log(kstep)))
            kstep = nmaxiter_new / 20.0 * 20.0
            nmaxiter_new = int(abs(log(1.0e-5)) * abs(log(abs(xlim_min)) / log(kstep)))
            print(f"    nmaxiter is changed to {nmaxiter_new}, kstep is changed to {kstep}\n")
        xmin, S2_min = self.search_xlim_m(variance, ivar, xstep, xlim_min, S2_0, S2_target, nmaxiter_new, S2_eps, kstep)
        if print_level:
            print(f"    calculate in xrange: {xmin:g} - {x0:g} at S2={S2_min:g}")
        if scale == 'log':
            xstep = (log(x0) - log(xmin)) / (ncal - 1)
        else:
            xstep = (x0 - xmin) / (ncal - 1)
        for i in range(1, ncal):
            if scale == 'log':
                x = x0 * exp(-i * xstep)
            else:
                x = x0 - i * xstep
            _S2 = self.S2_with_x_ivar(x, ivar, variance)
            x_list.append(x)
            p_list.append(exp(-_S2 / 2.0) / p0)

#zipを使って、x_listとp_listをx_listの値準にソートする
        x_list, p_list = zip(*sorted(zip(x_list, p_list)))
#       print("x  y")
#       for x, p in zip(x_list, p_list):
#           print(f"{x:g}\t{p:g}")

        return x_list, p_list

    def get_xrange_list(self, itarget, nx, varname = None, xmin = None, xmax = None, scale = None, print_level = 1):
        if varname is None:
            varname = self.varname[itarget]
        if xmin is None:
            xmin = self.kmin[itarget]
        if xmax is None:
            xmax = self.kmax[itarget]
        if scale is None:
            scale = self.pk_scale[itarget]

        if scale == 'log':
            if print_level:
                print(f"   range {xmin} - {xmax} for {nx} points in {scale} scale")
            logx0 = log(xmin)
            logx1 = log(xmax)
            logxstep = (logx1 - logx0) / (nx - 1)
            xrange_list = [exp(logx0 + i * logxstep) for i in range(nx)]
        else:
            if print_level:
                print(f"   range {xmin} - {xmax} for {nx} points in linear scale")
            xstep = (xmax - xmin) / (nx - 1)
            xrange_list = [xmin + i * xstep for i in range(nx)]

        if xmax < xrange_list[nx - 1]:
            xrange_list[nx - 1] = xmax
            if print_level:
                print(f"    xrange_list for {varname}:", [f"{x:10.4g}" for x in xrange_list])

        return xrange_list
