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 import cv as cv
  11 import cv2 as cv2
  12 import datetime as dt
  13
  14 errorq  = 0.1
  15 errorD  = 1E-1
  16 epsilon = 1E-10
  17 vrate = 0.8
  18 p = 0.7
  19 coteCarre = 50
  20 distanceEmmissionMax = 30
  21 nbiter = 1000
  22 POS = 1
  23 POS_NUL = 2
  24 POSINF1 = 3
  25 init = []
  26 fichier_init="config_initiale_default.txt"
  27
  28
  29
  30 # construction du  graphe
  31 n = 10
  32 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)]
  33
  34 #n=3
  35 #lg= [(0,1,23),(1,0,15),(1,2,45)]
  36 sink = n-1
  37
  38 def distance(d1,d2):
  39     return mt.sqrt(sum([(d1[t]-d2[t])**2 for t in d1]))
  40
  41
  42 def genereGraph():
  43     test = False
  44     while not test :
  45         G = nx.DiGraph()
  46         l = [(random()*coteCarre, random()*coteCarre) for _  in range(n)]
  47         for io in range(len(l)) :
  48             for ie in range(len(l)) :
  49                 if ie != io  :
  50                     dist = mt.sqrt((l[io][0]-l[ie][0])**2 + (l[io][1]-l[ie][1])**2)
  51                     if dist <= distanceEmmissionMax :
  52                         G.add_edge(io,ie,weight=dist)
  53                         G.add_edge(ie,io,weight=dist)
  54         test = not(any([ not(nx.has_path(G,o,sink)) for o in  G.nodes() if sink in G.nodes() and o != sink]))
  55     return (G,l)
  56
  57
  58 def afficheGraph(G,l,tx,ty,sink):
  59     r = 20
  60     img = cv.CreateImage ((tx, ty), 32, 3)
  61     cv.Rectangle(img, (0,0),(tx,ty), cv.Scalar(255,255,255), thickness=-1)
  62     def px((x,y)):
  63         u = float(tx)/(coteCarre + 2*distanceEmmissionMax)
  64         return(int(distanceEmmissionMax*u + x * u),int(distanceEmmissionMax*u + y * u))
  65     for i in set(range(len(l)))-set([sink]):
  66         (x,y) = l[i]
  67         pix,piy = px((x,y))
  68         demx = distanceEmmissionMax*tx/(coteCarre+2*distanceEmmissionMax)
  69         cv.Circle(img, (pix,piy),demx, cv.Scalar(125,125,125))
  70
  71     for i in set(range(len(l)))-set([sink]):
  72         (x,y) = l[i]
  73         pix,piy = px((x,y))
  74         cv.Circle(img, (pix,piy),r, cv.Scalar(125,125,125),thickness=-1)
  75
  76     #sink
  77     (x,y) = l[sink]
  78     pix,piy = px((x,y))
  79
  80     cv.Rectangle(img, (pix-r/2,piy-r/2),(pix+r/2,piy+r/2), cv.Scalar(125,125,125), thickness=-1)
  81
  82     for i in range(len(l)):
  83         for j in range(len(l)):
  84
  85             if np.linalg.norm(np.array(l[i])-np.array(l[j])) < distanceEmmissionMax :
  86                 (xi,yi) = l[i]
  87                 pixi,piyi = px((xi,yi))
  88                 (xj,yj) = l[j]
  89                 pixj,piyj = px((xj,yj))
  90                 cv.Line(img, (pixi,piyi), (pixj,piyj), cv.Scalar(125,125,125))
  91
  92
  93     """
  94     for i in range(len(l)):
  95         (x,y) = l[i]
  96         pix,piy = px((x,y))
  97         print i
  98         textColor = (0, 0, 255)  # red
  99         font = cv2.FONT_HERSHEY_SIMPLEX
 100         imgp =
 101         cv2.putText(img, str(i), (pix-r/4,piy-r/2),font, 3.0, textColor)#,thickn    """
 102     cv.SaveImage("SensorNetwork.png",img)
 103
 104 G = nx.DiGraph()
 105 G.add_weighted_edges_from(lg)
 106 #print nx.is_strongly_connected(G)
 107
 108
 109 #nx.draw(G)
 110 #pb.show()
 111
 112 #(G,l) = genereGraph()
 113 N = G.nodes()
 114 #V = list(set(sample(N,int(len(N)*vrate)))-set([sink]))
 115 V = list(set(N)-set([sink]))
 116 source = V
 117 print "source",source
 118
 119
 120 #afficheGraph(G,l,500,500,sink)
 121 #nx.draw(G)
 122 #pb.show()
 123
 124
 125
 126
 127
 128
 129
 130 #print G.edges(data=True)
 131 #TODO afficher le graphe et etre sur qu'il est connexe
 132
 133
 134
 135
 136
 137 L = range(len(G.edges()))
 138 d = [di['weight'] for (_,_,di) in G.edges(data=True)]
 139
 140 #print "L",L
 141 #print "d",d
 142
 143
 144 def ail(i,l):
 145     assert  l in L, " pb de dimennsion de l: "+str(l)+ " "+ str(L)
 146     r = 0
 147     (o,e) = G.edges()[l]
 148     if o == i :
 149         r = 1
 150     elif e == i :
 151         r = -1
 152     return r
 153
 154
 155
 156 # Constantes
 157 a = [[ail(i,l) for l in L ] for i in xrange(n)]
 158 aplus  = [[1 if  ail(i,l) > 0 else 0  for l in L ] for i in xrange(n)]
 159 amoins  = [[1 if  ail(i,l) < 0 else 0  for l in L ] for i in xrange(n)]
 160
 161
 162
 163
 164
 165 alpha = 0.5
 166 beta = 1.3E-8
 167 gamma = 55.54
 168 delta = 0.2
 169 amplifieur = 1
 170 sigma2 = 3500
 171 Bi = 5
 172 omega = 0.15
 173 D = 100
 174 path_loss_exp = 4
 175 cr = 0.5
 176 cs = [alpha + beta*(di**path_loss_exp) for di in d]
 177 #print "cs",cs
 178
 179
 180 #TODO
 181 def etahi(h,i,R):
 182     ret = 0
 183     if i == h :
 184         ret = R[h]
 185     elif i == sink  :
 186         ret = -R[h]
 187     return ret
 188
 189
 190 def psi(Ps,i):
 191     if i not in Ps:
 192         return 0
 193     else :
 194         return Ps[i]
 195
 196
 197
 198 def cmpPair(x1,x2):
 199     return cmp(x1[1],x2[1])
 200
 201 def aminp(l):
 202     l.sort(cmp=cmpPair)
 203     return l[0][0]
 204
 205
 206
 207
 208
 209 def AfficheVariation (up,vp,lap,wp,thetap,etap,qp,Psp,Rhp,xp,valeurFonctionDualep):
 210     global u, v, la, w, theta , q, Ps, Rh, eta, x,valeurFonctionDuale
 211
 212     print "du=",distance(u,up),
 213     print "dv=",distance(v,vp),
 214     print "dlambda=",distance(la,lap),
 215     print "dw=",distance(w,wp),
 216     print "dtheta=",abs(theta-thetap),
 217     print "deta=",distance(eta,etap),
 218     print "dq=",distance(q,qp),
 219     print "dPs=",distance(Ps,Psp),
 220     print "dRh=",distance(Rh,Rhp),
 221     print "dx=",distance(x,xp),
 222     print "dL=",abs(valeurFonctionDuale-valeurFonctionDualep),"\n"
 223
 224
 225 valeurFonctionDuale = 0
 226
 227 def entre0et1(x):
 228     r = x if (x >0 and x <= 1) else -1
 229     #print "ds entre0et1 x et r", x,r
 230     return r
 231
 232
 233 def xpos(x):
 234     r = x if x >0 else -1
 235     #print "ds xpos x et r", x,r
 236     return r
 237
 238 def xposounul(x):
 239     return x if x >= 0 else -1
 240
 241 #f_q,q[i],POS,[i]
 242 def armin(f,xini,xr,args):
 243     #xpos = lambda x : x if x > 0 else -1
 244     r = 0
 245     if xr == POS :
 246         r= opt.fmin_cobyla(f,xini,cons=[xpos],args=args,consargs=(),rhoend=errorD/100,iprint=0,maxfun=nbiter)
 247     elif xr == POS_NUL :
 248         r = opt.fmin_cobyla(f,xini,[xposounul],args,consargs=(),rhoend=errorD/100,iprint=0,maxfun=nbiter)
 249     elif xr == POSINF1:
 250         r = opt.fmin_cobyla(f,xini,[entre0et1],args,consargs=(),rhoend=errorD/100,iprint=0,maxfun=nbiter)
 251
 252
 253     return r
 254
 255
 256
 257
 258
 259
 260 def maj_theta(k):
 261     return omega/(mt.pow(k,0.5))
 262
 263
 264
 265 def maj_theta(k):
 266     return omega/mt.sqrt(k)
 267
 268
 269
 270
 271 def maj(k,maj_theta,mxg,idxexp):
 272     # initialisation des operateurs lagrangiens
 273     global u, v, la, w, theta , q,  Ps, Rh,  eta, x, valeurFonctionDuale
 274
 275     smax = 0
 276     #print "iteration",k
 277
 278     up = {}
 279     for h in V:
 280         for i in N:
 281             if not ASYNC or  random() < taux_succes:
 282                 s = eta[(h,i)]-sum([a[i][l]*x[(h,l)] for l in L])
 283                 if abs(s) > mxg :
 284                     print "ds calcul u",abs(s),idxexp
 285                     mxg = abs(s)
 286                 smax = max(smax,abs(s))
 287                 up[(h,i)] = u[(h,i)]-theta*s
 288             else :
 289                 up[(h,i)] = u[(h,i)]
 290
 291
 292
 293     vp = {}
 294     for h in V:
 295         if not ASYNC or  random() < taux_succes:
 296             s = Rh[h]- mt.log(float(sigma2)/D)/(gamma*mt.pow(Ps[h],float(2)/3))
 297             if abs(s) > mxg :
 298                 print "ds calcul v",abs(s),idxexp
 299                 mxg = abs(s)
 300             smax = max(smax,abs(s))
 301             vp[h] = max(0,v[h]-theta*s)
 302         else :
 303             vp[h] = v[h]
 304
 305
 306
 307     lap={}
 308     for i in N:
 309         if not ASYNC or  random() < taux_succes:
 310             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)
 311             if abs(s) > mxg :
 312                 print "ds calcul la",abs(s),idxexp,i
 313                 mxg = abs(s)
 314             smax = max(smax,abs(s))
 315             resla = la[i]-theta*s
 316             lap[i] = max(0,resla)
 317         else :
 318             lap[i] = la[i]
 319
 320
 321
 322
 323     wp={}
 324     for l in L:
 325         if not ASYNC or  random() < taux_succes:
 326             s = sum([a[i][l]*q[i] for i in N])
 327             if abs(s) > mxg :
 328                 print "ds calcul w",abs(s),idxexp
 329                 mxg = abs(s)
 330             smax = max(smax,abs(s))
 331             wp[l] = w[l] + theta*s
 332         else :
 333             wp[l] = w[l]
 334
 335
 336
 337     thetap = maj_theta(k)
 338
 339
 340     qp={}
 341     def f_q(qi,*args):
 342         fa = sum([a[i][l]*w[l] for l in L]) - la[i]*Bi
 343         return qi*qi + qi*fa
 344
 345     tq,tqa = 0,0
 346     for i in N:
 347         if not ASYNC or  random() < taux_succes:
 348             n1 = dt.datetime.now()
 349             c = -float(sum([a[i][l]*w[l] for l in L]) - la[i]*Bi)/(2*amplifieur)
 350             n2 = dt.datetime.now()
 351             cp = armin(f_q,q[i],POS,tuple([i]))
 352             n3 = dt.datetime.now()
 353             tq += (n2-n1).microseconds
 354             tqa += (n3-n2).microseconds
 355             rep = epsilon if c <= 0 else c
 356             qp[i] = rep
 357         else :
 358             qp[i] = q[i]
 359
 360
 361
 362     tp,tpa = 0,0
 363     Psp={}
 364     def f_Ps(psh,h):
 365         try:
 366             #print psh,(gamma*mt.pow(psh,float(2)/3))
 367             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)
 368             return res
 369         except :
 370             return 0
 371     for h in V:
 372         if not ASYNC or  random() < taux_succes:
 373             n1 = dt.datetime.now()
 374             t= float(-3*la[h]+mt.sqrt(9*(la[h]**2)+64*delta*v[h]*mt.log(float(sigma2)/D)/gamma))/(16*delta)
 375             rep = mt.pow(t,float(3)/5)
 376             n2 = dt.datetime.now()
 377             cp = armin(f_Ps,Ps[h],POS,tuple([h]))
 378             n3 = dt.datetime.now()
 379             tp += (n2-n1).microseconds
 380             tpa += (n3-n2).microseconds
 381             Psp[h]=rep
 382         else :
 383             Psp[h] = Ps[h]
 384
 385
 386
 387     etap={}
 388
 389     for h in V:
 390         for i in N :
 391             etap[(h,i)] = etahi(h,i,Rh)
 392
 393
 394
 395     Rhp={}
 396     tr,tra=0,0
 397     def f_Rh(rh,h):
 398         return delta*rh*rh-v[h]*rh-sum([u[(h,i)]*eta[(h,i)] for i in N])
 399     for h in V:
 400         if not ASYNC or  random() < taux_succes:
 401             rep = float(v[h])/(2*delta)
 402             n1 = dt.datetime.now()
 403             rep = float(v[h])/(2*delta)
 404             n2 = dt.datetime.now()
 405             cp = armin(f_Rh,Rh[h],POS_NUL,tuple([h]))
 406             n3 = dt.datetime.now()
 407             tr += (n2-n1).microseconds
 408             tra += (n3-n2).microseconds
 409             Rhp[h] = 0 if rep < 0 else rep
 410         else :
 411             Rhp[h] = Rh[h]
 412
 413
 414
 415
 416     xp={}
 417     tx,txa=0,0
 418     def f_x(xhl,h,l):
 419         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]))
 420         return r
 421
 422     for h in V:
 423         for l in L:
 424             if not ASYNC or  random() < taux_succes:
 425                 n1 = dt.datetime.now()
 426                 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)
 427                 n2 = dt.datetime.now()
 428                 cp = armin(f_x,x[(h,l)],POS_NUL,tuple([h,l]))
 429                 n3 = dt.datetime.now()
 430                 tx += (n2-n1).microseconds
 431                 txa += (n3-n2).microseconds
 432
 433                 xp[(h,l)] = 0 if rep < 0 else rep
 434             else :
 435                 xp[(h,l)] = x[(h,l)]
 436
 437
 438
 439
 440
 441     valeurFonctionDualep = 0
 442     valeurFonctionDualep += sum([amplifieur*q[i]*q[i] for i in N])
 443     valeurFonctionDualep += sum([sum([delta*(x[(h,l)]**2) for l in L]) for h in V])
 444     valeurFonctionDualep += sum([delta*(Rh[h]**2) for h in V])
 445     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])
 446     valeurFonctionDualep += sum([v[h]*(mt.log(float(sigma2)/D)/(gamma*mt.pow(Ps[h],float(2)/3)) - Rh[h]) for h in V])
 447     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 ])
 448     valeurFonctionDualep += sum([w[l]*sum([a[i][l]*q[i] for i in N]) for l in L])
 449
 450
 451     #AfficheVariation(up,vp,lap,wp,thetap,etap,qp,Psp,Rhp,xp,valeurFonctionDualep)
 452
 453     arret = abs(valeurFonctionDuale-valeurFonctionDualep)
 454
 455     return (up,vp,lap,wp,thetap,etap,qp,Psp,Rhp,xp,valeurFonctionDualep,arret,mxg,smax,tq+tp+tr+tx,tqa+tpa+tra+txa)
 456
 457
 458
 459
 460 """
 461 print "u",u
 462 print "v",v
 463 print "lambda",la
 464 print "w",w
 465 print "theta",theta
 466 print "eta", eta
 467 print "q",q
 468 print "Ps",Ps
 469 print "Rh",Rh
 470 print "x",x
 471
 472 """
 473
 474
 475
 476 def initialisation():
 477     global u, v, la, w, theta , q,  Ps, Rh,  eta, x,init
 478
 479     theta = omega
 480
 481     q = {}
 482     for i in N :
 483         q[i] = 0.15 + random()*0.05
 484
 485
 486     Ps= {}
 487     for h in V :
 488         Ps[h] =  0.2+random()*0.3
 489
 490     Rh = {}
 491     for vi in V:
 492         Rh[vi] = 0.1 + random()*0.1
 493
 494     eta = {}
 495     for h in V :
 496         for i in N:
 497             eta[(h,i)]= etahi(h,i,Rh)
 498
 499     x={}
 500     for h in V :
 501         for l in L:
 502             x[(h,l)]=0
 503
 504
 505
 506
 507     # initialisation des operateurs lagrangiens
 508     u={}
 509     for h in V :
 510         for i in N:
 511             u[(h,i)]= random()
 512
 513     v = {}
 514     for h in V :
 515         v[h] = Rh[h]
 516
 517     la = {}
 518     for i in N:
 519         la[i] = random()
 520
 521     w={}
 522     for l in L:
 523         w[l] =  random()
 524
 525     init = [deepcopy(q),deepcopy(Ps),deepcopy(Rh),deepcopy(eta),
 526             deepcopy(x),deepcopy(u),deepcopy(v),deepcopy(la),deepcopy(w)]
 527
 528
 529
 530 def initialisation_():
 531     global u, v, la, w, theta , q,  Ps, Rh,  eta, x,init
 532     fd = open(fichier_init,"r")
 533     l= fd.readline()
 534     init_p = eval(l)
 535     print init_p
 536     theta = omega
 537     (q,Ps,Rh,eta,x,u,v,la,w) = tuple([deepcopy(x) for x in init_p])
 538     init = [deepcopy(q),deepcopy(Ps),deepcopy(Rh),deepcopy(eta),
 539             deepcopy(x),deepcopy(u),deepcopy(v),deepcopy(la),deepcopy(w)]
 540
 541
 542
 543 def __evalue_maj_theta__(nbexp,out=False):
 544     global u, v, la, w, theta , q,  Ps, Rh,  eta, x, valeurFonctionDuale
 545     res = {}
 546     m = []
 547     itermax = 100000
 548
 549     def __maj_theta(k):
 550         mem = []
 551         om = omega/(mt.pow(k,0.75))
 552         return om
 553     liste_arret=[]
 554     for idxexp in range(nbexp):
 555         mxg = 0
 556         if not(out):
 557             initialisation()
 558         else :
 559             initialisation_()
 560
 561         k = 1
 562         arret = False
 563         sm = 0
 564         err, ar = 0,0
 565         dur,dura=0,0
 566         while k < itermax  and not arret :
 567             print "ici"
 568             (u,v,la,w,theta,eta,q,Ps,Rh,x,valeurFonctionDuale,ar,mxg,smax,ct,cta)=maj(k,__maj_theta,mxg,idxexp)
 569             err =  (max(q.values()) - min(q.values()))/min(q.values())
 570             dur += ct
 571             dura += cta
 572             arret = err <  errorq or ar < errorD
 573             k+=1
 574             variation = "+" if smax > sm else "-"
 575             if out : print variation,
 576             if k%10 ==0 :
 577                 print "k:", k,
 578                 "erreur sur q",
 579                 errorq, "erreur sur F", ar, "et q:", q
 580
 581                 if out : print "maxg=", mxg
 582                 mem = [deepcopy(q),deepcopy(Ps),deepcopy(Rh),deepcopy(eta),
 583                        deepcopy(x),deepcopy(u),deepcopy(v),deepcopy(la),deepcopy(w)]
 584             """
 585             if k%4500 == 0 :
 586
 587                 #print "#########\n",mem,"\#########\n"
 588             if k%4600 == 0 :
 589                 #print "#########m\n",mem,"\#########\n"
 590             """
 591
 592
 593             if smax - sm  > 500:
 594                 print "variation trop grande"
 595                 print "init"
 596                 print init
 597                 exit
 598             sm = smax
 599
 600         if k == itermax:
 601             print "nbre d'iteration trop grand"
 602             print "init"
 603             print init
 604             sy.exit(1)
 605         else :
 606             liste_arret += [(err, ar,k,errorD)]
 607
 608         print "###############"
 609         print k
 610         print "###############"
 611         m += [k]
 612
 613     #print liste_arret
 614     #print (min(m),max(m),float(sum(m))/nbexp,m),m
 615     return (liste_arret,dur,dura)
 616
 617
 618 def __une_seule_exp__(fichier_donees):
 619     global u, v, la, w, theta , q,  Ps, Rh,  eta, x, valeurFonctionDuale
 620     initialisation()
 621
 622
 623     fichier = open(fichier_donees, "r")
 624     instructions ={}
 625     for line in fichier:
 626          l = line.split("=")
 627          instructions[l[0]] = eval(l[1])
 628     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']
 629     nbexp = 1
 630     res = {}
 631     m = []
 632     itermax = 100000
 633
 634     def __maj_theta(k):
 635         om = omega/(mt.pow(k,0.75))
 636         return om
 637     for idxexp in range(nbexp):
 638         mxg = 0
 639         k = 1
 640         arret = False
 641         sm = 0
 642         while k < itermax  and not arret :
 643             (u,v,la,w,theta,eta,q,Ps,Rh,x,valeurFonctionDuale,ar,mxg,smax)=maj(k,__maj_theta,mxg,idxexp)
 644             errorq =  (max(q.values()) - min(q.values()))/min(q.values())
 645             arret = errorq <  error
 646             k+=1
 647             variation = "+" if smax > sm else "-"
 648             print variation,
 649             if k%100 ==0 :
 650                 print "k:",k,"erreur sur q", errorq, "et q:",q
 651                 print "maxg=", mxg
 652             if smax - sm  > 500:
 653                 print "variation trop grande"
 654                 print "init"
 655                 print init
 656                 exit
 657             sm = smax
 658
 659         if k == itermax:
 660             print "nbre d'iteration trop grand"
 661             print "init"
 662             print init
 663             exit
 664
 665         print "###############"
 666         print k
 667         print "###############"
 668         m += [k]
 669
 670     print (min(m),max(m),float(sum(m))/nbexp,m),m
 671
 672
 673
 674
 675 ASYNC = False
 676 #__une_seule_exp__("config_initiale.py")
 677 #pour evaluerle test d'arret
 678 """
 679 listederes=[]
 680 for k in list(np.logspace(-6, -3, num=15)):
 681     errorD = k
 682     listederes += __evalue_maj_theta__(5)
 683
 684 print listederes
 685 """
 686
 687
 688
 689 # pour evaluer argmin
 690 listederes=[]
 691 k=0
 692 for k in list(np.logspace(-5, -3, num=10)):
 693     print "k",k
 694     errorD = k
 695     listederes += [(k,__evalue_maj_theta__(3))]
 696     k+=1
 697 print listederes
 698
 699
 700
 701 print"#######################"
 702
 703
 704 #ASYNC = True
 705 #taux_succes = 0.9
 706 #__evalue_maj_theta__()
 707
 708
 709
 710
 711 """
 712 print "u",u
 713 print "v",v
 714 print "lambda",la
 715 print "w",w
 716 print "theta",theta
 717 print "eta", eta
 718 print "q",q
 719 print "Ps",Ps
 720 print "Rh",Rh
 721 print "x",x
 722 print "L",valeurFonctionDuale
 723
 724 """
 725
 726 # relation ente les variables primaires et secondaires ?
 727