# -*- coding: utf-8 -*-
#!python
#cython: boundscheck=False
#cython: cdivision=True
#cython: embedsignature=True
"""
Created on Tue May  6 13:24:42 2014

@author: emmanuel


Algorithms similar to Till Kretschmar's master thesis (2013).



This module implements path builder algorithms to detect
feed forward structure.
"""
import scipy as sp
import numpy as np
from util.urng cimport URNG
from util.max_seed import max_seed
from libc.stdlib cimport malloc, free
from scipy.stats.mstats import gmean

ctypedef unsigned long long IntCounter


def array_to_list(arr):
    """
    converts an adjacency matrix arr into a list of lists
    where the k-th entry of the list contains another list
    which entries are the labels of the outgoing vertices
    of vertex k
    """
    return [[k for k,l in enumerate(out)  if l] for out in arr.T ]



def extension_ratio_from_nr_paths_sampled(paths):
    first_zero = paths.shape[0]
    for k in range(paths.shape[0]):
        if paths[k] == 0:
            first_zero = k
            break
    
    paths = paths.astype(np.float)
    if first_zero >= 2:
        paths = paths[:first_zero]
        ratio_arr = paths[1:-1] / paths[2:]
        ratio = gmean(ratio_arr)
    else:
        ratio = 1
    return ratio



def ideal_probabilities_from_branching(arr, max_length, ratio_paths_to_nodes):
    """
    probabilities[0] = 1
    probabilities[1] = (nr nodes) / (nr edges)
    probabiliteis[2:] = calculated from tentative builder
    """
    MAX_LENGTH_FOR_PROBABILITY_CALCULATION = 10
    NR_REPETITIONS_FOR_PROBABILITY_CALCULATION = 8

    E = (arr != 0).astype(float)
    for k in range(E.shape[0]):
        E[k,k] = 0
        
    probabilities = sp.zeros(max_length)
    probabilities[0] = 1
    probabilities[1] = float(E.shape[0]) / E.sum()

    # first a heuristic initialization of the probabilities
    A = sp.dot(E,E)
    for k in range(A.shape[0]):
        A[k,k] = 0
    probabilities[2:] =  E.sum() / A.sum()

    # update heuristic probabilities with sampled ones
    max_length_for_calc = min(MAX_LENGTH_FOR_PROBABILITY_CALCULATION, E.shape[0])
    graph = array_to_list(E)
    paths, cycles = _run_pes(graph, max_length_for_calc, probabilities,
                             NR_REPETITIONS_FOR_PROBABILITY_CALCULATION)
                             
                                
    
    probabilities[2:] *= extension_ratio_from_nr_paths_sampled(paths)
    probabilities *= ratio_paths_to_nodes**(1. / max_length)
    return probabilities
    


cdef IntCounter node_in_path(int node, int * path, int path_length):
    """
    Returns
    -------
    
    Returns 1 of node is in the path of length path_length
    and 0 otherwise.
    
    In pure Python that would be sth. like "node in path"
    """
    cdef int k
    for k in range(path_length):
        if node == path[k]:
            return 1
    return 0


cdef void extend_path(int *path, int **c_graph, int *out_degrees,
                      int max_length,
                      int current_length,
                      float *c_probabilities,
                      IntCounter *pathlist,
                      IntCounter *cyclelist,
                      URNG rng):
    """
    Arguments
    ---------
    
        path:
            list of integers which make up the path to be extended
        c_graph : A list of lists of the type [ [1, 2], [0], [1] ]
            graph[k] is the set of vertices vertex k projects to
        out_degrees : 
            the out degrees of the vertices
        max_length:
            max path length to sample
        current_length:
            length of current path
        c_probabilities:
            probabilitiy with which to explore a possible path further
        pathlist: 
            pathlist[n] : nr of paths with n+1 vertices
        cyclelist:
            cyclelist[n] : nr of cycles with n+1 vertices
    """
    # since we have a path, increas number of paths of current_length by 1
    pathlist[current_length-1] += 1
    
    #check if the current path is a cycle
    # if it is a cycle increase cycle counter by 1
    cdef int k
    cdef int first_node = path[0]
    cdef int last_node = path[current_length-1]
    for k in range(out_degrees[last_node]): 
        if c_graph[last_node][k] == first_node:
            cyclelist[current_length-1] += 1
            break
        
    # explore outgoing edges of last node of the current path
    # extend the path if possible with a certain probability
    cdef int * new_path
    if current_length < max_length:
        for k in range(out_degrees[last_node]):
            if not node_in_path(c_graph[last_node][k], path, current_length ):
                if rng.rvs() < c_probabilities[current_length]:
                    path[current_length] = c_graph[last_node][k]
                    extend_path(path, c_graph, out_degrees, max_length,
                                current_length + 1,
                                c_probabilities, pathlist, cyclelist, rng)
  







