from collections.abc import Mapping
from ctypes import c_int, c_int32, c_size_t, c_double, c_char_p, c_bool, POINTER
from weakref import WeakValueDictionary

import numpy as np
from numpy.ctypeslib import as_array
import scipy.stats

from openmc.exceptions import AllocationError, InvalidIDError
from openmc.data.reaction import REACTION_NAME
from . import _dll, Nuclide
from .core import _FortranObjectWithID
from .error import _error_handler
from .filter import _get_filter


__all__ = ['Tally', 'tallies', 'global_tallies', 'num_realizations']

# Tally functions
_dll.openmc_extend_tallies.argtypes = [c_int32, POINTER(c_int32), POINTER(c_int32)]
_dll.openmc_extend_tallies.restype = c_int
_dll.openmc_extend_tallies.errcheck = _error_handler
_dll.openmc_get_tally_index.argtypes = [c_int32, POINTER(c_int32)]
_dll.openmc_get_tally_index.restype = c_int
_dll.openmc_get_tally_index.errcheck = _error_handler
_dll.openmc_global_tallies.argtypes = [POINTER(POINTER(c_double))]
_dll.openmc_global_tallies.restype = c_int
_dll.openmc_global_tallies.errcheck = _error_handler
_dll.openmc_tally_get_active.argtypes = [c_int32, POINTER(c_bool)]
_dll.openmc_tally_get_active.restype = c_int
_dll.openmc_tally_get_active.errcheck = _error_handler
_dll.openmc_tally_get_estimator.argtypes = [c_int32, POINTER(c_int)]
_dll.openmc_tally_get_estimator.restype = c_int
_dll.openmc_tally_get_estimator.errcheck = _error_handler
_dll.openmc_tally_get_id.argtypes = [c_int32, POINTER(c_int32)]
_dll.openmc_tally_get_id.restype = c_int
_dll.openmc_tally_get_id.errcheck = _error_handler
_dll.openmc_tally_get_filters.argtypes = [
    c_int32, POINTER(POINTER(c_int32)), POINTER(c_size_t)]
_dll.openmc_tally_get_filters.restype = c_int
_dll.openmc_tally_get_filters.errcheck = _error_handler
_dll.openmc_tally_get_multiply_density.argtypes = [c_int32, POINTER(c_bool)]
_dll.openmc_tally_get_multiply_density.restype = c_int
_dll.openmc_tally_get_multiply_density.errcheck = _error_handler
_dll.openmc_tally_get_n_realizations.argtypes = [c_int32, POINTER(c_int32)]
_dll.openmc_tally_get_n_realizations.restype = c_int
_dll.openmc_tally_get_n_realizations.errcheck = _error_handler
_dll.openmc_tally_get_nuclides.argtypes = [
    c_int32, POINTER(POINTER(c_int)), POINTER(c_int)]
_dll.openmc_tally_get_nuclides.restype = c_int
_dll.openmc_tally_get_nuclides.errcheck = _error_handler
_dll.openmc_tally_get_scores.argtypes = [
    c_int32, POINTER(POINTER(c_int)), POINTER(c_int)]
_dll.openmc_tally_get_scores.restype = c_int
_dll.openmc_tally_get_scores.errcheck = _error_handler
_dll.openmc_tally_get_type.argtypes = [c_int32, POINTER(c_int32)]
_dll.openmc_tally_get_type.restype = c_int
_dll.openmc_tally_get_type.errcheck = _error_handler
_dll.openmc_tally_get_writable.argtypes = [c_int32, POINTER(c_bool)]
_dll.openmc_tally_get_writable.restype = c_int
_dll.openmc_tally_get_writable.errcheck = _error_handler
_dll.openmc_tally_reset.argtypes = [c_int32]
_dll.openmc_tally_reset.restype = c_int
_dll.openmc_tally_reset.errcheck = _error_handler
_dll.openmc_tally_results.argtypes = [
    c_int32, POINTER(POINTER(c_double)), POINTER(c_size_t*3)]
