import math as mt
from random import *
import networkx as nx 
import numpy as np
import pylab as pb
from itertools import *
from scipy import optimize as opt
from copy import deepcopy 
import sys as sy


error  = 0.1
epsilon = 1E-10
vrate = 0.8
p = 0.7
coteCarre = 50
distanceEmmissionMax = 30
nbiter = 1000
POS = 1
POS_NUL = 2
POSINF1 = 3
init = []
fichier_init="config_initiale_default.txt"



# construction du  graphe 
n = 10
lg = [(0, 1, 22.004323820359151), (0, 2, 28.750632705280324), (0, 3, 29.68069293796183), (0, 4, 8.547146256271331), (0, 5, 28.079672647730469), (0, 7, 23.017867703525138), (0, 8, 6.1268526078857208), (0, 9, 24.573433868296771), (1, 0, 22.004323820359151), (1, 2, 18.807277287689722), (1, 3, 18.982897767602783), (1, 4, 16.848855991756174), (1, 5, 17.042671653231526), (1, 6, 16.410544777532913), (1, 7, 25.598808236367063), (1, 8, 20.175759189503321), (1, 9, 12.843763853990932), (2, 0, 28.750632705280324), (2, 1, 18.807277287689722), (2, 3, 1.0957062702237066), (2, 4, 29.159997765424084), (2, 5, 1.8557839901886808), (2, 6, 23.122260476726876), (2, 9, 6.052562826627808), (3, 0, 29.68069293796183), (3, 1, 18.982897767602783), (3, 2, 1.0957062702237066), (3, 4, 29.884008054261855), (3, 5, 1.9922790489539697), (3, 6, 22.479228556182363), (3, 9, 6.4359869969688059), (4, 0, 8.547146256271331), (4, 1, 16.848855991756174), (4, 2, 29.159997765424084), (4, 3, 29.884008054261855), (4, 5, 28.006189408396626), (4, 7, 15.774839848636024), (4, 8, 3.6206480052249144), (4, 9, 23.804744370383144), (5, 0, 28.079672647730469), (5, 1, 17.042671653231526), (5, 2, 1.8557839901886808), (5, 3, 1.9922790489539697), (5, 4, 28.006189408396626), (5, 6, 21.492976178079076), (5, 8, 29.977996181215822), (5, 9, 4.4452006140146185), (6, 1, 16.410544777532913), (6, 2, 23.122260476726876), (6, 3, 22.479228556182363), (6, 5, 21.492976178079076), (6, 9, 20.04488615603487), (7, 0, 23.017867703525138), (7, 1, 25.598808236367063), (7, 4, 15.774839848636024), (7, 8, 16.915923579829141), (8, 0, 6.1268526078857208), (8, 1, 20.175759189503321), (8, 4, 3.6206480052249144), (8, 5, 29.977996181215822), (8, 7, 16.915923579829141), (8, 9, 25.962918470558208), (9, 0, 24.573433868296771), (9, 1, 12.843763853990932), (9, 2, 6.052562826627808), (9, 3, 6.4359869969688059), (9, 4, 23.804744370383144), (9, 5, 4.4452006140146185), (9, 6, 20.04488615603487), (9, 8, 25.962918470558208)]

#n=3
#lg= [(0,1,23),(1,0,15),(1,2,45)]
sink = n-1


def genereGraph():
    test = False 
    while not test :
        G = nx.DiGraph()
        l = [(random()*coteCarre, random()*coteCarre) for _  in range(n)]
        for io in range(len(l)) :
            for ie in range(len(l)) :
                if ie != io  : 
                    dist = mt.sqrt((l[io][0]-l[ie][0])**2 + (l[io][1]-l[ie][1])**2)
                    if dist <= distanceEmmissionMax :
                        G.add_edge(io,ie,weight=dist)
                        G.add_edge(ie,io,weight=dist)
        test = not(any([ not(nx.has_path(G,o,sink)) for o in  G.nodes() if sink in G.nodes() and o != sink]))
    return G

