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