import collections
import time

import numpy
import scipy.optimize
import scipy.stats
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from matplotlib.colors import BoundaryNorm

from threeML.io.rich_display import display
import numpy as np

from threeML.minimizer import minimization
from threeML.exceptions import custom_exceptions
from threeML.io.table import Table
from threeML.parallel.parallel_client import ParallelClient
from threeML.io.progress_bar import ProgressBar
from threeML.config.config import threeML_config
from threeML.exceptions.custom_exceptions import custom_warnings


class ReducingNumberOfThreads(Warning):
    pass


class ReducingNumberOfSteps(Warning):
    pass


class NotANumberInLikelihood(Warning):
    pass


class JointLikelihood(object):
    def __init__(self, likelihood_model, data_list, verbose=False):
        """
        Implement a joint likelihood analysis.

        :param likelihood_model: the model for the likelihood analysis
        :param data_list: the list of data sets (plugin instances) to be used in this analysis
        :param verbose: (True or False) print every step in the -log likelihood minimization
        :return:
        """

        # Process optional keyword parameters
        self.verbose = verbose

        self._likelihood_model = likelihood_model

        self._data_list = data_list

        for dataset in self._data_list.values():

            dataset.set_model(self._likelihood_model)

        # This is to keep track of the number of calls to the likelihood
        # function
        self.ncalls = 0

        self.set_minimizer('MINUIT')

        # Initial set of free parameters

        self._free_parameters = self._likelihood_model.free_parameters

        # Initially set the value of _current_minimum to None, it will be change by the fit() method

        self._current_minimum = None

        self._minimizer = None

    @property
    def likelihood_model(self):
        """
        :return: likelihood model for this analysis
        """

        return self._likelihood_model

    @property
    def data_list(self):
        """
        :return: data list for this analysis
        """

        return self._data_list

    @property
    def current_minimum(self):
        """
        :return: current minimum of the joint likelihood (available only after the fit() method)
        """
        return self._current_minimum

    @property
    def minimizer(self):
        """
        :return: an instance of the minimizer used in the fit (available only after the fit() method)
        """
        return self._minimizer

    def _update_free_parameters(self):

        """Update the dictionary of free parameters"""

        self._free_parameters = self._likelihood_model.free_parameters

    def fit(self):
        """
        Perform a fit of the current likelihood model on the datasets

        :param pre_fit: (True or False) If True, perform a pre-fit with only normalizations free (experimental)
        :return: a dictionary with the results on the parameters, and the values of the likelihood at the minimum
                 for each dataset and the total one.
        """

        # Update the list of free parameters, to be safe against changes the user might do between
        # the creation of this class and the calling of this method

        self._update_free_parameters()

        # Now check and fix if needed all the deltas of the parameters
        # to 10% of their value (otherwise the fit will be super-slow)

        for k, v in self._free_parameters.iteritems():

            if abs(v.delta) < abs(v.value) * 0.1:

                v.delta = abs(v.value) * 0.1

        # Instance the minimizer

        self._minimizer = self.Minimizer(self.minus_log_like_profile,
                                         self._free_parameters)

        # Perform the fit

        xs, log_likelihood_minimum = self._minimizer.minimize()

        # Store the current minimum for the -log likelihood

        self._current_minimum = float(log_likelihood_minimum)

        # Print results

        print("Best fit values:\n")

        self._minimizer.print_fit_results()

        print("Nuisance parameters:\n")

        nuisance_parameters = collections.OrderedDict()

        for dataset in self._data_list.values():

            for pName in dataset.get_nuisance_parameters():

                nuisance_parameters[(dataset.get_name(), pName)] = dataset.nuisanceParameters[pName]

        nuisance_parameters_table = self._get_table_of_parameters(nuisance_parameters)

        display(nuisance_parameters_table)

        print("\nCorrelation matrix:\n")

        self._minimizer.print_correlation_matrix()

        print("\nMinimum of -logLikelihood is: %s\n" % log_likelihood_minimum)

        print("Contributions to the -logLikelihood at the minimum:\n")

        # Fill the dictionary with the values of the -log likelihood (total, and dataset by dataset)

        minus_log_likelihood_values = collections.OrderedDict()

        minus_log_likelihood_values['total'] = log_likelihood_minimum

        data = []
        name_length = 0

        for dataset in self._data_list.values():

            ml = dataset.inner_fit() * (-1)

            minus_log_likelihood_values[dataset.get_name()] = ml

            name_length = max(name_length, len(dataset.get_name()) + 1)
            data.append([dataset.get_name(), ml])

        log_like_values_table = Table(rows=data, names=["Detector", "-LogL"], dtype=('S%i' % name_length, float))

        log_like_values_table['-LogL'].format = '2.2f'

        display(log_like_values_table)

        # Prepare the dictionary with the results

        results = collections.OrderedDict()

        results['parameters'] = xs
        results['minusLogLike'] = minus_log_likelihood_values

        return results

    def get_errors(self):
        """
        Compute the errors on the parameters using an accurate algorithm.

        :return: a dictionary containing the asymmetric errors for each parameter.
        """

        # Check that the user performed a fit first

        assert self._current_minimum is not None, "You have to run the .fit method before calling errors."

        return self._minimizer.get_errors()

    def get_contours(self, param_1, param_1_minimum, param_1_maximum, param_1_n_steps,
                     param_2=None, param_2_minimum=None, param_2_maximum=None, param_2_n_steps=None,
                     progress=True, **options):
        """
        Generate confidence contours for the given parameters by stepping for the given number of steps between
        the given boundaries. Call it specifying only source_1, param_1, param_1_minimum and param_1_maximum to
        generate the profile of the likelihood for parameter 1. Specify all parameters to obtain instead a 2d
        contour of param_1 vs param_2.

        NOTE: if using parallel computation, param_1_n_steps must be an integer multiple of the number of running
        engines. If that is not the case, the code will reduce the number of steps to match that requirement, and
        issue a warning

        :param param_1: fully qualified name of the first parameter or parameter instance
        :param param_1_minimum: lower bound for the range for the first parameter
        :param param_1_maximum: upper bound for the range for the first parameter
        :param param_1_n_steps: number of steps for the first parameter
        :param param_2: fully qualified name of the second parameter or parameter instance
        :param param_2_minimum: lower bound for the range for the second parameter
        :param param_2_maximum: upper bound for the range for the second parameter
        :param param_2_n_steps: number of steps for the second parameter
        :param progress: (True or False) whether to display progress or not
        :param log: by default the steps are taken linearly. With this optional parameter you can provide a tuple of
                    booleans which specify whether the steps are to be taken logarithmically. For example,
                    'log=(True,False)' specify that the steps for the first parameter are to be taken logarithmically,
                    while they are linear for the second parameter. If you are generating the profile for only one
                    parameter, you can specify 'log=(True,)' or 'log=(False,)' (optional)
        :return: a tuple containing an array corresponding to the steps for the first parameter, an array corresponding
                 to the steps for the second parameter (or None if stepping only in one direction), a matrix of size
                 param_1_steps x param_2_steps containing the value of the function at the corresponding points in the
                 grid. If param_2_steps is None (only one parameter), then this reduces to an array of
                 size param_1_steps.
        """

        if hasattr(param_1,"value"):

            # Substitute with the name
            param_1 = param_1.path

        if hasattr(param_2,'value'):

            param_2 = param_2.path

        # Check that the parameters exist
        assert param_1 in self._likelihood_model.free_parameters, "Parameter %s is not a free parameters of the " \
                                                                 "current model" % param_1

        if param_2 is not None:
            assert param_2 in self._likelihood_model.free_parameters, "Parameter %s is not a free parameters of the " \
                                                                      "current model" % param_2


        # Check that we have a valid fit

        assert self._current_minimum is not None, "You have to run the .fit method before calling get_contours."

        # Then restore the best fit

        self._minimizer._restore_best_fit()

        # Check minimal assumptions about the procedure

        assert not (param_1 == param_2), "You have to specify two different parameters"

        assert param_1_minimum < param_1_maximum, "Minimum larger than maximum for parameter 1"

        if param_2 is not None:
            assert param_2_minimum < param_2_maximum, "Minimum larger than maximum for parameter 2"

        # Check whether we are parallelizing or not

        if not threeML_config['parallel']['use-parallel']:

            a, b, cc = self.minimizer.contours(param_1, param_1_minimum, param_1_maximum, param_1_n_steps,
                                               param_2, param_2_minimum, param_2_maximum, param_2_n_steps,
                                               progress, **options)

            # Collapse the second dimension of the results if we are doing a 1d contour

            if param_2 is None:
                cc = cc[:, 0]

        else:

            # With parallel computation

            # In order to distribute fairly the computation, the strategy is to parallelize the computation
            # by assigning to the engines one "line" of the grid at the time

            # Connect to the engines

            client = ParallelClient(**options)

            # Get the number of engines

            n_engines = client.get_number_of_engines()

            # Check whether the number of threads is larger than the number of steps in the first direction

            if n_engines > param_1_n_steps:

                n_engines = int(param_1_n_steps)

                custom_warnings.warn("The number of engines is larger than the number of steps. Using only %s engines."
                                     % n_engines, ReducingNumberOfThreads)

            # Check if the number of steps is divisible by the number
            # of threads, otherwise issue a warning and make it so

            if float(param_1_n_steps) % n_engines != 0:
                # Set the number of steps to an integer multiple of the engines
                # (note that // is the floor division, also called integer division)

                param_1_n_steps = (param_1_n_steps // n_engines) * n_engines

                custom_warnings.warn("Number of steps is not a multiple of the number of threads. Reducing steps to %s"
                                     % param_1_n_steps, ReducingNumberOfSteps)

            # Compute the number of splits, i.e., how many lines in the grid for each engine.
            # (note that this is guaranteed to be an integer number after the previous checks)

            p1_split_steps = param_1_n_steps // n_engines

            # Prepare arrays for results

            if param_2 is None:

                # One array
                pcc = numpy.zeros(param_1_n_steps)

                pa = numpy.linspace(param_1_minimum, param_1_maximum, param_1_n_steps)
                pb = None

            else:

                pcc = numpy.zeros((param_1_n_steps, param_2_n_steps))

                # Prepare the two axes of the parameter space
                pa = numpy.linspace(param_1_minimum, param_1_maximum, param_1_n_steps)
                pb = numpy.linspace(param_2_minimum, param_2_maximum, param_2_n_steps)

            # Define the parallel worker which will go through the computation

            # NOTE: I only divide
            # on the first parameter axis so that the different
            # threads are more or less well mixed for points close and
            # far from the best fit

            def worker(start_index):

                # Re-create the minimizer

                # backup_freeParameters = copy.deepcopy(self.freeParameters)

                this_minimizer = self.Minimizer(self.minus_log_like_profile,
                                                self._free_parameters)

                this_p1min = pa[start_index * p1_split_steps]
                this_p1max = pa[(start_index + 1) * p1_split_steps - 1]