totot

[desynchronisation-controle.git] / exp_controle_asynchrone / simulMWSN.py

diff --git a/exp_controle_asynchrone/simulMWSN.py b/exp_controle_asynchrone/simulMWSN.py
index 4844593705e7776e1bbc936dfd621a089fbb6a28..8e286a33e631f1e1832c862cae9079f5091ffc00 100644 (file)
--- a/exp_controle_asynchrone/simulMWSN.py
+++ b/exp_controle_asynchrone/simulMWSN.py
@@ -7,9 +7,12 @@ from itertools import *
 from scipy import optimize as opt
 from copy import deepcopy 
 import sys as sy
 from scipy import optimize as opt
 from copy import deepcopy 
 import sys as sy

+import cv as cv
+import cv2 as cv2 
+import datetime as dt

 
 

-
-error  = 0.1
+errorq  = 0.1
+errorD  = 1E-2

 epsilon = 1E-10
 vrate = 0.8
 p = 0.7
 epsilon = 1E-10
 vrate = 0.8
 p = 0.7


@@ -32,6 +35,9 @@ lg = [(0, 1, 22.004323820359151), (0, 2, 28.750632705280324), (0, 3, 29.68069293
 #lg= [(0,1,23),(1,0,15),(1,2,45)]
 sink = n-1
 
 #lg= [(0,1,23),(1,0,15),(1,2,45)]
 sink = n-1
 

+def distance(d1,d2):
+    return mt.sqrt(sum([(d1[t]-d2[t])**2 for t in d1]))
+

 
 def genereGraph():
     test = False 
 
 def genereGraph():
     test = False 


@@ -46,7 +52,54 @@ def genereGraph():
                         G.add_edge(io,ie,weight=dist)
                         G.add_edge(ie,io,weight=dist)
         test = not(any([ not(nx.has_path(G,o,sink)) for o in  G.nodes() if sink in G.nodes() and o != sink]))
                         G.add_edge(io,ie,weight=dist)
                         G.add_edge(ie,io,weight=dist)
         test = not(any([ not(nx.has_path(G,o,sink)) for o in  G.nodes() if sink in G.nodes() and o != sink]))

-    return G
+    return (G,l)
+
+
+def afficheGraph(G,l,tx,ty,sink):
+    r = 20
+    img = cv.CreateImage ((tx, ty), 32, 3)
+    cv.Rectangle(img, (0,0),(tx,ty), cv.Scalar(255,255,255), thickness=-1)
+    def px((x,y)): 
+        u = float(tx)/(coteCarre + 2*distanceEmmissionMax)
+        return(int(distanceEmmissionMax*u + x * u),int(distanceEmmissionMax*u + y * u))
+    for i in set(range(len(l)))-set([sink]):
+        (x,y) = l[i]
+        pix,piy = px((x,y))
+        demx = distanceEmmissionMax*tx/(coteCarre+2*distanceEmmissionMax)
+        cv.Circle(img, (pix,piy),demx, cv.Scalar(125,125,125))     
+    
+    for i in set(range(len(l)))-set([sink]):        
+        (x,y) = l[i]
+        pix,piy = px((x,y))
+        cv.Circle(img, (pix,piy),r, cv.Scalar(125,125,125),thickness=-1)     
+    
+    #sink
+    (x,y) = l[sink]
+    pix,piy = px((x,y))
+
+    cv.Rectangle(img, (pix-r/2,piy-r/2),(pix+r/2,piy+r/2), cv.Scalar(125,125,125), thickness=-1)
+
+    for i in range(len(l)):
+        for j in range(len(l)):
+            
+            if np.linalg.norm(np.array(l[i])-np.array(l[j])) < distanceEmmissionMax :
+                (xi,yi) = l[i]
+                pixi,piyi = px((xi,yi))
+                (xj,yj) = l[j]
+                pixj,piyj = px((xj,yj))
+                cv.Line(img, (pixi,piyi), (pixj,piyj), cv.Scalar(125,125,125))
+    
+
+    """
+    for i in range(len(l)):
+        (x,y) = l[i]
+        pix,piy = px((x,y))
+        print i
+        textColor = (0, 0, 255)  # red
+        font = cv2.FONT_HERSHEY_SIMPLEX
+        imgp = 
+        cv2.putText(img, str(i), (pix-r/4,piy-r/2),font, 3.0, textColor)#,thickn    """
+    cv.SaveImage("SensorNetwork.png",img)

 
 G = nx.DiGraph()
 G.add_weighted_edges_from(lg)
 
 G = nx.DiGraph()
 G.add_weighted_edges_from(lg)


@@ -56,6 +109,19 @@ G.add_weighted_edges_from(lg)
 #nx.draw(G)
 #pb.show()
 
 #nx.draw(G)
 #pb.show()
 

+#(G,l) = genereGraph()
+N = G.nodes()
+#V = list(set(sample(N,int(len(N)*vrate)))-set([sink]))
+V = list(set(N)-set([sink]))
+source = V
+print "source",source
+
+
+#afficheGraph(G,l,500,500,sink)
+#nx.draw(G)
+#pb.show()
+
+    

 
 
 
 
 
 


@@ -64,15 +130,8 @@ G.add_weighted_edges_from(lg)
 #print G.edges(data=True)
 #TODO afficher le graphe et etre sur qu'il est connexe
 
 #print G.edges(data=True)
 #TODO afficher le graphe et etre sur qu'il est connexe
 

-N = G.nodes()
-
-
-#V = list(set(sample(N,int(len(N)*vrate)))-set([sink]))
-V = list(set(N)-set([sink]))

 
 
 
 

-source = V
-print "source",source

 
 
 L = range(len(G.edges()))
 
 
 L = range(len(G.edges()))


@@ -107,7 +166,6 @@ alpha = 0.5
 beta = 1.3E-8
 gamma = 55.54
 delta = 0.2
 beta = 1.3E-8
 gamma = 55.54
 delta = 0.2

-zeta = 0.1

 amplifieur = 1
 sigma2 = 3500
 Bi = 5
 amplifieur = 1
 sigma2 = 3500
 Bi = 5


@@ -146,8 +204,6 @@ def aminp(l):
 
 
 
 
 
 

-def distance(d1,d2):
-    return mt.sqrt(sum([(d1[t]-d2[t])**2 for t in d1]))

 
 
 def AfficheVariation (up,vp,lap,wp,thetap,etap,qp,Psp,Rhp,xp,valeurFonctionDualep):
 
 
 def AfficheVariation (up,vp,lap,wp,thetap,etap,qp,Psp,Rhp,xp,valeurFonctionDualep):


@@ -187,19 +243,11 @@ def armin(f,xini,xr,args):
     #xpos = lambda x : x if x > 0 else -1 
     r = 0
     if xr == POS :
     #xpos = lambda x : x if x > 0 else -1 
     r = 0
     if xr == POS :

-        #print "strictement pos"
-        #print "parametres passes a cobyla",xini,xpos,args,"consargs=(),rhoend=1E-5,iprint=0,maxfun=1000"
-        r= opt.fmin_cobyla(f,xini,cons=[xpos],args=args,consargs=(),rhoend=1E-3,iprint=0,maxfun=nbiter)
-        #print "le min str pos est",r 
-        #print "le min pos   est", r,"avec xini",xini
+        r= opt.fmin_cobyla(f,xini,cons=[xpos],args=args,consargs=(),rhoend=errorD/100,iprint=0,maxfun=nbiter)

     elif xr == POS_NUL :
     elif xr == POS_NUL :

-        #print "pos ou nul"
-        r = opt.fmin_cobyla(f,xini,[xposounul],args,consargs=(),rhoend=1E-3,iprint=0,maxfun=nbiter)
-        # print "le min pos est", r
-        #print "le min pos null  est", r,"avec xini",xini
+        r = opt.fmin_cobyla(f,xini,[xposounul],args,consargs=(),rhoend=errorD/100,iprint=0,maxfun=nbiter)

     elif xr == POSINF1:
     elif xr == POSINF1:

-        r = opt.fmin_cobyla(f,xini,[entre0et1],args,consargs=(),rhoend=1E-3,iprint=0,maxfun=nbiter)
-        #print "le min pos inf 1 est", r,"avec xini",xini    
+        r = opt.fmin_cobyla(f,xini,[entre0et1],args,consargs=(),rhoend=errorD/100,iprint=0,maxfun=nbiter)

 
     
     return r
 
     
     return r


@@ -220,7 +268,7 @@ def maj_theta(k):
 
 
 
 
 
 

-def maj(k,maj_theta,mxg,idxexp):
+def maj(k,maj_theta,mxg,idxexp,comppsh=False):

     # initialisation des operateurs lagrangiens
     global u, v, la, w, theta , q,  Ps, Rh,  eta, x, valeurFonctionDuale
     
     # initialisation des operateurs lagrangiens
     global u, v, la, w, theta , q,  Ps, Rh,  eta, x, valeurFonctionDuale
     


@@ -290,39 +338,59 @@ def maj(k,maj_theta,mxg,idxexp):
 
 
     qp={}
 
 
     qp={}

-    def f_q(qi,i):
+    def f_q(qi,*args):

         fa = sum([a[i][l]*w[l] for l in L]) - la[i]*Bi 
         return qi*qi + qi*fa    
 
         fa = sum([a[i][l]*w[l] for l in L]) - la[i]*Bi 
         return qi*qi + qi*fa    
 

+    tq,tqa = 0,0

     for i in N:
         if not ASYNC or  random() < taux_succes:
     for i in N:
         if not ASYNC or  random() < taux_succes:

+            n1 = dt.datetime.now()

             c = -float(sum([a[i][l]*w[l] for l in L]) - la[i]*Bi)/(2*amplifieur)
             c = -float(sum([a[i][l]*w[l] for l in L]) - la[i]*Bi)/(2*amplifieur)

+            n2 = dt.datetime.now()
+            #cp = armin(f_q,q[i],POS,tuple([i]))
+            n3 = dt.datetime.now()
+            tq += (n2-n1).microseconds
+            tqa += (n3-n2).microseconds

             rep = epsilon if c <= 0 else c
             qp[i] = rep
         else :
             qp[i] = q[i]
 
 
             rep = epsilon if c <= 0 else c
             qp[i] = rep
         else :
             qp[i] = q[i]
 
 

-        

 
 

- 
+    tp,tpa = 0,0

     Psp={}
     Psp={}

-    #print "maj des des Psh" 

     def f_Ps(psh,h):
     def f_Ps(psh,h):

-        #print "ds f_ps",psh, v[h]* mt.log(float(sigma2)/D)/(gamma*((psh**2)**(float(2)/3))) +la[h]*psh
-         return v[h]* mt.log(float(sigma2)/D)/(gamma*mt.pow(float(2)/3)) +la[h]*psh
-     for h in V:
-         if not ASYNC or  random() < taux_succes:
-             """
-             lah = 0.05 if la[h] == 0 else  la[h]
-             rep = mt.pow(float(2*v[h]*mt.log(float(sigma2)/D))/(3*gamma*lah),float(3)/5)
-             Psp[h] = epsilon if rep <= 0 else rep
-             """
-             t= float(-3*la[h]+mt.sqrt(9*(la[h]**2)+64*zeta*v[h]*mt.log(float(sigma2)/D)/gamma))/(16*zeta)
-             #print t
-             rep = mt.pow(t,float(3)/5)
-             Psp[h]=rep
-         else :
+        try:
+            #print psh,(gamma*mt.pow(psh,float(2)/3))
+            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)
+            return res
+        except : 
+            return 0
+    for h in V:
+        if not ASYNC or  random() < taux_succes:
+            n1 = dt.datetime.now()
+            #calcul nouvelle solution 
+            if comppsh :
+                if la[h] != 0 :
+                    t= float(2*v[h]*mt.log(float(sigma2)/D))/(3*gamma*la[h])
+                    rep = mt.pow(t,float(3)/5)                
+                else :
+                    rep = 1000
+            else :                
+                t= float(-3*la[h]+mt.sqrt(9*(la[h]**2)+64*delta*v[h]*mt.log(float(sigma2)/D)/gamma))/(16*delta)
+                rep = mt.pow(t,float(3)/5)
+
+            
+
+            n2 = dt.datetime.now()
+            #cp = armin(f_Ps,Ps[h],POS,tuple([h]))
+            n3 = dt.datetime.now()
+            tp += (n2-n1).microseconds
+            tpa += (n3-n2).microseconds
+            Psp[h]=rep
+        else :

             Psp[h] = Ps[h]
 
 
             Psp[h] = Ps[h]
 
 


@@ -336,26 +404,43 @@ def maj(k,maj_theta,mxg,idxexp):
 
 
     Rhp={}
 
 
     Rhp={}

+    tr,tra=0,0

     def f_Rh(rh,h):
         return delta*rh*rh-v[h]*rh-sum([u[(h,i)]*eta[(h,i)] for i in N])
     def f_Rh(rh,h):
         return delta*rh*rh-v[h]*rh-sum([u[(h,i)]*eta[(h,i)] for i in N])

-

     for h in V:
         if not ASYNC or  random() < taux_succes:
             rep = float(v[h])/(2*delta)
     for h in V:
         if not ASYNC or  random() < taux_succes:
             rep = float(v[h])/(2*delta)

+            n1 = dt.datetime.now()
+            rep = float(v[h])/(2*delta)
+            n2 = dt.datetime.now()
+            #cp = armin(f_Rh,Rh[h],POS_NUL,tuple([h]))
+            n3 = dt.datetime.now()
+            tr += (n2-n1).microseconds
+            tra += (n3-n2).microseconds

             Rhp[h] = 0 if rep < 0 else rep
         else : 
             Rhp[h] = Rh[h]
           
             Rhp[h] = 0 if rep < 0 else rep
         else : 
             Rhp[h] = Rh[h]
           

+
+
+

     xp={}
     xp={}

+    tx,txa=0,0

     def f_x(xhl,h,l):
         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]))
         return r
     def f_x(xhl,h,l):
         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]))
         return r

