From: couchot Date: Thu, 24 Oct 2013 12:39:00 +0000 (+0200) Subject: ajout des fichiers python d'expériementation X-Git-Url: https://bilbo.iut-bm.univ-fcomte.fr/and/gitweb/desynchronisation-controle.git/commitdiff_plain/21fd71af1979f07838bd6c2d85a03d23ed733f93 ajout des fichiers python d'expériementation --- diff --git a/exp_controle_asynchrone/simulMWSN.py b/exp_controle_asynchrone/simulMWSN.py new file mode 100644 index 0000000..fbd48aa --- /dev/null +++ b/exp_controle_asynchrone/simulMWSN.py @@ -0,0 +1,550 @@ +import math as mt +from random import * +import networkx as nx +import numpy as np +import pylab as pb +from itertools import * +from scipy import optimize as opt +from copy import deepcopy +error = 0.1 +epsilon = 1E-10 +vrate = 0.8 +p = 0.7 +coteCarre = 50 +distanceEmmissionMax = 30 +nbiter = 1000 +POS = 1 +POS_NUL = 2 +POSINF1 = 3 +init = [] + + + +# construction du graphe +n = 10 +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)] + +#n=3 +#lg= [(0,1,23),(1,0,15),(1,2,45)] +sink = n-1 + + +def genereGraph(): + test = False + while not test : + G = nx.DiGraph() + l = [(random()*coteCarre, random()*coteCarre) for _ in range(n)] + for io in range(len(l)) : + for ie in range(len(l)) : + if ie != io : + dist = mt.sqrt((l[io][0]-l[ie][0])**2 + (l[io][1]-l[ie][1])**2) + if dist <= distanceEmmissionMax : + 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 + +G = nx.DiGraph() +G.add_weighted_edges_from(lg) +#print nx.is_strongly_connected(G) + + +#nx.draw(G) +#pb.show() + + + + + + +#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())) +d = [di['weight'] for (_,_,di) in G.edges(data=True)] + +#print "L",L +#print "d",d + + +def ail(i,l): + assert l in L, " pb de dimennsion de l: "+str(l)+ " "+ str(L) + r = 0 + (o,e) = G.edges()[l] + if o == i : + r = 1 + elif e == i : + r = -1 + return r + + + +# Constantes +a = [[ail(i,l) for l in L ] for i in xrange(n)] +aplus = [[1 if ail(i,l) > 0 else 0 for l in L ] for i in xrange(n)] +amoins = [[1 if ail(i,l) < 0 else 0 for l in L ] for i in xrange(n)] + + + + + +alpha = 0.5 +beta = 1.3E-8 +gamma = 55.54 +delta = 0.2 +amplifieur = 1 +sigma2 = 3500 +Bi = 5 +omega = 0.15 +D = 100 +path_loss_exp = 4 +cr = 0.5 +cs = [alpha + beta*(di**path_loss_exp) for di in d] +#print "cs",cs + + +#TODO +def etahi(h,i,R): + ret = 0 + if i == h : + ret = R[h] + elif i == sink : + ret = -R[h] + return ret + + +def psi(Ps,i): + if i not in Ps: + return 0 + else : + return Ps[i] + + + +def cmpPair(x1,x2): + return cmp(x1[1],x2[1]) + +def aminp(l): + l.sort(cmp=cmpPair) + return l[0][0] + + + +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): + global u, v, la, w, theta , q, Ps, Rh, eta, x,valeurFonctionDuale + + print "du=",distance(u,up), + print "dv=",distance(v,vp), + print "dlambda=",distance(la,lap), + print "dw=",distance(w,wp), + print "dtheta=",abs(theta-thetap), + print "deta=",distance(eta,etap), + print "dq=",distance(q,qp), + print "dPs=",distance(Ps,Psp), + print "dRh=",distance(Rh,Rhp), + print "dx=",distance(x,xp), + print "dL=",abs(valeurFonctionDuale-valeurFonctionDualep),"\n" + + +valeurFonctionDuale = 0 + +def entre0et1(x): + r = x if (x >0 and x <= 1) else -1 + #print "ds entre0et1 x et r", x,r + return r + + +def xpos(x): + r = x if x >0 else -1 + #print "ds xpos x et r", x,r + return r + +def xposounul(x): + return x if x >= 0 else -1 + +#f_q,q[i],POS,[i] +def armin(f,xini,xr,args): + #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 + 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 + 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 + + + return r + + + + + + +def maj_theta(k): + return omega/(mt.pow(k,0.5)) + + + +def maj_theta(k): + return omega/mt.sqrt(k) + + + + +def maj(k,maj_theta,mxg,idxexp): + # initialisation des operateurs lagrangiens + global u, v, la, w, theta , q, Ps, Rh, eta, x, valeurFonctionDuale + + smax = 0 + #print "iteration",k + + up = {} + for h in V: + for i in N: + if not ASYNC or random() < taux_succes: + s = eta[(h,i)]-sum([a[i][l]*x[(h,l)] for l in L]) + if abs(s) > mxg : + print "ds calcul u",abs(s),idxexp + mxg = abs(s) + smax = max(smax,abs(s)) + up[(h,i)] = u[(h,i)]-theta*s + else : + up[(h,i)] = u[(h,i)] + + + + vp = {} + for h in V: + if not ASYNC or random() < taux_succes: + s = Rh[h]- mt.log(float(sigma2)/D)/(gamma*mt.pow(Ps[h],float(1)/3)) + if abs(s) > mxg : + print "ds calcul v",abs(s),idxexp + mxg = abs(s) + smax = max(smax,abs(s)) + vp[h] = max(0,v[h]-theta*s) + else : + vp[h] = v[h] + + + + lap={} + for i in N: + if not ASYNC or random() < taux_succes: + 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) + if abs(s) > mxg : + print "ds calcul la",abs(s),idxexp,i + mxg = abs(s) + smax = max(smax,abs(s)) + resla = la[i]-theta*s + lap[i] = max(0,resla) + else : + lap[i] = la[i] + + + + + wp={} + for l in L: + if not ASYNC or random() < taux_succes: + s = sum([a[i][l]*q[i] for i in N]) + if abs(s) > mxg : + print "ds calcul w",abs(s),idxexp + mxg = abs(s) + smax = max(smax,abs(s)) + wp[l] = w[l] + theta*s + else : + wp[l] = w[l] + + + + thetap = maj_theta(k) + + + qp={} + def f_q(qi,i): + fa = sum([a[i][l]*w[l] for l in L]) - la[i]*Bi + return qi*qi + qi*fa + + for i in N: + if not ASYNC or random() < taux_succes: + c = -float(sum([a[i][l]*w[l] for l in L]) - la[i]*Bi)/(2*amplifieur) + rep = epsilon if c <= 0 else c + qp[i] = rep + else : + qp[i] = q[i] + + + + + + Psp={} + #print "maj des des Psh" + def f_Ps(psh,h): + #print "ds f_ps",psh, v[h]* mt.log(float(sigma2)/D)/(gamma*((psh**2)**(float(1)/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 = (float(2*v[h]*mt.log(float(sigma2)/D))/mt.pow(3*gamma*lah,float(3)/5)) + Psp[h] = epsilon if rep <= 0 else rep + else : + Psp[h] = Ps[h] + + + + etap={} + + for h in V: + for i in N : + etap[(h,i)] = etahi(h,i,Rh) + + + + Rhp={} + 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) + Rhp[h] = 0 if rep < 0 else rep + else : + Rhp[h] = Rh[h] + + xp={} + 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: + 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) + xp[(h,l)] = 0 if rep < 0 else rep + else : + xp[(h,l)] = x[(h,l)] + + + + + + valeurFonctionDualep = 0 + valeurFonctionDualep += sum([amplifieur*q[i]*q[i] for i in N]) + valeurFonctionDualep += sum([sum([delta*(x[(h,l)]**2) for l in L]) for h in V]) + valeurFonctionDualep += sum([delta*(Rh[h]**2) for h in V]) + 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]) + valeurFonctionDualep += sum([v[h]*(mt.log(float(sigma2)/D)/(gamma*mt.pow(Ps[h],float(2)/3)) - Rh[h]) for h in V]) + 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 ]) + valeurFonctionDualep += sum([w[l]*sum([a[i][l]*q[i] for i in N]) for l in L]) + + + #AfficheVariation(up,vp,lap,wp,thetap,etap,qp,Psp,Rhp,xp,valeurFonctionDualep) + + arret = abs(valeurFonctionDuale-valeurFonctionDualep) < error + + return (up,vp,lap,wp,thetap,etap,qp,Psp,Rhp,xp,valeurFonctionDualep,arret,mxg,smax) + + + + +""" +print "u",u +print "v",v +print "lambda",la +print "w",w +print "theta",theta +print "eta", eta +print "q",q +print "Ps",Ps +print "Rh",Rh +print "x",x + +""" + + + +def initialisation(): + global u, v, la, w, theta , q, Ps, Rh, eta, x,init + + theta = omega + + q = {} + for i in N : + q[i] = 0.15 + random()*0.05 + + + Ps= {} + for h in V : + Ps[h] = 0.2+random()*0.3 + + Rh = {} + for vi in V: + Rh[vi] = 0.1 + random()*0.1 + + eta = {} + for h in V : + for i in N: + eta[(h,i)]= etahi(h,i,Rh) + + x={} + for h in V : + for l in L: + x[(h,l)]=0 + + + + + # initialisation des operateurs lagrangiens + u={} + for h in V : + for i in N: + u[(h,i)]= random() + + v = {} + for h in V : + v[h] = Rh[h] + + la = {} + for i in N: + la[i] = random() + + w={} + for l in L: + w[l] = random() + + init = [deepcopy(q),deepcopy(Ps),deepcopy(Rh),deepcopy(eta), + deepcopy(x),deepcopy(u),deepcopy(v),deepcopy(la),deepcopy(w)] + + + +def __evalue_maj_theta__(): + global u, v, la, w, theta , q, Ps, Rh, eta, x, valeurFonctionDuale + nbexp = 5 + res = {} + m = [] + itermax = 10000 + + def __maj_theta(k): + om = omega/(mt.pow(k,0.75)) + return om + for idxexp in range(nbexp): + mxg = 0 + initialisation() + k = 1 + arret = False + sm = 0 + while k < itermax and not arret : + (u,v,la,w,theta,eta,q,Ps,Rh,x,valeurFonctionDuale,ar,mxg,smax)=maj(k,__maj_theta,mxg,idxexp) + errorq = (max(q.values()) - min(q.values()))/min(q.values()) + 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 + if smax - sm > 500: + print "variation trop grande" + print "init" + print init + exit + 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 +__evalue_maj_theta__() +#ASYNC = True +#taux_succes = 0.9 +#__evalue_maj_theta__() + + +""" +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 + k+=1 + errorq = abs(min(q.values())-max(q.values())) + print "errorq",errorq + arret = errorq < error +""" + +""" + 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 +""" + + +""" + k +=1 + +""" + + + +print "u",u +print "v",v +print "lambda",la +print "w",w +print "theta",theta +print "eta", eta +print "q",q +print "Ps",Ps +print "Rh",Rh +print "x",x +print "L",valeurFonctionDuale + + + +# relation ente les variables primaires et secondaires ? + diff --git a/exp_controle_asynchrone/simulMWSN.py~ b/exp_controle_asynchrone/simulMWSN.py~ new file mode 100644 index 0000000..171c682 --- /dev/null +++ b/exp_controle_asynchrone/simulMWSN.py~ @@ -0,0 +1,528 @@ +import math as mt +from random import * +import networkx as nx +import numpy as np +import pylab as pb +from itertools import * +from scipy import optimize as opt + +error = 0.01 +epsilon = 1E-10 +vrate = 0.8 +p = 0.7 +coteCarre = 50 +distanceEmmissionMax = 30 +nbiter = 1000 +POS = 1 +POS_NUL = 2 +POSINF1 = 3 + + + + +# construction du graphe +n = 10 +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)] + +#n=3 +#lg= [(0,1,23),(1,0,15),(1,2,45)] +sink = n-1 + + +def genereGraph(): + test = False + while not test : + G = nx.DiGraph() + l = [(random()*coteCarre, random()*coteCarre) for _ in range(n)] + for io in range(len(l)) : + for ie in range(len(l)) : + if ie != io : + dist = mt.sqrt((l[io][0]-l[ie][0])**2 + (l[io][1]-l[ie][1])**2) + if dist <= distanceEmmissionMax : + 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 + +G = nx.DiGraph() +G.add_weighted_edges_from(lg) +#print nx.is_strongly_connected(G) + + +#nx.draw(G) +#pb.show() + + + + + + +#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())) +d = [di['weight'] for (_,_,di) in G.edges(data=True)] + +#print "L",L +#print "d",d + + +def ail(i,l): + assert l in L, " pb de dimennsion de l: "+str(l)+ " "+ str(L) + r = 0 + (o,e) = G.edges()[l] + if o == i : + r = 1 + elif e == i : + r = -1 + return r + + + +# Constantes +a = [[ail(i,l) for l in L ] for i in xrange(n)] +aplus = [[1 if ail(i,l) > 0 else 0 for l in L ] for i in xrange(n)] +amoins = [[1 if ail(i,l) < 0 else 0 for l in L ] for i in xrange(n)] + + + + + +alpha = 0.5 +beta = 1.3E-8 +gamma = 55.54 +delta = 0.2 +amplifieur = 1 +sigma2 = 3500 +Bi = 5 +omega = 0.15 +D = 100 +path_loss_exp = 4 +cr = 0.5 +cs = [alpha + beta*(di**path_loss_exp) for di in d] +#print "cs",cs + + +#TODO +def etahi(h,i,R): + ret = 0 + if i == h : + ret = R[h] + elif i == sink : + ret = -R[h] + return ret + + +def psi(Ps,i): + if i not in Ps: + return 0 + else : + return Ps[i] + + + +def cmpPair(x1,x2): + return cmp(x1[1],x2[1]) + +def aminp(l): + l.sort(cmp=cmpPair) + return l[0][0] + + + +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): + global u, v, la, w, theta , q, Ps, Rh, eta, x,valeurFonctionDuale + + print "du=",distance(u,up), + print "dv=",distance(v,vp), + print "dlambda=",distance(la,lap), + print "dw=",distance(w,wp), + print "dtheta=",abs(theta-thetap), + print "deta=",distance(eta,etap), + print "dq=",distance(q,qp), + print "dPs=",distance(Ps,Psp), + print "dRh=",distance(Rh,Rhp), + print "dx=",distance(x,xp), + print "dL=",abs(valeurFonctionDuale-valeurFonctionDualep),"\n" + + +valeurFonctionDuale = 0 + +def entre0et1(x): + r = x if (x >0 and x <= 1) else -1 + #print "ds entre0et1 x et r", x,r + return r + + +def xpos(x): + r = x if x >0 else -1 + #print "ds xpos x et r", x,r + return r + +def xposounul(x): + return x if x >= 0 else -1 + +#f_q,q[i],POS,[i] +def armin(f,xini,xr,args): + #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 + 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 + 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 + + + return r + + + + + + +def maj_theta(k): + return omega/(mt.pow(k,0.5)) + + + +def maj_theta(k): + return omega/mt.sqrt(k) + + + + +def maj(k,maj_theta): + # initialisation des operateurs lagrangiens + global u, v, la, w, theta , q, Ps, Rh, eta, x, valeurFonctionDuale + + + #print "iteration",k + + up = {} + for h in V: + for i in N: + if not ASYNC or random() < taux_perte: + up[(h,i)] = u[(h,i)]-theta*(eta[(h,i)]-sum([a[i][l]*x[(h,l)] for l in L])) + else : + up[(h,i)] = u[(h,i)] + + + + vp = {} + for h in V: + #print "vp",gamma*mt.pow(Ps[h],float(1)/3) + #vp[h] = max(0,v[h]-theta*(Rh[h]- mt.log(float(sigma2)/D)/(gamma*mt.pow(Ps[h],float(1)/3)))) + if not ASYNC or random() < taux_perte: + vp[h] = max(0,v[h]-theta*(Rh[h]- mt.log(float(sigma2)/D)/(gamma*mt.pow(Ps[h],float(1)/3)))) + else : + vp[h] = v[h] + + + + lap={} + for i in N: + if not ASYNC or random() < taux_perte: + resla = la[i]-theta*((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))) + lap[i] = max(0,resla) + else : + lap[i] = la[i] + + + + + wp={} + for l in L: + if not ASYNC or random() < taux_perte: + wp[l] = w[l] + theta*(sum([a[i][l]*q[i] for i in N])) + else : + wp[l] = w[l] + + + + thetap = maj_theta(k) + + #print "k",k,"maj theta", thetap + + + qp={} + def f_q(qi,i): + fa = sum([a[i][l]*w[l] for l in L]) - la[i]*Bi + return qi*qi + qi*fa + + for i in N: + if not ASYNC or random() < taux_perte: + c = -float(sum([a[i][l]*w[l] for l in L]) - la[i]*Bi)/(2*amplifieur) + rep = epsilon if c <= 0 else c + qp[i] = rep + else : + qp[i] = q[i] + + + + + + Psp={} + #print "maj des des Psh" + def f_Ps(psh,h): + #print "ds f_ps",psh, v[h]* mt.log(float(sigma2)/D)/(gamma*((psh**2)**(float(1)/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_perte: + lah = 0.05 if la[h] == 0 else la[h] + rep = (float(2*v[h]*mt.log(float(sigma2)/D))/mt.pow(3*gamma*lah,float(3)/5)) + Psp[h] = epsilon if rep <= 0 else rep + else : + Psp[h] = Ps[h] + + + + etap={} + + for h in V: + for i in N : + etap[(h,i)] = etahi(h,i,Rh) + + + + Rhp={} + 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_perte: + rep = float(v[h])/(2*delta) + Rhp[h] = 0 if rep < 0 else rep + else : + Rhp[h] = Rh[h] + + xp={} + 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_perte: + 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) + xp[(h,l)] = 0 if rep < 0 else rep + else : + xp[(h,l)] = x[(h,l)] + + + + + + valeurFonctionDualep = 0 + valeurFonctionDualep += sum([amplifieur*q[i]*q[i] for i in N]) + valeurFonctionDualep += sum([sum([delta*(x[(h,l)]**2) for l in L]) for h in V]) + valeurFonctionDualep += sum([delta*(Rh[h]**2) for h in V]) + 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]) + valeurFonctionDualep += sum([v[h]*(mt.log(float(sigma2)/D)/(gamma*mt.pow(Ps[h],float(2)/3)) - Rh[h]) for h in V]) + 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 ]) + valeurFonctionDualep += sum([w[l]*sum([a[i][l]*q[i] for i in N]) for l in L]) + + + #AfficheVariation(up,vp,lap,wp,thetap,etap,qp,Psp,Rhp,xp,valeurFonctionDualep) + + arret = abs(valeurFonctionDuale-valeurFonctionDualep) < error + + return (up,vp,lap,wp,thetap,etap,qp,Psp,Rhp,xp,valeurFonctionDualep,arret) + + + + +""" +print "u",u +print "v",v +print "lambda",la +print "w",w +print "theta",theta +print "eta", eta +print "q",q +print "Ps",Ps +print "Rh",Rh +print "x",x + +""" + + + +def initialisation(): + global u, v, la, w, theta , q, Ps, Rh, eta, x, valeurFonctionDuale + + theta = omega + + + # TODO Ahmed : initialiser avec qi >0 et somme (ail.qi) = 0 Ei/i = Ti + #q = [0 for i in N] + q = {} + for i in N : + q[i] = 0.15 + random()*0.05 + + + + + + #TODO Ahmed doit donner les valeurs initiales pour Ps + Ps= {} + for h in V : + Ps[h] = 0.2+random()*0.3 + + + + + #TODO Ahmed doit donner les valeurs initiales pour Rh + Rh = {} + for vi in V: + Rh[vi] = 0.1 + random()*0.1 + + + + eta = {} + for h in V : + for i in N: + eta[(h,i)]= etahi(h,i,Rh) + + + + # Ok pour 0 + x={} + for h in V : + for l in L: + x[(h,l)]=0#random() + + + + + # initialisation des operateurs lagrangiens + u={} + for h in V : + for i in N: + u[(h,i)]= random() + + v = {} + for h in V : + v[h] = Rh[h] + + la = {} + for i in N: + la[i] = random() + + w={} + for l in L: + w[l] = random() + + + +def __evalue_maj_theta__(): + global u, v, la, w, theta , q, Ps, Rh, eta, x, valeurFonctionDuale + nbexp = 50 + res = {} + m = [] + def __maj_theta(k): + om = omega/(mt.pow(k,0.75)) + return om + for idxexp in range(nbexp): + + initialisation() + k = 1 + arret = False + while k < 1000000 and not arret : + (u,v,la,w,theta,eta,q,Ps,Rh,x,valeurFonctionDuale,ar)=maj(k,__maj_theta) + errorq = abs(min(q.values())-max(q.values())) + arret = errorq < error + k+=1 + if k%5000 ==0 : + print "k:",k,"erreur sur q", errorq, "et q:",q + print "###############" + print k + print "###############" + m += [k] + print (min(m),max(m),float(sum(m))/nbexp,m),m + + + + + +#ASYNC = False +#__evalue_maj_theta() +ASYNC = True +taux_perte = 0.9 +__evalue_maj_theta__() + + +""" +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 + k+=1 + errorq = abs(min(q.values())-max(q.values())) + print "errorq",errorq + arret = errorq < error +""" + +""" + 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 +""" + + +""" + k +=1 + +""" + + + +print "u",u +print "v",v +print "lambda",la +print "w",w +print "theta",theta +print "eta", eta +print "q",q +print "Ps",Ps +print "Rh",Rh +print "x",x +print "L",valeurFonctionDuale + + + +# relation ente les variables primaires et secondaires ? + diff --git a/exp_controle_asynchrone/simulMWSNASYNC.py b/exp_controle_asynchrone/simulMWSNASYNC.py new file mode 100644 index 0000000..3376a8e --- /dev/null +++ b/exp_controle_asynchrone/simulMWSNASYNC.py @@ -0,0 +1,509 @@ +import math as mt +from random import * +import networkx as nx +import numpy as np +import pylab as pb +from itertools import * +from scipy import optimize as opt + +error = 1E-5 +epsilon = 1E-10 +vrate = 0.8 +p = 0.7 +coteCarre = 50 +distanceEmmissionMax = 30 +nbiter = 1000 +POS = 1 +POS_NUL = 2 +POSINF1 = 3 + +ASYNC = False +taux_perte = 0.99 + +# construction du graphe +n = 10 +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)] + +#n=3 +#lg= [(0,1,23),(1,0,15),(1,2,45)] +sink = n-1 + + +def genereGraph(): + test = False + while not test : + G = nx.DiGraph() + l = [(random()*coteCarre, random()*coteCarre) for _ in range(n)] + for io in range(len(l)) : + for ie in range(len(l)) : + if ie != io : + dist = mt.sqrt((l[io][0]-l[ie][0])**2 + (l[io][1]-l[ie][1])**2) + if dist <= distanceEmmissionMax : + 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 + +G = nx.DiGraph() +G.add_weighted_edges_from(lg) +#print nx.is_strongly_connected(G) + + +#nx.draw(G) +#pb.show() + + + + + + +#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())) +d = [di['weight'] for (_,_,di) in G.edges(data=True)] + + + +def ail(i,l): + assert l in L, " pb de dimennsion de l: "+str(l)+ " "+ str(L) + r = 0 + (o,e) = G.edges()[l] + if o == i : + r = 1 + elif e == i : + r = -1 + return r + + + +# Constantes +a = [[ail(i,l) for l in L ] for i in xrange(n)] +aplus = [[1 if ail(i,l) > 0 else 0 for l in L ] for i in xrange(n)] +amoins = [[1 if ail(i,l) < 0 else 0 for l in L ] for i in xrange(n)] + + + + + +alpha = 0.5 +beta = 1.3E-8 +gamma = 55.54 +delta = 0.2 +amplifieur = 1 +sigma2 = 3500 +Bi = 5 +omega = 0.15 +D = 100 +path_loss_exp = 4 +cr = 0.5 +cs = [alpha + beta*(di**path_loss_exp) for di in d] +#print "cs",cs + + +#TODO +def etahi(h,i,R): + ret = 0 + if i == h : + ret = R[h] + elif i == sink : + ret = -R[h] + return ret + + +def psi(Ps,i): + if i not in Ps: + return 0 + else : + return Ps[i] + +# separer le calcul des operateurs lagrangien pour avoir un pas d'iteration de decallage + + +def cmpPair(x1,x2): + return cmp(x1[1],x2[1]) + +def aminp(l): + l.sort(cmp=cmpPair) + return l[0][0] + + + +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): + global u, v, la, w, theta , q, Ps, Rh, eta, x,valeurFonctionDuale + + print "du=",distance(u,up), + print "dv=",distance(v,vp), + print "dlambda=",distance(la,lap), + print "dw=",distance(w,wp), + print "dtheta=",abs(theta-thetap), + print "deta=",distance(eta,etap), + print "dq=",distance(q,qp), + print "dPs=",distance(Ps,Psp), + print "dRh=",distance(Rh,Rhp), + print "dx=",distance(x,xp), + print "dL=",abs(valeurFonctionDuale-valeurFonctionDualep),"\n" + + +valeurFonctionDuale = 0 + +def entre0et1(x): + r = x if (x >0 and x <= 1) else -1 + #print "ds entre0et1 x et r", x,r + return r + + +def xpos(x): + r = x if x >0 else -1 + #print "ds xpos x et r", x,r + return r + +def xposounul(x): + return x if x >= 0 else -1 + +#f_q,q[i],POS,[i] +def armin(f,xini,xr,args): + #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 + 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 + 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 + + + return r + + + + + + +def maj_theta(k): + return omega/(mt.pow(k,0.5)) + + + +def maj_theta(k): + return omega/mt.sqrt(k) + + + + + +def maj(k,maj_theta): + # initialisation des operateurs lagrangiens + global u, v, la, w, theta , q, Ps, Rh, eta, x, valeurFonctionDuale + + + #print "iteration",k + + up = {} + for h in V: + for i in N: + if not ASYNC or random() < taux_perte: + up[(h,i)] = u[(h,i)]-theta*(eta[(h,i)]-sum([a[i][l]*x[(h,l)] for l in L])) + else : + up[(h,i)] = u[(h,i)] + + + + + vp = {} + for h in V: + #print "vp",gamma*mt.pow(Ps[h],float(1)/3) + if not ASYNC or random() < taux_perte: + vp[h] = max(0,v[h]-theta*(Rh[h]- mt.log(float(sigma2)/D)/(gamma*mt.pow(Ps[h],float(1)/3)))) + else : + vp[h] = v[h] + + + lap={} + for i in N: + if not ASYNC or random() < taux_perte: + resla = la[i]-theta*((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))) + lap[i] = max(0,resla) + else : + lap[i] = la[i] + + wp={} + for l in L: + if not ASYNC or random() < taux_perte: + wp[l] = w[l] + theta*(sum([a[i][l]*q[i] for i in N])) + else : + wp[l] = w[l] + + thetap = maj_theta(k) + + + qp={} + def f_q(qi,i): + fa = sum([a[i][l]*w[l] for l in L]) - la[i]*Bi + return qi*qi + qi*fa + + for i in N: + if not ASYNC or random() < taux_perte: + c = -float(sum([a[i][l]*w[l] for l in L]) - la[i]*Bi)/(2*amplifieur) + rep = epsilon if c <= 0 else c + qp[i] = rep + else : + qp[i] = q[i] + + + + + + + + Psp={} + def f_Ps(psh,h): + 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_perte: + lah = 0.05 if la[h] == 0 else la[h] + rep = (float(2*v[h]*mt.log(float(sigma2)/D))/mt.pow(3*gamma*lah,float(3)/5)) + Psp[h] = epsilon if rep <= 0 else rep + else : + Psp[h] = Ps[h] + + + + + etap={} + for h in V: + for i in N : + etap[(h,i)] = etahi(h,i,Rh) + + + + Rhp={} + 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_perte: + rep = float(v[h])/(2*delta) + Rhp[h] = 0 if rep < 0 else rep + else : + Rhp[h] = Rh[h] + + xp={} + 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_perte: + 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) + xp[(h,l)] = 0 if rep < 0 else rep + else : + xp[(h,l)] = x[(h,l)] + + valeurFonctionDualep = 0 + valeurFonctionDualep += sum([amplifieur*q[i]*q[i] for i in N]) + valeurFonctionDualep += sum([sum([delta*mt.pow(x[(h,l)],2) for l in L]) for h in V]) + valeurFonctionDualep += sum([delta*(Rh[h]**2) for h in V]) + 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]) + valeurFonctionDualep += sum([v[h]*(mt.log(float(sigma2)/D)/(gamma*mt.pow(Ps[h],float(2)/3)) - Rh[h]) for h in V]) + 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 ]) + valeurFonctionDualep += sum([w[l]*sum([a[i][l]*q[i] for i in N]) for l in L]) + + + AfficheVariation(up,vp,lap,wp,thetap,etap,qp,Psp,Rhp,xp,valeurFonctionDualep) + + arret = abs(valeurFonctionDuale-valeurFonctionDualep) < error + + return (up,vp,lap,wp,thetap,etap,qp,Psp,Rhp,xp,valeurFonctionDualep,arret) + + + + +""" +print "u",u +print "v",v +print "lambda",la +print "w",w +print "theta",theta +print "eta", eta +print "q",q +print "Ps",Ps +print "Rh",Rh +print "x",x + +""" + + + +def initialisation(): + global u, v, la, w, theta , q, Ps, Rh, eta, x, valeurFonctionDuale + + theta = omega + + + # TODO Ahmed : initialiser avec qi >0 et somme (ail.qi) = 0 Ei/i = Ti + #q = [0 for i in N] + q = {} + for i in N : + q[i] = 0.15 + random()*0.05 + + + + + + #TODO Ahmed doit donner les valeurs initiales pour Ps + Ps= {} + for h in V : + Ps[h] = 0.2+random()*0.3 + + + + + #TODO Ahmed doit donner les valeurs initiales pour Rh + Rh = {} + for vi in V: + Rh[vi] = 0.1 + random()*0.1 + + + + eta = {} + for h in V : + for i in N: + eta[(h,i)]= etahi(h,i,Rh) + + + + # Ok pour 0 + x={} + for h in V : + for l in L: + x[(h,l)]=0#random() + + + + + # initialisation des operateurs lagrangiens + u={} + for h in V : + for i in N: + u[(h,i)]= random() + + v = {} + for h in V : + v[h] = Rh[h] + + la = {} + for i in N: + la[i] = random() + + w={} + for l in L: + w[l] = random() + + + + + +def __evalue_maj_theta(): + global u, v, la, w, theta , q, Ps, Rh, eta, x, valeurFonctionDuale + nbexp = 1 + res = {} + for l in range(1,11): + m = 0 + def __maj_theta(k): + om = omega/(mt.pow(k,float(l)/10)) + #om = omega/(mt.pow(k,0.5)) + #print "ici",k,om + return om + for _ in range(nbexp): + 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 + k+=1 + m += k + res[l]= m/nbexp + print l, res[l] + print res + + +__evalue_maj_theta() +#ASYNC = True +#__evalue_maj_theta() + + +""" +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 + k+=1 +""" + + +""" + 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 +""" + + +""" + k +=1 + +""" + +""" +print "nb it",k +print "u",u +print "v",v +print "lambda",la +print "w",w +print "theta",theta +print "eta", eta +print "q",q +print "Ps",Ps +print "Rh",Rh +print "x",x +print "L",valeurFonctionDuale +""" + + +# relation ente les variables primaires et secondaires ? +