import numpy as np
from scipy.special import gammaln, psi

from pybasicbayes.abstractions import GibbsSampling, MeanField, MeanFieldSVI
from pyhawkes.internals.distributions import Dirichlet

class DirichletImpulseResponses(GibbsSampling, MeanField, MeanFieldSVI):
    """
    Encapsulates the impulse response vector distribution. In the
    discrete time Hawkes model this is a set of Dirichlet-distributed
    vectors of length B for each pair of processes, k and k', which
    we denote $\bbeta^{(k,k')}. This class contains all K^2 vectors.
    """
    def __init__(self, model, gamma=None):
        """
        Initialize a set of Dirichlet weight vectors.
        :param K:     The number of processes in the model.
        :param B:     The number of basis functions in the model.
        :param gamma: The Dirichlet prior parameter. If none it will be set
                      to a symmetric prior with parameter 1.
        """
        # assert isinstance(model, DiscreteTimeNetworkHawkesModel), \
        #        "model must be a DiscreteTimeNetworkHawkesModel"
        self.model = model
        self.K = model.K
        self.B = model.B

        if gamma is not None:
            assert np.isscalar(gamma) or \
                   (isinstance(gamma, np.ndarray) and
                    gamma.shape == (self.B,)), \
                "gamma must be a scalar or a length B vector"

            if np.isscalar(gamma):
                self.gamma = gamma * np.ones(self.B)
            else:
                self.gamma = gamma
        else:
            self.gamma = np.ones(self.B)

        # Initialize with a draw from the prior
        self.g = np.empty((self.K, self.K, self.B))
        self.resample()

        # Initialize mean field parameters
        self.mf_gamma = self.gamma[None, None, :] * np.ones((self.K, self.K, self.B))

    @property
    def impulses(self):
        basis = self.model.basis.basis
        return np.tensordot(basis, self.g, axes=[1,2])

    def rvs(self, size=[]):
        """
        Sample random variables from the Dirichlet impulse response distribution.
        :param size:
        :return:
        """
        pass

    def log_likelihood(self, x):
        '''
        log likelihood (either log probability mass function or log probability
        density function) of x, which has the same type as the output of rvs()
        '''
        assert isinstance(x, np.ndarray) and x.shape == (self.K,self.K,self.B), \
            "x must be a KxKxB array of impulse responses"

        gamma = self.gamma
        # Compute the normalization constant
        Z = gammaln(gamma).sum() - gammaln(gamma.sum())
        # Add the likelihood of x
        return self.K**2 * Z + ((gamma-1.0)[None,None,:] * np.log(x)).sum()

    def log_probability(self):
        return self.log_likelihood(self.g)

    def resample(self, data=[]):
        """
        Resample the
        """
        ss = np.zeros((self.K, self.K, self.B)) + \
             sum([d.compute_ir_ss() for d in data])

        for k1 in range(self.K):
            for k2 in range(self.K):
                alpha_post = self.gamma + ss[k1, k2, :]
                self.g[k1,k2,:] = np.random.dirichlet(alpha_post)

    def expected_g(self):
        # \sum_{b} \gamma_b
        trm2 = self.mf_gamma.sum(axis=2)
        E_g = self.mf_gamma / trm2[:,:,None]
        return E_g

    def expected_log_g(self):
        E_lng = np.zeros_like(self.mf_gamma)

        # \psi(\sum_{b} \gamma_b)
        trm2 = psi(self.mf_gamma.sum(axis=2))
        for b in range(self.B):
            E_lng[:,:,b] = psi(self.mf_gamma[:,:,b]) - trm2

        return E_lng

    def mf_update_gamma(self, data, minibatchfrac=1.0, stepsize=1.0):
        """
        Update gamma given E[Z]
        :return:
        """
        exp_ss = sum([d.compute_exp_ir_ss() for d in data])
        gamma_hat = self.gamma + exp_ss / minibatchfrac
        self.mf_gamma = (1.0 - stepsize) * self.mf_gamma + stepsize * gamma_hat

    def expected_log_likelihood(self,x):
        pass

    def meanfieldupdate(self, data=[]):
        self.mf_update_gamma(data)

    def meanfield_sgdstep(self, data, minibatchfrac, stepsize):
        self.mf_update_gamma(data, minibatchfrac=minibatchfrac, stepsize=stepsize)

    def get_vlb(self):
        """
        Variational lower bound for \lambda_k^0
        E[LN p(g | \gamma)] -
        E[LN q(g | \tilde{\gamma})]
        :return:
        """
        vlb = 0

        # First term
        # E[LN p(g | \gamma)]
        E_ln_g = self.expected_log_g()
        vlb += Dirichlet(self.gamma[None, None, :]).negentropy(E_ln_g=E_ln_g).sum()

        # Second term
        # E[LN q(g | \tilde{gamma})]
        vlb -= Dirichlet(self.mf_gamma).negentropy().sum()

        return vlb

    def resample_from_mf(self):
        """
        Resample from the mean field distribution
        :return:
        """
        self.g = np.zeros((self.K, self.K, self.B))
        for k1 in range(self.K):
            for k2 in range(self.K):
                self.g[k1,k2,:] = np.random.dirichlet(self.mf_gamma[k1,k2,:])