-

     
     for h in V:
         for l in L:
             if not ASYNC or  random() < taux_succes:
     
     for h in V:
         for l in L:
             if not ASYNC or  random() < taux_succes:

+                n1 = dt.datetime.now()

                 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)
                 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)

+                n2 = dt.datetime.now()
+                #cp = armin(f_x,x[(h,l)],POS_NUL,tuple([h,l]))
+                n3 = dt.datetime.now()
+                tx += (n2-n1).microseconds
+                txa += (n3-n2).microseconds
+                

                 xp[(h,l)] = 0 if rep < 0 else rep 
             else :
                 xp[(h,l)] = x[(h,l)] 
                 xp[(h,l)] = 0 if rep < 0 else rep 
             else :
                 xp[(h,l)] = x[(h,l)] 


@@ -376,9 +461,9 @@ def maj(k,maj_theta,mxg,idxexp):
 
     #AfficheVariation(up,vp,lap,wp,thetap,etap,qp,Psp,Rhp,xp,valeurFonctionDualep)
 
 
     #AfficheVariation(up,vp,lap,wp,thetap,etap,qp,Psp,Rhp,xp,valeurFonctionDualep)
 

-    arret = abs(valeurFonctionDuale-valeurFonctionDualep) < error
+    arret = abs(valeurFonctionDuale-valeurFonctionDualep) 

 
 

