modification de refs

[desynchronisation-controle.git] / exp_controle_asynchrone / simulMWSN.py

diff --git a/exp_controle_asynchrone/simulMWSN.py b/exp_controle_asynchrone/simulMWSN.py
index 218ecb7bc27746b049b3d9b3db3edca1caf48a53..7ca4bd5d22977de781f11709bac26c1c130e6253 100644 (file)
--- a/exp_controle_asynchrone/simulMWSN.py
+++ b/exp_controle_asynchrone/simulMWSN.py
@@ -9,14 +9,15 @@ from copy import deepcopy
 import sys as sy
 import cv as cv
 import cv2 as cv2 
 import sys as sy
 import cv as cv
 import cv2 as cv2 

+import datetime as dt

 
 errorq  = 0.1
 
 errorq  = 0.1

-errorD  = 1E-5
+errorD  = 1E-1

 epsilon = 1E-10
 vrate = 0.8
 p = 0.7
 coteCarre = 50
 epsilon = 1E-10
 vrate = 0.8
 p = 0.7
 coteCarre = 50

-distanceEmmissionMax = 20
+distanceEmmissionMax = 30

 nbiter = 1000
 POS = 1
 POS_NUL = 2
 nbiter = 1000
 POS = 1
 POS_NUL = 2


@@ -108,13 +109,15 @@ G.add_weighted_edges_from(lg)
 #nx.draw(G)
 #pb.show()
 
 #nx.draw(G)
 #pb.show()
 

-(G,l) = genereGraph()
+#(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
 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)
+
+
+#afficheGraph(G,l,500,500,sink)

 #nx.draw(G)
 #pb.show()
 
 #nx.draw(G)
 #pb.show()
 


@@ -122,7 +125,7 @@ afficheGraph(G,l,500,500,sink)
 
 
 
 
 
 

-exit()
+

     
 #print G.edges(data=True)
 #TODO afficher le graphe et etre sur qu'il est connexe
     
 #print G.edges(data=True)
 #TODO afficher le graphe et etre sur qu'il est connexe


@@ -163,7 +166,6 @@ alpha = 0.5
 beta = 1.3E-8
 gamma = 55.54
 delta = 0.2
 beta = 1.3E-8
 gamma = 55.54
 delta = 0.2

-zeta = 0.1

 amplifieur = 1
 sigma2 = 3500
 Bi = 5
 amplifieur = 1
 sigma2 = 3500
 Bi = 5


@@ -241,19 +243,11 @@ def armin(f,xini,xr,args):
     #xpos = lambda x : x if x > 0 else -1 
     r = 0
     if xr == POS :
     #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 :
     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:
     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
 
     
     return r


@@ -344,39 +338,48 @@ def maj(k,maj_theta,mxg,idxexp):
 
 
     qp={}
 
 
     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    
 
         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:
     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)
             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]
 
 
             rep = epsilon if c <= 0 else c
             qp[i] = rep
         else :
             qp[i] = q[i]
 
 

-        

 
 

- 
+    tp,tpa = 0,0

     Psp={}
     Psp={}

-    #print "maj des des Psh" 

     def f_Ps(psh,h):
     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
+        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:
     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 :
+        if not ASYNC or  random() < taux_succes:
+            n1 = dt.datetime.now()
+            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]
 
 
             Psp[h] = Ps[h]
 
 


@@ -390,26 +393,43 @@ def maj(k,maj_theta,mxg,idxexp):
 
 
     Rhp={}
 
 
     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])
     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)
     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]
           
             Rhp[h] = 0 if rep < 0 else rep
         else : 
             Rhp[h] = Rh[h]
           

+
+
+

     xp={}
     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
     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:
     
     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)
                 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)] 
                 xp[(h,l)] = 0 if rep < 0 else rep 
             else :
                 xp[(h,l)] = x[(h,l)] 


@@ -432,7 +452,7 @@ def maj(k,maj_theta,mxg,idxexp):
 
     arret = abs(valeurFonctionDuale-valeurFonctionDualep) 
 
 
     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)

 
 
 
 
 
 


@@ -525,7 +545,7 @@ def __evalue_maj_theta__(nbexp,out=False):
     res = {}
     m = []
     itermax = 100000
     res = {}
     m = []
     itermax = 100000

- 
+

     def __maj_theta(k):
         mem = []
         om = omega/(mt.pow(k,0.75))
     def __maj_theta(k):
         mem = []
         om = omega/(mt.pow(k,0.75))


@@ -542,15 +562,18 @@ def __evalue_maj_theta__(nbexp,out=False):
         arret = False
         sm = 0
         err, ar = 0,0
         arret = False
         sm = 0
         err, ar = 0,0

+        dur,dura=0,0

         while k < itermax  and not arret : 
         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)
+            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())
             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,
             arret = err <  errorq or ar < errorD
             k+=1
             variation = "+" if smax > sm else "-"
             if out : print variation,

-            if k%500 ==0 :
+            if k%10 ==0 :

                 print "k:", k, 
                 "erreur sur q", 
                 errorq, "erreur sur F", ar, "et q:", q
                 print "k:", k, 
                 "erreur sur q", 
                 errorq, "erreur sur F", ar, "et q:", q


@@ -563,7 +586,7 @@ def __evalue_maj_theta__(nbexp,out=False):
                 
                 #print "#########\n",mem,"\#########\n"
             if k%4600 == 0 :
                 
                 #print "#########\n",mem,"\#########\n"
             if k%4600 == 0 :

-                #print "#########\n",mem,"\#########\n"
+                #print "#########m\n",mem,"\#########\n"

             """
                 
 
             """
                 
 


@@ -589,7 +612,7 @@ def __evalue_maj_theta__(nbexp,out=False):
  
     #print liste_arret    
     #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
+    return (liste_arret,dur,dura)

 
 
 def __une_seule_exp__(fichier_donees):
 
 
 def __une_seule_exp__(fichier_donees):


@@ -651,12 +674,33 @@ def __une_seule_exp__(fichier_donees):
 
 ASYNC = False
 #__une_seule_exp__("config_initiale.py")
 
 ASYNC = False
 #__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
 listederes=[]
 for k in list(np.logspace(-6, -3, num=15)):
     errorD = k
     listederes += __evalue_maj_theta__(5)
     
 print listederes

+"""
+
+
+
+# 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
+
+
+
+print"#######################"
+
+

 #ASYNC = True
 #taux_succes = 0.9
 #__evalue_maj_theta__()
 #ASYNC = True
 #taux_succes = 0.9
 #__evalue_maj_theta__()