def pes(arr, int max_length, probabilities=None, repeat=1, ratio_paths_to_nodes=1):
    """
    Sample strong edge paths uniformly using the path enumeration sampling algorithm [#wernicke_2006_pes]_.

    Parameters
    ----------
        arr : 2d array
            Directed adjacency matrix.
            
        max_length: int
            Sample paths up to length max_length.
            
        probabilities : 1d array of length max_length
            Probabilities[n] is the probability of
            exploring a path of length n further.
            
        repeat: int
            Repeat the builder procedure ``repeat`` times and sum up the so
            obtained paths/cycles.

        ratio_paths_to_nodes : float
            Try to sample ratio_paths_to_nodes * nr_nodes paths of length max_length.

            
    Returns
    -------
    
        paths : list
           List of nr of paths of given length:

               * paths[0] = Nr of nodes
               * paths[1] = Nr of paths with 1 edge
               * paths[2] = Nr of paths consisting of 2 edges
               * etc. ....
            
        cycles : list
         List of nr of cycles of given length.


    .. note::
        The algorithm is rather unstable for paths of length larger than about 50;
        in these cases the number of returned paths can both exponentially increase or decay
        (i.e. no path is sampled).
        For shorter paths (length shorter than about 20), this algorithm is stable.
        Importantly, in this algorithm the paths are sampled uniformly, unlike in
        random-walk based sampling algorithms which do not yield uniform samples.



    .. [#wernicke_2006_pes] Wernicke, Sebastian. “Efficient Detection of Network Motifs.”
                        IEEE/ACM Trans. Comput. Biol. Bioinformatics 3, no. 4 (October 2006):
                        347–359. doi:10.1109/TCBB.2006.51.

    """

    graph  =  array_to_list(arr != 0)
    if probabilities is None:
        probabilities = ideal_probabilities_from_branching(arr, max_length, ratio_paths_to_nodes)
    
    
    paths, cycles = _run_pes(graph, max_length, probabilities, repeat)
    return paths, cycles


def _run_pes(graph,int max_length, probabilities, repeat):
    """
    Parameters
    ----------
        graph : A list of lists of the type [ [1, 2], [0], [1] ]
            graph[k] is the set of vertices vertex k projects to
            
        max_length: int
            sample paths up to length max_length
            
        probabilities : array of length max_length
            probabilities[n] is the probability of
            exploring a path of length n further
            
        repeat: int
            nr of times to repeat the builder
            
    Returns
    -------
        paths : list of nr of paths of given length
            paths[0] = Nr of nodes
            paths[1] = Nr of edges
            paths[2] = Nr of paths consisting of 2 edges
            etc. ....
        cycles L list of nr of cycles of given length
    """
    # init path and cycle list where to store nr of paths and cycles
    # of a given length
    cdef:
        IntCounter *pathlist = <IntCounter *>malloc(max_length * sizeof(IntCounter))
        IntCounter *cyclelist = <IntCounter *>malloc(max_length * sizeof(IntCounter))
        int k
        URNG rng = URNG(sp.random.randint(max_seed), 0., 1.)
    for k in range(max_length):
        pathlist[k] = 0
        cyclelist[k] = 0
    
    # convert Python list of probabilities to C array of float
    cdef float * c_probabilities = <float *> malloc(len(probabilities)*sizeof(float) )
    for k in range(len(probabilities)):
        c_probabilities[k] = float(probabilities[k])
    
    #convert graph list of list to C array type 
    cdef:
        int nr_nodes = len(graph)
        int *out_degrees = <int *>malloc(nr_nodes * sizeof(int))
        int **c_graph = <int **>malloc(nr_nodes * sizeof(int *))
        int curr_out_deg
    for k in range(nr_nodes):
        curr_out_deg = len(graph[k])
        out_degrees[k] = curr_out_deg
        c_graph[k] = <int *>malloc(curr_out_deg * sizeof(int))
        for l in range(curr_out_deg):
            c_graph[k][l] = graph[k][l]
    
    #make a stack of paths for depth first search
    cdef int *path = <int *>malloc(max_length * sizeof(int))
    
    for _ in range(repeat):    
        for k in range(nr_nodes):
            if rng.rvs() < c_probabilities[0]: #  : #c_probabilities[0] * RAND_MAX:
                # create initial path consisting of 1 single vertex
                path[0] = k
                extend_path(path, c_graph, out_degrees,max_length, 1,
                            c_probabilities, pathlist, cyclelist, rng)
        
    

                    
    # convert for output back to Python list types
    path_ret = [ pathlist[k] for k in range(max_length) ]
    cycle_ret = [ cyclelist[k] for k in range(max_length) ]
    
    
    #clean up
    free(path)
    free(c_probabilities)
    free(pathlist)
    free(cyclelist)
    for k in range(nr_nodes):
        free(c_graph[k])
    free(c_graph)
    free(out_degrees)
    
    return sp.asarray(path_ret), sp.asarray(cycle_ret)