-    return (up,vp,lap,wp,thetap,etap,qp,Psp,Rhp,xp,valeurFonctionDualep,arret,mxg,smax)
+    return (up,vp,lap,wp,thetap,etap,qp,Psp,Rhp,xp,valeurFonctionDualep,arret,mxg,smax,tq+tp+tr+tx,tqa+tpa+tra+txa)

 
 
 
 
 
 


@@ -471,11 +556,12 @@ def __evalue_maj_theta__(nbexp,out=False):
     res = {}
     m = []
     itermax = 100000
     res = {}
     m = []
     itermax = 100000

- 
+

     def __maj_theta(k):
         mem = []
         om = omega/(mt.pow(k,0.75))
         return om
     def __maj_theta(k):
         mem = []
         om = omega/(mt.pow(k,0.75))
         return om

+    liste_arret=[]

     for idxexp in range(nbexp):
         mxg = 0
         if not(out):
     for idxexp in range(nbexp):
         mxg = 0
         if not(out):


@@ -483,6 +569,85 @@ def __evalue_maj_theta__(nbexp,out=False):
         else :
             initialisation_()
             
         else :
             initialisation_()
             

+        k = 1
+        arret = False
+        sm = 0
+        err, ar = 0,0
+        dur,dura=0,0
+        while k < itermax  and not arret : 
+            print "ici"
+            (u,v,la,w,theta,eta,q,Ps,Rh,x,valeurFonctionDuale,ar,mxg,smax,ct,cta)=maj(k,__maj_theta,mxg,idxexp)
+            err =  (max(q.values()) - min(q.values()))/min(q.values())
+            dur += ct
+            dura += cta
+            arret = err <  errorq or ar < errorD
+            k+=1
+            variation = "+" if smax > sm else "-"
+            if out : print variation,
+            if k%10 ==0 :
+                print "k:", k, 
+                "erreur sur q", 
+                errorq, "erreur sur F", ar, "et q:", q
+
+                if out : print "maxg=", mxg
+                mem = [deepcopy(q),deepcopy(Ps),deepcopy(Rh),deepcopy(eta),
+                       deepcopy(x),deepcopy(u),deepcopy(v),deepcopy(la),deepcopy(w)]
+            """
+            if k%4500 == 0 :
+                
+                #print "#########\n",mem,"\#########\n"
+            if k%4600 == 0 :
+                #print "#########m\n",mem,"\#########\n"
+            """
+                
+
+            if smax - sm  > 500:
+                print "variation trop grande"
+                print "init"
+                print init
+                exit 
+            sm = smax
+        
+        if k == itermax:
+            print "nbre d'iteration trop grand"
+            print "init"
+            print init
+            sy.exit(1)
+        else :
+            liste_arret += [(err, ar,k,errorD)]
+
+        print "###############"
+        print k
+        print "###############"
+        m += [k]
+ 
+    #print liste_arret    
+    #print (min(m),max(m),float(sum(m))/nbexp,m),m
+    #return (liste_arret,dur,dura)
+    return (liste_arret,dur)
+
+
+def __une_seule_exp__(fichier_donees):
+    global u, v, la, w, theta , q,  Ps, Rh,  eta, x, valeurFonctionDuale 
+    initialisation()
+        
+
+    fichier = open(fichier_donees, "r")
+    instructions ={}
+    for line in fichier:    
+         l = line.split("=")
+         instructions[l[0]] = eval(l[1])
+    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']
+    nbexp = 1
+    res = {}
+    m = []
+    itermax = 100000
+ 
+    def __maj_theta(k):
+        om = omega/(mt.pow(k,0.75))
+        return om
+    for idxexp in range(nbexp):
+        mxg = 0

         k = 1
         arret = False
         sm = 0
         k = 1
         arret = False
         sm = 0


