In this challenge our goal is to order neurons from the fly connectome in maximally "feedforward" way. The high level motivation is gaining insights into direction of information flow in the brain.
The FlyWire connectome is modeled as a directed graph G=(V,E) where V - nodes - represent neurons, and E - weighted edges - represent synaptic connections. In this challenge, the goal is to "feed forward" sort the nodes of this directed graph so that the sum of edge weights pointing backward is minimized (or in other words, forward pointing edge weights maximized). Specifically, your goal is to produce an order O(V) -> {1, ..., N} that assigns a unique number to each node in V, so that the sum of edge weights u -> v where O(u) < U(v) is maximized.

Data
* Connectome - directed graph (download here)
(Columns: ['Source Node  ID', 'Target Node ID', 'Edge Weight']
First few rows of data:
      Source Node  ID      Target Node ID  Edge Weight
0  720575940629970489  720575940631267655            8
1  720575940629970489  720575940620005542           16
2  720575940629970489  720575940626081319            7
3  720575940629970489  720575940652258806            5
4  720575940629970489  720575940624508563           21)
* Benchmark solution (download here)
(Columns: ['Node ID', 'Order']
First few rows of data:
              Node ID  Order
0  720575940626087322      0
1  720575940611746258      1
2  720575940610412386      2
3  720575940624724919      3
4  720575940626365768      4)
Reference code
The benchmark solution provided above was produced using this library (sfas) (see source code here - self-contained, less than 250 lines including comments). You are welcome to alter the algorithm therein or implement your own.
This is the benchmark solution code:

import random
import pandas as pd


def nodes_from_connections_table(connections_table):
    return set([r[0] for r in connections_table]).union(set([r[1] for r in connections_table]))


def adjacency_lists(connections_table):
    node_list = nodes_from_connections_table(connections_table)
    input_lists = {n: [] for n in node_list}
    output_lists = {n: [] for n in node_list}
    for r in connections_table:
        input_lists[r[1]].append(r)
        output_lists[r[0]].append(r)
    return input_lists, output_lists


def compute_order(connections, verbosity=0, random_seed=0):
    total_edge_weight = sum([r[2] for r in connections])

    input_rows, output_rows = adjacency_lists(connections)

    random.seed(random_seed)
    node_ids = list(input_rows.keys())
    random.shuffle(node_ids)

    nodes_compacted = {k: i for i, k in enumerate(node_ids)}

    def compact_dict(rows, nidx):
        cd = {}
        for r in rows:
            cn = nodes_compacted[r[nidx]]
            cd[cn] = cd.get(cn, 0) + r[2]
        return cd

    input_dict = {nodes_compacted[k]: compact_dict(v, 0) for k, v in input_rows.items()}
    output_dict = {nodes_compacted[k]: compact_dict(v, 1) for k, v in output_rows.items()}

    feedback_val = 0

    anti_parallel_edges = []
    for nd, indict in input_dict.items():
        for ind, w in indict.items():
            if nd < ind:
                if nd in input_dict[ind]:
                    anti_parallel_edges.append((ind, nd))
    for (n1, n2) in anti_parallel_edges:
        w12 = output_dict[n1][n2]
        w21 = output_dict[n2][n1]
        if w12 > w21:
            del output_dict[n2][n1]
            del input_dict[n1][n2]
            output_dict[n1][n2] -= w21
            input_dict[n2][n1] -= w21
            feedback_val += w21
        else:
            del output_dict[n1][n2]
            del input_dict[n2][n1]
            output_dict[n2][n1] -= w12
            input_dict[n1][n2] -= w12
            feedback_val += w12

    def compute_bucket_key_for_node(n):
        if not output_dict.get(n):
            return 'sinks'
        elif not input_dict.get(n):
            return 'sources'
        else:
            return sum(output_dict[n].values()) - sum(input_dict[n].values())

    buckets = {}

    def insert_to_bucket(node, key):
        if key in buckets:
            buckets[key].add(node)
        else:
            buckets[key] = {node}

    for n in range(len(nodes_compacted)):
        insert_to_bucket(node=n, key=compute_bucket_key_for_node(n))

    left_list = []
    right_list = []

    num_nodes = len(nodes_compacted)

    while len(left_list) + len(right_list) < num_nodes:
        # Ensure there are nodes in 'sinks' or 'sources' before popping
        if buckets.get('sinks') and len(buckets['sinks']) > 0:
            n = buckets['sinks'].pop()
            right_list.append(n)
        elif buckets.get('sources') and len(buckets['sources']) > 0:
            n = buckets['sources'].pop()
            left_list.append(n)
        else:
            # Get the non-empty max delta bucket
            non_empty_deltas = {key for key in buckets.keys() - {'sinks', 'sources'} if len(buckets[key]) > 0}
            if non_empty_deltas:
                max_delta = max(non_empty_deltas)
                n = buckets[max_delta].pop()
                left_list.append(n)
                feedback_val += sum(input_dict[n].values())
            else:
                raise ValueError("No valid nodes to process; check for graph inconsistency")

    right_list.reverse()

    nodes_compacted_inverse = {v: k for k, v in nodes_compacted.items()}
    res = [nodes_compacted_inverse[n] for n in left_list + right_list]

    return res


def score_order(order, connections):
    node_to_order = {node: idx for idx, node in enumerate(order)}
    forward_sum = sum(weight for source, target, weight in connections
                      if node_to_order[source] < node_to_order[target])
    return forward_sum


# Read the data and compute the order
connectome_df = pd.read_csv('../connectome.csv')
connections = list(connectome_df.itertuples(index=False, name=None))

# Compute the order
order = compute_order(connections)

# Calculate the score
forward_sum = score_order(order, connections)

# Calculate total edge weight
total_edge_weight = sum(weight for _, _, weight in connections)

# Calculate percentage of forward weight
percentage_forward = (forward_sum / total_edge_weight) * 100

print(f"Score (sum of forward edge weights): {forward_sum}")
print(f"Total edge weight: {total_edge_weight}")
print(f"Percentage of forward edge weight: {percentage_forward:.2f}%")

the code above is a benchmark, the default code. the code solution that you come up with must be better than this. Currently, the benchmark is 70.4% and 29353623 sum of weights. Create a solution better than this. The solution doesn't have to follow the same style or method as the benchmark.