# cython: boundscheck=False, cdivision=True
import numpy as np
cimport numpy as np
from libc.math cimport log, exp

DTYPE = np.float64
ctypedef np.float_t DTYPE_t

cdef DTYPE_t plogis(DTYPE_t weight, DTYPE_t delta):
    return 1./(1. + exp(-weight/delta))

cdef DTYPE_t weight_success(int dist, int dt, DTYPE_t alpha, DTYPE_t delta,
                            DTYPE_t w1, DTYPE_t w2, DTYPE_t w3):
    """weight for successful infection, exponential time model"""
    cdef DTYPE_t structural, temporal, result
    structural = dist * log(delta)
    # structural = plogis(w1,delta) * plogis(w2,delta) * plogis(w3,delta)
    temporal = log(exp(alpha)-1.) - alpha*dt
    # temporal = 1 - exp(-alpha*dt)
    # if exp(-alpha*dt)==0.: print 'UNDERFLOW ERROR'
    # temporal = 1. / (1. + (dt - 1.)/alpha)**0.01 - 1. / (1. + dt/alpha)**0.01
    result = structural + temporal
    # print 'st', structural, temporal
    return result

cdef DTYPE_t weight_failure(int dist, int dt, DTYPE_t alpha, DTYPE_t delta,
                            DTYPE_t w1, DTYPE_t w2, DTYPE_t w3):
    """weight for failed infection, exponential time model"""
    cdef DTYPE_t structural, temporal, result
    structural = delta ** dist
    # structural = plogis(w1,delta) * plogis(w2,delta) * plogis(w3,delta)
    temporal = exp(-alpha * dt)
    # temporal = 1.  - 1. / (1. + dt/alpha)**0.01
    result = log(1. - structural + structural * temporal)
    # print 'stnv', structural, temporal
    return result

def ml(dict root_victims, dict victims, dict non_victims, DTYPE_t age,
       DTYPE_t alpha, DTYPE_t delta):
    cdef:
        int n_roots, n_victims, n_nodes, roots, i, dist, dt, t, l
        DTYPE_t beta, ll
        list parents, failures, successes
    n_roots, n_victims = len(root_victims), len(victims)
    n_nodes = 148152
    cdef:
        np.ndarray[DTYPE_t] probs = np.zeros(n_victims, dtype=DTYPE)
        np.ndarray[DTYPE_t] probs_fail = np.zeros(n_victims, dtype=DTYPE)
        np.ndarray[DTYPE_t] probs_nv = np.zeros(len(non_victims), dtype=DTYPE)

    # loop through victims
    for i, parents in enumerate(victims.itervalues()):
        # for each victim node i, compute the probability that all its parents
        # fail to infect it, also computes the probability that its most
        # likely parent infects it
        failures = [weight_failure(dist, dt, alpha, delta, w1, w2, w3)
                    for (dist, dt, w1, w2, w3) in parents]
        probs_fail[i] = sum(failures)
        successes = [weight_success(dist, dt, alpha, delta, w1, w2, w3)
                     for (dist, dt, w1, w2, w3) in parents]
        # find parent that maximizes log(p) - log(\tilde{p})
        probs[i] = max(s - failures[l] for l, s in enumerate(successes)) 

    # loop through non-victims
    for i, parents in enumerate(non_victims.itervalues()):
        # for each non victim node, compute the probability that all its
        # parents fail to infect it
        failures = [weight_failure(dist, dt, alpha, delta, w1, w2, w3)
                    for (dist, dt, w1, w2, w3) in parents]
        probs_nv[i] = sum(failures)

    # print successes
    # print failures
    # print probs

    # calculate log likelihood
    # probs.sort(); probs = probs[::-1] # sort probs in descending order
    # cdef:
    #     np.ndarray[DTYPE_t] cums = probs.cumsum()
    ll =  probs_fail.sum() # add probability that all edges to victims fail
    ll += probs_nv.sum() # add probability that all edges to non_victims fail

    # print 'probs', probs
    max_beta_add = float('-inf')
    # iterate over all victim nodes to find the optimal threshold
    for beta in np.arange(0.001, .2, .002):
        thresh = log(beta/(3012*(1.-beta)))
        # print 'beta:', beta, 'thresh:', thresh, 'infected:', len(probs[probs>=thresh])
        roots = n_roots + len(probs[probs<thresh])

        # add probability for realized edges and subtract probability these edges fail
        beta_add = (probs[probs>=thresh]).sum()
        # add probability for the seeds and non-seeds
        beta_add += roots * log(beta) + (n_nodes-roots) * log(1. - beta)

        if beta_add > max_beta_add:
            max_beta = beta
            max_roots = roots
            max_beta_add = beta_add
            # print 'beta:', max_beta, 'add:', max_beta_add, 'roots:', max_roots

    ll += max_beta_add
    roots = max_roots
    beta = max_beta
    # print n_nodes, n_roots, n_victims, max_i, roots
    return (beta, roots, ll)