@@ -492,86 +657,222 @@ def __evalue_maj_theta__(nbexp,out=False):
             arret = errorq <  error
             k+=1
             variation = "+" if smax > sm else "-"
             arret = errorq <  error
             k+=1
             variation = "+" if smax > sm else "-"

-            sm = smax

             print variation,
             if k%100 ==0 :
                 print "k:",k,"erreur sur q", errorq, "et q:",q
                 print "maxg=", mxg
             print variation,
             if k%100 ==0 :
                 print "k:",k,"erreur sur q", errorq, "et q:",q
                 print "maxg=", mxg

+            if smax - sm  > 500:
+                print "variation trop grande"
+                print "init"
+                print init
+                exit 
+            sm = smax
+
+        if k == itermax:
+            print "nbre d'iteration trop grand"
+            print "init"
+            print init
+            exit 
+
+        print "###############"
+        print k
+        print "###############"
+        m += [k]
+
+    print (min(m),max(m),float(sum(m))/nbexp,m),m
+
+
+
+
+ASYNC = False
+#__une_seule_exp__("config_initiale.py")
+#pour evaluerle test d'arret
+"""
+listederes=[]
+for k in list(np.logspace(-6, -3, num=15)):
+    errorD = k
+    listederes += __evalue_maj_theta__(5)
+    
+print listederes
+"""
+
+
+
+
+def __comparePsh_et_Psh_avec8_3__(nbexp,out=False):
+    print "tttttttttttt"
+    global u, v, la, w, theta , q,  Ps, Rh,  eta, x, valeurFonctionDuale 
+    res = {}
+    m,mp = [],[]
+    itermax = 100000
+
+    def __maj_theta(k):
+        mem = []
+        om = omega/(mt.pow(k,0.75))
+        return om
+    liste_arret=[]
+    for idxexp in range(nbexp):
+        mxg = 0
+        if not(out):
+            initialisation()
+        else :
+            initialisation_()
+            
+        k = 1
+        arret = False
+        sm = 0
+        err, ar = 0,0
+        dur,dura=0,0
+
+        # duppliquee 
+        while k < itermax  and not arret : 
+            (u,v,la,w,theta,eta,q,Ps,Rh,x,valeurFonctionDuale,ar,mxg,smax,ct,cta)=maj(k,__maj_theta,mxg,idxexp)
+            err =  (max(q.values()) - min(q.values()))/min(q.values())
+            dur += ct
+            dura += cta
+            arret = err <  errorq or ar < errorD
+            k+=1
+            variation = "+" if smax > sm else "-"
+            if out : print variation,
+            if k%100 ==0 :
+                print "k:", k, 
+                "erreur sur q", 
+                errorq, "erreur sur F", ar, "et q:", q
+
+                if out : print "maxg=", mxg

                 mem = [deepcopy(q),deepcopy(Ps),deepcopy(Rh),deepcopy(eta),
                        deepcopy(x),deepcopy(u),deepcopy(v),deepcopy(la),deepcopy(w)]
                 mem = [deepcopy(q),deepcopy(Ps),deepcopy(Rh),deepcopy(eta),
                        deepcopy(x),deepcopy(u),deepcopy(v),deepcopy(la),deepcopy(w)]