class SBMDirichletImpulseResponses(GibbsSampling):
    """
    A impulse response vector model with a set of Dirichlet-distributed
    vectors of length B for each pair of blocks, c and c', which
    we denote $\bbeta^{(c,c')}. This class contains all C^2 vectors.
    """
    def __init__(self, C, K, B, gamma=None):
        """
        Initialize a set of Dirichlet weight vectors.
        :param K:     The number of processes in the model.
        :param B:     The number of basis functions in the model.
        :param gamma: The Dirichlet prior parameter. If none it will be set
                      to a symmetric prior with parameter 1.
        """
        self.C = C
        self.K = K
        self.B = B

        if gamma is not None:
            assert np.isscalar(gamma) or \
                   (isinstance(gamma, np.ndarray) and
                    gamma.shape == (B,)), \
                "gamma must be a scalar or a length B vector"

            if np.isscalar(gamma):
                self.gamma = gamma * np.ones(B)
            else:
                self.gamma = gamma
        else:
            self.gamma = np.ones(self.B)

        # Initialize with a draw from the prior
        self.blockg = np.empty((self.C, self.C, self.B))
        self.resample()

        # Initialize mean field parameters
        self.mf_gamma = self.gamma[None, None, :] * np.ones((self.C, self.C, self.B))

    def rvs(self, size=[]):
        """
        Sample random variables from the Dirichlet impulse response distribution.
        :param size:
        :return:
        """
        raise NotImplementedError()

    def log_likelihood(self, x):
        '''
        log likelihood (either log probability mass function or log probability
        density function) of x, which has the same type as the output of rvs()
        '''
        raise NotImplementedError()
        assert isinstance(x, np.ndarray) and x.shape == (self.K,self.K,self.B), \
            "x must be a KxKxB array of impulse responses"

        gamma = self.gamma
        # Compute the normalization constant
        Z = gammaln(gamma).sum() - gammaln(gamma.sum())
        # Add the likelihood of x
        return self.K**2 * Z + ((gamma-1.0)[None,None,:] * np.log(x)).sum()

    def log_probability(self):
        return self.log_likelihood(self.g)

    def _get_suff_statistics(self, data):
        """
        Compute the sufficient statistics from the data set.
        :param data: a TxK array of event counts assigned to the background process
        :return:
        """
        raise NotImplementedError()
        # The only sufficient statistic is the KxKxB array of event counts assigned
        # to each of the basis functions
        if data is not None:
            ss = data.sum(axis=0)
        else:
            ss = np.zeros((self.K, self.K, self.B))

        return ss

    def resample(self, data=None):
        """
        Resample the
        :param data: a TxKxKxB array of parents. T time bins, K processes,
                     K parent processes, and B bases for each parent process.
        """
        raise NotImplementedError()
        assert data is None or \
               (isinstance(data, np.ndarray) and
                data.ndim == 4 and
                data.shape[1] == data.shape[2] == self.K
                and data.shape[3] == self.B), \
            "Data must be a TxKxKxB array of parents"


        ss = self._get_suff_statistics(data)
        for k1 in range(self.K):
            for k2 in range(self.K):
                alpha_post = self.gamma + ss[k1, k2, :]
                self.g[k1,k2,:] = np.random.dirichlet(alpha_post)

    def expected_g(self):
        """
        Compute the expected impulse response vector wrt c and mf_gamma
        :return:
        """
        raise NotImplementedError()
        # \sum_{b} \gamma_b
        trm2 = self.mf_gamma.sum(axis=2)
        E_g = self.mf_gamma / trm2[:,:,None]
        return E_g

    def expected_log_g(self):
        """
        Compute the expected log impulse response vector wrt c and mf_gamma
        :return:
        """
        raise NotImplementedError()
        E_lng = np.zeros_like(self.mf_gamma)

        # \psi(\sum_{b} \gamma_b)
        trm2 = psi(self.mf_gamma.sum(axis=2))
        for b in range(self.B):
            E_lng[:,:,b] = psi(self.mf_gamma[:,:,b]) - trm2

        return E_lng

    def mf_update_gamma(self, EZ):
        """
        Update gamma given E[Z]
        :return:
        """
        raise NotImplementedError()
        self.mf_gamma = self.gamma + EZ.sum(axis=0)

    def expected_log_likelihood(self,x):
        raise NotImplementedError()
        pass

    def meanfieldupdate(self, EZ):
        raise NotImplementedError()
        self.mf_update_gamma(EZ)

    def get_vlb(self):
        """
        Variational lower bound for \lambda_k^0
        E[LN p(g | \gamma)] -
        E[LN q(g | \tilde{\gamma})]
        :return:
        """
        raise NotImplementedError()
        vlb = 0

        # First term
        # E[LN p(g | \gamma)]
        E_ln_g = self.expected_log_g()
        vlb += Dirichlet(self.gamma[None, None, :]).negentropy(E_ln_g=E_ln_g).sum()

        # Second term
        # E[LN q(g | \tilde{gamma})]
        vlb -= Dirichlet(self.mf_gamma).negentropy().sum()

        return vlb

    def resample_from_mf(self):
        """
        Resample from the mean field distribution
        :return:
        """
        raise NotImplementedError()
        self.g = np.zeros((self.K, self.K, self.B))
        for k1 in range(self.K):
            for k2 in range(self.K):
                self.g[k1,k2,:] = np.random.dirichlet(self.mf_gamma[k1,k2,:])


