from scipy import optimize as opt
from copy import deepcopy
import sys as sy
+import cv as cv
+import cv2 as cv2
+import datetime as dt
-
-error = 0.1
+errorq = 0.1
+errorD = 1E-2
epsilon = 1E-10
vrate = 0.8
p = 0.7
#lg= [(0,1,23),(1,0,15),(1,2,45)]
sink = n-1
+def distance(d1,d2):
+ return mt.sqrt(sum([(d1[t]-d2[t])**2 for t in d1]))
+
def genereGraph():
test = False
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
+ return (G,l)
+
+
+def afficheGraph(G,l,tx,ty,sink):
+ r = 20
+ img = cv.CreateImage ((tx, ty), 32, 3)
+ cv.Rectangle(img, (0,0),(tx,ty), cv.Scalar(255,255,255), thickness=-1)
+ def px((x,y)):
+ u = float(tx)/(coteCarre + 2*distanceEmmissionMax)
+ return(int(distanceEmmissionMax*u + x * u),int(distanceEmmissionMax*u + y * u))
+ for i in set(range(len(l)))-set([sink]):
+ (x,y) = l[i]
+ pix,piy = px((x,y))
+ demx = distanceEmmissionMax*tx/(coteCarre+2*distanceEmmissionMax)
+ cv.Circle(img, (pix,piy),demx, cv.Scalar(125,125,125))
+
+ for i in set(range(len(l)))-set([sink]):
+ (x,y) = l[i]
+ pix,piy = px((x,y))
+ cv.Circle(img, (pix,piy),r, cv.Scalar(125,125,125),thickness=-1)
+
+ #sink
+ (x,y) = l[sink]
+ pix,piy = px((x,y))
+
+ cv.Rectangle(img, (pix-r/2,piy-r/2),(pix+r/2,piy+r/2), cv.Scalar(125,125,125), thickness=-1)
+
+ for i in range(len(l)):
+ for j in range(len(l)):
+
+ if np.linalg.norm(np.array(l[i])-np.array(l[j])) < distanceEmmissionMax :
+ (xi,yi) = l[i]
+ pixi,piyi = px((xi,yi))
+ (xj,yj) = l[j]
+ pixj,piyj = px((xj,yj))
+ cv.Line(img, (pixi,piyi), (pixj,piyj), cv.Scalar(125,125,125))
+
+
+ """
+ for i in range(len(l)):
+ (x,y) = l[i]
+ pix,piy = px((x,y))
+ print i
+ textColor = (0, 0, 255) # red
+ font = cv2.FONT_HERSHEY_SIMPLEX
+ imgp =
+ cv2.putText(img, str(i), (pix-r/4,piy-r/2),font, 3.0, textColor)#,thickn """
+ cv.SaveImage("SensorNetwork.png",img)
G = nx.DiGraph()
G.add_weighted_edges_from(lg)
#nx.draw(G)
#pb.show()
+#(G,l) = genereGraph()
+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
+
+
+#afficheGraph(G,l,500,500,sink)
+#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)]
beta = 1.3E-8
gamma = 55.54
delta = 0.2
-zeta = 0.1
amplifieur = 1
sigma2 = 3500
Bi = 5
-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):
#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
+ r= opt.fmin_cobyla(f,xini,cons=[xpos],args=args,consargs=(),rhoend=errorD/100,iprint=0,maxfun=nbiter)
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
+ r = opt.fmin_cobyla(f,xini,[xposounul],args,consargs=(),rhoend=errorD/100,iprint=0,maxfun=nbiter)
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
+ r = opt.fmin_cobyla(f,xini,[entre0et1],args,consargs=(),rhoend=errorD/100,iprint=0,maxfun=nbiter)
return r
-def maj(k,maj_theta,mxg,idxexp):
+def maj(k,maj_theta,mxg,idxexp,comppsh=False):
# initialisation des operateurs lagrangiens
global u, v, la, w, theta , q, Ps, Rh, eta, x, valeurFonctionDuale
qp={}
- def f_q(qi,i):
+ def f_q(qi,*args):
fa = sum([a[i][l]*w[l] for l in L]) - la[i]*Bi
return qi*qi + qi*fa
+ tq,tqa = 0,0
for i in N:
if not ASYNC or random() < taux_succes:
+ n1 = dt.datetime.now()
c = -float(sum([a[i][l]*w[l] for l in L]) - la[i]*Bi)/(2*amplifieur)
+ n2 = dt.datetime.now()
+ #cp = armin(f_q,q[i],POS,tuple([i]))
+ n3 = dt.datetime.now()
+ tq += (n2-n1).microseconds
+ tqa += (n3-n2).microseconds
rep = epsilon if c <= 0 else c
qp[i] = rep
else :
qp[i] = q[i]
-
-
+ tp,tpa = 0,0
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(2)/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))/(16*zeta)
- #print t
- rep = mt.pow(t,float(3)/5)
- Psp[h]=rep
- else :
+ try:
+ #print psh,(gamma*mt.pow(psh,float(2)/3))
+ res= v[h]* mt.log(float(sigma2)/D)/(gamma*mt.pow(psh,float(2)/3)) +la[h]*psh + delta*mt.pow(psh,float(8)/3)
+ return res
+ except :
+ return 0
+ for h in V:
+ if not ASYNC or random() < taux_succes:
+ n1 = dt.datetime.now()
+ #calcul nouvelle solution
+ if comppsh :
+ if la[h] != 0 :
+ t= float(2*v[h]*mt.log(float(sigma2)/D))/(3*gamma*la[h])
+ rep = mt.pow(t,float(3)/5)
+ else :
+ rep = 1000
+ else :
+ t= float(-3*la[h]+mt.sqrt(9*(la[h]**2)+64*delta*v[h]*mt.log(float(sigma2)/D)/gamma))/(16*delta)
+ rep = mt.pow(t,float(3)/5)
+
+
+
+ n2 = dt.datetime.now()
+ #cp = armin(f_Ps,Ps[h],POS,tuple([h]))
+ n3 = dt.datetime.now()
+ tp += (n2-n1).microseconds
+ tpa += (n3-n2).microseconds
+ Psp[h]=rep
+ else :
Psp[h] = Ps[h]
Rhp={}
+ tr,tra=0,0
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)
+ n1 = dt.datetime.now()
+ rep = float(v[h])/(2*delta)
+ n2 = dt.datetime.now()
+ #cp = armin(f_Rh,Rh[h],POS_NUL,tuple([h]))
+ n3 = dt.datetime.now()
+ tr += (n2-n1).microseconds
+ tra += (n3-n2).microseconds
Rhp[h] = 0 if rep < 0 else rep
else :
Rhp[h] = Rh[h]
+
+
+
xp={}
+ tx,txa=0,0
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:
+ n1 = dt.datetime.now()
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)
+ n2 = dt.datetime.now()
+ #cp = armin(f_x,x[(h,l)],POS_NUL,tuple([h,l]))
+ n3 = dt.datetime.now()
+ tx += (n2-n1).microseconds
+ txa += (n3-n2).microseconds
+
xp[(h,l)] = 0 if rep < 0 else rep
else :
xp[(h,l)] = x[(h,l)]
#AfficheVariation(up,vp,lap,wp,thetap,etap,qp,Psp,Rhp,xp,valeurFonctionDualep)
- arret = abs(valeurFonctionDuale-valeurFonctionDualep) < error
+ arret = abs(valeurFonctionDuale-valeurFonctionDualep)
- return (up,vp,lap,wp,thetap,etap,qp,Psp,Rhp,xp,valeurFonctionDualep,arret,mxg,smax)
+ return (up,vp,lap,wp,thetap,etap,qp,Psp,Rhp,xp,valeurFonctionDualep,arret,mxg,smax,tq+tp+tr+tx,tqa+tpa+tra+txa)
def __evalue_maj_theta__(nbexp,out=False):
global u, v, la, w, theta , q, Ps, Rh, eta, x, valeurFonctionDuale
- nbexp = 10
res = {}
m = []
itermax = 100000
-
+
def __maj_theta(k):
mem = []
om = omega/(mt.pow(k,0.75))
return om
+ liste_arret=[]
for idxexp in range(nbexp):
mxg = 0
if not(out):
k = 1
arret = False
sm = 0
+ err, ar = 0,0
+ dur,dura=0,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
+ print "ici"
+ (u,v,la,w,theta,eta,q,Ps,Rh,x,valeurFonctionDuale,ar,mxg,smax,ct,cta)=maj(k,__maj_theta,mxg,idxexp)
+ err = (max(q.values()) - min(q.values()))/min(q.values())
+ dur += ct
+ dura += cta
+ arret = err < errorq or ar < errorD
k+=1
variation = "+" if smax > sm else "-"
- print variation,
- if k%100 ==0 :
- print "k:",k,"erreur sur q", errorq, "et q:",q
- print "maxg=", mxg
+ if out : print variation,
+ if k%10 ==0 :
+ print "k:", k,
+ "erreur sur q",
+ errorq, "erreur sur F", ar, "et q:", q
+
+ if out : 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"
+
+ #print "#########\n",mem,"\#########\n"
if k%4600 == 0 :
- print "#########\n",mem,"\#########\n"
-
-
+ #print "#########m\n",mem,"\#########\n"
+ """
+
if smax - sm > 500:
print "variation trop grande"
print init
exit
sm = smax
-
+
if k == itermax:
print "nbre d'iteration trop grand"
print "init"
print init
sy.exit(1)
+ else :
+ liste_arret += [(err, ar,k,errorD)]
print "###############"
print k
print "###############"
m += [k]
-
- print (min(m),max(m),float(sum(m))/nbexp,m),m
-
+
+ #print liste_arret
+ #print (min(m),max(m),float(sum(m))/nbexp,m),m
+ #return (liste_arret,dur,dura)
+ return (liste_arret,dur)
def __une_seule_exp__(fichier_donees):
ASYNC = False
-__une_seule_exp__("config_initiale.py")
-#__evalue_maj_theta__()
+#__une_seule_exp__("config_initiale.py")
+#pour evaluerle test d'arret
+"""
+listederes=[]
+for k in list(np.logspace(-6, -3, num=15)):
+ errorD = k
+ listederes += __evalue_maj_theta__(5)
+
+print listederes
+"""
+
+
+
+
+def __comparePsh_et_Psh_avec8_3__(nbexp,out=False):
+ print "tttttttttttt"
+ global u, v, la, w, theta , q, Ps, Rh, eta, x, valeurFonctionDuale
+ res = {}
+ m,mp = [],[]
+ itermax = 100000
+
+ def __maj_theta(k):
+ mem = []
+ om = omega/(mt.pow(k,0.75))
+ return om
+ liste_arret=[]
+ for idxexp in range(nbexp):
+ mxg = 0
+ if not(out):
+ initialisation()
+ else :
+ initialisation_()
+
+ k = 1
+ arret = False
+ sm = 0
+ err, ar = 0,0
+ dur,dura=0,0
+
+ # duppliquee
+ while k < itermax and not arret :
+ (u,v,la,w,theta,eta,q,Ps,Rh,x,valeurFonctionDuale,ar,mxg,smax,ct,cta)=maj(k,__maj_theta,mxg,idxexp)
+ err = (max(q.values()) - min(q.values()))/min(q.values())
+ dur += ct
+ dura += cta
+ arret = err < errorq or ar < errorD
+ k+=1
+ variation = "+" if smax > sm else "-"
+ if out : print variation,
+ if k%100 ==0 :
+ print "k:", k,
+ "erreur sur q",
+ errorq, "erreur sur F", ar, "et q:", q
+
+ if out : 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 "#########m\n",mem,"\#########\n"
+ """
+
+
+ if smax - sm > 500:
+ print "variation trop grande"
+ print "init"
+ print init
+ exit
+ sm = smax
+
+ if k == itermax:
+ print "nbre d'iteration trop grand"
+ print "init"
+ print init
+ sy.exit(1)
+ else :
+ liste_arret += [(err, ar,k,errorD)]
+
+ print "###############"
+ print k," avec optym"
+ print "###############"
+ m += [k]
+ ##fin de duplication
+
+
+
+ #remise a zero
+ q,Ps,Rh,eta,x,u,v,la,w=init[0],init[1],init[2],init[3],init[4],init[5],init[6],init[7],init[8]
+ k = 1
+ arret = False
+ sm = 0
+ err, ar = 0,0
+ dur,dura=0,0
+ # duppliquee
+ while k < itermax and not arret :
+ (u,v,la,w,theta,eta,q,Ps,Rh,x,valeurFonctionDuale,ar,mxg,smax,ct,cta)=maj(k,__maj_theta,mxg,idxexp,True)
+ err = (max(q.values()) - min(q.values()))/min(q.values())
+ dur += ct
+ dura += cta
+ arret = err < errorq or ar < errorD
+ k+=1
+ variation = "+" if smax > sm else "-"
+ if out : print variation,
+ if k%100 ==0 :
+ print "k:", k,
+ "erreur sur q",
+ errorq, "erreur sur F", ar, "et q:", q
+
+ if out : 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 "#########m\n",mem,"\#########\n"
+ """
+
+
+ if smax - sm > 500:
+ print "variation trop grande"
+ print "init"
+ print init
+ exit
+ sm = smax
+
+ if k == itermax:
+ print "nbre d'iteration trop grand"
+ print "init"
+ print init
+ sy.exit(1)
+ else :
+ liste_arret += [(err, ar,k,errorD)]
+
+ print "###############"
+ print k,"sans optym"
+ print "###############"
+ mp += [k]
+
+
+ #print liste_arret
+ #print (min(m),max(m),float(sum(m))/nbexp,m),m
+ #return (liste_arret,dur,dura)
+ return (m,mp)
+
+
+
+# pour evaluer argmin
+"""
+listederes=[]
+k=0
+for k in list(np.logspace(-5, -3, num=10)):
+ print "k",k
+ errorD = k
+ listederes += [(k,__evalue_maj_theta__(3))]
+ k+=1
+print listederes
+"""
+
+# pour evaluer psh
+listederes=[]
+for k in list(np.logspace(-5, -3, num=10)):
+ print "Precision",k
+ errorD = k
+ listederes += [(k,__comparePsh_et_Psh_avec8_3__(10))]
+ k+=1
+print listederes
+
+
+
+
+
+
+print"#######################"
+
+
#ASYNC = True
#taux_succes = 0.9
#__evalue_maj_theta__()
-
+"""
print "u",u
print "v",v
print "lambda",la
print "x",x
print "L",valeurFonctionDuale
-
+"""
# relation ente les variables primaires et secondaires ?