+            """

             if k%4500 == 0 :
             if k%4500 == 0 :

-                print "#########\n",mem,"\#########\n"
+                
+                #print "#########\n",mem,"\#########\n"

             if k%4600 == 0 :
             if k%4600 == 0 :

-                print "#########\n",mem,"\#########\n"
-
-
+                #print "#########m\n",mem,"\#########\n"
+            """
+                

 
             if smax - sm  > 500:
                 print "variation trop grande"
                 print "init"
                 print init
 
             if smax - sm  > 500:
                 print "variation trop grande"
                 print "init"
                 print init

-                sy.exit(0)
+                exit 
+            sm = smax
+        

         if k == itermax:
             print "nbre d'iteration trop grand"
             print "init"
             print init
             sy.exit(1)
         if k == itermax:
             print "nbre d'iteration trop grand"
             print "init"
             print init
             sy.exit(1)

+        else :
+            liste_arret += [(err, ar,k,errorD)]

 
         print "###############"
 
         print "###############"

-        print k
+        print k," avec optym"

         print "###############"
         m += [k]
         print "###############"
         m += [k]

-
-    print (min(m),max(m),float(sum(m))/nbexp,m),m
+        ##fin de duplication 

 
 
 
 
 
 

+        #remise a zero 
+        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]
+        k = 1
+        arret = False
+        sm = 0
+        err, ar = 0,0
+        dur,dura=0,0
+        # duppliquee 
+        while k < itermax  and not arret : 
+            (u,v,la,w,theta,eta,q,Ps,Rh,x,valeurFonctionDuale,ar,mxg,smax,ct,cta)=maj(k,__maj_theta,mxg,idxexp,True)
+            err =  (max(q.values()) - min(q.values()))/min(q.values())
+            dur += ct
+            dura += cta
+            arret = err <  errorq or ar < errorD
+            k+=1
+            variation = "+" if smax > sm else "-"
+            if out : print variation,
+            if k%100 ==0 :
+                print "k:", k, 
+                "erreur sur q", 
+                errorq, "erreur sur F", ar, "et q:", q

 
 

