diff --git a/repair/apps/asmfa/graphs/graph.py b/repair/apps/asmfa/graphs/graph.py index 09066c247..e251701db 100644 --- a/repair/apps/asmfa/graphs/graph.py +++ b/repair/apps/asmfa/graphs/graph.py @@ -724,6 +724,7 @@ def _get_affected_flows(self, solution_part): origin__activity = af.origin_activity, destination__activity = af.destination_activity ) + kwargs = { 'strategy_material': af.material.id } diff --git a/repair/apps/asmfa/graphs/graphwalker.py b/repair/apps/asmfa/graphs/graphwalker.py index 1ad106d7e..261f9b4cf 100644 --- a/repair/apps/asmfa/graphs/graphwalker.py +++ b/repair/apps/asmfa/graphs/graphwalker.py @@ -1,48 +1,79 @@ try: import graph_tool as gt - from graph_tool import util, search + from graph_tool import search from graph_tool.search import BFSVisitor except ModuleNotFoundError: class BFSVisitor: pass import copy -import numpy as np -from django.db.models import Q -import time - -from repair.apps.asmfa.models.flows import FractionFlow -from repair.apps.asmfa.models import Actor class NodeVisitor(BFSVisitor): - def __init__(self, name, amount, visited, change, - balance_factor): + def __init__(self, name, amount, change, + balance_factor, forward=True): self.id = name self.amount = amount - self.visited = visited self.change = change self.balance_factor = balance_factor + self.forward = forward + + def examine_vertex(self, u): + vertex_id = int(u) + out_degree = u.out_degree() + if not out_degree: + return + + bf = self.balance_factor[vertex_id] + sum_in_deltas = u.in_degree(weight=self.change) + balanced_delta = sum_in_deltas * bf + sum_out_f = u.out_degree(weight=self.amount) + if sum_out_f: + amount_factor = balanced_delta / sum_out_f + else: + amount_factor = balanced_delta / out_degree - def examine_edge(self, e): - """Compute the amount change on each inflow for the vertex + for e_out in u.out_edges(): + amount_delta = self.amount[e_out] * amount_factor + if self.forward: + self.change[e_out] += amount_delta + else: + if abs(self.change[e_out]) < abs(amount_delta): + self.change[e_out] = amount_delta + else: + self.change[e_out] += amount_delta + + +class NodeVisitorBalanceDeltas(BFSVisitor): + + def __init__(self, name, amount, change, + balance_factor): + self.id = name + self.amount = amount + self.change = change + self.balance_factor = balance_factor - This function is invoked on a edge as it is popped from the queue. - """ - u = e.target() + def examine_vertex(self, u): vertex_id = int(u) + out_degree = u.out_degree() + if not out_degree: + return - balanced_delta = self.change[e] * self.balance_factor[vertex_id] - edges_out = list(u.out_edges()) - sum_out_f = sum(self.amount[out_f] for out_f in edges_out) + sum_in_deltas = u.in_degree(self.change) + sum_out_deltas = u.out_degree(self.change) + bf = self.balance_factor[vertex_id] + balanced_out_deltas = sum_out_deltas / bf + balanced_delta = sum_in_deltas - balanced_out_deltas + if abs(balanced_delta) < 0.0000001: + return + sum_out_f = u.out_degree(weight=self.amount) if sum_out_f: amount_factor = balanced_delta / sum_out_f - elif edges_out: - amount_factor = balanced_delta / len(edges_out) - for e_out in edges_out: - if not (self.visited[e_out] and self.visited[e]): - self.change[e_out] += self.amount[e_out] * amount_factor - self.visited[e_out] = True + else: + amount_factor = balanced_delta / out_degree + for e_out in u.out_edges(): + amount_delta = self.amount[e_out] * amount_factor + self.change[e_out] += amount_delta def traverse_graph(g, edge, delta, upstream=True): @@ -61,47 +92,173 @@ def traverse_graph(g, edge, delta, upstream=True): Edge ProperyMap (float) The signed change on the edges """ - # Property map for keeping track of the visited edge. Once an edge has - # been visited it won't be processed anymore. + plot = False amount = g.ep.amount - visited = g.new_edge_property("bool", val=False) change = g.new_edge_property("float", val=0.0) - change[edge] = delta - visited[edge] = True + total_change = g.new_edge_property("float", val=0.0) + + if plot: + # prepare plotting of intermediate results + from repair.apps.asmfa.tests import flowmodeltestdata + flowmodeltestdata.plot_materials(g, file='materials.png') + flowmodeltestdata.plot_amounts(g,'amounts.png', 'amount') + g.ep.change = change # We are only interested in the edges that define the solution g.set_edge_filter(g.ep.include) + MAX_ITERATIONS = 20 + balance_factor = g.vp.downstream_balance_factor.a + + # make a first run with the given changes to the implementation edge # By default we go upstream first, because 'demand dictates supply' if upstream: + node = edge.source() g.set_reversed(True) - balance_factor = 1 / g.vp.downstream_balance_factor.a - node = edge.target() + balance_factor = 1 / balance_factor else: + node = edge.target() g.set_reversed(False) - balance_factor = g.vp.downstream_balance_factor.a - node = edge.source() - node_visitor = NodeVisitor(g.vp["id"], amount, visited, change, + # initialize the node-visitors + node_visitor = NodeVisitor(g.vp["id"], amount, change, + balance_factor) + node_visitor2 = NodeVisitorBalanceDeltas(g.vp["id"], amount, change, balance_factor) + + node_visitor.forward = True + total_change.a[:] = 0 + new_delta = delta + i = 0 + change[edge] = new_delta + # start in one direction search.bfs_search(g, node, node_visitor) + change[edge] = new_delta - # now go downstream, if we started upstream - # (or upstream, if we started downstream) + if plot: + ## Plot changes after forward run + g.ep.change.a[:] = change.a + flowmodeltestdata.plot_amounts(g,'plastic_deltas.png', 'change') - g.set_reversed(not g.is_reversed()) - node = edge.target() if g.is_reversed() else edge.source() - # reverse the balancing factors - node_visitor.balance_factor = 1 / node_visitor.balance_factor - # print("\nTraversing in 2. direction") + node = reverse_graph(g, node_visitor, node_visitor2, edge) search.bfs_search(g, node, node_visitor) + change[edge] = new_delta + + if plot: + ## Plot changes after backward run + g.ep.change.a[:] = change.a + flowmodeltestdata.plot_amounts(g,'plastic_deltas.png', 'change') + + # balance out the changes + search.bfs_search(g, node, node_visitor2) + change[edge] = new_delta + + # add up the total changes + total_change.a += change.a + + if plot: + ## Plot total changes + g.ep.change.a[:] = total_change.a + flowmodeltestdata.plot_amounts(g,f'plastic_deltas_{i}.png', 'change') + + node = reverse_graph(g, node_visitor, node_visitor2, edge) + + if upstream: + if node.in_degree(): + sum_f = node.in_degree(weight=total_change) + new_delta = delta - sum_f + else: + new_delta = 0 + else: + if node.out_degree(): + sum_f = node.out_degree(weight=total_change) + new_delta = delta - sum_f + else: + new_delta = 0 + i += 1 + + + while i < MAX_ITERATIONS and abs(new_delta) > 0.00001: + change.a[:] = 0 + change[edge] = new_delta + + # start in one direction + + search.bfs_search(g, node, node_visitor) + change[edge] = 0 + + if plot: + ## Plot changes after forward run + g.ep.change.a[:] = change.a + flowmodeltestdata.plot_amounts(g,'plastic_deltas.png', 'change') + + + # now go downstream, if we started upstream + # (or upstream, if we started downstream) + node = reverse_graph(g, node_visitor, node_visitor2, edge) + if upstream: + sum_f = node.out_degree(weight=total_change) + \ + node.out_degree(weight=change) + else: + sum_f = node.in_degree(weight=total_change) + \ + node.in_degree(weight=change) + new_delta = delta - sum_f + change[edge] = new_delta + search.bfs_search(g, node, node_visitor) + + + if plot: + ## Plot changes after backward run + g.ep.change.a[:] = change.a + flowmodeltestdata.plot_amounts(g,'plastic_deltas.png', 'change') + + # balance out the changes + search.bfs_search(g, node, node_visitor2) + change[edge] = 0 + + if plot: + ## Plot changes after balancing + g.ep.change.a[:] = change.a + flowmodeltestdata.plot_amounts(g,'plastic_deltas.png', 'change') + + # add up the total changes + total_change.a += change.a + + node = reverse_graph(g, node_visitor, node_visitor2, edge) + + if plot: + ## Plot total changes + g.ep.change.a[:] = total_change.a + flowmodeltestdata.plot_amounts(g,f'plastic_deltas_{i}.png', 'change') + + if upstream: + if node.in_degree(): + sum_f = node.in_degree(weight=total_change) + new_delta = delta - sum_f + else: + new_delta = 0 + else: + if node.out_degree(): + sum_f = node.out_degree(weight=total_change) + new_delta = delta - sum_f + else: + new_delta = 0 + i += 1 # finally clean up - del visited g.set_reversed(False) g.clear_filters() - return node_visitor.change + return total_change + + +def reverse_graph(g, node_visitor: NodeVisitor, node_visitor2, edge): + g.set_reversed(not g.is_reversed()) + node_visitor.balance_factor = 1 / node_visitor.balance_factor + node = edge.target() if not g.is_reversed() else edge.source() + node_visitor.forward = not node_visitor.forward + node_visitor2.balance_factor = 1 / node_visitor2.balance_factor + return node class GraphWalker: diff --git a/repair/apps/asmfa/tests/flowmodeltestdata.py b/repair/apps/asmfa/tests/flowmodeltestdata.py index 127b6019f..9ce2c87ff 100644 --- a/repair/apps/asmfa/tests/flowmodeltestdata.py +++ b/repair/apps/asmfa/tests/flowmodeltestdata.py @@ -106,8 +106,7 @@ def plastic_package_graph(): G.vp.id[10] = 'Stock 2' G.vp.id[11] = 'Waste 2' flow = G.new_edge_property("object") - eid = G.new_edge_property( - "int") # need a persistent edge id, because graph-tool can reindex the edges + eid = G.new_edge_property("int") # need a persistent edge id, because graph-tool can reindex the edges G.edge_properties["flow"] = flow G.edge_properties["eid"] = eid e = G.add_edge(G.vertex(0), G.vertex(1)) @@ -153,12 +152,23 @@ def plastic_package_graph(): return split -def plot_amounts(g, file=None): - """Plots the graph with the 'amount' property on the edges into a file""" +def plot_amounts(g, file=None, prop='amount'): + """ + Plots the graph with the property defined by `prop` + on the edges into a file + + Parameters + ---------- + g : graph_tool.Graph + file : str + prop : str, optional (default='amount') + + """ + quantities = g.ep[prop] vertex_ids = [f'{int(v)}' for v in g.vertices()] vertex_text = g.new_vertex_property("string", vals=vertex_ids) mass_text = g.new_edge_property("string", - vals=[str(round(i, 2))for i in g.ep.amount]) + vals=[str(round(i, 3)) for i in quantities]) gt.draw.graph_draw(g, vertex_size=20, vertex_text=vertex_text, vprops={"text_position": -1, "font_size": 10}, diff --git a/repair/apps/asmfa/tests/test_graph.py b/repair/apps/asmfa/tests/test_graph.py index 6da888492..f83e29c22 100644 --- a/repair/apps/asmfa/tests/test_graph.py +++ b/repair/apps/asmfa/tests/test_graph.py @@ -75,7 +75,6 @@ def test_plastic_packaging(self): Consumption --> Recycling -0.12 Recycling --> Production -0.06 """ - return plastic = flowmodeltestdata.plastic_package_graph() gw = GraphWalker(plastic) change = gw.graph.new_edge_property('float') @@ -100,47 +99,86 @@ def test_plastic_packaging(self): else: gw.graph.ep.include[e] = False result = gw.calculate(implementation_edges, deltas) + + # for each flow, test if the results are approximately + # the expected value, taking some delta due to the balancing + # into account for i, e in enumerate(result.edges()): - print(f"{result.vp.id[e.source()]} --> {result.vp.id[e.target()]} / {result.ep.material[e]}: {result.ep.amount[e]}") + print(f"{result.vp.id[e.source()]} --> {result.vp.id[e.target()]} " + f"/ {result.ep.material[e]}: {result.ep.amount[e]}") if result.vp.id[e.source()] == 'Packaging' \ and result.vp.id[e.target()] == 'Consumption' \ and result.ep.material[e] == 'plastic': expected = 5.0 - 0.3 - self.assertAlmostEqual(result.ep.amount[e], expected, 2, 'Packaging->Consumption') + self.assertAlmostEqual( + result.ep.amount[e], expected, delta=0.2, + msg='Packaging->Consumption') elif result.vp.id[e.source()] == 'Oil rig' \ and result.vp.id[e.target()] == 'Oil refinery' \ and result.ep.material[e] == 'crude oil': expected = 20.0 - 0.3 - self.assertAlmostEqual(result.ep.amount[e], expected, 2, 'Oil rig->Oil refinery') + self.assertAlmostEqual( + result.ep.amount[e], expected, delta=0.2, + msg='Oil rig->Oil refinery') elif result.vp.id[e.source()] == 'Oil refinery' \ and result.vp.id[e.target()] == 'Production' \ and result.ep.material[e] == 'plastic': expected = 4.0 - 0.24 - self.assertAlmostEqual(result.ep.amount[e], expected, 2, 'Oil refinery->Production') + self.assertAlmostEqual( + result.ep.amount[e], expected, delta=0.2, + msg='Oil refinery->Production') elif result.vp.id[e.source()] == 'Production' \ and result.vp.id[e.target()] == 'Packaging' \ and result.ep.material[e] == 'plastic': expected = 5.0 - 0.3 - self.assertAlmostEqual(result.ep.amount[e], expected, 2, 'Production->Packaging') + self.assertAlmostEqual( + result.ep.amount[e], expected, delta=0.2, + msg='Production->Packaging') elif result.vp.id[e.source()] == 'Consumption' \ and result.vp.id[e.target()] == 'Burn' \ and result.ep.material[e] == 'plastic': expected = 3.0 - 0.18 - self.assertAlmostEqual(result.ep.amount[e], expected, 2, 'Consumption->Burn') + self.assertAlmostEqual( + result.ep.amount[e], expected, delta=0.1, + msg='Consumption->Burn') elif result.vp.id[e.source()] == 'Consumption' \ and result.vp.id[e.target()] == 'Recycling' \ and result.ep.material[e] == 'plastic': expected = 2.0 - 0.12 - self.assertAlmostEqual(result.ep.amount[e], expected, 2, 'Consumption->Recycling') + self.assertAlmostEqual( + result.ep.amount[e], expected, delta=0.1, + msg='Consumption->Recycling') elif result.vp.id[e.source()] == 'Recycling' \ and result.vp.id[e.target()] == 'Production' \ and result.ep.material[e] == 'plastic': expected = 1.0 - 0.06 - self.assertAlmostEqual(result.ep.amount[e], expected, 2, 'Recycling->Production') + self.assertAlmostEqual( + result.ep.amount[e], expected, delta=0.06, + msg='Recycling->Production') else: self.assertAlmostEqual(result.ep.amount[e], - gw.graph.ep.amount[e], places=2) - + gw.graph.ep.amount[e], + delta=result.ep.amount[e]/10) + + # test if the changes in all nodes are balanced + result.ep.change.a = result.ep.amount.a - gw.graph.ep.amount.a + + for u in result.vertices(): + # dangling nodes with no in- or outflows can be ignored + if not (u.in_degree() and u.out_degree()): + continue + # for the rest, the sum of the in-deltas should equal to the + # sum of the out-deltas, adjusted with the balancing factor + sum_in_deltas = u.in_degree(result.ep.change) + sum_out_deltas = u.out_degree(result.ep.change) + balanced_out_deltas = sum_out_deltas / bf[u] + balanced_delta = sum_in_deltas - balanced_out_deltas + self.assertAlmostEqual( + balanced_delta, 0, places=4, + msg=f'Node {int(u)} not balanced, deltas in: {sum_in_deltas}, ' + f'deltas out: {sum_out_deltas}, bf: {bf[u]}, ' + f'balanced deltas out: {balanced_out_deltas}' + ) def test_milk_production(self): """Reduce milk production between Farm->Packaging @@ -1251,6 +1289,7 @@ def affect_biofuel_chain(solpart): strat_digest_factor = strat_out_digester_sum / strat_in_digester_sum self.assertAlmostEqual(sq_digest_factor, strat_digest_factor, + places=1, msg=f'the factor at actor {biodigester} in ' 'strategy is not the same as in status quo') @@ -1281,7 +1320,7 @@ def assert_balance_factor(activity): sf_out = out_flows.aggregate(amount=Sum('strategy_amount'))['amount'] sq_factor = (sq_out / sq_in) if sq_out and sq_in else 1 sf_factor = (sf_out / sf_in) if sf_out and sf_in else 1 - self.assertAlmostEqual(sq_factor, sf_factor, + self.assertAlmostEqual(sq_factor, sf_factor, 1, msg='the balance factor at actor ' f'{actor} in strategy is not the ' 'same as in status quo')