X-Git-Url: https://bilbo.iut-bm.univ-fcomte.fr/and/gitweb/desynchronisation-controle.git/blobdiff_plain/90e98a89726a12ad0e1d68d39c993aa73b3c4230..e6cd9df3d469f7916d513610d3bc22ab055f790a:/exp_controle_asynchrone/simulMWSN.py diff --git a/exp_controle_asynchrone/simulMWSN.py b/exp_controle_asynchrone/simulMWSN.py index 4844593..218ecb7 100644 --- a/exp_controle_asynchrone/simulMWSN.py +++ b/exp_controle_asynchrone/simulMWSN.py @@ -7,14 +7,16 @@ from itertools import * from scipy import optimize as opt from copy import deepcopy import sys as sy +import cv as cv +import cv2 as cv2 - -error = 0.1 +errorq = 0.1 +errorD = 1E-5 epsilon = 1E-10 vrate = 0.8 p = 0.7 coteCarre = 50 -distanceEmmissionMax = 30 +distanceEmmissionMax = 20 nbiter = 1000 POS = 1 POS_NUL = 2 @@ -32,6 +34,9 @@ lg = [(0, 1, 22.004323820359151), (0, 2, 28.750632705280324), (0, 3, 29.68069293 #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 @@ -46,7 +51,54 @@ def genereGraph(): 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) @@ -56,23 +108,27 @@ 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() + +exit() #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())) @@ -146,8 +202,6 @@ def aminp(l): -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): @@ -311,7 +365,7 @@ def maj(k,maj_theta,mxg,idxexp): 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: + for h in V: if not ASYNC or random() < taux_succes: """ lah = 0.05 if la[h] == 0 else la[h] @@ -376,7 +430,7 @@ def maj(k,maj_theta,mxg,idxexp): #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) @@ -476,6 +530,7 @@ def __evalue_maj_theta__(nbexp,out=False): mem = [] om = omega/(mt.pow(k,0.75)) return om + liste_arret=[] for idxexp in range(nbexp): mxg = 0 if not(out): @@ -486,92 +541,130 @@ def __evalue_maj_theta__(nbexp,out=False): k = 1 arret = False sm = 0 + err, ar = 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 + err = (max(q.values()) - min(q.values()))/min(q.values()) + + arret = err < errorq or ar < errorD 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 out : print variation, + if k%500 ==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 "#########\n",mem,"\#########\n" + """ + if smax - sm > 500: print "variation trop grande" print "init" print init - sy.exit(0) + 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 liste_arret + #print (min(m),max(m),float(sum(m))/nbexp,m),m + return liste_arret - print (min(m),max(m),float(sum(m))/nbexp,m),m +def __une_seule_exp__(fichier_donees): + global u, v, la, w, theta , q, Ps, Rh, eta, x, valeurFonctionDuale + initialisation() + + fichier = open(fichier_donees, "r") + instructions ={} + for line in fichier: + l = line.split("=") + instructions[l[0]] = eval(l[1]) + u, v, la, w, q, Ps, Rh, eta, x, = instructions['u'], instructions['v'], instructions['la'], instructions['w'], instructions['q'], instructions['Ps'], instructions['Rh'], instructions['eta'], instructions['x'] + nbexp = 1 + res = {} + m = [] + itermax = 100000 + + def __maj_theta(k): + om = omega/(mt.pow(k,0.75)) + return om + for idxexp in range(nbexp): + mxg = 0 + 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 "-" + 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 + sm = smax + 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__(1,True) +ASYNC = False +#__une_seule_exp__("config_initiale.py") +listederes=[] +for k in list(np.logspace(-6, -3, num=15)): + errorD = k + listederes += __evalue_maj_theta__(5) + +print listederes #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 @@ -584,7 +677,7 @@ print "Rh",Rh print "x",x print "L",valeurFonctionDuale - +""" # relation ente les variables primaires et secondaires ?