]> AND Private Git Repository - desynchronisation-controle.git/blobdiff - exp_controle_asynchrone/simulMWSN.py
Logo AND Algorithmique Numérique Distribuée

Private GIT Repository
beaucoup decourbes
[desynchronisation-controle.git] / exp_controle_asynchrone / simulMWSN.py
index 8bbeddb23f5f6772df28d180efbc2f9706a179c7..218ecb7bc27746b049b3d9b3db3edca1caf48a53 100644 (file)
@@ -6,17 +6,23 @@ import pylab as pb
 from itertools import *
 from scipy import optimize as opt
 from copy import deepcopy 
-error  = 0.1
+import sys as sy
+import cv as cv
+import cv2 as cv2 
+
+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
 POSINF1 = 3
 init = []
+fichier_init="config_initiale_default.txt"
 
 
 
@@ -28,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 
@@ -42,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)
@@ -52,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()))
@@ -103,6 +163,7 @@ alpha = 0.5
 beta = 1.3E-8
 gamma = 55.54
 delta = 0.2
+zeta = 0.1
 amplifieur = 1
 sigma2 = 3500
 Bi = 5
@@ -141,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):
@@ -240,7 +299,7 @@ def maj(k,maj_theta,mxg,idxexp):
     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))
+            s = Rh[h]- mt.log(float(sigma2)/D)/(gamma*mt.pow(Ps[h],float(2)/3))
             if abs(s) > mxg :
                 print "ds calcul v",abs(s),idxexp
                 mxg = abs(s) 
@@ -304,14 +363,20 @@ def maj(k,maj_theta,mxg,idxexp):
     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
+        #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 = (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 :
+         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 :
             Psp[h] = Ps[h]
 
 
@@ -365,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)
 
@@ -442,52 +507,89 @@ def initialisation():
 
 
 
-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 = 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
-        initialisation()
+        if not(out):
+            initialisation()
+        else :
+            initialisation_()
+            
         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 "-"
-            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"
+            if k%4600 == 0 :
+                #print "#########\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
-            exit 
+            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
 
 
 def __une_seule_exp__(fichier_donees):
@@ -548,8 +650,13 @@ def __une_seule_exp__(fichier_donees):
 
 
 ASYNC = False
-__une_seule_exp__("config_initiale.py")
-#__evalue_maj_theta__()
+#__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__()
@@ -557,7 +664,7 @@ __une_seule_exp__("config_initiale.py")
 
 
 
-
+"""
 print "u",u
 print "v",v
 print "lambda",la
@@ -570,7 +677,7 @@ print "Rh",Rh
 print "x",x
 print "L",valeurFonctionDuale
 
-
+"""
 
 # relation ente les variables primaires et secondaires ?