G = nx.DiGraph()
G.add_weighted_edges_from(lg)
#print nx.is_strongly_connected(G)


#nx.draw(G)
#pb.show()





    
#print G.edges(data=True)
#TODO afficher le graphe et etre sur qu'il est connexe

N = G.nodes()


#V = list(set(sample(N,int(len(N)*vrate)))-set([sink]))
V = list(set(N)-set([sink]))


source = V
print "source",source


L = range(len(G.edges()))
d = [di['weight'] for (_,_,di) in G.edges(data=True)]

#print "L",L
#print "d",d


def ail(i,l):
    assert  l in L, " pb de dimennsion de l: "+str(l)+ " "+ str(L)
    r = 0
    (o,e) = G.edges()[l]
    if o == i :
        r = 1 
    elif e == i :
        r = -1
    return r



# Constantes 
a = [[ail(i,l) for l in L ] for i in xrange(n)]
aplus  = [[1 if  ail(i,l) > 0 else 0  for l in L ] for i in xrange(n)]
amoins  = [[1 if  ail(i,l) < 0 else 0  for l in L ] for i in xrange(n)]





alpha = 0.5
beta = 1.3E-8
gamma = 55.54
delta = 0.2
zeta = 0.1
amplifieur = 1
sigma2 = 3500
Bi = 5
omega = 0.15
D = 100
path_loss_exp = 4
cr = 0.5
cs = [alpha + beta*(di**path_loss_exp) for di in d]
#print "cs",cs


#TODO 
def etahi(h,i,R):
    ret = 0
    if i == h :
        ret = R[h]
    elif i == sink  :
        ret = -R[h]
    return ret


def psi(Ps,i):
    if i not in Ps:
        return 0
    else :
        return Ps[i]
    


def cmpPair(x1,x2):
    return cmp(x1[1],x2[1])

def aminp(l):
    l.sort(cmp=cmpPair)
    return l[0][0]



def distance(d1,d2):
    return mt.sqrt(sum([(d1[t]-d2[t])**2 for t in d1]))


def AfficheVariation (up,vp,lap,wp,thetap,etap,qp,Psp,Rhp,xp,valeurFonctionDualep):
    global u, v, la, w, theta , q, Ps, Rh, eta, x,valeurFonctionDuale
    
    print "du=",distance(u,up),
    print "dv=",distance(v,vp),
    print "dlambda=",distance(la,lap),
    print "dw=",distance(w,wp),
    print "dtheta=",abs(theta-thetap),
    print "deta=",distance(eta,etap),
    print "dq=",distance(q,qp),
    print "dPs=",distance(Ps,Psp),
    print "dRh=",distance(Rh,Rhp),
    print "dx=",distance(x,xp),
    print "dL=",abs(valeurFonctionDuale-valeurFonctionDualep),"\n"


valeurFonctionDuale = 0    

def entre0et1(x):
    r = x if (x >0 and x <= 1) else -1
    #print "ds entre0et1 x et r", x,r 
    return r

    
def xpos(x):
    r = x if x >0 else -1
    #print "ds xpos x et r", x,r 
    return r

def xposounul(x):
    return x if x >= 0 else -1

#f_q,q[i],POS,[i]
def armin(f,xini,xr,args):
    #xpos = lambda x : x if x > 0 else -1 
    r = 0
    if xr == POS :
        #print "strictement pos"
        #print "parametres passes a cobyla",xini,xpos,args,"consargs=(),rhoend=1E-5,iprint=0,maxfun=1000"
        r= opt.fmin_cobyla(f,xini,cons=[xpos],args=args,consargs=(),rhoend=1E-3,iprint=0,maxfun=nbiter)
        #print "le min str pos est",r 
        #print "le min pos   est", r,"avec xini",xini
    elif xr == POS_NUL :
        #print "pos ou nul"
        r = opt.fmin_cobyla(f,xini,[xposounul],args,consargs=(),rhoend=1E-3,iprint=0,maxfun=nbiter)
        # print "le min pos est", r
        #print "le min pos null  est", r,"avec xini",xini
    elif xr == POSINF1:
        r = opt.fmin_cobyla(f,xini,[entre0et1],args,consargs=(),rhoend=1E-3,iprint=0,maxfun=nbiter)
        #print "le min pos inf 1 est", r,"avec xini",xini    

    
    return r






