#!/usr/bin/env python
# -*- coding: iso-8859-15 -*-
"""
ReTrace - find branching pathways in metabolic networks
Copyright (C) 2009 Esa Pitkдnen
See file COPYING for license notice.

Implementation of Yen's k shortest simple (loopless) paths.

Running time is O(kn(m + n log n)) in a graph with n nodes
and m edges.

"""

#try:
#    import psyco # try to enable psyco JIT
#    psyco.full()
#except:
#    pass

import sys, os, random, time, math

from graph import randomGraph, graphGenerators, graph
from eppstein import dijkstra, priodict

class Path(list):
    def __hash__(self):
        h = ""
        for e in self.__iter__():
            h += "%s*" % (e)
        return hash(h)

def getPath(s, t, P, G):
    cost = 0.0
    n = t
    p = Path()

    if len(P) == 0 or t not in P:
        return p, 0

    while n != s:
        cost += G.E[P[n]][n]
        p.append(n)
        n = P[n]
    p.append(n)
    p.reverse()
    return p, cost
import heapq
def kspSimpleYen(G, K, s, t):
    assert(K > 0)
    X = priodict.priorityDictionary()  # paths
    X = []
    removals = {}

    found = set()
    pathways = []

    (D, P) = dijkstra.Dijkstra(G.E, s, t)
    p, cost = getPath(s, t, P, G)
    p.cost = cost
    p.dix = 0
    #X[p] = cost
    heapq.heappush(X,(p,cost))
    removals[p] = []
    found.add(p)

    k = 0
    I = 0
    #for p in X:
    while X:
        p, _cost = heapq.heappop(X)
        R = removals.pop(p)

        pathways.append(p)

        k += 1
        if k == K:
            break

        # Add to removals every edge from nodes v_{l-1}, l = |p|
        # on the path p, starting from the end of the path

        for (u, v, c) in R:
            del G.E[u][v]

        R2 = []
        i = p.dix
        while i < len(p) - 1:
            I += 1
            # Delete the next edge on the path
            pi = p[i]
            pii = p[i + 1]
            R2.append((pi, pii, G.E[pi][pii]))
            del G.E[pi][pii]

            (D, P) = dijkstra.Dijkstra(G.E, pi, t)
            p2, cost = getPath(pi, t, P, G)

            if len(p2) > 0:
                p3 = Path(p[0:i])
                p3.extend(p2)

                if p3 not in found:
                    p3.cost = p.cost + cost
                    p3.dix = i
                    X[p3] = p3.cost
                    removals[p3] = list(R2)
                    found.add(p3)

            R3 = []
            for u in G.E[pi]:
                if u != pii:
                    R3.append((pi, u, G.E[pi][u]))

            for (u, v, c) in R3:
                del G.E[u][v]
            R2.extend(R3)
            i += 1

        # Restore graph

        for (u, v, c) in R2:
            G.E[u][v] = c

        for (u, v, c) in R:
            G.E[u][v] = c

    print 'I =',I
    return pathways

def testKSP():
    G = graphGenerators.martinsPascoalExample()
    s = "N1"
    t = "N5"
    P = kspSimpleYen(G, 100, s, t)
    for p in P:
        print p

def meansd(L):
    mean = 1.0 * reduce(lambda x, y: x + y, L) / len(L)
    sum = 0.0
    for x in L:
        sum += (x - mean) * (x - mean)
    sd = math.sqrt(sum / len(L))
    return mean, sd


def testKSP2():
    edgetonoderatio = 1.5

    krange = range(50, 55, 5)
    nrange = range(200, 1000, 10)
    rounds = 1
    rrange = range(rounds)

    o = sys.stdout

    s = 1
    t = 2

    o.flush()
    for k in krange:
        for n in nrange:

            o.write("%d\t%d\t" % (k, rounds))
            o.flush()
            m = int(edgetonoderatio * n)
            G = randomGraph.randomGraph(n, m)
            o.write("%d\t%d\t" % (G.numNodes(), G.numEdges()))
            o.flush()
            times = []

            for r in rrange:
                V = list(G.V)
                s = random.choice(V)
                t = random.choice(V)

                st = time.time()

                paths = kspSimpleYen(G, k, s, t)

                et = time.time()

                times.append(et - st)

            mean, sd = meansd(times)
            o.write("%f\t%f\t%d\n" % (mean, sd, len(paths)))
            o.flush()

def removeCyclicPaths(P):
    Q = []
    for p in P:
        nodes = set()
        cyclic = False
        for n in p:
            if n in nodes:
                cyclic = True
                break
            else:
                nodes.add(n)
        if not cyclic:
            Q.append(p)
    return Q


G =  {0: {1: 4, 4: 4, 5: 254, },
1: {2: 4, 4: 254, 5: 4, },
2: {3: 4, 5: 254, 6: 5, 7: 254, },
3: {6: 254, 7: 6, },
4: {5: 4, 8: 4, 9: 254, },
5: {6: 4, 8: 254, 9: 4, },
6: {7: 4, 9: 254, 10: 4, },
7: {11: 5, },
8: {9: 4, 12: 4, 13: 254, },
9: {10: 4, 12: 254, 13: 4, },
10: {11: 4, 13: 254, },
11: {15: 5, },
12: {13: 4, },
13: {},
15: {},
}
def cost(G, p):
    return sum([G1[x[0]][x[1]] for x in zip(p,p[1:])])
def triangular2usual(G):
   G1={}
   for v1 in G.iterkeys():
      for (v2, cost2) in G[v1].iteritems():
         if not G1.has_key(v1): G1[v1]={}
         if not G1.has_key(v2): G1[v2]={}
         G1[v1][v2]=cost2
         G1[v2][v1]=cost2
   return G1
import itertools
def path_overlay(p):
   allp = set(reduce(lambda x, y: x+y, p))
   cnt = 0.0;
   for v in allp:
      for pp in p:
         if pp.count(v): cnt += 1
   return cnt/len(allp)

def overlays(G, paths, n):
   for pp in paths:
      #pp = [int(x) for x in p.split(',')]
      s=0
      for v in range(len(pp)-1):
         s += G[pp[v]][pp[v+1]]
      print pp,s
   l_overlay = [(path_overlay(c),sum([cost(G,p) for p in c]),c) for c in itertools.combinations(paths,n)]
   m1,m2 = max([x[0] for x in l_overlay]), float(max([x[1] for x in l_overlay]))
   l_overlay = [(x[0]/m1, x[1]/m2, x[2]) for x in l_overlay]
   l_overlay.sort(lambda x,y: cmp(x[0]**2+x[1]**2,y[0]**2+y[1]**2))
   print 'paths:', len(paths)
   print 'overlays:(',len(l_overlay),')', l_overlay[0]

if __name__ == "__main__":
    #testKSP2()
    G1 = triangular2usual(G)
    G2 = graph.Graph()
    for v1 in G1:
        G2.addNode(v1)
    for v1 in G1:
        for v2 in G1[v1]:
            G2.addEdge(v1,v2,G1[v1][v2])
    paths = kspSimpleYen(G2,20,0,15)
    sp = []
    for p in paths:
        vp = zip(p,p[1:])
        s = cost(G,p)
        sp.append((s,tuple(p)))
    #sp.sort(lambda x,y: cmp(x[0],y[0]))
    for _sp in sp:
        print _sp
    overlays(G1,paths,4)