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-2
  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,comppsh=False):
 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             #calcul nouvelle solution
 375             if comppsh :
 376                 if la[h] != 0 :
 377                     t= float(2*v[h]*mt.log(float(sigma2)/D))/(3*gamma*la[h])
 378                     rep = mt.pow(t,float(3)/5)
 379                 else :
 380                     rep = 1000
 381             else :
 382                 t= float(-3*la[h]+mt.sqrt(9*(la[h]**2)+64*delta*v[h]*mt.log(float(sigma2)/D)/gamma))/(16*delta)
 383                 rep = mt.pow(t,float(3)/5)
 384
 385
 386
 387             n2 = dt.datetime.now()
 388             #cp = armin(f_Ps,Ps[h],POS,tuple([h]))
 389             n3 = dt.datetime.now()
 390             tp += (n2-n1).microseconds
 391             tpa += (n3-n2).microseconds
 392             Psp[h]=rep
 393         else :
 394             Psp[h] = Ps[h]
 395
 396
 397
 398     etap={}
 399
 400     for h in V:
 401         for i in N :
 402             etap[(h,i)] = etahi(h,i,Rh)
 403
 404
 405
 406     Rhp={}
 407     tr,tra=0,0
 408     def f_Rh(rh,h):
 409         return delta*rh*rh-v[h]*rh-sum([u[(h,i)]*eta[(h,i)] for i in N])
 410     for h in V:
 411         if not ASYNC or  random() < taux_succes:
 412             rep = float(v[h])/(2*delta)
 413             n1 = dt.datetime.now()
 414             rep = float(v[h])/(2*delta)
 415             n2 = dt.datetime.now()
 416             #cp = armin(f_Rh,Rh[h],POS_NUL,tuple([h]))
 417             n3 = dt.datetime.now()
 418             tr += (n2-n1).microseconds
 419             tra += (n3-n2).microseconds
 420             Rhp[h] = 0 if rep < 0 else rep
 421         else :
 422             Rhp[h] = Rh[h]
 423
 424
 425
 426
 427     xp={}
 428     tx,txa=0,0
 429     def f_x(xhl,h,l):
 430         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]))
 431         return r
 432
 433     for h in V:
 434         for l in L:
 435             if not ASYNC or  random() < taux_succes:
 436                 n1 = dt.datetime.now()
 437                 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)
 438                 n2 = dt.datetime.now()
 439                 #cp = armin(f_x,x[(h,l)],POS_NUL,tuple([h,l]))
 440                 n3 = dt.datetime.now()
 441                 tx += (n2-n1).microseconds
 442                 txa += (n3-n2).microseconds
 443
 444                 xp[(h,l)] = 0 if rep < 0 else rep
 445             else :
 446                 xp[(h,l)] = x[(h,l)]
 447
 448
 449
 450
 451
 452     valeurFonctionDualep = 0
 453     valeurFonctionDualep += sum([amplifieur*q[i]*q[i] for i in N])
 454     valeurFonctionDualep += sum([sum([delta*(x[(h,l)]**2) for l in L]) for h in V])
 455     valeurFonctionDualep += sum([delta*(Rh[h]**2) for h in V])
 456     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])
 457     valeurFonctionDualep += sum([v[h]*(mt.log(float(sigma2)/D)/(gamma*mt.pow(Ps[h],float(2)/3)) - Rh[h]) for h in V])
 458     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 ])
 459     valeurFonctionDualep += sum([w[l]*sum([a[i][l]*q[i] for i in N]) for l in L])
 460
 461
 462     #AfficheVariation(up,vp,lap,wp,thetap,etap,qp,Psp,Rhp,xp,valeurFonctionDualep)
 463
 464     arret = abs(valeurFonctionDuale-valeurFonctionDualep)
 465
 466     return (up,vp,lap,wp,thetap,etap,qp,Psp,Rhp,xp,valeurFonctionDualep,arret,mxg,smax,tq+tp+tr+tx,tqa+tpa+tra+txa)
 467
 468
 469
 470
 471 """
 472 print "u",u
 473 print "v",v
 474 print "lambda",la
 475 print "w",w
 476 print "theta",theta
 477 print "eta", eta
 478 print "q",q
 479 print "Ps",Ps
 480 print "Rh",Rh
 481 print "x",x
 482
 483 """
 484
 485
 486
 487 def initialisation():
 488     global u, v, la, w, theta , q,  Ps, Rh,  eta, x,init
 489
 490     theta = omega
 491
 492     q = {}
 493     for i in N :
 494         q[i] = 0.15 + random()*0.05
 495
 496
 497     Ps= {}
 498     for h in V :
 499         Ps[h] =  0.2+random()*0.3
 500
 501     Rh = {}
 502     for vi in V:
 503         Rh[vi] = 0.1 + random()*0.1
 504
 505     eta = {}
 506     for h in V :
 507         for i in N:
 508             eta[(h,i)]= etahi(h,i,Rh)
 509
 510     x={}
 511     for h in V :
 512         for l in L:
 513             x[(h,l)]=0
 514
 515
 516
 517
 518     # initialisation des operateurs lagrangiens
 519     u={}
 520     for h in V :
 521         for i in N:
 522             u[(h,i)]= random()
 523
 524     v = {}
 525     for h in V :
 526         v[h] = Rh[h]
 527
 528     la = {}
 529     for i in N:
 530         la[i] = random()
 531
 532     w={}
 533     for l in L:
 534         w[l] =  random()
 535
 536     init = [deepcopy(q),deepcopy(Ps),deepcopy(Rh),deepcopy(eta),
 537             deepcopy(x),deepcopy(u),deepcopy(v),deepcopy(la),deepcopy(w)]
 538
 539
 540
 541 def initialisation_():
 542     global u, v, la, w, theta , q,  Ps, Rh,  eta, x,init
 543     fd = open(fichier_init,"r")
 544     l= fd.readline()
 545     init_p = eval(l)
 546     print init_p
 547     theta = omega
 548     (q,Ps,Rh,eta,x,u,v,la,w) = tuple([deepcopy(x) for x in init_p])
 549     init = [deepcopy(q),deepcopy(Ps),deepcopy(Rh),deepcopy(eta),
 550             deepcopy(x),deepcopy(u),deepcopy(v),deepcopy(la),deepcopy(w)]
 551
 552
 553
 554 def __evalue_maj_theta__(nbexp,out=False):
 555     global u, v, la, w, theta , q,  Ps, Rh,  eta, x, valeurFonctionDuale
 556     res = {}
 557     m = []
 558     itermax = 100000
 559
 560     def __maj_theta(k):
 561         mem = []
 562         om = omega/(mt.pow(k,0.75))
 563         return om
 564     liste_arret=[]
 565     for idxexp in range(nbexp):
 566         mxg = 0
 567         if not(out):
 568             initialisation()
 569         else :
 570             initialisation_()
 571
 572         k = 1
 573         arret = False
 574         sm = 0
 575         err, ar = 0,0
 576         dur,dura=0,0
 577         while k < itermax  and not arret :
 578             print "ici"
 579             (u,v,la,w,theta,eta,q,Ps,Rh,x,valeurFonctionDuale,ar,mxg,smax,ct,cta)=maj(k,__maj_theta,mxg,idxexp)
 580             err =  (max(q.values()) - min(q.values()))/min(q.values())
 581             dur += ct
 582             dura += cta
 583             arret = err <  errorq or ar < errorD
 584             k+=1
 585             variation = "+" if smax > sm else "-"
 586             if out : print variation,
 587             if k%10 ==0 :
 588                 print "k:", k,
 589                 "erreur sur q",
 590                 errorq, "erreur sur F", ar, "et q:", q
 591
 592                 if out : print "maxg=", mxg
 593                 mem = [deepcopy(q),deepcopy(Ps),deepcopy(Rh),deepcopy(eta),
 594                        deepcopy(x),deepcopy(u),deepcopy(v),deepcopy(la),deepcopy(w)]
 595             """
 596             if k%4500 == 0 :
 597
 598                 #print "#########\n",mem,"\#########\n"
 599             if k%4600 == 0 :
 600                 #print "#########m\n",mem,"\#########\n"
 601             """
 602
 603
 604             if smax - sm  > 500:
 605                 print "variation trop grande"
 606                 print "init"
 607                 print init
 608                 exit
 609             sm = smax
 610
 611         if k == itermax:
 612             print "nbre d'iteration trop grand"
 613             print "init"
 614             print init
 615             sy.exit(1)
 616         else :
 617             liste_arret += [(err, ar,k,errorD)]
 618
 619         print "###############"
 620         print k
 621         print "###############"
 622         m += [k]
 623
 624     #print liste_arret
 625     #print (min(m),max(m),float(sum(m))/nbexp,m),m
 626     #return (liste_arret,dur,dura)
 627     return (liste_arret,dur)
 628
 629
 630 def __une_seule_exp__(fichier_donees):
 631     global u, v, la, w, theta , q,  Ps, Rh,  eta, x, valeurFonctionDuale
 632     initialisation()
 633
 634
 635     fichier = open(fichier_donees, "r")
 636     instructions ={}
 637     for line in fichier:
 638          l = line.split("=")
 639          instructions[l[0]] = eval(l[1])
 640     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']
 641     nbexp = 1
 642     res = {}
 643     m = []
 644     itermax = 100000
 645
 646     def __maj_theta(k):
 647         om = omega/(mt.pow(k,0.75))
 648         return om
 649     for idxexp in range(nbexp):
 650         mxg = 0
 651         k = 1
 652         arret = False
 653         sm = 0
 654         while k < itermax  and not arret :
 655             (u,v,la,w,theta,eta,q,Ps,Rh,x,valeurFonctionDuale,ar,mxg,smax)=maj(k,__maj_theta,mxg,idxexp)
 656             errorq =  (max(q.values()) - min(q.values()))/min(q.values())
 657             arret = errorq <  error
 658             k+=1
 659             variation = "+" if smax > sm else "-"
 660             print variation,
 661             if k%100 ==0 :
 662                 print "k:",k,"erreur sur q", errorq, "et q:",q
 663                 print "maxg=", mxg
 664             if smax - sm  > 500:
 665                 print "variation trop grande"
 666                 print "init"
 667                 print init
 668                 exit
 669             sm = smax
 670
 671         if k == itermax:
 672             print "nbre d'iteration trop grand"
 673             print "init"
 674             print init
 675             exit
 676
 677         print "###############"
 678         print k
 679         print "###############"
 680         m += [k]
 681
 682     print (min(m),max(m),float(sum(m))/nbexp,m),m
 683
 684
 685
 686
 687 ASYNC = False
 688 #__une_seule_exp__("config_initiale.py")
 689 #pour evaluerle test d'arret
 690 """
 691 listederes=[]
 692 for k in list(np.logspace(-6, -3, num=15)):
 693     errorD = k
 694     listederes += __evalue_maj_theta__(5)
 695
 696 print listederes
 697 """
 698
 699
 700
 701
 702 def __comparePsh_et_Psh_avec8_3__(nbexp,out=False):
 703     print "tttttttttttt"
 704     global u, v, la, w, theta , q,  Ps, Rh,  eta, x, valeurFonctionDuale
 705     res = {}
 706     m,mp = [],[]
 707     itermax = 100000
 708
 709     def __maj_theta(k):
 710         mem = []
 711         om = omega/(mt.pow(k,0.75))
 712         return om
 713     liste_arret=[]
 714     for idxexp in range(nbexp):
 715         mxg = 0
 716         if not(out):
 717             initialisation()
 718         else :
 719             initialisation_()
 720
 721         k = 1
 722         arret = False
 723         sm = 0
 724         err, ar = 0,0
 725         dur,dura=0,0
 726
 727         # duppliquee
 728         while k < itermax  and not arret :
 729             (u,v,la,w,theta,eta,q,Ps,Rh,x,valeurFonctionDuale,ar,mxg,smax,ct,cta)=maj(k,__maj_theta,mxg,idxexp)
 730             err =  (max(q.values()) - min(q.values()))/min(q.values())
 731             dur += ct
 732             dura += cta
 733             arret = err <  errorq or ar < errorD
 734             k+=1
 735             variation = "+" if smax > sm else "-"
 736             if out : print variation,
 737             if k%100 ==0 :
 738                 print "k:", k,
 739                 "erreur sur q",
 740                 errorq, "erreur sur F", ar, "et q:", q
 741
 742                 if out : print "maxg=", mxg
 743                 mem = [deepcopy(q),deepcopy(Ps),deepcopy(Rh),deepcopy(eta),
 744                        deepcopy(x),deepcopy(u),deepcopy(v),deepcopy(la),deepcopy(w)]
 745             """
 746             if k%4500 == 0 :
 747
 748                 #print "#########\n",mem,"\#########\n"
 749             if k%4600 == 0 :
 750                 #print "#########m\n",mem,"\#########\n"
 751             """
 752
 753
 754             if smax - sm  > 500:
 755                 print "variation trop grande"
 756                 print "init"
 757                 print init
 758                 exit
 759             sm = smax
 760
 761         if k == itermax:
 762             print "nbre d'iteration trop grand"
 763             print "init"
 764             print init
 765             sy.exit(1)
 766         else :
 767             liste_arret += [(err, ar,k,errorD)]
 768
 769         print "###############"
 770         print k," avec optym"
 771         print "###############"
 772         m += [k]
 773         ##fin de duplication
 774
 775
 776
 777         #remise a zero
 778         q,Ps,Rh,eta,x,u,v,la,w=init[0],init[1],init[2],init[3],init[4],init[5],init[6],init[7],init[8]
 779         k = 1
 780         arret = False
 781         sm = 0
 782         err, ar = 0,0
 783         dur,dura=0,0
 784         # duppliquee
 785         while k < itermax  and not arret :
 786             (u,v,la,w,theta,eta,q,Ps,Rh,x,valeurFonctionDuale,ar,mxg,smax,ct,cta)=maj(k,__maj_theta,mxg,idxexp,True)
 787             err =  (max(q.values()) - min(q.values()))/min(q.values())
 788             dur += ct
 789             dura += cta
 790             arret = err <  errorq or ar < errorD
 791             k+=1
 792             variation = "+" if smax > sm else "-"
 793             if out : print variation,
 794             if k%100 ==0 :
 795                 print "k:", k,
 796                 "erreur sur q",
 797                 errorq, "erreur sur F", ar, "et q:", q
 798
 799                 if out : print "maxg=", mxg
 800                 mem = [deepcopy(q),deepcopy(Ps),deepcopy(Rh),deepcopy(eta),
 801                        deepcopy(x),deepcopy(u),deepcopy(v),deepcopy(la),deepcopy(w)]
 802             """
 803             if k%4500 == 0 :
 804
 805                 #print "#########\n",mem,"\#########\n"
 806             if k%4600 == 0 :
 807                 #print "#########m\n",mem,"\#########\n"
 808             """
 809
 810
 811             if smax - sm  > 500:
 812                 print "variation trop grande"
 813                 print "init"
 814                 print init
 815                 exit
 816             sm = smax
 817
 818         if k == itermax:
 819             print "nbre d'iteration trop grand"
 820             print "init"
 821             print init
 822             sy.exit(1)
 823         else :
 824             liste_arret += [(err, ar,k,errorD)]
 825
 826         print "###############"
 827         print k,"sans optym"
 828         print "###############"
 829         mp += [k]
 830
 831
 832     #print liste_arret
 833     #print (min(m),max(m),float(sum(m))/nbexp,m),m
 834     #return (liste_arret,dur,dura)
 835     return (m,mp)
 836
 837
 838
 839 # pour evaluer argmin
 840 """
 841 listederes=[]
 842 k=0
 843 for k in list(np.logspace(-5, -3, num=10)):
 844     print "k",k
 845     errorD = k
 846     listederes += [(k,__evalue_maj_theta__(3))]
 847     k+=1
 848 print listederes
 849 """
 850
 851 # pour evaluer psh
 852 listederes=[]
 853 for k in list(np.logspace(-5, -3, num=10)):
 854     print "Precision",k
 855     errorD = k
 856     listederes += [(k,__comparePsh_et_Psh_avec8_3__(10))]
 857     k+=1
 858 print listederes
 859
 860
 861
 862
 863
 864
 865 print"#######################"
 866
 867
 868 #ASYNC = True
 869 #taux_succes = 0.9
 870 #__evalue_maj_theta__()
 871
 872
 873
 874
 875 """
 876 print "u",u
 877 print "v",v
 878 print "lambda",la
 879 print "w",w
 880 print "theta",theta
 881 print "eta", eta
 882 print "q",q
 883 print "Ps",Ps
 884 print "Rh",Rh
 885 print "x",x
 886 print "L",valeurFonctionDuale
 887
 888 """
 889
 890 # relation ente les variables primaires et secondaires ?
 891