def maj_theta(k):
    return omega/(mt.pow(k,0.5))



def maj_theta(k):
    return omega/mt.sqrt(k)




def maj(k,maj_theta,mxg,idxexp):
    # initialisation des operateurs lagrangiens
    global u, v, la, w, theta , q,  Ps, Rh,  eta, x, valeurFonctionDuale
    
    smax = 0
    #print "iteration",k
    
    up = {}        
    for h in V:
        for i in N:
            if not ASYNC or  random() < taux_succes:
                s = eta[(h,i)]-sum([a[i][l]*x[(h,l)] for l in L])
                if abs(s) > mxg :
                    print "ds calcul u",abs(s),idxexp
                    mxg = abs(s)                 
                smax = max(smax,abs(s))
                up[(h,i)] = u[(h,i)]-theta*s
            else :
                up[(h,i)] = u[(h,i)]


    
    vp = {}
    for h in V:
        if not ASYNC or  random() < taux_succes:
            s = Rh[h]- mt.log(float(sigma2)/D)/(gamma*mt.pow(Ps[h],float(1)/3))
            if abs(s) > mxg :
                print "ds calcul v",abs(s),idxexp
                mxg = abs(s) 
            smax = max(smax,abs(s))
            vp[h] = max(0,v[h]-theta*s)
        else :
            vp[h] = v[h]



    lap={}
    for i in N:
        if not ASYNC or  random() < taux_succes:
            s = q[i]*Bi -sum([aplus[i][l]*cs[l]*sum([x[(h,l)] for h in V]) for l in L])-cr*sum([amoins[i][l]*sum([x[(h,l)] for h in V]) for l in L])-psi(Ps,i)
            if abs(s) > mxg :
                print "ds calcul la",abs(s),idxexp,i
                mxg = abs(s) 
            smax = max(smax,abs(s))
            resla = la[i]-theta*s
            lap[i] = max(0,resla) 
        else :
            lap[i] = la[i]



                        
    wp={}
    for l in L:
        if not ASYNC or  random() < taux_succes:
            s = sum([a[i][l]*q[i] for i in N])
            if abs(s) > mxg :
                print "ds calcul w",abs(s),idxexp
                mxg = abs(s) 
            smax = max(smax,abs(s))
            wp[l] = w[l] + theta*s 
        else : 
            wp[l] = w[l]

        

    thetap = maj_theta(k)


    qp={}
    def f_q(qi,i):
        fa = sum([a[i][l]*w[l] for l in L]) - la[i]*Bi 
        return qi*qi + qi*fa    

    for i in N:
        if not ASYNC or  random() < taux_succes:
            c = -float(sum([a[i][l]*w[l] for l in L]) - la[i]*Bi)/(2*amplifieur)
            rep = epsilon if c <= 0 else c
            qp[i] = rep
        else :
            qp[i] = q[i]


        

 
    Psp={}
    #print "maj des des Psh" 
    def f_Ps(psh,h):
        #print "ds f_ps",psh, v[h]* mt.log(float(sigma2)/D)/(gamma*((psh**2)**(float(1)/3))) +la[h]*psh
        return v[h]* mt.log(float(sigma2)/D)/(gamma*mt.pow(float(2)/3)) +la[h]*psh
    for h in V:
        if not ASYNC or  random() < taux_succes:
            
            lah = 0.05 if la[h] == 0 else  la[h]
            rep = mt.pow(float(2*v[h]*mt.log(float(sigma2)/D))/(3*gamma*lah),float(3)/5)
            Psp[h] = epsilon if rep <= 0 else rep
            """
            t= float(-3*la[h]+mt.sqrt(9*la[h]**2+64*zeta*v[h]*mt.log(float(sigma2)/D)/gamma))/float(16*zeta)
            #print t
            rep = mt.pow(t,float(3)/5)
            Psp[h]=rep
            """