_dll.openmc_tally_results.restype = c_int
_dll.openmc_tally_results.errcheck = _error_handler
_dll.openmc_tally_set_active.argtypes = [c_int32, c_bool]
_dll.openmc_tally_set_active.restype = c_int
_dll.openmc_tally_set_active.errcheck = _error_handler
_dll.openmc_tally_set_filters.argtypes = [c_int32, c_size_t, POINTER(c_int32)]
_dll.openmc_tally_set_filters.restype = c_int
_dll.openmc_tally_set_filters.errcheck = _error_handler
_dll.openmc_tally_set_estimator.argtypes = [c_int32, c_char_p]
_dll.openmc_tally_set_estimator.restype = c_int
_dll.openmc_tally_set_estimator.errcheck = _error_handler
_dll.openmc_tally_set_id.argtypes = [c_int32, c_int32]
_dll.openmc_tally_set_id.restype = c_int
_dll.openmc_tally_set_id.errcheck = _error_handler
_dll.openmc_tally_set_multiply_density.argtypes = [c_int32, c_bool]
_dll.openmc_tally_set_multiply_density.restype = c_int
_dll.openmc_tally_set_multiply_density.errcheck = _error_handler
_dll.openmc_tally_set_nuclides.argtypes = [c_int32, c_int, POINTER(c_char_p)]
_dll.openmc_tally_set_nuclides.restype = c_int
_dll.openmc_tally_set_nuclides.errcheck = _error_handler
_dll.openmc_tally_set_scores.argtypes = [c_int32, c_int, POINTER(c_char_p)]
_dll.openmc_tally_set_scores.restype = c_int
_dll.openmc_tally_set_scores.errcheck = _error_handler
_dll.openmc_tally_set_type.argtypes = [c_int32, c_char_p]
_dll.openmc_tally_set_type.restype = c_int
_dll.openmc_tally_set_type.errcheck = _error_handler
_dll.openmc_tally_set_writable.argtypes = [c_int32, c_bool]
_dll.openmc_tally_set_writable.restype = c_int
_dll.openmc_tally_set_writable.errcheck = _error_handler
_dll.openmc_remove_tally.argtypes = [c_int32]
_dll.openmc_remove_tally.restype = c_int
_dll.openmc_remove_tally.errcheck = _error_handler
_dll.tallies_size.restype = c_size_t


_SCORES = {
    -1: 'flux', -2: 'total', -3: 'scatter', -4: 'nu-scatter',
    -5: 'absorption', -6: 'fission', -7: 'nu-fission', -8: 'kappa-fission',
    -9: 'current', -10: 'events', -11: 'delayed-nu-fission',
    -12: 'prompt-nu-fission', -13: 'inverse-velocity', -14: 'fission-q-prompt',
    -15: 'fission-q-recoverable', -16: 'decay-rate', -17: 'pulse-height'
}
_ESTIMATORS = {
    0: 'analog', 1: 'tracklength', 2: 'collision'
}
_TALLY_TYPES = {
    0: 'volume', 1: 'mesh-surface', 2: 'surface', 3: 'pulse-height'
}


def global_tallies():
    """Mean and standard deviation of the mean for each global tally.

    Returns
    -------
    list of tuple
        For each global tally, a tuple of (mean, standard deviation)

    """
    ptr = POINTER(c_double)()
    _dll.openmc_global_tallies(ptr)
    array = as_array(ptr, (4, 3))

    # Get sum, sum-of-squares, and number of realizations
    sum_ = array[:, 1]
    sum_sq = array[:, 2]
    n = num_realizations()

    # Determine mean
    if n > 0:
        mean = sum_ / n
    else:
        mean = sum_.copy()

    # Determine standard deviation
    nonzero = np.abs(mean) > 0
    stdev = np.empty_like(mean)
    stdev.fill(np.inf)
    if n > 1:
        stdev[nonzero] = np.sqrt((sum_sq[nonzero]/n - mean[nonzero]**2)/(n - 1))

    return list(zip(mean, stdev))


def num_realizations():
    """Number of realizations of global tallies."""
    return c_int32.in_dll(_dll, 'n_realizations').value


