Implement model selection for temporal graphs

[IngoScholtes]

In the current implementation, the function MultiOrderModel.estimate_order requires an instance of PathData, i.e. it is not possible to directly do the model selection / order detection for time-respecting paths.

[IngoScholtes]

In the branch trp-model-selection, I Implemented a proposed generalizaton of the likelihood ratio test to walks up to a max. length in an event graph, which is one possible generalization of the model selection to time-respecting paths (or walks)

[M-Lampert]

I checked your implementation for the path extraction. With the `for'-loop-based DFS implementation, it will probably be very slow for larger graphs as the warning says:

pathpyG/src/pathpyG/algorithms/temporal.py

Lines 61 to 138 in f6c18ff

	def extract_causal_paths(g: TemporalGraph, delta: float | int = 1, max_paths: Optional[int] = None) -> PathData:
	"""Enumerate maximal time-respecting paths in a temporal graph into a [PathData][pathpyG.PathData] object.
	A maximal time-respecting path is a sequence of temporal edges `(e_1, ..., e_L)` such that
	each `e_{i+1}` continues `e_i` within `delta` (as defined by [lift_order_temporal][pathpyG.algorithms.lift_order_temporal]),
	`e_1` has no time-respecting predecessor edge, and `e_L` has no time-respecting successor edge.
	If an edge has more than one valid predecessor or successor, it is part of more than one
	maximal path, so the returned [PathData][pathpyG.PathData] can contain overlapping paths
	that share individual edges.
	Warning:
	The number of maximal time-respecting paths can grow exponentially in the branching
	factor of the temporal graph (an edge with `b` valid continuations multiplies the
	number of paths passing through it by `b`). This function enumerates paths explicitly
	and is intended as a **reference implementation for small graphs** - e.g. to validate
	higher-order model likelihoods computed directly from the event graph - and not as a
	scalable path extraction routine. For large or densely connected temporal graphs, build
	a [MultiOrderModel][pathpyG.MultiOrderModel] directly via
	[MultiOrderModel.from_temporal_graph][pathpyG.MultiOrderModel.from_temporal_graph] instead,
	which never materializes individual paths. Use `max_paths` to fail fast rather than
	exhausting memory on inputs that are too large.
	Args:
	g: The temporal graph to extract paths from.
	delta: The maximum time difference between two consecutive edges in a path.
	max_paths: If set, raise a `RuntimeError` once more than this many maximal paths have
	been found.
	Returns:
	PathData: A [PathData][pathpyG.PathData] object containing all maximal time-respecting
	paths, each stored with weight 1.0.
	Examples:
	>>> import pathpyG as pp
	>>> g = pp.TemporalGraph.from_edge_list([("a", "b", 1), ("b", "c", 5), ("c", "d", 9), ("c", "e", 9)])
	>>> paths = pp.algorithms.extract_causal_paths(g, delta=4)
	>>> sorted(paths.get_walk(i) for i in range(paths.num_paths))
	[('a', 'b', 'c', 'd'), ('a', 'b', 'c', 'e')]
	"""
	num_events = g.data.edge_index.size(1)
	paths = PathData(mapping=g.mapping, device=g.data.edge_index.device)
	if num_events == 0:
	return paths
	ho_index = lift_order_temporal(g, delta)
	src = g.data.edge_index[0].tolist()
	dst = g.data.edge_index[1].tolist()
	successors: list[list[int]] = [[] for _ in range(num_events)]
	has_predecessor = [False] * num_events
	for s, d in ho_index.t().tolist():
	successors[s].append(d)
	has_predecessor[d] = True
	roots = [i for i in range(num_events) if not has_predecessor[i]]
	node_seqs: list[list[int]] = []
	# DFS over the event graph from each root to each reachable leaf. The event graph is
	# acyclic because continuations strictly increase in time, so this always terminates.
	stack: list[tuple[int, list[int]]] = [(root, [src[root], dst[root]]) for root in roots]
	while stack:
	event, seq = stack.pop()
	succs = successors[event]
	if not succs:
	node_seqs.append(seq)
	if max_paths is not None and len(node_seqs) > max_paths:
	raise RuntimeError(
	f"Number of maximal time-respecting paths exceeds max_paths={max_paths}. "
	"extract_causal_paths does not scale to graphs with a high branching factor; "
	"consider building a MultiOrderModel directly from the temporal graph instead."
	)
	else:
	for succ in succs:
	stack.append((succ, seq + [dst[succ]]))
	node_id_seqs = [g.mapping.to_ids(seq) for seq in node_seqs]
	paths.append_walks(node_id_seqs, weights=[1.0] * len(node_seqs))
	return paths

Note that I already implemented a Tensor-based causal-path extraction at some point. You can find it in the implementation concepts notebook: https://www.pathpy.net/0.2.0-dev/tutorial/implementation_concepts/#longest-path-extraction
It should be considerably faster, since it parallelises the walks from all possible roots. It can get prohibitive in terms of memory, though.

We should probably think about a sampling-based/random-walk-based approach in that case. We could use the random-walk-implementation from torch-cluster for that: https://github.com/rusty1s/pytorch_cluster#randomwalk-sampling

[IngoScholtes]

Indeed, this method is slow but it is also not what is used for the model selection :-)

Let me explain:

The extract causal paths method can be used to extract a PathData object that contains all maximal time-respecting paths, which is mainly for testing and for illustrative examples in small toy example temporal graphs. In principle, this object could be passed to the estimate_order function to perform the model selection for a temporal graph.

However, the new estimate_order method in the class MultiOrderModel now takes PathData object as optional argument. only, to perform the model selection for a given set of paths or walks. By default, the statistics of paths up to a max. length of k is used for this purpose, which is calculated in the same shot as the edge weight for higher-order models. However, this requires to calculate the MultiOrderModel from time-respecting walks up to a given max. length to be consistent. This is done with the new function MultiOrderModel.from_time_respecting_walks which internally invokes algorithms.temporal.walk_counts based on the lift_order function. However, this is done without for loops again, so this should be efficient.

For MultiOrderModels created from PathData, there is a new function attach_path_statistics that computes the same statistic that can be used by the test, so for both data sources the same method can be used.

It is possible that this code can be further simplified and it definitely requires more testing, but at least the results look reasonable and all previous tests suceed. With the help of CLaude, I also found a bug in the previous order estimation that did not properly consider path weights, so this should be an improvement anyway.