X-Git-Url: https://bilbo.iut-bm.univ-fcomte.fr/and/gitweb/desynchronisation-controle.git/blobdiff_plain/b6e0e34a9afa3b4060edb6098e6fcc397979ffe8..fe438aa764b5d1ad7ea61512cffbb3761eff99b7:/exp_controle_asynchrone/simulMWSN.py?ds=inline diff --git a/exp_controle_asynchrone/simulMWSN.py b/exp_controle_asynchrone/simulMWSN.py index 43083c4..45326b6 100644 --- a/exp_controle_asynchrone/simulMWSN.py +++ b/exp_controle_asynchrone/simulMWSN.py @@ -9,9 +9,10 @@ from copy import deepcopy import sys as sy import cv as cv import cv2 as cv2 +import datetime as dt errorq = 0.1 -errorD = 1E-5 +errorD = 1E-2 epsilon = 1E-10 vrate = 0.8 p = 0.7 @@ -59,11 +60,12 @@ def afficheGraph(G,l,tx,ty,sink): 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)): - return(int(tx*x/coteCarre),ty-int(ty*y/coteCarre)) + 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 + 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]): @@ -113,6 +115,8 @@ N = G.nodes() V = list(set(N)-set([sink])) source = V print "source",source + + #afficheGraph(G,l,500,500,sink) #nx.draw(G) #pb.show() @@ -122,7 +126,6 @@ print "source",source - #print G.edges(data=True) #TODO afficher le graphe et etre sur qu'il est connexe @@ -163,7 +166,6 @@ alpha = 0.5 beta = 1.3E-8 gamma = 55.54 delta = 0.2 -zeta = 0.1 amplifieur = 1 sigma2 = 3500 Bi = 5 @@ -241,19 +243,11 @@ 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 + r= opt.fmin_cobyla(f,xini,cons=[xpos],args=args,consargs=(),rhoend=errorD/100,iprint=0,maxfun=nbiter) 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: - 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 @@ -274,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 @@ -344,39 +338,59 @@ def maj(k,maj_theta,mxg,idxexp): 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 + tq,tqa = 0,0 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) + 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] - - + tp,tpa = 0,0 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(2)/3))) +la[h]*psh - return v[h]* mt.log(float(sigma2)/D)/(gamma*mt.pow(float(2)/3)) +la[h]*psh + 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: - """ - 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 : + 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 = 10 + 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] @@ -390,26 +404,43 @@ def maj(k,maj_theta,mxg,idxexp): 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]) - 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] + + + 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 - 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) + 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)] @@ -432,7 +463,7 @@ def maj(k,maj_theta,mxg,idxexp): 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) @@ -525,7 +556,7 @@ def __evalue_maj_theta__(nbexp,out=False): res = {} m = [] itermax = 100000 - + def __maj_theta(k): mem = [] om = omega/(mt.pow(k,0.75)) @@ -542,15 +573,18 @@ def __evalue_maj_theta__(nbexp,out=False): arret = False sm = 0 err, ar = 0,0 + dur,dura=0,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) + 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%500 ==0 : + if k%10 ==0 : print "k:", k, "erreur sur q", errorq, "erreur sur F", ar, "et q:", q @@ -563,7 +597,7 @@ def __evalue_maj_theta__(nbexp,out=False): #print "#########\n",mem,"\#########\n" if k%4600 == 0 : - #print "#########\n",mem,"\#########\n" + #print "#########m\n",mem,"\#########\n" """ @@ -589,7 +623,8 @@ def __evalue_maj_theta__(nbexp,out=False): #print liste_arret #print (min(m),max(m),float(sum(m))/nbexp,m),m - return liste_arret + #return (liste_arret,dur,dura) + return (liste_arret,dur) def __une_seule_exp__(fichier_donees): @@ -651,12 +686,185 @@ def __une_seule_exp__(fichier_donees): 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)] + """ + 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," avec optym" + print "###############" + m += [k] + ##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" + """ + + + 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) + + + +# pour evaluer argmin +""" +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 +print listederes +""" + +# 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 + + + + + + +print"#######################" + + #ASYNC = True #taux_succes = 0.9 #__evalue_maj_theta__()