+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
+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 = []
+
+
+
+# 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
+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 = (float(2*v[h]*mt.log(float(sigma2)/D))/mt.pow(3*gamma*lah,float(3)/5))
+ Psp[h] = epsilon if rep <= 0 else rep
+ 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 __evalue_maj_theta__():
+ global u, v, la, w, theta , q, Ps, Rh, eta, x, valeurFonctionDuale
+ nbexp = 5
+ res = {}
+ m = []
+ itermax = 10000
+
+ def __maj_theta(k):
+ om = omega/(mt.pow(k,0.75))
+ return om
+ for idxexp in range(nbexp):
+ mxg = 0
+ 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
+ if smax - sm > 500:
+ print "variation trop grande"
+ print "init"
+ print init
+ exit
+ if k == itermax:
+ print "nbre d'iteration trop grand"
+ print "init"
+ print init
+ exit
+
+ print "###############"
+ print k
+ print "###############"
+ m += [k]
+
+ print (min(m),max(m),float(sum(m))/nbexp,m),m
+
+
+
+
+
+ASYNC = False
+__evalue_maj_theta__()
+#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 ?
+