class ContinuousTimeImpulseResponses(GibbsSampling):
    """
    Continuous time impulse response model with logistic normal
    impulse response functions.
    """
    def __init__(self, model, mu_0=0., lmbda_0=1., alpha_0=1., beta_0=1.):
        self.model = model
        self.K = model.K
        self.dt_max = model.dt_max

        self.mu_0 = mu_0
        self.lmbda_0 = lmbda_0
        self.alpha_0 = alpha_0
        self.beta_0 = beta_0

        from pyhawkes.utils.utils import sample_nig
        self.mu, self.tau = \
            sample_nig(self.mu_0 * np.ones((self.K, self.K)),
                       self.lmbda_0 * np.ones((self.K, self.K)),
                       self.alpha_0 * np.ones((self.K, self.K)),
                       self.beta_0 * np.ones((self.K, self.K)))

    @property
    def impulses(self):
        N_pts = 50
        t = np.linspace(0, self.dt_max, N_pts)
        ir = np.zeros((N_pts, self.K, self.K))
        for k1 in range(self.K):
            for k2 in range(self.K):
                ir[:,k1,k2] = self.impulse(t, k1, k2)
        return t, ir

    # TODO: Rename this
    def impulse(self, dt, k1, k2):
        """
        Impulse response induced by an event on process k1 on
        the rate of process k2 at lag dt
        """
        from pyhawkes.utils.utils import logit
        mu, tau, dt_max = self.mu[k1,k2], self.tau[k1,k2], self.dt_max
        Z = dt * (dt_max - dt)/dt_max * np.sqrt(2*np.pi/tau)
        return 1./Z * np.exp(-tau/2. * (logit(dt/dt_max) - mu)**2)

    def rvs(self, size=[]):
        """
        Sample random variables from the Dirichlet impulse response distribution.
        :param size:
        :return:
        """
        pass

    def log_likelihood(self, x):
        '''
        log likelihood (either log probability mass function or log probability
        density function) of x, which has the same type as the output of rvs()
        '''
        return 0

    def log_probability(self):
        return self.log_likelihood((self.mu, self.tau))

    def resample(self, data=[]):
        """
        Resample the
        :param data: a TxKxKxB array of parents. T time bins, K processes,
                     K parent processes, and B bases for each parent process.
        """
        mu_0, lmbda_0, alpha_0, beta_0 = self.mu_0, self.lmbda_0, self.alpha_0, self.beta_0
        assert data is None or isinstance(data, list)

        # 0: count, # 1: Sum of scaled dt, #2: Sum of sq scaled dt
        ss = np.zeros((3, self.K, self.K))
        for d in data:
            ss += d.compute_imp_suff_stats()

        n = ss[0]
        xbar = np.nan_to_num(ss[1] / n)
        xvar = ss[2]

        alpha_post = alpha_0 + n / 2.
        # beta_post = beta_0 + 0.5 * xvar
        # beta_post += 0.5 * lmbda_0 * n / (lmbda_0 + n) * (xbar-mu_0)**2
        beta_post = beta_0 + 0.5 * xvar + 0.5 * lmbda_0 * n * (xbar-mu_0)**2 / (lmbda_0+n)

        lmbda_post = lmbda_0 + n
        mu_post = (lmbda_0 * mu_0 + n * xbar) / (lmbda_0 + n)

        from pyhawkes.utils.utils import sample_nig
        self.mu, self.tau = \
            sample_nig(mu_post, lmbda_post, alpha_post, beta_post)

        assert np.isfinite(self.mu).all()
        assert np.isfinite(self.tau).all()