############reprendre ici

        else :
            Psp[h] = Ps[h]



    etap={}

    for h in V:
        for i in N :
            etap[(h,i)] = etahi(h,i,Rh)
            


    Rhp={}
    def f_Rh(rh,h):
        return delta*rh*rh-v[h]*rh-sum([u[(h,i)]*eta[(h,i)] for i in N])

    for h in V:
        if not ASYNC or  random() < taux_succes:
            rep = float(v[h])/(2*delta)
            Rhp[h] = 0 if rep < 0 else rep
        else : 
            Rhp[h] = Rh[h]
          
    xp={}
    def f_x(xhl,h,l):
        r =  delta*xhl*xhl + xhl*(cs[l]*sum([la[i]*aplus[i][l] for i in N]) +cr*sum([la[i]*amoins[i][l] for i in N])+sum([u[(h,i)]*a[i][l] for i in N]))
        return r

    
    for h in V:
        for l in L:
            if not ASYNC or  random() < taux_succes:
                rep = -float(cs[l]*sum([la[i]*aplus[i][l] for i in N]) +cr*sum([la[i]*amoins[i][l] for i in N])+sum([u[(h,i)]*a[i][l] for i in N]))/(2*delta)
                xp[(h,l)] = 0 if rep < 0 else rep 
            else :
                xp[(h,l)] = x[(h,l)] 





    valeurFonctionDualep = 0
    valeurFonctionDualep += sum([amplifieur*q[i]*q[i] for i in N])  
    valeurFonctionDualep += sum([sum([delta*(x[(h,l)]**2) for l in L]) for h in V]) 
    valeurFonctionDualep += sum([delta*(Rh[h]**2) for h in V])  
    valeurFonctionDualep += sum([sum([u[(h,i)]*(sum([ a[i][l]*x[(h,l)] for l in L])- eta[(h,i)]) for i in N]) for h in V])
    valeurFonctionDualep += sum([v[h]*(mt.log(float(sigma2)/D)/(gamma*mt.pow(Ps[h],float(2)/3)) - Rh[h]) for h in V]) 
    valeurFonctionDualep += sum([la[i]*(psi(Ps,i) +sum([aplus[i][l]*cs[l]*sum([x[(h,l)] for h in V]) for l in L])+ cr*sum([ amoins[i][l]*sum([x[(h,l)] for h in V]) for l in L]) -q[i]*Bi)  for i in N ]) 
    valeurFonctionDualep += sum([w[l]*sum([a[i][l]*q[i] for i in N]) for l in L])
    

    #AfficheVariation(up,vp,lap,wp,thetap,etap,qp,Psp,Rhp,xp,valeurFonctionDualep)

    arret = abs(valeurFonctionDuale-valeurFonctionDualep) < error

    return (up,vp,lap,wp,thetap,etap,qp,Psp,Rhp,xp,valeurFonctionDualep,arret,mxg,smax)




"""
print "u",u
print "v",v
print "lambda",la
print "w",w
print "theta",theta
print "eta", eta
print "q",q
print "Ps",Ps
print "Rh",Rh
print "x",x

"""