+                if out : print "maxg=", mxg
+                mem = [deepcopy(q),deepcopy(Ps),deepcopy(Rh),deepcopy(eta),
+                       deepcopy(x),deepcopy(u),deepcopy(v),deepcopy(la),deepcopy(w)]
+            """
+            if k%4500 == 0 :
+                
+                #print "#########\n",mem,"\#########\n"
+            if k%4600 == 0 :
+                #print "#########m\n",mem,"\#########\n"
+            """
+                

 
 

-ASYNC = False
-__evalue_maj_theta__(1,True)
+            if smax - sm  > 500:
+                print "variation trop grande"
+                print "init"
+                print init
+                exit 
+            sm = smax
+        
+        if k == itermax:
+            print "nbre d'iteration trop grand"
+            print "init"
+            print init
+            sy.exit(1)
+        else :
+            liste_arret += [(err, ar,k,errorD)]

 
 

+        print "###############"
+        print k,"sans optym"
+        print "###############"
+        mp += [k]

 
 

+ 
+    #print liste_arret    
+    #print (min(m),max(m),float(sum(m))/nbexp,m),m
+    #return (liste_arret,dur,dura)
+    return (m,mp)

 
 

-#ASYNC = True
-#taux_succes = 0.9
-#__evalue_maj_theta__()

 
 
 
 

