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