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     res = {}
 472     m = []
 473     itermax = 100000
 474
 475     def __maj_theta(k):
 476         mem = []
 477         om = omega/(mt.pow(k,0.75))
 478         return om
 479     for idxexp in range(nbexp):
 480         mxg = 0
 481         if not(out):
 482             initialisation()
 483         else :
 484             initialisation_()
 485
 486         k = 1
 487         arret = False
 488         sm = 0
 489         while k < itermax  and not arret :
 490             (u,v,la,w,theta,eta,q,Ps,Rh,x,valeurFonctionDuale,ar,mxg,smax)=maj(k,__maj_theta,mxg,idxexp)
 491             errorq =  (max(q.values()) - min(q.values()))/min(q.values())
 492             arret = errorq <  error
 493             k+=1
 494             variation = "+" if smax > sm else "-"
 495             sm = smax
 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                 sy.exit(0)
 514         if k == itermax:
 515             print "nbre d'iteration trop grand"
 516             print "init"
 517             print init
 518             sy.exit(1)
 519
 520         print "###############"
 521         print k
 522         print "###############"
 523         m += [k]
 524
 525     print (min(m),max(m),float(sum(m))/nbexp,m),m
 526
 527
 528
 529
 530
 531 ASYNC = False
 532 __evalue_maj_theta__(1,True)
 533
 534
 535
 536 #ASYNC = True
 537 #taux_succes = 0.9
 538 #__evalue_maj_theta__()
 539
 540
 541 """
 542 initialisation()
 543 k = 1
 544 arret = False
 545 while k < 10000  and not arret :
 546     (u,v,la,w,theta,eta,q,Ps,Rh,x,valeurFonctionDuale,ar)=maj(k,maj_theta)
 547     arret = ar
 548     k+=1
 549     errorq =  abs(min(q.values())-max(q.values()))
 550     print "errorq",errorq
 551     arret = errorq <  error
 552 """
 553
 554 """
 555     print "u",u
 556     print "w",w
 557     print "theta",theta
 558     print "eta", eta
 559     print "q",q
 560     print "v",v
 561     print "lambda",la
 562     print "Ps",Ps
 563     print "Rh",Rh
 564     print "x",x
 565 """
 566
 567
 568 """
 569     k +=1
 570
 571 """
 572
 573
 574
 575 print "u",u
 576 print "v",v
 577 print "lambda",la
 578 print "w",w
 579 print "theta",theta
 580 print "eta", eta
 581 print "q",q
 582 print "Ps",Ps
 583 print "Rh",Rh
 584 print "x",x
 585 print "L",valeurFonctionDuale
 586
 587
 588
 589 # relation ente les variables primaires et secondaires ?
 590