A little functional thinking might nevertheless help. Like the view point that
mode directed tabling builds a map from some tuples to a lattice. m : T^n → L,
and a non-registered tuple has the bottom value ⊥ of that lattice.
Otherwise, U R welcome and I have to pass the ball to Jan W., and
thank him for bringing mode directed tabling to SWI-Prolog. Its really a
enjoyable feature and matches amazingly AoC 2023.
Sorry to dwell on part 1 day 21 @j4n_bur53, which by now has been solved 5 different ways on this thread, but I think that the fastest possible solution would be the following which takes advantage of the odd even counting:
I think everything we’ve discussed so far will have higher time complexity than that for part 1.
EDIT: Technically it doesn’t even have to be shortest paths – it can be any path at all between S and the rest of the nodes.
Here comes a 6-th way to solve it and maybe some insight?
You can also use CHR to implement mode directed tabling. Here is a
code that solves part 1 day 21. Internally I named it solutionchr.p, since
I derived it from my mode directed tabling solution.p:
:- use_module(library(chr)).
:- chr_constraint go/1, dist/4, neigh/2.
go(N) ==> garden(X, Y, 0'S),
neigh(0,1), neigh(1,0), neigh(0,-1), neigh(-1,0),
dist(X, Y, N, 0).
dist(X,Y,N,D1) \ dist(X,Y,N,D2) <=> D1 =< D2 | true.
dist(H, J, N, L), neigh(DX, DY) ==> L < N, X is H+DX, Y is J+DY,
garden(X, Y, C), C \= 0'# | R is L+1, dist(X, Y, N, R).
count(N, X) :-
go(N),
aggregate_all(count,
(current_chr_constraint(dist(_,_,N,K)), (N+K) mod 2 =:= 0), X).
But I am not a frequent CHR user, so maybe this is not the
best way to write it with CHR. But it seems to work fine.
Here the result for the tester 11x11 from the part 1 description:
?- count(6, X).
X = 16
Edit 27.12.2023
Credits go to this SWI-Prolog discourse thread, where it was made the
impression that for the Dijkstra algorithm some CHR with priorities
would be needed, or alternatively that we would want to have
a fibonacci heap. Well this is even not wrong, if one wants to
remodel some common strategy how Dijkstra algorithm proceeds,
that allows early stopping for the shortes path. But if you compute
distances everywhere you don’t aim at early stopping, you anyway
need to exhaust. You can inspect my code above, I basically used the
same rules as presented Dijkstra CHR thread but without any priorities,
and since the min operator on an integer lattice it is used, I tend
to claim that any operation order saturation will give the same.
I tried the above on the whole input – seems to take ages and then blow up on account of memory. Did you get it working?
Also this 15 LOC solution completes in about 1.3M inferences for the 64 step problem which I think might be faster than the min based tabling solution above, even though it doesn’t use any optimisation. It suggests that maybe a bit too much is happening behind the scenes.
Yeah my CHR solution doesn’t scale so well:
?- set_prolog_flag(stack_limit, 2_147_483_648).
true.
?- time(count(64,X)), write(X), nl, fail.
% 901,974,873 inferences, 84.953 CPU in 85.473 seconds (99% CPU, 10617324 Lips)
Compare to the mode directed tabling:
But I see the same accepted answer for both. So both seem to be correct.
@j4n_bur53 I did a quick hard-coded to 64 steps solution on this basis to see what I’d get. Would you mind comparing it to your min tabling solution for time (not inferences, they seem to completely not correspond to time when comparing tabled solutions) ?
With the 15 LOC implementation I get 0.7-1.5 sec to calculate for 64 steps on the whole input. With this walk + odd/even I get 0.16-0.26 sec: a massive improvement even in the worse case scenario.
grid(_,_,[]) --> [].
grid(X0,_,Xs) --> "\n", {X is X0+1}, !, grid(X,1,Xs).
grid(X0,Y0,[p(X0-Y0,V)|Xs]) --> [V],{Y is Y0+1}, !, grid(X0,Y,Xs).
grid(X0,Y0,Xs) --> `.`, {Y is Y0+1}, !, grid(X0,Y,Xs).
next(X0-Y, X-Y) :- X is 1+X0, p(X-Y,_).
next(X0-Y, X-Y) :- X is -1+X0, p(X-Y,_).
next(X-Y0, X-Y) :- Y is 1+Y0, p(X-Y,_).
next(X-Y0, X-Y) :- Y is -1+Y0, p(X-Y,_).
tourist(_,[],Acc,Acc) :- !.
tourist(Vs,[_-N|T],A,C) :- N > 64, !, tourist(Vs,T,A,C).
tourist(Vs,[H-_|T],A,C) :- ht_get(Vs, H, _), !, tourist(Vs,T,A,C).
tourist(Vs,[H-N0|T],A,C) :-
(N0 mod 2 =:= 0-> A1 is A+1; A1=A),
ht_put(Vs, H, N0), N is N0+1,
findall(X-N, (next(H,X), \+ ht_get(Vs,X,_)), Next),
append(T,Next,Stack),
!, tourist(Vs, Stack,A1,C).
solve(File, Part1) :-
phrase_from_file(grid(1, 1, Ps0), File),
include([p(_,V)]>>(V\=0'#), Ps0, Ps),
retractall(p(_,_)), maplist(asserta, Ps),
ht_new(Vs), tourist(Vs, [66-66-0], 0, Part1).
Its not as fast as mode directed tabling. I guess your append/3
call slows it down. Here again the mode directed tabling:
And this is what I get for your code on the same machine:
?- time(solve(X)).
% 2,893,213 inferences, 0.125 CPU in 0.128 seconds (97% CPU, 23145704 Lips)
BTW your solution gives the same result, so I guess its correct.
The real benchmark is part 2. What is the speed there, can you
do the extrapolation with your new take?
Ok, I could break my own speed record for day 21, which was
based on tabling. Lets say tabling means using the SWI-Prolog tabling
directive. Now I have a solution that uses memoization, i.e. assertz/1 and
retract/1, and BFS flood fill, without computing shortest distances at all.
Here are some timing results, the BFS flood fill beats the
shortest distances in part 1, didn’t try part 2 yet:
/* tabling and shortest distance */
?- time(count(64,X)).
% 435,474 inferences, 0.063 CPU in 0.058 seconds (108% CPU, 6967584 Lips)
/* memoization and BFS flood fill */
?- time(count(64,X)).
% 181,776 inferences, 0.016 CPU in 0.037 seconds (42% CPU, 11633664 Lips)
The results agree. Here is the straight forward code, it makes use
of logical view of retract/1 predicate as available in most Prolog systems,
to enumerate a current frontier, and replace it by a new frontier,
without interferring in the enumeration of the frontier:
:- dynamic flooded/2.
:- dynamic frontier/2.
neigh(0, 1).
neigh(1, 0).
neigh(0, -1).
neigh(-1, 0).
init :-
retractall(frontier(_,_)),
retractall(flooded(_,_)),
garden(Y, X, 0'S),
assertz(frontier(Y, X)),
assertz(flooded(Y, X)),
fail.
init.
step :-
retract(frontier(Y, X)),
neigh(DY, DX),
H is Y+DY, J is X+DX,
garden(H, J, C), C \= 0'#,
\+ flooded(H, J),
assertz(frontier(H, J)),
assertz(flooded(H, J)),
fail.
step.
solve(N) :-
init,
between(1, N, _),
step,
fail.
solve(_).
count(N, C) :-
solve(N),
aggregate_all(count, (flooded(X,Y), (N+X+Y) mod 2 =:= 0), C).
Edit 30.12.2023
Credits go to @chansey97 for grilling me over the last days, so that flood
fill almost came out of my ears. The (N+X+Y) mod 2 discovery is already
a few days old though, predates a little bit the heated fight.
Of course I want also thank @emiruz for the challenge, and @jan for not
banning me despite using strong language here and then.
Day 23 was straightforward, here is an implementation in 46 LOC. I feel these later day problems are getting a bit too “just-so”. I feel like I learned more transferable things from the earlier problems.
Here is day 24 in 29 LOC. Part1 was easy using either the simplex or clpq library. Part 2 is easy with Python (due to sympy or z3) but hard in Prolog.
The problem: you have objects in 3D space defined by their positions X and their velocity V – both given as 3D vectors. Define a new object given by \bar{X},\bar{V}, which collides with all other objects at some time. I.e. X_i + V_iT_i = \bar{X} + \bar{V}T_i where i indexes all the other objects.
For languages with a symbolic or SMT solver, you can just transcribe this and get a solution.
In Prolog, under clpq, the non-linear constraints remain passive, so it produces no solutions. It so happens that the position and velocity are specified to be integers, but under clpfd, the problem requires domains to be specified. All the numbers in the question have massive ranges other than velocity, which I only know because I solved the problem by other means, so in my answer I specify is the domain for velocity: deeply unsatisfying.
Then there are the idiosyncrasies of clpfd. The domain for velocity is -900 to 900, but if I change it to -1000 to 1000, I get “variables are uninstantiated” related errors. It was surprising that the range of a domain could cause errors like those in the first place!
@BarrensZeppelin, am I right in saying you had the same issues?
Is there a neat way to tackle this in Prolog?
I experienced the same issues, yes. I used an even smaller domain for the velocities.
It turns out that you can turn the problem into solving a system of linear equations. See here. But I didn’t explore that idea.
This is an excellent observation – I’ll look to refactor into it, but I still feel there is a SMT solver shaped hole in the Prolog toolbox 
For posterity the crux of the linearisation idea is as follows. Note, it is simple but far from obvious to the uninitiated.
Starting from \bar{X} + \bar{V}T_i = X_i + V_iT_i, we rearrange to \bar{X} - X_i = T_i(V_i-\bar{V}). At this point you have to know that it is the case that any vector multiplied by its cross product is zero. I.e. (\bar{X} - X_i)\cdot (\bar{V} - V_i) = T_i(V_i-\bar{V})\cdot (\bar{V} - V_i) = 0, leaving us with (\bar{X} - X_i)\cdot (\bar{V} - V_i) = 0: a so-called “bilinear” form. NB: the operation is cross-product NOT multiplication.
The non-linear term \bar{X}\bar{V} can be eliminated by reproducing the equations for different values of i and equating them. NB: E.g. pick i=0 and i=1 and equate the two resulting equations. Doing this twice (e.g. with 0 and 2 next time) creates two vector equations, which when flattened out (since each vector has 3 components) gives 6 equations for 6 unknowns, which can then be solved in the usual Ax=b way: another stumbling block for Prolog here since matrix ops are not in the standard library, but hopefully clpq will do.
Haven’t been following this in detail, but looking at the original question there are a few issues:
clpfd can be used to find solutions.T (time) is unbounded at the high end. That means constraints like A-X #= (I-D)*T result in unbounded X. No matter what clp you’re using, the more “explicit” constaints you can define, e.g., time and space, the better to minimize the search space that labelling predicates have to explore.clpq can be used for continuous domain problems, but as you pointed out it may not help given its passive support for non-linear constraints. Possible solutions to explore: linearize the problem, as you’ve shown in subsequent posts, or find a continuous domain clp that supports active non-linear constraints, like add-on pack clpBNR (full disclosure: I maintain it.).I wanted to but clpBNR seemed to have problems even with part1. The operative part is this:
cross(Min, Max, h(Px1-Py1-_,Vx1-Vy1-_), h(Px2-Py2-_, Vx2-Vy2-_)) :-
{ (Y - Py1) / Vy1 = (X - Px1) / Vx1 },
{ (Y - Py2) / Vy2 = (X - Px2) / Vx2 },
{ X >= Min }, { X =< Max },
{ Y >= Min }, { Y =< Max }.
I changed “=” to “==” to make it clpBNR compatible and ran the same code which works with clpq. It takes clpq maybe 3 seconds to complete. I stopped the clpBNR version after about 2 minutes.
Full code is here.
Input file is here in case you want to replicate it.
So here is a bilinear form solution in Prolog and Python – they should be the same solution. It works in Python but I just can’t make it work in Prolog under clpq. I don’t see why because its linear and not under specified. @BarrensZeppelin , can you see at a glance what I’ve done wrong? It looks like a clpq issue to me…
Put this data into sample.txt:
19, 13, 30 @ -2, 1, -2
18, 19, 22 @ -1, -1, -2
20, 25, 34 @ -2, -2, -4
12, 31, 28 @ -1, -2, -1
20, 19, 15 @ 1, -5, -3
Prolog version (call solve("sample.txt", Part2)):
:- use_module(library(dcg/basics)).
:- use_module(library(clpq)).
line(h(A-B-C,D-E-F)) -->
number(A), ", ", number(B), ", ", number(C), " @",
blanks, number(D), ",", blanks, number(E), ",", blanks, number(F).
lines([]) --> [].
lines([L|Ls]) --> line(L), "\n", lines(Ls).
lines([L]) --> line(L).
cp(A-B-C, D-E-F, X-Y-Z) :- X=B*F-C*E, Y=C*D-A*F, Z=A*E-B*D.
sub(A-B-C, E-F-G, X-Y-Z) :- X=A-E, Y=B-F, Z=C-G.
bilinear(h(X1,V1), h(X2,V2), X-Y-Z, I-J-K) :-
% X_hat x (V1-V2) + (X1-X2) x V_hat = X1 x V1 - X2 x V2
% where X_hat = X,Y,Z, V_hat = I,J,K.
sub(V1,V2,V_sub), sub(X1,X2,X_sub),
cp(X-Y-Z, V_sub, X_cp), cp(X_sub, I-J-K, V_cp),
cp(X1,V1,X1_cp), cp(X2,V2,X2_cp),
{ X_cp + V_cp = X1_cp - X2_cp }.
solve(File, Part2) :-
phrase_from_file(lines(Ls), File),
Ls=[A,B,C|_],
bilinear(A, B, X-Y-Z, I-J-K),
bilinear(B, C, X-Y-Z, I-J-K),
Part2 is X+Y+Z.
Python:
from sympy import symbols, Eq, solve, Matrix
# Define symbols for the unknowns
X, Y, Z, I, J, K = symbols('X Y Z I J K')
# Sample data
sample_data = """
19, 13, 30 @ -2, 1, -2
18, 19, 22 @ -1, -1, -2
20, 25, 34 @ -2, -2, -4
12, 31, 28 @ -1, -2, -1
20, 19, 15 @ 1, -5, -3
"""
sample_objects = [
([int(n) for n in pos.split(',')],
[int(n) for n in vel.split(',')])
for line in sample_data.strip().split('\n')
for pos, vel in [line.split('@')]
]
# Function to calculate the cross product of two vectors
def cross_product(v1, v2):
return Matrix([v1[1]*v2[2] - v1[2]*v2[1], v1[2]*v2[0] -\
v1[0]*v2[2], v1[0]*v2[1] - v1[1]*v2[0]])
cross_product_equations = []
# For the first 3 pairs of objects, create an equation based on
# the cross product
for i in range(3):
for j in range(i + 1, len(sample_objects)):
(Xi, Yi, Zi), (Ii, Ji, Ki) = sample_objects[i]
(Xj, Yj, Zj), (Ij, Jj, Kj) = sample_objects[j]
# Convert to Matrix form
Pi = Matrix([Xi, Yi, Zi])
Vi = Matrix([Ii, Ji, Ki])
Pj = Matrix([Xj, Yj, Zj])
Vj = Matrix([Ij, Jj, Kj])
# Create equation for the difference in cross products
eq = Eq(cross_product(Matrix([X, Y, Z]), Vi - Vj) + \
cross_product(Matrix([I, J, K]), Pj - Pi),
cross_product(Pi, Vi) - cross_product(Pj, Vj))
cross_product_equations.append(eq)
# Solve the system of cross product equations
cross_product_solution = \
solve(cross_product_equations, (X, Y, Z, I, J, K))
print(cross_product_solution)
EDIT: Works now; wasn’t a clpq problem.
Just starting with the example data in the original question (converted to rule):
lines([ h(19-13-30, -2- 1- -2),
h(18-19-22, -1- -1- -2),
h(20-25-34, -2- -2- -4),
h(12-31-28, -1- -2- -1),
h(20-19-15, 1- -5- -3)
]).
and the clpBNR form of the cross predicate:
cross(Min, Max, h(Px1-Py1-_,Vx1-Vy1-_), h(Px2-Py2-_, Vx2-Vy2-_)) :-
{ (Y - Py1) / Vy1 == (X - Px1) / Vx1 },
{ (Y - Py2) / Vy2 == (X - Px2) / Vx2 },
{ X >= Min }, { X =< Max },
{ Y >= Min }, { Y =< Max }.
yields:
?- lines(Ls), member(A,Ls), member(B,Ls), A @> B, cross(7,27,A,B).
Ls = [h(19-13-30, -2-1- -2), h(18-19-22, -1- -1- -2), h(20-25-34, -2- -2- -4), h(12-31-28, -1- -2- -1), h(20-19-15, 1- -5- -3)],
A = h(19-13-30, -2-1- -2),
B = h(18-19-22, -1- -1- -2) ;
Ls = [h(19-13-30, -2-1- -2), h(18-19-22, -1- -1- -2), h(20-25-34, -2- -2- -4), h(12-31-28, -1- -2- -1), h(20-19-15, 1- -5- -3)],
A = h(20-25-34, -2- -2- -4),
B = h(19-13-30, -2-1- -2) ;
Ls = [h(19-13-30, -2-1- -2), h(18-19-22, -1- -1- -2), h(20-25-34, -2- -2- -4), h(12-31-28, -1- -2- -1), h(20-19-15, 1- -5- -3)],
A = h(20-25-34, -2- -2- -4),
B = h(20-19-15, 1- -5- -3) ;
Ls = [h(19-13-30, -2-1- -2), h(18-19-22, -1- -1- -2), h(20-25-34, -2- -2- -4), h(12-31-28, -1- -2- -1), h(20-19-15, 1- -5- -3)],
A = h(20-19-15, 1- -5- -3),
B = h(19-13-30, -2-1- -2) ;
Ls = [h(19-13-30, -2-1- -2), h(18-19-22, -1- -1- -2), h(20-25-34, -2- -2- -4), h(12-31-28, -1- -2- -1), h(20-19-15, 1- -5- -3)],
A = h(20-19-15, 1- -5- -3),
B = h(18-19-22, -1- -1- -2) ;
false.
This is a bit hard to check against the example output, so I added [X,Y]:
cross0(Min, Max, h(Px1-Py1-_,Vx1-Vy1-_), h(Px2-Py2-_, Vx2-Vy2-_),[X,Y]) :-
{ (Y - Py1) / Vy1 == (X - Px1) / Vx1 },
{ (Y - Py2) / Vy2 == (X - Px2) / Vx2 },
{ X >= Min }, { X =< Max },
{ Y >= Min }, { Y =< Max }.
producing:
?- lines(Ls), member(A,Ls), member(B,Ls), A @> B, cross0(7,27,A,B,[X,Y]).
Ls = [h(19-13-30, -2-1- -2), h(18-19-22, -1- -1- -2), h(20-25-34, -2- -2- -4), h(12-31-28, -1- -2- -1), h(20-19-15, 1- -5- -3)],
A = h(19-13-30, -2-1- -2),
B = h(18-19-22, -1- -1- -2),
Y:: 15.3333333333333339...,
X:: 14.3333333333333339... ;
Ls = [h(19-13-30, -2-1- -2), h(18-19-22, -1- -1- -2), h(20-25-34, -2- -2- -4), h(12-31-28, -1- -2- -1), h(20-19-15, 1- -5- -3)],
A = h(20-25-34, -2- -2- -4),
B = h(19-13-30, -2-1- -2),
Y:: 16.6666666666666679...,
X:: 11.6666666666666661... ;
Ls = [h(19-13-30, -2-1- -2), h(18-19-22, -1- -1- -2), h(20-25-34, -2- -2- -4), h(12-31-28, -1- -2- -1), h(20-19-15, 1- -5- -3)],
A = h(20-25-34, -2- -2- -4),
B = h(20-19-15, 1- -5- -3),
Y:: 24.0000000000000000...,
X:: 19.0000000000000000... ;
Ls = [h(19-13-30, -2-1- -2), h(18-19-22, -1- -1- -2), h(20-25-34, -2- -2- -4), h(12-31-28, -1- -2- -1), h(20-19-15, 1- -5- -3)],
A = h(20-19-15, 1- -5- -3),
B = h(19-13-30, -2-1- -2),
Y:: 11.7777777777777786...,
X:: 21.4444444444444429... ;
Ls = [h(19-13-30, -2-1- -2), h(18-19-22, -1- -1- -2), h(20-25-34, -2- -2- -4), h(12-31-28, -1- -2- -1), h(20-19-15, 1- -5- -3)],
A = h(20-19-15, 1- -5- -3),
B = h(18-19-22, -1- -1- -2),
Y:: 20.6666666666666679...,
X:: 19.6666666666666679... ;
false.
Since there is no constraint on the time, this produces all the solutions including those for T<0. This is consistent with the 5 solutions documented in the original question and takes about 0.2 sec. on my hardware.
The easiest way to omit solutions for T<0 is to note that each of the constraint LHS (or RHS) defines implicitly the value of T so adding that (plus taking advantage of clpBNR declarations to simplify a bit):
cross1(Min, Max, h(Px1-Py1-_,Vx1-Vy1-_), h(Px2-Py2-_, Vx2-Vy2-_),[X,Y]) :-
[X,Y]::real(Min,Max),
{ (Y - Py1) / Vy1 == (X - Px1) / Vx1,
(Y - Py2) / Vy2 == (X - Px2) / Vx2,
(Y - Py1) / Vy1 >= 0, (Y - Py2) / Vy2 >= 0 % constrain implicit T
}.
produces:
?- lines(Ls), member(A,Ls), member(B,Ls), A @> B, cross1(7,27,A,B,[X,Y]).
Ls = [h(19-13-30, -2-1- -2), h(18-19-22, -1- -1- -2), h(20-25-34, -2- -2- -4), h(12-31-28, -1- -2- -1), h(20-19-15, 1- -5- -3)],
A = h(19-13-30, -2-1- -2),
B = h(18-19-22, -1- -1- -2),
X:: 14.3333333333333339...,
Y:: 15.3333333333333339... ;
Ls = [h(19-13-30, -2-1- -2), h(18-19-22, -1- -1- -2), h(20-25-34, -2- -2- -4), h(12-31-28, -1- -2- -1), h(20-19-15, 1- -5- -3)],
A = h(20-25-34, -2- -2- -4),
B = h(19-13-30, -2-1- -2),
X:: 11.6666666666666661...,
Y:: 16.6666666666666679... ;
false.
which is the two “positive” answers in the original question (and saves a little time). Also not that this problem only requires constraint propagation, i.e., no labeling step is necessary.
So I have some confidence that the code is correct. But the next question is how does it scale from 5 data points in the question example to 300 points(“input.txt”). That’s a huge increase in complexity C(300,2) = 44850 vs. C(5,2) = 10 (according to an online calculator) - a factor of 4485. So if the time for N=5 is ~0.2 secs., I would expect the time for N=300 to be around 15 minutes.
The reason clpq is that quick is that that’s the time taken just to process the constraints - there’s no propagation time since non-linear constraints remain passive.
The problem can be rephrased to be purely linear if you wish but I found that it ran fastest in the form presented. Sorry if this is naive, but I thought that processing more constraints should imply fewer possibilities considered and therefore quicker processing time, else what is the point?
T<0 are not solutions, and should not be returned. clpq does not return them.
Also, please see part 2 in the post above. I can’t make it work in clpq but I can’t seem to convert it to clpBNR either 
I don’t understand what you’re trying to say. More constraints should prune the search space, but “processing” the constraints, i.e., adding them to the constraint network in the first place will actually take longer - there’s more constraints to process. Don’t confuse this initial “compilation” step with constraint propagation (consistency checking) and subsequent labelling (search).
As previously noted, clpq makes no attempt to propagate non-linear constraints and just defers until they become linear. (OTOH, clpBNR has no internal solver for systems of linear constraints, but that can be done externally, e.g., using library(simplex). So there are tradeoffs.)
Then you have to add constraints such that no solutions are generated when T<0; otherwise you haven’t fully constrained the problem.
With a little further experimentation, I also see that the question example is more than a little confusing. It does find path crossings but I don’t think they’re collisions because the particles arrive at the crossing point at different times. If I add a predicate that calculates the time each particle reaches the crossing, you see that they’re different:
time_at(h(Px-Py-_,Vx-Vy-_),[X,Y],T) :- % calculate time T of arrival at [X,Y]
{T is (X-Px)/Vx, T is (Y-Py)/Vy}. % just use interval arithmetic to find T
Now:
?- lines(Ls), member(A,Ls), member(B,Ls), A @> B, cross1(7,27,A,B,[X,Y]), time_at(A,[X,Y],Ta), time_at(B,[X,Y],Tb).
Ls = [h(19-13-30, -2-1- -2), h(18-19-22, -1- -1- -2), h(20-25-34, -2- -2- -4), h(12-31-28, -1- -2- -1), h(20-19-15, 1- -5- -3)],
A = h(19-13-30, -2-1- -2),
B = h(18-19-22, -1- -1- -2),
X:: 14.3333333333333339...,
Y:: 15.3333333333333339...,
Ta:: 2.3333333333333335...,
Tb:: 3.6666666666666665... ;
Ls = [h(19-13-30, -2-1- -2), h(18-19-22, -1- -1- -2), h(20-25-34, -2- -2- -4), h(12-31-28, -1- -2- -1), h(20-19-15, 1- -5- -3)],
A = h(20-25-34, -2- -2- -4),
B = h(19-13-30, -2-1- -2),
X:: 11.6666666666666661...,
Y:: 16.6666666666666679...,
Ta:: 4.1666666666666670...,
Tb:: 3.6666666666666665... ;
false.
Note that Ta and Tb, the times the particles arrive at [X,Y] are different, so how is that a collision? Maybe I’m missing something. (I also note that the Z dimension isn’t considered at all, so not really a realistic problem statement IMO.)
So what does clpq return?
If by part 2 you mean the linear Python combo, I’ll defer to others as I don’t really understand what you’re doing. I think I posted a simple clpBNR solution, based on your original code, that agrees with the original question example (at least to the extent documented). It’s going to take a while (~15 min. ?) to run it on the full data set, and I’m not sure how to check the results in any case.