class Tally(_FortranObjectWithID):
    """Tally stored internally.

    This class exposes a tally that is stored internally in the OpenMC
    library. To obtain a view of a tally with a given ID, use the
    :data:`openmc.lib.tallies` mapping.

    Parameters
    ----------
    uid : int or None
        Unique ID of the tally
    new : bool
        When `index` is None, this argument controls whether a new object is
        created or a view of an existing object is returned.
    index : int or None
        Index in the `tallies` array.

    Attributes
    ----------
    id : int
        ID of the tally
    estimator: str
        Estimator type of tally (analog, tracklength, collision)
    filters : list
        List of tally filters
    mean : numpy.ndarray
        An array containing the sample mean for each bin
    multiply_density : bool
        Whether reaction rates should be multiplied by atom density

        .. versionadded:: 0.14.0
    nuclides : list of str
        List of nuclides to score results for
    num_realizations : int
        Number of realizations
    results : numpy.ndarray
        Array of tally results
    std_dev : numpy.ndarray
        An array containing the sample standard deviation for each bin
    type : str
        Type of tally (volume, mesh_surface, surface)

    """
    __instances = WeakValueDictionary()

    def __new__(cls, uid=None, new=True, index=None):
        mapping = tallies
        if index is None:
            if new:
                # Determine ID to assign
                if uid is None:
                    uid = max(mapping, default=0) + 1
                else:
                    if uid in mapping:
                        raise AllocationError('A tally with ID={} has already '
                                              'been allocated.'.format(uid))

                index = c_int32()
                _dll.openmc_extend_tallies(1, index, None)
                index = index.value
            else:
                index = mapping[uid]._index

        if index not in cls.__instances:
            instance = super().__new__(cls)
            instance._index = index
            if uid is not None:
                instance.id = uid
            cls.__instances[index] = instance

        return cls.__instances[index]

    @property
    def active(self):
        active = c_bool()
        _dll.openmc_tally_get_active(self._index, active)
        return active.value

    @active.setter
    def active(self, active):
        _dll.openmc_tally_set_active(self._index, active)

    @property
    def type(self):
        type = c_int32()
        _dll.openmc_tally_get_type(self._index, type)
        return _TALLY_TYPES[type.value]

    @type.setter
    def type(self, type):
        _dll.openmc_tally_set_type(self._index, type.encode())

    @property
    def estimator(self):
        estimator = c_int()
        _dll.openmc_tally_get_estimator(self._index, estimator)
        return _ESTIMATORS[estimator.value]

    @estimator.setter
    def estimator(self, estimator):
        _dll.openmc_tally_set_estimator(self._index, estimator.encode())

    @property
    def id(self):
        tally_id = c_int32()
        _dll.openmc_tally_get_id(self._index, tally_id)
        return tally_id.value

    @id.setter
    def id(self, tally_id):
        _dll.openmc_tally_set_id(self._index, tally_id)

    @property
    def filters(self):
        filt_idx = POINTER(c_int32)()
        n = c_size_t()
        _dll.openmc_tally_get_filters(self._index, filt_idx, n)
        return [_get_filter(filt_idx[i]) for i in range(n.value)]

    @filters.setter
    def filters(self, filters):
        # Get filter indices as int32_t[]
        n = len(filters)
        indices = (c_int32*n)(*(f._index for f in filters))

        _dll.openmc_tally_set_filters(self._index, n, indices)

    def find_filter(self, filter_type):
        """
        Returns the first instance of a filter matching the specified type

        Parameters
        ----------
        filter_type : subclass of openmc.lib.Filter
            The filter type to match when retrieving a filter instance

        Returns
        -------
        filter : openmc.lib.Filter
            The filter instance matching the input filter type

        Raises
        ------
        ValueError if a filter instance matching the input filter type cannot be found.
        """
        for filter in self.filters:
            if isinstance(filter, filter_type):
                return filter

        raise ValueError(f'No filter of type {filter_type} on tally {self.id}')

    @property
    def mean(self):
        n = self.num_realizations
        sum_ = self.results[:, :, 1]
        if n > 0:
            return sum_ / n
        else:
            return sum_.copy()

    @property
    def nuclides(self):
        nucs = POINTER(c_int)()
        n = c_int()
        _dll.openmc_tally_get_nuclides(self._index, nucs, n)
        return [Nuclide(nucs[i]).name if nucs[i] >= 0 else 'total'
                for i in range(n.value)]

    @nuclides.setter
    def nuclides(self, nuclides):
        nucs = (c_char_p * len(nuclides))()
        nucs[:] = [x.encode() for x in nuclides]
        _dll.openmc_tally_set_nuclides(self._index, len(nuclides), nucs)

    @property
    def num_realizations(self):
        n = c_int32()
        _dll.openmc_tally_get_n_realizations(self._index, n)
        return n.value

    @property
    def results(self):
        data = POINTER(c_double)()
        shape = (c_size_t*3)()
        _dll.openmc_tally_results(self._index, data, shape)
        return as_array(data, tuple(shape))

    @property
    def scores(self):
        scores_as_int = POINTER(c_int)()
        n = c_int()
        try:
            _dll.openmc_tally_get_scores(self._index, scores_as_int, n)
        except AllocationError:
            return []
        else:
            scores = []
            for i in range(n.value):
                if scores_as_int[i] in _SCORES:
                    scores.append(_SCORES[scores_as_int[i]])
                elif scores_as_int[i] in REACTION_NAME:
                    scores.append(REACTION_NAME[scores_as_int[i]])
                else:
                    scores.append(str(scores_as_int[i]))
            return scores

    @scores.setter
    def scores(self, scores):
        scores_ = (c_char_p * len(scores))()
        scores_[:] = [x.encode() for x in scores]
        _dll.openmc_tally_set_scores(self._index, len(scores), scores_)

    @property
    def std_dev(self):
        results = self.results
        std_dev = np.empty(results.shape[:2])
        std_dev.fill(np.inf)

        n = self.num_realizations
        if n > 1:
            # Get sum and sum-of-squares from results
            sum_ = results[:, :, 1]
            sum_sq = results[:, :, 2]

            # Determine non-zero entries
            mean = sum_ / n
            nonzero = np.abs(mean) > 0

            # Calculate sample standard deviation of the mean
            std_dev[nonzero] = np.sqrt(
                (sum_sq[nonzero]/n - mean[nonzero]**2)/(n - 1))

        return std_dev

    @property
    def writable(self):
        writable = c_bool()
        _dll.openmc_tally_get_writable(self._index, writable)
        return writable.value

    @writable.setter
    def writable(self, writable):
        _dll.openmc_tally_set_writable(self._index, writable)

    @property
    def multiply_density(self):
        multiply_density = c_bool()
        _dll.openmc_tally_get_multiply_density(self._index, multiply_density)
        return multiply_density.value

    @multiply_density.setter
    def multiply_density(self, multiply_density):
        _dll.openmc_tally_set_multiply_density(self._index, multiply_density)

    def reset(self):
        """Reset results and num_realizations of tally"""
        _dll.openmc_tally_reset(self._index)

    def ci_width(self, alpha=0.05):
        """Confidence interval half-width based on a Student t distribution

        Parameters
        ----------
        alpha : float
            Significance level (one minus the confidence level!)

        Returns
        -------
        float
            Half-width of a two-sided (1 - :math:`alpha`) confidence interval

        """
        half_width = self.std_dev.copy()
        n = self.num_realizations
        if n > 1:
            half_width *= scipy.stats.t.ppf(1 - alpha/2, n - 1)
        return half_width


class _TallyMapping(Mapping):
    def __getitem__(self, key):
        index = c_int32()
        try:
            _dll.openmc_get_tally_index(key, index)
        except (AllocationError, InvalidIDError) as e:
            # __contains__ expects a KeyError to work correctly
            raise KeyError(str(e))
        return Tally(index=index.value)

    def __iter__(self):
        for i in range(len(self)):
            yield Tally(index=i).id

    def __len__(self):
        return _dll.tallies_size()

    def __repr__(self):
        return repr(dict(self))

    def __delitem__(self, key):
        """Delete a tally from tally vector and remove the ID,index pair from tally"""
        _dll.openmc_remove_tally(self[key]._index)

tallies = _TallyMapping()
