import pylab as pb
from itertools import *
from scipy import optimize as opt
+from copy import deepcopy
+import sys as sy
-error = 0.01
+
+error = 0.1
epsilon = 1E-10
vrate = 0.8
p = 0.7
POS = 1
POS_NUL = 2
POSINF1 = 3
-
+init = []
+fichier_init="config_initiale_default.txt"
beta = 1.3E-8
gamma = 55.54
delta = 0.2
+zeta = 0.1
amplifieur = 1
sigma2 = 3500
Bi = 5
-def maj(k,maj_theta):
+def maj(k,maj_theta,mxg,idxexp):
# initialisation des operateurs lagrangiens
- global u, v, la, w, theta , q, Ps, Rh, eta, x, valeurFonctionDuale
-
-
+ 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_perte:
- up[(h,i)] = u[(h,i)]-theta*(eta[(h,i)]-sum([a[i][l]*x[(h,l)] for l in L]))
+ 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:
- #print "vp",gamma*mt.pow(Ps[h],float(1)/3)
- #vp[h] = max(0,v[h]-theta*(Rh[h]- mt.log(float(sigma2)/D)/(gamma*mt.pow(Ps[h],float(1)/3))))
- if not ASYNC or random() < taux_perte:
- vp[h] = max(0,v[h]-theta*(Rh[h]- mt.log(float(sigma2)/D)/(gamma*mt.pow(Ps[h],float(1)/3))))
+ 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_perte:
- resla = la[i]-theta*((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 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_perte:
- wp[l] = w[l] + theta*(sum([a[i][l]*q[i] for i in N]))
+ 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)
- #print "k",k,"maj theta", thetap
-
qp={}
def f_q(qi,i):
return qi*qi + qi*fa
for i in N:
- if not ASYNC or random() < taux_perte:
+ 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
#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_perte:
+ 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))
+ 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]
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_perte:
+ if not ASYNC or random() < taux_succes:
rep = float(v[h])/(2*delta)
Rhp[h] = 0 if rep < 0 else rep
else :
for h in V:
for l in L:
- if not ASYNC or random() < taux_perte:
+ 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 :
arret = abs(valeurFonctionDuale-valeurFonctionDualep) < error
- return (up,vp,lap,wp,thetap,etap,qp,Psp,Rhp,xp,valeurFonctionDualep,arret)
+ return (up,vp,lap,wp,thetap,etap,qp,Psp,Rhp,xp,valeurFonctionDualep,arret,mxg,smax)
def initialisation():
- global u, v, la, w, theta , q, Ps, Rh, eta, x, valeurFonctionDuale
+ global u, v, la, w, theta , q, Ps, Rh, eta, x,init
theta = omega
-
- # TODO Ahmed : initialiser avec qi >0 et somme (ail.qi) = 0 Ei/i = Ti
- #q = [0 for i in N]
q = {}
for i in N :
q[i] = 0.15 + random()*0.05
-
-
-
- #TODO Ahmed doit donner les valeurs initiales pour Ps
Ps= {}
for h in V :
Ps[h] = 0.2+random()*0.3
-
-
-
- #TODO Ahmed doit donner les valeurs initiales pour Rh
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)
-
-
- # Ok pour 0
x={}
for h in V :
for l in L:
- x[(h,l)]=0#random()
+ x[(h,l)]=0
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__():
+
+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
- nbexp = 50
res = {}
m = []
+ itermax = 100000
+
def __maj_theta(k):
+ mem = []
om = omega/(mt.pow(k,0.75))
return om
for idxexp in range(nbexp):
-
- initialisation()
+ mxg = 0
+ if not(out):
+ initialisation()
+ else :
+ initialisation_()
+
k = 1
arret = False
- while k < 1000000 and not arret :
- (u,v,la,w,theta,eta,q,Ps,Rh,x,valeurFonctionDuale,ar)=maj(k,__maj_theta)
- errorq = abs(min(q.values())-max(q.values()))
+ 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
- if k%5000 ==0 :
+ 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()
-ASYNC = True
-taux_perte = 0.9
-__evalue_maj_theta__()
+ASYNC = False
+__evalue_maj_theta__(1,True)
+
+
+
+#ASYNC = True
+#taux_succes = 0.9
+#__evalue_maj_theta__()
"""