+# pour evaluer argmin

 """
 """

-initialisation()
-k = 1
-arret = False
-while k < 10000  and not arret : 
-    (u,v,la,w,theta,eta,q,Ps,Rh,x,valeurFonctionDuale,ar)=maj(k,maj_theta)
-    arret = ar
+listederes=[]
+k=0
+for k in list(np.logspace(-5, -3, num=10)):
+    print "k",k
+    errorD = k
+    listederes += [(k,__evalue_maj_theta__(3))]

     k+=1
     k+=1

-    errorq =  abs(min(q.values())-max(q.values()))
-    print "errorq",errorq
-    arret = errorq <  error
+print listederes

 """
 
 """
 

-"""
-    print "u",u
-    print "w",w
-    print "theta",theta
-    print "eta", eta
-    print "q",q
-    print "v",v
-    print "lambda",la
-    print "Ps",Ps
-    print "Rh",Rh
-    print "x",x
-"""
+# pour evaluer psh
+listederes=[]
+for k in list(np.logspace(-5, -3, num=10)):
+    print "Precision",k
+    errorD = k
+    listederes += [(k,__comparePsh_et_Psh_avec8_3__(10))]
+    k+=1
+print listederes

 
 
 
 

-"""
-    k +=1

 
 

-"""

 
 
 
 
 
 

+print"#######################"
+
+
+#ASYNC = True
+#taux_succes = 0.9
+#__evalue_maj_theta__()
+
+
+
+
+"""

 print "u",u
 print "v",v
 print "lambda",la
 print "u",u
 print "v",v
 print "lambda",la


@@ -584,7 +885,7 @@ print "Rh",Rh
 print "x",x
 print "L",valeurFonctionDuale
 
 print "x",x
 print "L",valeurFonctionDuale
 

-
+"""

 
 # relation ente les variables primaires et secondaires ?
 
 
 # relation ente les variables primaires et secondaires ?