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(1)/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(1)/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))/float(16*zeta)
 322             #print t
 323             rep = mt.pow(t,float(3)/5)
 324             Psp[h]=rep
 325             """
 326 ############reprendre ici
 327
 328         else :
 329             Psp[h] = Ps[h]
 330
 331
 332
 333     etap={}
 334
 335     for h in V:
 336         for i in N :
 337             etap[(h,i)] = etahi(h,i,Rh)
 338
 339
 340
 341     Rhp={}
 342     def f_Rh(rh,h):
 343         return delta*rh*rh-v[h]*rh-sum([u[(h,i)]*eta[(h,i)] for i in N])
 344
 345     for h in V:
 346         if not ASYNC or  random() < taux_succes:
 347             rep = float(v[h])/(2*delta)
 348             Rhp[h] = 0 if rep < 0 else rep
 349         else :
 350             Rhp[h] = Rh[h]
 351
 352     xp={}
 353     def f_x(xhl,h,l):
 354         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]))
 355         return r
 356
 357
 358     for h in V:
 359         for l in L:
 360             if not ASYNC or  random() < taux_succes:
 361                 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)
 362                 xp[(h,l)] = 0 if rep < 0 else rep
 363             else :
 364                 xp[(h,l)] = x[(h,l)]
 365
 366
 367
 368
 369
 370     valeurFonctionDualep = 0
 371     valeurFonctionDualep += sum([amplifieur*q[i]*q[i] for i in N])
 372     valeurFonctionDualep += sum([sum([delta*(x[(h,l)]**2) for l in L]) for h in V])
 373     valeurFonctionDualep += sum([delta*(Rh[h]**2) for h in V])
 374     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])
 375     valeurFonctionDualep += sum([v[h]*(mt.log(float(sigma2)/D)/(gamma*mt.pow(Ps[h],float(2)/3)) - Rh[h]) for h in V])
 376     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 ])
 377     valeurFonctionDualep += sum([w[l]*sum([a[i][l]*q[i] for i in N]) for l in L])
 378
 379
 380     #AfficheVariation(up,vp,lap,wp,thetap,etap,qp,Psp,Rhp,xp,valeurFonctionDualep)
 381
 382     arret = abs(valeurFonctionDuale-valeurFonctionDualep) < error
 383
 384     return (up,vp,lap,wp,thetap,etap,qp,Psp,Rhp,xp,valeurFonctionDualep,arret,mxg,smax)
 385
 386
 387
 388
 389 """
 390 print "u",u
 391 print "v",v
 392 print "lambda",la
 393 print "w",w
 394 print "theta",theta
 395 print "eta", eta
 396 print "q",q
 397 print "Ps",Ps
 398 print "Rh",Rh
 399 print "x",x
 400
 401 """
 402
 403
 404
 405 def initialisation():
 406     global u, v, la, w, theta , q,  Ps, Rh,  eta, x,init
 407
 408     theta = omega
 409
 410     q = {}
 411     for i in N :
 412         q[i] = 0.15 + random()*0.05
 413
 414
 415     Ps= {}
 416     for h in V :
 417         Ps[h] =  0.2+random()*0.3
 418
 419     Rh = {}
 420     for vi in V:
 421         Rh[vi] = 0.1 + random()*0.1
 422
 423     eta = {}
 424     for h in V :
 425         for i in N:
 426             eta[(h,i)]= etahi(h,i,Rh)
 427
 428     x={}
 429     for h in V :
 430         for l in L:
 431             x[(h,l)]=0
 432
 433
 434
 435
 436     # initialisation des operateurs lagrangiens
 437     u={}
 438     for h in V :
 439         for i in N:
 440             u[(h,i)]= random()
 441
 442     v = {}
 443     for h in V :
 444         v[h] = Rh[h]
 445
 446     la = {}
 447     for i in N:
 448         la[i] = random()
 449
 450     w={}
 451     for l in L:
 452         w[l] =  random()
 453
 454     init = [deepcopy(q),deepcopy(Ps),deepcopy(Rh),deepcopy(eta),
 455             deepcopy(x),deepcopy(u),deepcopy(v),deepcopy(la),deepcopy(w)]
 456
 457
 458
 459 def initialisation_():
 460     global u, v, la, w, theta , q,  Ps, Rh,  eta, x,init
 461     fd = open(fichier_init,"r")
 462     l= fd.readline()
 463     init_p = eval(l)
 464     print init_p
 465     theta = omega
 466     (q,Ps,Rh,eta,x,u,v,la,w) = tuple([deepcopy(x) for x in init_p])
 467     init = [deepcopy(q),deepcopy(Ps),deepcopy(Rh),deepcopy(eta),
 468             deepcopy(x),deepcopy(u),deepcopy(v),deepcopy(la),deepcopy(w)]
 469
 470
 471
 472 def __evalue_maj_theta__(nbexp,out=False):
 473     global u, v, la, w, theta , q,  Ps, Rh,  eta, x, valeurFonctionDuale
 474     res = {}
 475     m = []
 476     itermax = 100000
 477
 478     def __maj_theta(k):
 479         mem = []
 480         om = omega/(mt.pow(k,0.75))
 481         return om
 482     for idxexp in range(nbexp):
 483         mxg = 0
 484         if not(out):
 485             initialisation()
 486         else :
 487             initialisation_()
 488
 489         k = 1
 490         arret = False
 491         sm = 0
 492         while k < itermax  and not arret :
 493             (u,v,la,w,theta,eta,q,Ps,Rh,x,valeurFonctionDuale,ar,mxg,smax)=maj(k,__maj_theta,mxg,idxexp)
 494             errorq =  (max(q.values()) - min(q.values()))/min(q.values())
 495             arret = errorq <  error
 496             k+=1
 497             variation = "+" if smax > sm else "-"
 498             sm = smax
 499             print variation,
 500             if k%100 ==0 :
 501                 print "k:",k,"erreur sur q", errorq, "et q:",q
 502                 print "maxg=", mxg
 503                 mem = [deepcopy(q),deepcopy(Ps),deepcopy(Rh),deepcopy(eta),
 504                        deepcopy(x),deepcopy(u),deepcopy(v),deepcopy(la),deepcopy(w)]
 505             if k%4500 == 0 :
 506                 print "#########\n",mem,"\#########\n"
 507             if k%4600 == 0 :
 508                 print "#########\n",mem,"\#########\n"
 509
 510
 511
 512             if smax - sm  > 500:
 513                 print "variation trop grande"
 514                 print "init"
 515                 print init
 516                 sy.exit(0)
 517         if k == itermax:
 518             print "nbre d'iteration trop grand"
 519             print "init"
 520             print init
 521             sy.exit(1)
 522
 523         print "###############"
 524         print k
 525         print "###############"
 526         m += [k]
 527
 528     print (min(m),max(m),float(sum(m))/nbexp,m),m
 529
 530
 531
 532
 533
 534 ASYNC = False
 535 __evalue_maj_theta__(1,True)
 536
 537
 538
 539 #ASYNC = True
 540 #taux_succes = 0.9
 541 #__evalue_maj_theta__()
 542
 543
 544 """
 545 initialisation()
 546 k = 1
 547 arret = False
 548 while k < 10000  and not arret :
 549     (u,v,la,w,theta,eta,q,Ps,Rh,x,valeurFonctionDuale,ar)=maj(k,maj_theta)
 550     arret = ar
 551     k+=1
 552     errorq =  abs(min(q.values())-max(q.values()))
 553     print "errorq",errorq
 554     arret = errorq <  error
 555 """
 556
 557 """
 558     print "u",u
 559     print "w",w
 560     print "theta",theta
 561     print "eta", eta
 562     print "q",q
 563     print "v",v
 564     print "lambda",la
 565     print "Ps",Ps
 566     print "Rh",Rh
 567     print "x",x
 568 """
 569
 570
 571 """
 572     k +=1
 573
 574 """
 575
 576
 577
 578 print "u",u
 579 print "v",v
 580 print "lambda",la
 581 print "w",w
 582 print "theta",theta
 583 print "eta", eta
 584 print "q",q
 585 print "Ps",Ps
 586 print "Rh",Rh
 587 print "x",x
 588 print "L",valeurFonctionDuale
 589
 590
 591
 592 # relation ente les variables primaires et secondaires ?
 593