def initialisation():
    global u, v, la, w, theta , q,  Ps, Rh,  eta, x,init 

    theta = omega

    q = {}
    for i in N :
        q[i] = 0.15 + random()*0.05


    Ps= {}
    for h in V :
        Ps[h] =  0.2+random()*0.3  

    Rh = {}
    for vi in V:
        Rh[vi] = 0.1 + random()*0.1

    eta = {}
    for h in V :
        for i in N:
            eta[(h,i)]= etahi(h,i,Rh)

    x={}
    for h in V :
        for l in L:
            x[(h,l)]=0


 

    # initialisation des operateurs lagrangiens
    u={}
    for h in V :
        for i in N:
            u[(h,i)]= random()

    v = {}
    for h in V :
        v[h] = Rh[h]

    la = {}
    for i in N:
        la[i] = random()

    w={}
    for l in L:
        w[l] =  random()

    init = [deepcopy(q),deepcopy(Ps),deepcopy(Rh),deepcopy(eta),
            deepcopy(x),deepcopy(u),deepcopy(v),deepcopy(la),deepcopy(w)]



def initialisation_():
    global u, v, la, w, theta , q,  Ps, Rh,  eta, x,init 
    fd = open(fichier_init,"r")
    l= fd.readline()
    init_p = eval(l)
    print init_p
    theta = omega
    (q,Ps,Rh,eta,x,u,v,la,w) = tuple([deepcopy(x) for x in init_p])
    init = [deepcopy(q),deepcopy(Ps),deepcopy(Rh),deepcopy(eta),
            deepcopy(x),deepcopy(u),deepcopy(v),deepcopy(la),deepcopy(w)]



def __evalue_maj_theta__(nbexp,out=False):
    global u, v, la, w, theta , q,  Ps, Rh,  eta, x, valeurFonctionDuale 
    res = {}
    m = []
    itermax = 100000
 
    def __maj_theta(k):
        mem = []
        om = omega/(mt.pow(k,0.75))
        return om
    for idxexp in range(nbexp):
        mxg = 0
        if not(out):
            initialisation()
        else :
            initialisation_()
            
        k = 1
        arret = False
        sm = 0
        while k < itermax  and not arret : 
            (u,v,la,w,theta,eta,q,Ps,Rh,x,valeurFonctionDuale,ar,mxg,smax)=maj(k,__maj_theta,mxg,idxexp)
            errorq =  (max(q.values()) - min(q.values()))/min(q.values())
            arret = errorq <  error
            k+=1
            variation = "+" if smax > sm else "-"
            sm = smax
            print variation,
            if k%100 ==0 :
                print "k:",k,"erreur sur q", errorq, "et q:",q
                print "maxg=", mxg
                mem = [deepcopy(q),deepcopy(Ps),deepcopy(Rh),deepcopy(eta),
                       deepcopy(x),deepcopy(u),deepcopy(v),deepcopy(la),deepcopy(w)]
            if k%4500 == 0 :
                print "#########\n",mem,"\#########\n"
            if k%4600 == 0 :
                print "#########\n",mem,"\#########\n"



            if smax - sm  > 500:
                print "variation trop grande"
                print "init"
                print init
                sy.exit(0)
        if k == itermax:
            print "nbre d'iteration trop grand"
            print "init"
            print init
            sy.exit(1)

        print "###############"
        print k
        print "###############"
        m += [k]

    print (min(m),max(m),float(sum(m))/nbexp,m),m





ASYNC = False
__evalue_maj_theta__(1,True)



#ASYNC = True
#taux_succes = 0.9
#__evalue_maj_theta__()


"""
initialisation()
k = 1
arret = False
while k < 10000  and not arret : 
    (u,v,la,w,theta,eta,q,Ps,Rh,x,valeurFonctionDuale,ar)=maj(k,maj_theta)
    arret = ar
    k+=1
    errorq =  abs(min(q.values())-max(q.values()))
    print "errorq",errorq
    arret = errorq <  error
"""

"""
    print "u",u
    print "w",w
    print "theta",theta
    print "eta", eta
    print "q",q
    print "v",v
    print "lambda",la
    print "Ps",Ps
    print "Rh",Rh
    print "x",x
"""


"""
    k +=1

"""



print "u",u
print "v",v
print "lambda",la
print "w",w
print "theta",theta
print "eta", eta
print "q",q
print "Ps",Ps
print "Rh",Rh
print "x",x
print "L",valeurFonctionDuale



# relation ente les variables primaires et secondaires ?