from .baseanalysis import BaseAnalysis
from modelflow import Input
import numpy as np
import numpy.linalg as la
import sys


class RelativeCycleAnalysis(BaseAnalysis):
    """
    Compares relative to ER the expected number of cycles which is approximated here
    by :math:`(p*n)^{length}`, where ``p`` is the connectivity, ``n`` the number of excitatory nodes
    and ``length`` the cycle length.

    Note that this approximation ignores the non-independent nature of the Bernoulli
    product on the diagonal of the matrix power. This approximation fails for very small
    and very sparse networks but is precise for large and dense networks.

    Configuration parameter: length, integer

    Output keys: relative_cycles_<length>
    """
    length = Input()

    def analyse(self):
        self._verify_length()
        p_exc = np.count_nonzero(self.network["E", "E"]) / self.network.nr_exc ** 2
        self._verify_connectivity_range(p_exc)
        if p_exc == 0:
            relative_cycles = 0
        else:
            measured_cycles = np.trace(la.matrix_power((self.network["E", "E"] != 0).astype(float), self.length))
            expected = self._expected_cycles(p_exc)
            relative_cycles = measured_cycles / expected if expected != 0 else 0
        return {self._result_key(): relative_cycles}

    def _expected_cycles(self, p_exc):
        return (self.network.nr_exc * p_exc) ** self.length

    def _verify_length(self):
        if self.length < 3 or self.length > 8:
            print("Cycle analysis only tested for length in [3,8]. Got {length}".format(length=self.length), file=sys.stderr)

    def _verify_connectivity_range(self, p):
        if p < .07 or p > .4:
            print("Warning Cycle analysis only tested for p in [0.07, 0.4]. Got {p}".format(p=p), file=sys.stderr)

    def _result_key(self):
        return "relative_cycles_" + str(self.length)
