exp_controle_asynchrone/simulMWSN.py

   1 import math as mt
   2 from random import *
   3 import networkx as nx
   4 import numpy as np
   5 import pylab as pb
   6 from itertools import *
   7 from scipy import optimize as opt
   8 from copy import deepcopy
   9 import sys as sy
  10
  11
  12 error  = 0.1
  13 epsilon = 1E-10
  14 vrate = 0.8
  15 p = 0.7
  16 coteCarre = 50
  17 distanceEmmissionMax = 30
  18 nbiter = 1000
  19 POS = 1
  20 POS_NUL = 2
  21 POSINF1 = 3
  22 init = []
  23 fichier_init="config_initiale_default.txt"
  24
  25
  26
  27 # construction du  graphe
  28 n = 10
  29 lg = [(0, 1, 22.004323820359151), (0, 2, 28.750632705280324), (0, 3, 29.68069293796183), (0, 4, 8.547146256271331), (0, 5, 28.079672647730469), (0, 7, 23.017867703525138), (0, 8, 6.1268526078857208), (0, 9, 24.573433868296771), (1, 0, 22.004323820359151), (1, 2, 18.807277287689722), (1, 3, 18.982897767602783), (1, 4, 16.848855991756174), (1, 5, 17.042671653231526), (1, 6, 16.410544777532913), (1, 7, 25.598808236367063), (1, 8, 20.175759189503321), (1, 9, 12.843763853990932), (2, 0, 28.750632705280324), (2, 1, 18.807277287689722), (2, 3, 1.0957062702237066), (2, 4, 29.159997765424084), (2, 5, 1.8557839901886808), (2, 6, 23.122260476726876), (2, 9, 6.052562826627808), (3, 0, 29.68069293796183), (3, 1, 18.982897767602783), (3, 2, 1.0957062702237066), (3, 4, 29.884008054261855), (3, 5, 1.9922790489539697), (3, 6, 22.479228556182363), (3, 9, 6.4359869969688059), (4, 0, 8.547146256271331), (4, 1, 16.848855991756174), (4, 2, 29.159997765424084), (4, 3, 29.884008054261855), (4, 5, 28.006189408396626), (4, 7, 15.774839848636024), (4, 8, 3.6206480052249144), (4, 9, 23.804744370383144), (5, 0, 28.079672647730469), (5, 1, 17.042671653231526), (5, 2, 1.8557839901886808), (5, 3, 1.9922790489539697), (5, 4, 28.006189408396626), (5, 6, 21.492976178079076), (5, 8, 29.977996181215822), (5, 9, 4.4452006140146185), (6, 1, 16.410544777532913), (6, 2, 23.122260476726876), (6, 3, 22.479228556182363), (6, 5, 21.492976178079076), (6, 9, 20.04488615603487), (7, 0, 23.017867703525138), (7, 1, 25.598808236367063), (7, 4, 15.774839848636024), (7, 8, 16.915923579829141), (8, 0, 6.1268526078857208), (8, 1, 20.175759189503321), (8, 4, 3.6206480052249144), (8, 5, 29.977996181215822), (8, 7, 16.915923579829141), (8, 9, 25.962918470558208), (9, 0, 24.573433868296771), (9, 1, 12.843763853990932), (9, 2, 6.052562826627808), (9, 3, 6.4359869969688059), (9, 4, 23.804744370383144), (9, 5, 4.4452006140146185), (9, 6, 20.04488615603487), (9, 8, 25.962918470558208)]
  30
  31 #n=3
  32 #lg= [(0,1,23),(1,0,15),(1,2,45)]
  33 sink = n-1
  34
  35
  36 def genereGraph():
  37     test = False
  38     while not test :
  39         G = nx.DiGraph()
  40         l = [(random()*coteCarre, random()*coteCarre) for _  in range(n)]
  41         for io in range(len(l)) :
  42             for ie in range(len(l)) :
  43                 if ie != io  :
  44                     dist = mt.sqrt((l[io][0]-l[ie][0])**2 + (l[io][1]-l[ie][1])**2)
  45                     if dist <= distanceEmmissionMax :
  46                         G.add_edge(io,ie,weight=dist)
  47                         G.add_edge(ie,io,weight=dist)
  48         test = not(any([ not(nx.has_path(G,o,sink)) for o in  G.nodes() if sink in G.nodes() and o != sink]))
  49     return G
  50
  51 G = nx.DiGraph()
  52 G.add_weighted_edges_from(lg)
  53 #print nx.is_strongly_connected(G)
  54
  55
  56 #nx.draw(G)
  57 #pb.show()
  58
  59
  60
  61
  62
  63
  64 #print G.edges(data=True)
  65 #TODO afficher le graphe et etre sur qu'il est connexe
  66
  67 N = G.nodes()
  68
  69
  70 #V = list(set(sample(N,int(len(N)*vrate)))-set([sink]))
  71 V = list(set(N)-set([sink]))
  72
  73
  74 source = V
  75 print "source",source
  76
  77
  78 L = range(len(G.edges()))
  79 d = [di['weight'] for (_,_,di) in G.edges(data=True)]
  80
  81 #print "L",L
  82 #print "d",d
  83
  84
  85 def ail(i,l):
  86     assert  l in L, " pb de dimennsion de l: "+str(l)+ " "+ str(L)
  87     r = 0
  88     (o,e) = G.edges()[l]
  89     if o == i :
  90         r = 1
  91     elif e == i :
  92         r = -1
  93     return r
  94
  95
  96
  97 # Constantes
  98 a = [[ail(i,l) for l in L ] for i in xrange(n)]
  99 aplus  = [[1 if  ail(i,l) > 0 else 0  for l in L ] for i in xrange(n)]
 100 amoins  = [[1 if  ail(i,l) < 0 else 0  for l in L ] for i in xrange(n)]
 101
 102
 103
 104
 105
 106 alpha = 0.5
 107 beta = 1.3E-8
 108 gamma = 55.54
 109 delta = 0.2
 110 zeta = 0.1
 111 amplifieur = 1
 112 sigma2 = 3500
 113 Bi = 5
 114 omega = 0.15
 115 D = 100
 116 path_loss_exp = 4
 117 cr = 0.5
 118 cs = [alpha + beta*(di**path_loss_exp) for di in d]
 119 #print "cs",cs
 120
 121
 122 #TODO
 123 def etahi(h,i,R):
 124     ret = 0
 125     if i == h :
 126         ret = R[h]
 127     elif i == sink  :
 128         ret = -R[h]
 129     return ret
 130
 131
 132 def psi(Ps,i):
 133     if i not in Ps:
 134         return 0
 135     else :
 136         return Ps[i]
 137
 138
 139
 140 def cmpPair(x1,x2):
 141     return cmp(x1[1],x2[1])
 142
 143 def aminp(l):
 144     l.sort(cmp=cmpPair)
 145     return l[0][0]
 146
 147
 148
 149 def distance(d1,d2):
 150     return mt.sqrt(sum([(d1[t]-d2[t])**2 for t in d1]))
 151
 152
 153 def AfficheVariation (up,vp,lap,wp,thetap,etap,qp,Psp,Rhp,xp,valeurFonctionDualep):
 154     global u, v, la, w, theta , q, Ps, Rh, eta, x,valeurFonctionDuale
 155
 156     print "du=",distance(u,up),
 157     print "dv=",distance(v,vp),
 158     print "dlambda=",distance(la,lap),
 159     print "dw=",distance(w,wp),
 160     print "dtheta=",abs(theta-thetap),
 161     print "deta=",distance(eta,etap),
 162     print "dq=",distance(q,qp),
 163     print "dPs=",distance(Ps,Psp),
 164     print "dRh=",distance(Rh,Rhp),
 165     print "dx=",distance(x,xp),
 166     print "dL=",abs(valeurFonctionDuale-valeurFonctionDualep),"\n"
 167
 168
 169 valeurFonctionDuale = 0
 170
 171 def entre0et1(x):
 172     r = x if (x >0 and x <= 1) else -1
 173     #print "ds entre0et1 x et r", x,r
 174     return r
 175
 176
 177 def xpos(x):
 178     r = x if x >0 else -1
 179     #print "ds xpos x et r", x,r
 180     return r
 181
 182 def xposounul(x):
 183     return x if x >= 0 else -1
 184
 185 #f_q,q[i],POS,[i]
 186 def armin(f,xini,xr,args):
 187     #xpos = lambda x : x if x > 0 else -1
 188     r = 0
 189     if xr == POS :
 190         #print "strictement pos"
 191         #print "parametres passes a cobyla",xini,xpos,args,"consargs=(),rhoend=1E-5,iprint=0,maxfun=1000"
 192         r= opt.fmin_cobyla(f,xini,cons=[xpos],args=args,consargs=(),rhoend=1E-3,iprint=0,maxfun=nbiter)
 193         #print "le min str pos est",r
 194         #print "le min pos   est", r,"avec xini",xini
 195     elif xr == POS_NUL :
 196         #print "pos ou nul"
 197         r = opt.fmin_cobyla(f,xini,[xposounul],args,consargs=(),rhoend=1E-3,iprint=0,maxfun=nbiter)
 198         # print "le min pos est", r
 199         #print "le min pos null  est", r,"avec xini",xini
 200     elif xr == POSINF1:
 201         r = opt.fmin_cobyla(f,xini,[entre0et1],args,consargs=(),rhoend=1E-3,iprint=0,maxfun=nbiter)
 202         #print "le min pos inf 1 est", r,"avec xini",xini
 203
 204
 205     return r
 206
 207
 208
 209
 210
 211
 212 def maj_theta(k):
 213     return omega/(mt.pow(k,0.5))
 214
 215
 216
 217 def maj_theta(k):
 218     return omega/mt.sqrt(k)
 219
 220
 221
 222
 223 def maj(k,maj_theta,mxg,idxexp):
 224     # initialisation des operateurs lagrangiens
 225     global u, v, la, w, theta , q,  Ps, Rh,  eta, x, valeurFonctionDuale
 226
 227     smax = 0
 228     #print "iteration",k
 229
 230     up = {}
 231     for h in V:
 232         for i in N:
 233             if not ASYNC or  random() < taux_succes:
 234                 s = eta[(h,i)]-sum([a[i][l]*x[(h,l)] for l in L])
 235                 if abs(s) > mxg :
 236                     print "ds calcul u",abs(s),idxexp
 237                     mxg = abs(s)
 238                 smax = max(smax,abs(s))
 239                 up[(h,i)] = u[(h,i)]-theta*s
 240             else :
 241                 up[(h,i)] = u[(h,i)]
 242
 243
 244
 245     vp = {}
 246     for h in V:
 247         if not ASYNC or  random() < taux_succes:
 248             s = Rh[h]- mt.log(float(sigma2)/D)/(gamma*mt.pow(Ps[h],float(2)/3))
 249             if abs(s) > mxg :
 250                 print "ds calcul v",abs(s),idxexp
 251                 mxg = abs(s)
 252             smax = max(smax,abs(s))
 253             vp[h] = max(0,v[h]-theta*s)
 254         else :
 255             vp[h] = v[h]
 256
 257
 258
 259     lap={}
 260     for i in N:
 261         if not ASYNC or  random() < taux_succes:
 262             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)
 263             if abs(s) > mxg :
 264                 print "ds calcul la",abs(s),idxexp,i
 265                 mxg = abs(s)
 266             smax = max(smax,abs(s))
 267             resla = la[i]-theta*s
 268             lap[i] = max(0,resla)
 269         else :
 270             lap[i] = la[i]
 271
 272
 273
 274
 275     wp={}
 276     for l in L:
 277         if not ASYNC or  random() < taux_succes:
 278             s = sum([a[i][l]*q[i] for i in N])
 279             if abs(s) > mxg :
 280                 print "ds calcul w",abs(s),idxexp
 281                 mxg = abs(s)
 282             smax = max(smax,abs(s))
 283             wp[l] = w[l] + theta*s
 284         else :
 285             wp[l] = w[l]
 286
 287
 288
 289     thetap = maj_theta(k)
 290
 291
 292     qp={}
 293     def f_q(qi,i):
 294         fa = sum([a[i][l]*w[l] for l in L]) - la[i]*Bi
 295         return qi*qi + qi*fa
 296
 297     for i in N:
 298         if not ASYNC or  random() < taux_succes:
 299             c = -float(sum([a[i][l]*w[l] for l in L]) - la[i]*Bi)/(2*amplifieur)
 300             rep = epsilon if c <= 0 else c
 301             qp[i] = rep
 302         else :
 303             qp[i] = q[i]
 304
 305
 306
 307
 308
 309     Psp={}
 310     #print "maj des des Psh"
 311     def f_Ps(psh,h):
 312         #print "ds f_ps",psh, v[h]* mt.log(float(sigma2)/D)/(gamma*((psh**2)**(float(2)/3))) +la[h]*psh
 313          return v[h]* mt.log(float(sigma2)/D)/(gamma*mt.pow(float(2)/3)) +la[h]*psh
 314      for h in V:
 315          if not ASYNC or  random() < taux_succes:
 316              """
 317              lah = 0.05 if la[h] == 0 else  la[h]
 318              rep = mt.pow(float(2*v[h]*mt.log(float(sigma2)/D))/(3*gamma*lah),float(3)/5)
 319              Psp[h] = epsilon if rep <= 0 else rep
 320              """
 321              t= float(-3*la[h]+mt.sqrt(9*(la[h]**2)+64*zeta*v[h]*mt.log(float(sigma2)/D)/gamma))/(16*zeta)
 322              #print t
 323              rep = mt.pow(t,float(3)/5)
 324              Psp[h]=rep
 325          else :
 326             Psp[h] = Ps[h]
 327
 328
 329
 330     etap={}
 331
 332     for h in V:
 333         for i in N :
 334             etap[(h,i)] = etahi(h,i,Rh)
 335
 336
 337
 338     Rhp={}
 339     def f_Rh(rh,h):
 340         return delta*rh*rh-v[h]*rh-sum([u[(h,i)]*eta[(h,i)] for i in N])
 341
 342     for h in V:
 343         if not ASYNC or  random() < taux_succes:
 344             rep = float(v[h])/(2*delta)
 345             Rhp[h] = 0 if rep < 0 else rep
 346         else :
 347             Rhp[h] = Rh[h]
 348
 349     xp={}
 350     def f_x(xhl,h,l):
 351         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]))
 352         return r
 353
 354
 355     for h in V:
 356         for l in L:
 357             if not ASYNC or  random() < taux_succes:
 358                 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)
 359                 xp[(h,l)] = 0 if rep < 0 else rep
 360             else :
 361                 xp[(h,l)] = x[(h,l)]
 362
 363
 364
 365
 366
 367     valeurFonctionDualep = 0
 368     valeurFonctionDualep += sum([amplifieur*q[i]*q[i] for i in N])
 369     valeurFonctionDualep += sum([sum([delta*(x[(h,l)]**2) for l in L]) for h in V])
 370     valeurFonctionDualep += sum([delta*(Rh[h]**2) for h in V])
 371     valeurFonctionDualep += sum([sum([u[(h,i)]*(sum([ a[i][l]*x[(h,l)] for l in L])- eta[(h,i)]) for i in N]) for h in V])
 372     valeurFonctionDualep += sum([v[h]*(mt.log(float(sigma2)/D)/(gamma*mt.pow(Ps[h],float(2)/3)) - Rh[h]) for h in V])
 373     valeurFonctionDualep += sum([la[i]*(psi(Ps,i) +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]) -q[i]*Bi)  for i in N ])
 374     valeurFonctionDualep += sum([w[l]*sum([a[i][l]*q[i] for i in N]) for l in L])
 375
 376
 377     #AfficheVariation(up,vp,lap,wp,thetap,etap,qp,Psp,Rhp,xp,valeurFonctionDualep)
 378
 379     arret = abs(valeurFonctionDuale-valeurFonctionDualep) < error
 380
 381     return (up,vp,lap,wp,thetap,etap,qp,Psp,Rhp,xp,valeurFonctionDualep,arret,mxg,smax)
 382
 383
 384
 385
 386 """
 387 print "u",u
 388 print "v",v
 389 print "lambda",la
 390 print "w",w
 391 print "theta",theta
 392 print "eta", eta
 393 print "q",q
 394 print "Ps",Ps
 395 print "Rh",Rh
 396 print "x",x
 397
 398 """
 399
 400
 401
 402 def initialisation():
 403     global u, v, la, w, theta , q,  Ps, Rh,  eta, x,init
 404
 405     theta = omega
 406
 407     q = {}
 408     for i in N :
 409         q[i] = 0.15 + random()*0.05
 410
 411
 412     Ps= {}
 413     for h in V :
 414         Ps[h] =  0.2+random()*0.3
 415
 416     Rh = {}
 417     for vi in V:
 418         Rh[vi] = 0.1 + random()*0.1
 419
 420     eta = {}
 421     for h in V :
 422         for i in N:
 423             eta[(h,i)]= etahi(h,i,Rh)
 424
 425     x={}
 426     for h in V :
 427         for l in L:
 428             x[(h,l)]=0
 429
 430
 431
 432
 433     # initialisation des operateurs lagrangiens
 434     u={}
 435     for h in V :
 436         for i in N:
 437             u[(h,i)]= random()
 438
 439     v = {}
 440     for h in V :
 441         v[h] = Rh[h]
 442
 443     la = {}
 444     for i in N:
 445         la[i] = random()
 446
 447     w={}
 448     for l in L:
 449         w[l] =  random()
 450
 451     init = [deepcopy(q),deepcopy(Ps),deepcopy(Rh),deepcopy(eta),
 452             deepcopy(x),deepcopy(u),deepcopy(v),deepcopy(la),deepcopy(w)]
 453
 454
 455
 456 def initialisation_():
 457     global u, v, la, w, theta , q,  Ps, Rh,  eta, x,init
 458     fd = open(fichier_init,"r")
 459     l= fd.readline()
 460     init_p = eval(l)
 461     print init_p
 462     theta = omega
 463     (q,Ps,Rh,eta,x,u,v,la,w) = tuple([deepcopy(x) for x in init_p])
 464     init = [deepcopy(q),deepcopy(Ps),deepcopy(Rh),deepcopy(eta),
 465             deepcopy(x),deepcopy(u),deepcopy(v),deepcopy(la),deepcopy(w)]
 466
 467
 468
 469 def __evalue_maj_theta__(nbexp,out=False):
 470     global u, v, la, w, theta , q,  Ps, Rh,  eta, x, valeurFonctionDuale
 471     nbexp = 10
 472     res = {}
 473     m = []
 474     itermax = 100000
 475
 476     def __maj_theta(k):
 477         mem = []
 478         om = omega/(mt.pow(k,0.75))
 479         return om
 480     for idxexp in range(nbexp):
 481         mxg = 0
 482         if not(out):
 483             initialisation()
 484         else :
 485             initialisation_()
 486
 487         k = 1
 488         arret = False
 489         sm = 0
 490         while k < itermax  and not arret :
 491             (u,v,la,w,theta,eta,q,Ps,Rh,x,valeurFonctionDuale,ar,mxg,smax)=maj(k,__maj_theta,mxg,idxexp)
 492             errorq =  (max(q.values()) - min(q.values()))/min(q.values())
 493             arret = errorq <  error
 494             k+=1
 495             variation = "+" if smax > sm else "-"
 496             print variation,
 497             if k%100 ==0 :
 498                 print "k:",k,"erreur sur q", errorq, "et q:",q
 499                 print "maxg=", mxg
 500                 mem = [deepcopy(q),deepcopy(Ps),deepcopy(Rh),deepcopy(eta),
 501                        deepcopy(x),deepcopy(u),deepcopy(v),deepcopy(la),deepcopy(w)]
 502             if k%4500 == 0 :
 503                 print "#########\n",mem,"\#########\n"
 504             if k%4600 == 0 :
 505                 print "#########\n",mem,"\#########\n"
 506
 507
 508
 509             if smax - sm  > 500:
 510                 print "variation trop grande"
 511                 print "init"
 512                 print init
 513                 exit
 514             sm = smax
 515
 516         if k == itermax:
 517             print "nbre d'iteration trop grand"
 518             print "init"
 519             print init
 520             sy.exit(1)
 521
 522         print "###############"
 523         print k
 524         print "###############"
 525         m += [k]
 526
 527     print (min(m),max(m),float(sum(m))/nbexp,m),m
 528
 529
 530
 531 def __une_seule_exp__(fichier_donees):
 532     global u, v, la, w, theta , q,  Ps, Rh,  eta, x, valeurFonctionDuale
 533     initialisation()
 534
 535
 536     fichier = open(fichier_donees, "r")
 537     instructions ={}
 538     for line in fichier:
 539          l = line.split("=")
 540          instructions[l[0]] = eval(l[1])
 541     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']
 542     nbexp = 1
 543     res = {}
 544     m = []
 545     itermax = 100000
 546
 547     def __maj_theta(k):
 548         om = omega/(mt.pow(k,0.75))
 549         return om
 550     for idxexp in range(nbexp):
 551         mxg = 0
 552         k = 1
 553         arret = False
 554         sm = 0
 555         while k < itermax  and not arret :
 556             (u,v,la,w,theta,eta,q,Ps,Rh,x,valeurFonctionDuale,ar,mxg,smax)=maj(k,__maj_theta,mxg,idxexp)
 557             errorq =  (max(q.values()) - min(q.values()))/min(q.values())
 558             arret = errorq <  error
 559             k+=1
 560             variation = "+" if smax > sm else "-"
 561             print variation,
 562             if k%100 ==0 :
 563                 print "k:",k,"erreur sur q", errorq, "et q:",q
 564                 print "maxg=", mxg
 565             if smax - sm  > 500:
 566                 print "variation trop grande"
 567                 print "init"
 568                 print init
 569                 exit
 570             sm = smax
 571
 572         if k == itermax:
 573             print "nbre d'iteration trop grand"
 574             print "init"
 575             print init
 576             exit
 577
 578         print "###############"
 579         print k
 580         print "###############"
 581         m += [k]
 582
 583     print (min(m),max(m),float(sum(m))/nbexp,m),m
 584
 585
 586
 587
 588 ASYNC = False
 589 __une_seule_exp__("config_initiale.py")
 590 #__evalue_maj_theta__()
 591 #ASYNC = True
 592 #taux_succes = 0.9
 593 #__evalue_maj_theta__()
 594
 595
 596
 597
 598
 599 print "u",u
 600 print "v",v
 601 print "lambda",la
 602 print "w",w
 603 print "theta",theta
 604 print "eta", eta
 605 print "q",q
 606 print "Ps",Ps
 607 print "Rh",Rh
 608 print "x",x
 609 print "L",valeurFonctionDuale
 610
 611
 612
 613 # relation ente les variables primaires et secondaires ?
 614