From 3f7dd12e8d881d7c056ea6bbe05dc6d760ef4293 Mon Sep 17 00:00:00 2001 From: Alfred Galichon Date: Thu, 11 Jan 2024 18:39:33 +0100 Subject: [PATCH] update --- mec/__pycache__/nf.cpython-39.pyc | Bin 12633 -> 15314 bytes mec/nf.py | 116 ++++++++++++++++++++++++++---- 2 files changed, 104 insertions(+), 12 deletions(-) diff --git a/mec/__pycache__/nf.cpython-39.pyc b/mec/__pycache__/nf.cpython-39.pyc index a3268145113a9b300bcbb630e3d9ee9ec9f95658..83ccfe0350ed5e41db8547407c9a0cce794ceeb2 100644 GIT binary patch delta 5278 zcma)AeQ;b=6@T~b+qYlIX0vHG*?e`Iwn>+?p_Z?LP%Ix6Y6(SJD7CO`a$lNF)6Mqo zOIyNQ7n_zQNSxxW!#Iv*;|CfYkxHQqD1(Y1DuSQ_tHKC4j?OqE{-YE|{GIzY-6V)( zGVhmj&%Ni|d(S!d+;j74-;4XZ)3voh4L>h*mMpnm)83@U+#e4bAH$b-?&;pIkLxGy z@2g=mGEn<<;CE}6ufGbNoAHS#K>%<%+B>?Z&k~7|fkbwMTlvH#iPH@qO=Qn-)1sf4 zW!58ewTtC5;bwVp_+;o>(r_KYR{3uDR!p$7u01lt%HI57I%mx=(SlBy@+_A3)oo{% zyd`pV8@dCeuQ$EN%8NRZju6xXhA?82@hKAOjWo%pB472=aVm9MhFXTW*eB z%%E_j>!X5G@P)WUP4yXI9;e6SXON}j%Jt_~!uQ892 z#24i&@%E4+zu>l6b1XDS6x$V;XiM)qEZS8ZQ$^o^^RBps|S<)lEOaqjtQlfk~172RP zYWYk19n{-Lu)rwf=C8?fTF*c0)1bt6&l6Q0T-7&*j&IA$+HW}f5m4Blk|7o1 zLuGD9?*sCg_Mx>(#_w7-Niu3t2Z*|jfV7J&0j&bsx5@q`E7O0<7ANI%%i3k8qn7mya%4_Tf2q5eVe0*;R+1r>xG zP0z$zn5z4{(3BS~t>kljj(?hZv>`v_Yp&P0cD!vFzn{ET)$$? zN(8(*ICSQ!V^Zl$PS{`hU7@KFK8|>;aJ^@^ zekZ^UMB=CI0QcHKNI+ks#?3~}4i?cw^rPu(gc1a4s8oZmmiW^wMC}T_OUD5S+F`rS zl^JJ+&=GBbaX$|n@`LB4WfW=@FF21^dxT+Jh$8j$2Cy(m@aBJ z3@ws-Lk&eN0c!ODo!3n1rNz_QG+7`9UX0^>9M+xi0L~|Pjdl~bjYDx@{Ww{u*~ZwR zX0pqW&W$E*w|QPLos{`r1nE4mO69=cBWHl!@bZ%g_?)S%%|XJEi0J z`!G1+U@@ts6I!WbTEh{g;~DNnEwjE$Y=raNTc>&%(9QiX;zt2vE3dq6qnrx2%3ZY| zKVMw^A>-ef9j|ICp-QZ4h}R34Q#OTHiED{{9f7*&s`JvN_CE3F&fCA_e~#CftgG8x z6cp`q3J6k#oaCw!H&W+KfaEfVK*hf?g0p%=lh~O~6odib+CJ4yHDC1O$Q-TrVQpR~3U;cc6rkhAsL1yAN#$rWe zt>q{p!)y_8l(#I-ogd4+a!>rPLB6?9Y)2C4|u5kuYqJ*INzUZf*r1HE_5 zSdm|mM=?U467v`D9s^}L6jZWpRlG+nSFQAQjPItE8Xz>5U}Npo_3i0V+m?G*9q>Pb zZrroGql|r(6bsIY_uU@?p*O*v+WUCT0-4(q_WBtlu!pe#_ z&=j;RL5E|u*!vSP-Xi!Nfm&&Gxg@hrUej}4k=9Egc{#rPb^vUK@+PkB)RnDpc~vSd zvL~!LA@KEY<5t7pL0ZexxsYa>cKA5yGCUm--0%h~!Y<)G3 z!Z^4w#1OusKGr(dQ46K|$u_e;0T5?wbqtL;w(=AVv|z2c8+``sq3=V6!qw=}wSu=` zX3v^8f^4Z5G31ku$1y5H{IC#Fc?*8kGAq-7i;BWmKs{UV z?Y(6-jFTTQ^jtgV(s|o*EnVzExG7Vs2n?jN{Z`%{5qUS-$?7V$Dy3mL2-8=}*>Yek zH8RY_$cSBa12MDKeqxr*F*v~LSI^?IcPwiiuorqMT2yP9=cud5QB7oX7S?0PviNkr zPBmp1k9i#Z`A3U3Gf%-hpV6D#hmZW0-vrgjp?FJ7A~unBwIVXFu4=J4nU@|KCAe)G z4~w`~^5W@X*d{j|L)HF9{!Y9NJgCUFf{Qv$y$24u$rS+?h^sj*jS@kqAdClz8 z1SQa-|7S4>FW&}LBH;}K0VqI$7ywm^*Gzz3Ggx=HR7d;@9igS-K!+>uo4Q_EY+5W| z>-of$>T2EsGYh)iW}i3ewyh=}sz`Ri3vO}KHTUS@{}LiCH~3nL;wPl#X@X}6eoF8x z!OsZj7>b`0{DR;)0&;P(Ui_J8RivmYLB1_UNY5tiE0k|9czH$>quj+YE^pz!kqDc#&X| zfG$))ca>`?b#0~|s0m+$;a{vGU delta 2904 zcmaJ@eQaA-6@T}A*w0R!IC1L4jpI0JlI|^Snx!l0lBQy#?FQ{ksuF=ZTF1_PiJiE% zcb~m%-V4>%sVOThY^Njqfz$+zN$tp3)ge^a_*fYm>IPyUq|}xd>3&w&@+559v{FZnhAUem0+TLENzQ;VLF~KaN8NHG>yjdH}+Cu${BGTx>;PPD`k?M0S)C+E!IkDNR)BY#C+a3>LmnN0(a=Pr7F6R* zXU#0g!XKKR%jZKERcY!WNQ8H>lW;m*j-^Oxj9{FgkswIWgs=#))}_J!D{`7-DLPd> zoh#dxD$c-PS|`!^zP8SmTWr!h&O)$p1DtJ}m9BgpK59F%^&}aQSD!=5QmdX^DsAVS zs#Z)^*7-#GPs3v5e$N@wcp9pa*cM5il{_i4s&_g&<76wg6|QQPY*9%6?;?RT^lY^5 zK1O5#p{jF{nirkKCy5UcpUm5(`GS~xIBD}4+mUgu8rhP4%n=b1wc+f>JI7?8mx$j< zP$qr9oiEuUo1aOhs&(8hrbJ0H;yj5wNf07fgikhxBi|%aT9CG6)W3s!Wf*AxTkH`! zeSqMWAW~wUgB{UN-9C?0RhvtuFQ~P=50kizqtuZWY42HhCi;UGdCUJM_#zzb*d3#{ zi`5&tFAyblA42K^{HP-xoTEKyK1b9lT#AJZ2gz0VRV)gVv60sABWGzsN|K^jCs{)N zMUVV}6mK&qr7y$B@rh<>L6td{7V!Oc zS08&7zSecT;UzNZDHU>Ar-TQ;*%X7h?f`oY_I4Nhq^2h3rA)cVQ&=Arx(*+PU*ZD6 zMR>Qno4o-acRyOYN|!PWxtCKD_P+H{%+WQvnq4Wu~ft1%Zi9*|NqU|2o*6Zyy<F)EPGu%U7@>thZL?H zQY?Q$F!ZtlbG?JdaZD^BBU`5un#iXMB>-{uV@!x2tG0V8A zT+2xV432FM00`{31sRior3yHyb>n7f?)Y-5tZ96jQ53;(yz-TaMNldu@a~) zs$Luh0yFu0WSlJs3Nvxgni8X|f$6Hr{F29gN&E3U%39QVxYoCWodC9Z{NBe=vF6Y) zQL<-b`OC_@#c`^#+Ox=S@e09D2^y&t;rY!$b_cw?d3eum+_aU$j4Ov#-pF-sJi!iY zJjlJ=L`vsP+{gV$8N8X-@p_~@Jj4UM0Vyw)hSdU}zpZz+4h_}g`V1~-7eYx{W+m0d zg2Hm~m$YTw#pJuH>vwh6aIr94)AgdqzNg7mID+-nXreDZLSZiI`_i?l+_7Xe} zySMi2sH z9_2{%%i4mrpf4CDlNqkp2Hu7RAHXlMonV~c zMmfn}+W#o<{=g^fb$EYUci@+#pos<&dmXI7jyTPVs;`u>lQSY!%&b#>3r{{&V0v(T ze?{+gdL4@;EWe17ZU@0V1hT;1AZn68#?l~4_K3W1ZMH~B7v8!fV!cTML4wy1{*B8a zj*2~|der8+AQp9;RXXALU?+49&6zUEJ~0XBhPon0QBo>%FMdn<4uE`%Fk_y9n+jqd0 zfxrlLWC@2HrQ;M^W-H29m6E!wEh>v_TIDL&o-iq`wWMRc>xrt7v2(csJl3`WoSplJ zW$hiH3l1UR&#RmYl!bMnZ$dVPGfn>esSNEsxDU>EwZq-smiur-$;keJG$I831hTGg z{?NWddosY(-R)?~q3Wano2(LBGx!C|sEy8qpi?4*%b93*mRM6TlTKXEzqKN%;M$KUGr`3-Ru{yuUadk*d#?PD2uV00(D3g<@0 H*w+67?8Aj! diff --git a/mec/nf.py b/mec/nf.py index d172d86..5a2084d 100644 --- a/mec/nf.py +++ b/mec/nf.py @@ -209,6 +209,7 @@ def determine_departing_arc(self,entering_a, basis=None): def iterate(self, draw = False, verbose=0): entering_as = self.determine_entering_arcs() + print('entering = ',entering_as) if not entering_as: if verbose>0: print('Optimal solution found.\n=======================') @@ -219,6 +220,7 @@ def iterate(self, draw = False, verbose=0): else: entering_a=entering_as[0] departing_a = self.determine_departing_arc(entering_a) + print('entering_a=', entering_a,'departing_a=', departing_a) if departing_a is None: if verbose>0: print('Unbounded solution.') @@ -232,18 +234,50 @@ def iterate(self, draw = False, verbose=0): self.draw(p_z = p_z,mu_a=mu_a, gain_a = g_a, entering_a = entering_a, departing_a = departing_a) self.tableau_update(entering_a,departing_a) + self.tableau return(2) -class EQF_problem(Network_problem): +class EQF_problem: def __init__(self, nodesList, arcsList, galois_xy, q_z, active_basis=None, zero_node=0, pos=None, seed=777, verbose=0): - c_a = np.zeros(len(arcsList)) - Network_problem.__init__(self, nodesList, arcsList, c_a, q_z, active_basis, zero_node, pos, seed, verbose) + + self.zero_node = zero_node + self.nbz = len(nodesList) + self.nba = len(arcsList) + self.nodesList = nodesList + self.arcsList = arcsList + self.nodesDict = {node:node_ind for (node_ind,node) in enumerate(self.nodesList)} + self.arcsDict = {arc:arc_ind for (arc_ind,arc) in enumerate(self.arcsList)} + if verbose>1: + print('Number of nodes = '+str(self.nbz)+'; number of arcs = '+str(self.nba)+'.') + + data = np.concatenate([-np.ones(self.nba),np.ones(self.nba)]) + arcsIndices = list(range(self.nba)) + arcsOrigins = [self.nodesDict[o] for o,d in self.arcsList] + arcsDestinations = [self.nodesDict[d] for o,d in self.arcsList] + + znotzero = [i for i in range(self.nbz) if i != zero_node] + self.q_z = q_z + self.q0_z = q_z[znotzero] + + self.nabla_a_z = np.array(sp.csr_matrix((data, (arcsIndices+arcsIndices, arcsOrigins+arcsDestinations)), shape = (self.nba,self.nbz)).todense()) + self.nabla0_a_z = self.nabla_a_z[:,znotzero] + + self.basis = two_phase(self.nabla0_a_z.T,self.q0_z) + assert len(self.basis) == (self.nbz - 1) + + self.digraph = nx.DiGraph() + self.digraph.add_edges_from(arcsList) + if pos is None: + pos = nx.spring_layout(self.digraph, seed=seed) + self.pos = pos + + #arcsNames = [str(x)+str(y) for (x,y) in self.arcsList] self.galois_xy = galois_xy self.create_pricing_tree() - def create_pricing_tree(self, basis =None, verbose = False): + def create_pricing_tree(self, verbose = False): the_graph = nx.DiGraph() - the_graph.add_edges_from([self.arcsList[a] for a in self.basis(basis)]) + the_graph.add_edges_from([self.arcsList[a] for a in self.basis]) self.tree = {} @@ -260,10 +294,12 @@ def create_anytree(node, parent=None): create_anytree(self.nodesList[self.zero_node]) if verbose: - for pre, fill, node in RenderTree(self.tree[self.nodesList[self.zero_node]]): - print("%s%s" % (pre, node.name)) - + self.print_pricing_tree() + def print_pricing_tree(self): + for pre, fill, node in RenderTree(self.tree[self.nodesList[self.zero_node]]): + print("%s%s" % (pre, node.name)) + def psol_z(self, current_price=0): nodename = self.nodesList[self.zero_node] self.set_prices_r(nodename,current_price) @@ -278,9 +314,65 @@ def set_prices_r(self,nodename, current_price=0): for child in self.tree[nodename].children: self.set_prices_r(child.name,self.galois_xy[(child.name,nodename)](current_price)) - def cost_improvement_a(self, basis = None): - print('hello') - p_z = self.psol_z() - return np.array([self.galois_xy[(x,y)] (p_z[self.nodesDict[y]]) - p_z[self.nodesDict[x]] for (x,y) in self.arcsList]) + def cut_pricing_tree(self,a_exiting): # returns root of second connected component + x,y = self.arcsList[a_exiting] + print (x,y) + if (self.tree[y].parent == self.tree[x]): + return y + elif (self.tree[x].parent == self.tree[y]): + return x + else: + print('Error in pricing tree during cut phase.') + + def paste_pricing_tree(self,a_entering,z_oldroot): + x,y = self.arcsList[a_entering] + + if (self.tree[z_oldroot] in self.tree[y].ancestors): + z_newroot , z_prec = y,x + elif (self.tree[z_oldroot] in self.tree[x].ancestors): + z_newroot , z_prec = x,y + else: + print('Error in pricing tree during paste phase.') + z = z_newroot + while (z_prec != z_oldroot): + znext = self.tree[z].parent.name + self.tree[z].parent = self.tree[z_prec] + z_prec = z + z = znext + + def iterate(self, draw = False, verbose=0): + + p_z = self.psol_z() + cost_improvement_a = np.array([self.galois_xy[(x,y)] (p_z[self.nodesDict[y]]) - p_z[self.nodesDict[x]] for (x,y) in self.arcsList]) + entering_as = np.where(cost_improvement_a )[0].tolist() + print('entering = ',entering_as) + if not entering_as: + if verbose>0: + print('Optimal solution found.\n=======================') + if draw: + mu_a,p_z,g_a = self.musol_a(),self.p0sol_z(),self.gain_a() + self.draw(p_z = p_z,mu_a=mu_a) + return(0) + else: + entering_a=entering_as[0] + departing_a = self.determine_departing_arc(entering_a) + print('entering_a=', entering_a,'departing_a=', departing_a) + if departing_a is None: + if verbose>0: + print('Unbounded solution.') + return(1) + else: + if verbose>1: + print('entering=',entering_a) + print('departing=',departing_a) + if draw: + mu_a,p_z,g_a = self.musol_a(),self.p0sol_z(),self.gain_a() + self.draw(p_z = p_z,mu_a=mu_a, gain_a = g_a, entering_a = entering_a, departing_a = departing_a) + + z_oldroot = self.cut_pricing_tree(departing_a) + self.paste_pricing_tree(entering_a,z_oldroot) + self.basis.remove(departing_a) + self.basis.append(entering_a) + return(2) \ No newline at end of file