minimum spanning tree lua HPPL

[robmio]

Good morning Albert Chan, here is the program, which seems to work, but it gives wrong results. For instance:
Points: = {{3,4.5},{1,5},{5,6.5},{2.5,5},{7,3},{1.5,1},{4,2.5},{2,2},{8,4},{3,9},{4.5,0.5}​,{5,2}};
Edges: = {{3,4.5,1,5,3,4.5,5,6.5,3,4.5,2.5,5,3,4.5,7,3,3,4.5,1.5,1,3,4.5,4,2.5,3,4.5,2,2,​3,4.5,8,4,3,4.5,3,9,3,4.5,4.5,0.5,3,4.5,5,2},{1,5,5,6.5,1,5,2.5,5,1,5,7,3,1,5,1.​5,1,1,5,4,2.5,1,5,2,2,1,5,8,4,1,5,3,9,1,5,4.5,0.5,1,5,5,2},{5,6.5,2.5,5,5,6.5,7,​3,5,6.5,1.5,1,5,6.5,4,2.5,5,6.5,2,2,5,6.5,8,4,5,6.5,3,9,5,6.5,4.5,0.5,5,6.5,5,2}​,{2.5,5,7,3,2.5,5,1.5,1,2.5,5,4,2.5,2.5,5,2,2,2.5,5,8,4,2.5,5,3,9,2.5,5,4.5,0.5,​2.5,5,5,2},{7,3,1.5,1,7,3,4,2.5,7,3,2,2,7,3,8,4,7,3,3,9,7,3,4.5,0.5,7,3,5,2},{1.​5,1,4,2.5,1.5,1,2,2,1.5,1,8,4,1.5,1,3,9,1.5,1,4.5,0.5,1.5,1,5,2},{4,2.5,2,2,4,2.​5,8,4,4,2.5,3,9,4,2.5,4.5,0.5,4,2.5,5,2},{2,2,8,4,2,2,3,9,2,2,4.5,0.5,2,2,5,2},{​8,4,3,9,8,4,4.5,0.5,8,4,5,2},{3,9,4.5,0.5,3,9,5,2},{4.5,0.5,5,2}};


To compare the wrong result to the right result, I would like to send you an image, but I cannot find the "attach image" option. So I'm sending you a message on the forum.
Thanks for your availability, greetings, Roberto.

Code:


#cas
dist(x1,y1,x2,y2):=
BEGIN
RETURN sqrt((ABS(x2-x1))^2+(ABS(y2-y1))^2);
END;

same(x1,y1,x2,y2,edg):=
BEGIN
FOR jj FROM 1 TO length(edg) DO
IF (x1==edg(jj)(1) AND y1==edg(jj)(2) AND x2==edg(jj)(3) AND y2==edg(jj)(4)) OR (x2==edg(jj)(1) AND y2==edg(jj)(2) AND x1==edg(jj)(3) AND y1==edg(jj)(4))
THEN
RETURN true;
END;
END;
RETURN false;
END;

doesexist(lst,item):=
BEGIN
FOR uu FROM 1 TO length(lst) DO
IF lst(uu)==item THEN
RETURN true;
END;
END;
RETURN false;
END;

mstRob(points,edges):=
BEGIN
verticesi:={points(1)};
tree:={};
WHILE length(verticesi) != length(points) DO
  lnn:=10000;
  temp1:={};
  temp2:={};
  FOR j FROM 1 TO length(verticesi) DO
    FOR pp FROM 1 TO length(points) DO
      IF doesexist(verticesi,points(pp))==false THEN
        IF dist(points(pp)(1),points(pp)(2),verticesi(j)(1),verticesi(j)(2))<lnn AND same(points(pp)(1),points(pp)(2),verticesi(j)(1),verticesi(j)(2),edges) THEN
          lnn:=dist(points(pp)(1),points(pp)(2),verticesi(j)(1),verticesi(j)(2));
          temp1:=points(pp);
          temp2:=verticesi(j);
        END;
      END;
    END;
  END;
  verticesi:=append(verticesi,temp1);
  tree:=append(tree,{temp1(1),temp1(2),temp2(1),temp2(2)});
END;
RETURN tree;
END;
#end

[robmio]

I apologize, but I submitted the wrong attachment. See the two new attachments. It seems that the algorithm I presented does not take into account the minimum sum, but limits itself to joining the points in ascending order.
I can't figure out where I went wrong, translating the algorithm written in LUA (see https://github.com/asundheim/lua_minimum...er/mst.lua).

[Albert Chan]

It is more natural to use complex number to represent points: (x,y) ≡ x + y*i
List of valid edges has each edge with 2 points. (to make a segment)

With your example, we have:

Code:

pts := [ /* total 12 points */
(3,4.5),(1,5),(5,6.5),(2.5,5),(7,3),(1.5,1),
(4,2.5),(2,2),(8,4),(3,9),(4.5,0.5),(5,2)];

edgs := [ /* total 12*11/2 = 66 edges */
[(3,4.5),(1,5)],
[(3,4.5),(5,6.5)],
[(3,4.5),(2.5,5)],
[(3,4.5),(7,3)],
[(3,4.5),(1.5,1)],
[(3,4.5),(4,2.5)],
[(3,4.5),(2,2)],
[(3,4.5),(8,4)],
[(3,4.5),(3,9)],
[(3,4.5),(4.5,0.5)],
[(3,4.5),(5,2)],
[(1,5),(5,6.5)],
[(1,5),(2.5,5)],
[(1,5),(7,3)],
[(1,5),(1.5,1)],
[(1,5),(4,2.5)],
[(1,5),(2,2)],
[(1,5),(8,4)],
[(1,5),(3,9)],
[(1,5),(4.5,0.5)],
[(1,5),(5,2)],
[(5,6.5),(2.5,5)],
[(5,6.5),(7,3)],
[(5,6.5),(1.5,1)],
[(5,6.5),(4,2.5)],
[(5,6.5),(2,2)],
[(5,6.5),(8,4)],
[(5,6.5),(3,9)],
[(5,6.5),(4.5,0.5)],
[(5,6.5),(5,2)],
[(2.5,5),(7,3)],
[(2.5,5),(1.5,1)],
[(2.5,5),(4,2.5)],
[(2.5,5),(2,2)],
[(2.5,5),(8,4)],
[(2.5,5),(3,9)],
[(2.5,5),(4.5,0.5)],
[(2.5,5),(5,2)],
[(7,3),(1.5,1)],
[(7,3),(4,2.5)],
[(7,3),(2,2)],
[(7,3),(8,4)],
[(7,3),(3,9)],
[(7,3),(4.5,0.5)],
[(7,3),(5,2)],
[(1.5,1),(4,2.5)],
[(1.5,1),(2,2)],
[(1.5,1),(8,4)],
[(1.5,1),(3,9)],
[(1.5,1),(4.5,0.5)],
[(1.5,1),(5,2)],
[(4,2.5),(2,2)],
[(4,2.5),(8,4)],
[(4,2.5),(3,9)],
[(4,2.5),(4.5,0.5)],
[(4,2.5),(5,2)],
[(2,2),(8,4)],
[(2,2),(3,9)],
[(2,2),(4.5,0.5)],
[(2,2),(5,2)],
[(8,4),(3,9)],
[(8,4),(4.5,0.5)],
[(8,4),(5,2)],
[(3,9),(4.5,0.5)],
[(3,9),(5,2)],
[(4.5,0.5),(5,2)]];

This is the PPL code for mstTree(points, edges)

Code:

#cas
mstValid(edges, edge) := 
  contains(edges,edge) OR 
  contains(edges,reverse(edge));

mstTree(points,edges):=
BEGIN
LOCAL vertices, tree, j, k, p, v, d1, d2, pointsd, n, temp;
vertices, tree := {points[1]}, {}
n, pointsd := 1, SIZE(points);  // n = SIZE(vertices)
WHILE n < pointsd DO
  d1, temp := inf, (inf, inf)
  FOR j FROM 1 TO n DO
    v := vertices[j];
    FOR k FROM 1 TO pointsd DO
      p := points[k];
      IF contains(vertices, p) THEN continue END;      
      IF (d2 := abs(p-v)) >= d1 THEN continue END;
      IF NOT mstValid(edges, [p,v]) THEN continue END;
      d1, temp := d2, [p,v]
    END
  END
  tree[n] := temp;
  vertices[n+=1] := temp[1];
END
return tree;
END;
#end

Input pts, and edgs into HP Prime emulator, we get

CAS> r := mstTree(pts,edgs)

\(\left(\begin{array}{cc}
2.5+5.0*i & 3.0+4.5*i \\
1.0+5.0*i & 2.5+5.0*i \\
4.0+2.5*i & 3.0+4.5*i \\
5.0+2.0*i & 4.0+2.5*i \\
4.5+0.5*i & 5.0+2.0*i \\
2.0+2.0*i & 4.0+2.5*i \\
1.5+i & 2.0+2.0*i \\
7.0+3.0*i & 5.0+2.0*i \\
8.0+4.0*i & 7.0+3.0*i \\
5.0+6.5*i & 3.0+4.5*i \\
3.0+9.0*i & 5.0+6.5*i
\end{array}\right) \)

This matched the results of Lua code.
When plotted line segments, it matched your supplied image.

CAS> apply((x)->segment(x[0],x[1]),r)

Perhaps we should let the code make the valid edges ?
This way mstValid(edges, [p,v]) can be removed, with boost in speed.

[robmio]

(10-28-2021 04:39 PM)Albert Chan Wrote:  It is more natural to use complex number to represent points: (x,y) ≡ x + y*i
List of valid edges has each edge with 2 points. (to make a segment)

With your example, we have:

Code:

pts := [ /* total 12 points */
(3,4.5),(1,5),(5,6.5),(2.5,5),(7,3),(1.5,1),
(4,2.5),(2,2),(8,4),(3,9),(4.5,0.5),(5,2)];

edgs := [ /* total 12*11/2 = 66 edges */
[(3,4.5),(1,5)],
[(3,4.5),(5,6.5)],
[(3,4.5),(2.5,5)],
[(3,4.5),(7,3)],
[(3,4.5),(1.5,1)],
[(3,4.5),(4,2.5)],
[(3,4.5),(2,2)],
[(3,4.5),(8,4)],
[(3,4.5),(3,9)],
[(3,4.5),(4.5,0.5)],
[(3,4.5),(5,2)],
[(1,5),(5,6.5)],
[(1,5),(2.5,5)],
[(1,5),(7,3)],
[(1,5),(1.5,1)],
[(1,5),(4,2.5)],
[(1,5),(2,2)],
[(1,5),(8,4)],
[(1,5),(3,9)],
[(1,5),(4.5,0.5)],
[(1,5),(5,2)],
[(5,6.5),(2.5,5)],
[(5,6.5),(7,3)],
[(5,6.5),(1.5,1)],
[(5,6.5),(4,2.5)],
[(5,6.5),(2,2)],
[(5,6.5),(8,4)],
[(5,6.5),(3,9)],
[(5,6.5),(4.5,0.5)],
[(5,6.5),(5,2)],
[(2.5,5),(7,3)],
[(2.5,5),(1.5,1)],
[(2.5,5),(4,2.5)],
[(2.5,5),(2,2)],
[(2.5,5),(8,4)],
[(2.5,5),(3,9)],
[(2.5,5),(4.5,0.5)],
[(2.5,5),(5,2)],
[(7,3),(1.5,1)],
[(7,3),(4,2.5)],
[(7,3),(2,2)],
[(7,3),(8,4)],
[(7,3),(3,9)],
[(7,3),(4.5,0.5)],
[(7,3),(5,2)],
[(1.5,1),(4,2.5)],
[(1.5,1),(2,2)],
[(1.5,1),(8,4)],
[(1.5,1),(3,9)],
[(1.5,1),(4.5,0.5)],
[(1.5,1),(5,2)],
[(4,2.5),(2,2)],
[(4,2.5),(8,4)],
[(4,2.5),(3,9)],
[(4,2.5),(4.5,0.5)],
[(4,2.5),(5,2)],
[(2,2),(8,4)],
[(2,2),(3,9)],
[(2,2),(4.5,0.5)],
[(2,2),(5,2)],
[(8,4),(3,9)],
[(8,4),(4.5,0.5)],
[(8,4),(5,2)],
[(3,9),(4.5,0.5)],
[(3,9),(5,2)],
[(4.5,0.5),(5,2)]];

This is the PPL code for mstTree(points, edges)

Code:

#cas
mstValid(edges, edge) := 
  contains(edges,edge) OR 
  contains(edges,reverse(edge));

mstTree(points,edges):=
BEGIN
LOCAL vertices, tree, j, k, p, v, d1, d2, pointsd, n, temp;
vertices, tree := {points[1]}, {}
n, pointsd := 1, SIZE(points);  // n = SIZE(vertices)
WHILE n < pointsd DO
  d1, temp := inf, (inf, inf)
  FOR j FROM 1 TO n DO
    v := vertices[j];
    FOR k FROM 1 TO pointsd DO
      p := points[k];
      IF contains(vertices, p) THEN continue END;      
      IF (d2 := abs(p-v)) >= d1 THEN continue END;
      IF NOT mstValid(edges, [p,v]) THEN continue END;
      d1, temp := d2, [p,v]
    END
  END
  tree[n] := temp;
  vertices[n+=1] := temp[1];
END
return tree;
END;
#end

Input pts, and edgs into HP Prime emulator, we get

CAS> r := mstTree(pts,edgs)

\(\left(\begin{array}{cc}
2.5+5.0*i & 3.0+4.5*i \\
1.0+5.0*i & 2.5+5.0*i \\
4.0+2.5*i & 3.0+4.5*i \\
5.0+2.0*i & 4.0+2.5*i \\
4.5+0.5*i & 5.0+2.0*i \\
2.0+2.0*i & 4.0+2.5*i \\
1.5+i & 2.0+2.0*i \\
7.0+3.0*i & 5.0+2.0*i \\
8.0+4.0*i & 7.0+3.0*i \\
5.0+6.5*i & 3.0+4.5*i \\
3.0+9.0*i & 5.0+6.5*i
\end{array}\right) \)

This matched the results of Lua code.
When plotted line segments, it matched your supplied image.

CAS> apply((x)->segment(x[0],x[1]),r)

Perhaps we should let the code make the valid edges ?
This way mstValid(edges, [p,v]) can be removed, with boost in speed.

Congratulations, your solution is truly "elegant". However, I modified my code (see below), modifying the "same" statement, and the result obtained is like the one you proposed just now. I will save your code to compare it with my code.
Thanks a lot, Roberto.
Right code:

Code:


same(x1,y1,x2,y2,edg):=
BEGIN
FOR jj FROM 1 TO length(edg) DO
IF (x1≠edg(jj)(1) AND y1≠edg(jj)(2) AND x2≠edg(jj)(3) AND y2≠edg(jj)(4)) OR (x2≠edg(jj)(1) AND y2≠edg(jj)(2) AND x1≠edg(jj)(3) AND y1≠edg(jj)(4)) <----
THEN
RETURN true;
END;
END;
RETURN false;
END;

Wrong code:

Code:


same(x1,y1,x2,y2,edg):=
BEGIN
FOR jj FROM 1 TO length(edg) DO
IF (x1==edg(jj)(1) AND y1==edg(jj)(2) AND x2==edg(jj)(3) AND y2==edg(jj)(4)) OR (x2==edg(jj)(1) AND y2==edg(jj)(2) AND x1==edg(jj)(3) AND y1==edg(jj)(4)) <----
THEN
RETURN true;
END;
END;
RETURN false;
END;

[Albert Chan]

This version shuffled the list of points to vertices, and reduce loop counts by half.

Code:

#cas
mstValid(edges, edge) := 
  contains(edges,edge) OR 
  contains(edges,reverse(edge));

mstTree(points,edges):=
BEGIN
LOCAL tree, j, k, p, v, d1, d2, pointsd, n, temp;
tree := {}
n, pointsd := 1, SIZE(points);  // n = SIZE(vertices)
WHILE n < pointsd DO
  d1, temp := inf, (inf, inf)
  FOR j FROM 1 TO n DO
    v := points[j];             // filled vertices
    FOR k FROM n+1 TO pointsd DO
      p := points[k];
      IF (d2 := abs(p-v)) >= d1 THEN continue END;
      IF NOT mstValid(edges, [p,v]) THEN continue END;
      d1, temp := d2, [k,p,v]
    END
  END
  [k,p,v] := temp
  tree[n] := [p,v];
  points[k] := points[n+=1];    // make space
  points[n] := p;               // save vertice
END
return tree;
END;
#end

[robmio]

(10-28-2021 05:31 PM)Albert Chan Wrote:  This version shuffled the list of points to vertices, and reduce loop counts by half.

Code:

#cas
mstValid(edges, edge) := 
  contains(edges,edge) OR 
  contains(edges,reverse(edge));

mstTree(points,edges):=
BEGIN
LOCAL tree, j, k, p, v, d1, d2, pointsd, n, temp;
tree := {}
n, pointsd := 1, SIZE(points);  // n = SIZE(vertices)
WHILE n < pointsd DO
  d1, temp := inf, (inf, inf)
  FOR j FROM 1 TO n DO
    v := points[j];             // filled vertices
    FOR k FROM n+1 TO pointsd DO
      p := points[k];
      IF (d2 := abs(p-v)) >= d1 THEN continue END;
      IF NOT mstValid(edges, [p,v]) THEN continue END;
      d1, temp := d2, [k,p,v]
    END
  END
  [k,p,v] := temp
  tree[n] := [p,v];
  points[k] := points[n+=1];    // make space
  points[n] := p;               // save vertice
END
return tree;
END;
#end

Very good!
I ask a question: what if the points in the sample are not represented by a pair of numbers, but by three or four numbers, as in multivariate statistics? For instance:
points: = {{15,17,24,14}, {17,15,32,26}, {15,14,29,23}, {13,12,10,16}, {20,17,26,28}, {15,21,26,21}}.

edges:={{15,17,24,14,17,15,32,26,15,17,24,14,15,14,29,23,15,17,24,14,13,12,10,16​,15,17,24,14,20,17,26,28,15,17,24,14,15,21,26,21},{17,15,32,26,15,14,29,23,17,15​,32,26,13,12,10,16,17,15,32,26,20,17,26,28,17,15,32,26,15,21,26,21},{15,14,29,23​,13,12,10,16,15,14,29,23,20,17,26,28,15,14,29,23,15,21,26,21},{13,12,10,16,20,17​,26,28,13,12,10,16,15,21,26,21},{20,17,26,28,15,21,26,21}}.

With these points and edges, "my" code results in:

[[17,15,15,17],[15,14,17,15],[13,12,15,14],[20,17,17,15],[15,21,15,17]].

With these points and edges, what does your code provide, Albert?

[Albert Chan]

(10-28-2021 05:04 PM)robmio Wrote:  However, I modified my code (see below), modifying the "same" statement ...

No ! OP post were correct ! (except for minor leaking of variables, jj in same())

It was your inputs that were wrong.
Each point has 2 numbers, Edge have 2 points, thus 4 numbers.
Remember, each edge contains 4 values.

pts0 := [ [3,4.5], [1,5], ..., ,[5,2]]; /* 12 points */
edgs0 := [ [3,4.5 , 1,5], [3,4.5 , 5,6.5], ..., [4.5,0.5 , 5,2]]; /* 66 edges */

Running OP code, on my laptop, it finished in 0.273 sec, return these line segments.
(for comparison, my complex-number as points 2nd version, it finished in 0.006 sec)

CAS> mstRob(pts0, edgs0)
{{2.5,5,3,4.5},
{1,5,2.5,5},
{4,2.5,3,4.5},
{5,2,4,2.5},
{4.5,0.5,5,2},
{2,2,4,2.5},
{1.5,1,2,2},
{7,3,5,2},
{8,4,7,3},
{5,6.5,3,4.5},
{3,9,5,6.5}}

[Albert Chan]

Version 3, extending to multi-dimensional points.

List of possible edges may be huge, this version switched the logic, and disallowed "bad" edges.
Example, if edge [p1,p2] (and of course, [p2, p1]) is not allowed: badEdges = [[p1,p2], [p2,p1]]
Note that badEdges list length always even.

For comparing "distance", it uses norm2, squared, of difference of 2 points.

Without doing square roots, this version is even faster than my previous version.
Running same example, 12 points, 66 edges (badEdges = []) , it finished in 0.004 sec.

Code:

#cas  
mstTree(points, badEdges):=
BEGIN
LOCAL tree, j, k, p, v, d1, d2, pointsd, n, temp;
tree := {}
n, pointsd := 1, len(points);   // n = len(vertices)
WHILE n < pointsd DO
  d1, temp := inf, (inf, inf)
  FOR j FROM 1 TO n DO
    v := points[j];             // filled vertices
    FOR k FROM n+1 TO pointsd DO
      p := points[k];
      d2 = p .- v;
      IF (d2 *= d2) >= d1 THEN continue END;
      IF member([k,j], badEdges) THEN continue END;
      d1, temp := d2, [k,p,v];
    END
  END
  [k,p,v] := temp;
  tree[n] := [p,v];
  points[k] := points[n+=1];    // make space
  points[n] := p;               // save vertice
END
return tree;
END;
#end

With your example, 6 points, 6*5/2=15 edges (again, badEdges = []), this is the result:

CAS> pts := [[15,17,24,14], [17,15,32,26], [15,14,29,23], [13,12,10,16], [20,17,26,28], [15,21,26,21]]
CAS> mstTree(pts, [])

{[[15,21,26,21] , [15,17,24,14]],
  [[15,14,29,23] , [15,21,26,21]],
  [[17,15,32,26] , [15,14,29,23]],
  [[20,17,26,28] , [17,15,32,26]],
  [[13,12,10,16] , [15,17,24,14]]}

sum of line segments length = √69 + √62 + √23 + √53 + √229 ≈ 43.3893

[robmio]

(10-28-2021 09:07 PM)Albert Chan Wrote:  Version 3, extending to multi-dimensional points.

List of possible edges may be huge, this version switched the logic, and disallowed "bad" edges.
Example, if edge [p1,p2] (and of course, [p2, p1]) is not allowed: badEdges = [[p1,p2], [p2,p1]]
Note that badEdges list length always even.

For comparing "distance", it uses norm2, squared, of difference of 2 points.

Without doing square roots, this version is even faster than my previous version.
Running same example, 12 points, 66 edges (badEdges = []) , it finished in 0.004 sec.

Code:

#cas  
mstTree(points, badEdges):=
BEGIN
LOCAL tree, j, k, p, v, d1, d2, pointsd, n, temp;
tree := {}
n, pointsd := 1, len(points);   // n = len(vertices)
WHILE n < pointsd DO
  d1, temp := inf, (inf, inf)
  FOR j FROM 1 TO n DO
    v := points[j];             // filled vertices
    FOR k FROM n+1 TO pointsd DO
      p := points[k];
      d2 = p .- v;
      IF (d2 *= d2) >= d1 THEN continue END;
      IF member([k,j], badEdges) THEN continue END;
      d1, temp := d2, [k,p,v];
    END
  END
  [k,p,v] := temp;
  tree[n] := [p,v];
  points[k] := points[n+=1];    // make space
  points[n] := p;               // save vertice
END
return tree;
END;
#end

With your example, 6 points, 6*5/2=15 edges (again, badEdges = []), this is the result:

CAS> pts := [[15,17,24,14], [17,15,32,26], [15,14,29,23], [13,12,10,16], [20,17,26,28], [15,21,26,21]]
CAS> mstTree(pts, [])

{[[15,21,26,21] , [15,17,24,14]],
  [[15,14,29,23] , [15,21,26,21]],
  [[17,15,32,26] , [15,14,29,23]],
  [[20,17,26,28] , [17,15,32,26]],
  [[13,12,10,16] , [15,17,24,14]]}

sum of line segments length = √69 + √62 + √23 + √53 + √229 ≈ 43.3893

I understood: my mistake was to propose the edges to the code in this way:

{{3,4.5,1,5,3,4.5,5,6.5,3,4.5,2.5,5,3,4.5,7,3,3,4.5,1.5,1,3,4.5,4,2.5,3,4.5,2,2,​3,4.5,8,4,3,4.5,3,9,3,4.5,4.5,0.5,3,4.5,5,2},{1,5,5,6.5,1,5,2.5,5,1,5,7,3,1,5,1.​5,1,1,5,4,2.5,1,5,2,2,1,5,8,4,1,5,3,9,1,5,4.5,0.5,1,5,5,2},{5,6.5,2.5,5,5,6.5,7,​3,5,6.5,1.5,1,5,6.5,4,2.5,5,6.5,2,2,5,6.5,8,4,5,6.5,3,9,5,6.5,4.5,0.5,5,6.5,5,2}​,{2.5,5,7,3,2.5,5,1.5,1,2.5,5,4,2.5,2.5,5,2,2,2.5,5,8,4,2.5,5,3,9,2.5,5,4.5,0.5,​2.5,5,5,2},{7,3,1.5,1,7,3,4,2.5,7,3,2,2,7,3,8,4,7,3,3,9,7,3,4.5,0.5,7,3,5,2},{1.​5,1,4,2.5,1.5,1,2,2,1.5,1,8,4,1.5,1,3,9,1.5,1,4.5,0.5,1.5,1,5,2},{4,2.5,2,2,4,2.​5,8,4,4,2.5,3,9,4,2.5,4.5,0.5,4,2.5,5,2},{2,2,8,4,2,2,3,9,2,2,4.5,0.5,2,2,5,2},{​8,4,3,9,8,4,4.5,0.5,8,4,5,2},{3,9,4.5,0.5,3,9,5,2},{4.5,0.5,5,2}}

without making groups of 4 numbers into 4 numbers.
Thank you for having brilliantly solved the "multivariate" version as well. Are you a computer scientist?

Best wishes, Roberto.

[Albert Chan]

(10-29-2021 05:53 AM)robmio Wrote:  Thank you for having brilliantly solved the "multivariate" version as well. Are you a computer scientist?

No. I had never heard of minimum spanning tree

(10-28-2021 09:07 PM)Albert Chan Wrote:  Version 3, extending to multi-dimensional points
...

CAS> pts := [[15,17,24,14], [17,15,32,26], [15,14,29,23], [13,12,10,16], [20,17,26,28], [15,21,26,21]]
CAS> mstTree(pts, [])

{[[15,21,26,21] , [15,17,24,14]],
  [[15,14,29,23] , [15,21,26,21]],
  [[17,15,32,26] , [15,14,29,23]],
  [[20,17,26,28] , [17,15,32,26]],
  [[13,12,10,16] , [15,17,24,14]]}

sum of line segments length = √69 + √62 + √23 + √53 + √229 ≈ 43.3893

Instead of badEdges = [[p1,p2],[p2,p1], ...], we could simplify to [[1,2],[2,1], ...]
Same for tree = [[p(j1),p(k1)], [p(j2),p(k2)], ...], simplified to [[j1,k1],[j2,k2], ...]

But, this required pts not be shuffled (otherwise indexes meant nothing)
Version 4 solve the problems with another level of indirection (*)

Code:

#cas
mstTree(points, badEdges):=
BEGIN
LOCAL tree, j, k, v, d1, d2, pointsd, n, temp;
LOCAL m, mj, mk;                // mapped indexes
n, pointsd := 1, len(points);   // n = len(vertices)
tree, m := {}, range(1,pointsd+1);
WHILE n < pointsd DO
  d1 := temp := inf;
  FOR j FROM 1 TO n DO
    v := points[mj := m[j]];    // filled vertices
    FOR k FROM n+1 TO pointsd DO
      d2 = v .- points[mk := m[k]];
      IF (d2 *= d2) >= d1 THEN continue END;            
      IF member([mk,mj], badEdges) THEN continue END;
      d1, temp := d2, [k, [mk,mj]];
    END
  END
  k, tree[n] := temp;
  m[k] := m[n += 1];            // swap m[k], m[n+1]
  m[n] := temp[2][1];
END
return tree;
END;
#end

Redo the same example, with this new version

CAS> pts := [[15,17,24,14], [17,15,32,26], [15,14,29,23], [13,12,10,16], [20,17,26,28], [15,21,26,21]]
CAS> tree := mstTree(pts, [])

[[6,1],[3,6],[2,3],[5,2],[4,1]]

CAS> norm2(x) := sqrt(dot(x,x))
CAS> map(tree, x -> norm2(pts[x[1]] - pts[x[2]]))       → [√69,√62,√23,√53,√229]
CAS> sum(float(Ans))       → 43.3893190999

What if points #2, #3 line segment not allowed ?

CAS> tree := mstTree(pts, [[2,3],[3,2]])

[[6,1],[3,6],[5,3],[2,5],[4,1]]

CAS> map(tree, x -> norm2(pts[x[1]] - pts[x[2]]))       → [√69,√62,2*√17,√53,√229]
CAS> sum(float(Ans))       → 46.8396988279

(*)

[robmio]

Perfect, the substitution (in the output) of the coordinates with the numbers (referring to each node), will allow me to program the Friedman-Rafsky test.

Again thank you for solving the problem brilliantly.
Best regards, Roberto.

[robmio]

I address the question to Albert Chan: using the last algorithm you wrote, for a matrix of dimensions {100,4} the "physical" calculator PRIME_G2 takes about 13 seconds to calculate the "MST". Does your algorithm already include the implementation with the Fibonacci heap?

[Albert Chan]

Version 3 assumed all edges vaild, unless specified otherwise.
I think the speed is due to not spend time to validate edges.

By shuffling mapped indexes, it keep the points to examine low.
Points to examine reduced to about 1/3 (and, no code to validate vertices)

Note: test for member in list has O(n)

---

Points examined, for oringinal Lua code, n = total points, k = vertices examined

sum(n*k, k=1..n-1)
= n * n*(n-1)/2
= n^3/2 - n^2/2

Points examined, with shuffled indexes, splitting off vertices and non-vertices.

sum((n-k)*k, k=1..n-1)
= sum((n-1)*k - k*(k-1), k=1..n-1)       // see Funny Factorials and Slick Sums
= (n-1) * n*(n-1)/2 - n*(n-1)*(n-2)/3
= n*(n-1)*(n+1)/6
= n^3/6 - n/6

[robmio]

(11-02-2021 11:57 AM)Albert Chan Wrote:  Version 3 assumed all edges vaild, unless specified otherwise.
I think the speed is due to not spend time to validate edges.

By shuffling mapped indexes, it keep the points to examine low.
Points to examine reduced to about 1/3 (and, no code to validate vertices)

Note: test for member in list has O(n)

---

Points examined, for oringinal Lua code, n = total points, k = vertices examined

sum(n*k, k=1..n-1)
= n * n*(n-1)/2
= n^3/2 - n^2/2

Points examined, with shuffled indexes, splitting off vertices and non-vertices.

sum((n-k)*k, k=1..n-1)
= sum((n-1)*k - k*(k-1), k=1..n-1)       // see Funny Factorials and Slick Sums
= (n-1) * n*(n-1)/2 - n*(n-1)*(n-2)/3
= n*(n-1)*(n+1)/6
= n^3/6 - n/6

Good. I'll try to translate your algorithm into Python, to see if it gets even faster. Thanks for your help.

[Albert Chan]

(10-29-2021 10:52 AM)Albert Chan Wrote:  Instead of badEdges = [[p1,p2],[p2,p1], ...], we could simplify to [[1,2],[2,1], ...]

I thought it is trivial to supply 1 sided badEdges, then "doubled" them.
For python, joining two list is simply '+'

>>> bad = [[1,2],[3,6],[6,10]]
>>> bad2 = [ [k[1],k[0]] for k in bad]
>>> bad + bad2
[[1, 2], [3, 6], [6, 10], [2, 1], [6, 3], [10, 6]]

For HP Prime, close equivalent for '+' is extend.
But, instead of extending to end of list, it extend to items inside list.

Cas> bad := [[1,2],[3,6],[6,10]]
Cas> bad2 := map(bad, reverse)       → [[2,1],[6,3],[10,6]]
Cas> extend(bad, bad2)                     → [[1,2,2,1],[3,6,6,3],[6,10,10,6]]

To get what we wanted, we need to use transpose (3 times !)
This work, but ugly. Better solution welcome.

Cas> transpose(extend(transpose(bad), transpose(bad2))
[[1,2],[3,6],[6,10],[2,1],[6,3],[10,6]]


Version 5:

Code:

#cas
mstTree(points, badEdges):=
BEGIN
LOCAL tree, j, k, v, d1, d2, pointsd, n, temp;
LOCAL m, mj, mk;                // mapped indexes
n, pointsd := 1, len(points);   // n = len(vertices)
tree, m := {}, range(1,pointsd+1);
WHILE n < pointsd DO
  d1 := temp := inf;
  FOR j FROM 1 TO n DO
    v := points[mj := m[j]];    // filled vertices
    FOR k FROM n+1 TO pointsd DO
      d2 = v .- points[mk := m[k]];
      IF (d2 *= d2) >= d1 THEN continue END;            
      IF member([mj,mk], badEdges) THEN continue END;
      d1, temp := d2, [k, [mj,mk]];
    END
  END
  k, tree[n] =< temp;
  m[k] =< m[n += 1];            // swap m[k], m[n+1]
  m[n] =< temp[2][2];
END
return tree;
END;
#end

Since it is more natural to read from left to right, I flip the order.
Also, HP Prime always create a new list, even for just updating element.
We might as well update list element by reference, '=<' instead of ':='

OP example:

Cas> pts := [[3,4.5],[1,5],[5,6.5],[2.5,5],[7,3],[1.5,1],[4,2.5],[2,2],[8,4],[3,9],[4.5,0.5],[5,2]]
Cas> mstTree(pts, [])
[[1,4],[4,2],[1,7],[7,12],[12,11],[7,8],[8,6],[12,5],[5,9],[1,3],[3,10]]

It is easier to sketch the tree.

3 10
|
1 4 2
|
7 12 11
| |
| 5 9
|
8 6

[robmio]

(11-04-2021 03:47 PM)Albert Chan Wrote:  

(10-29-2021 10:52 AM)Albert Chan Wrote:  Instead of badEdges = [[p1,p2],[p2,p1], ...], we could simplify to [[1,2],[2,1], ...]

I thought it is trivial to supply 1 sided badEdges, then "doubled" them.
For python, joining two list is simply '+'

>>> bad = [[1,2],[3,6],[6,10]]
>>> bad2 = [ [k[1],k[0]] for k in bad]
>>> bad + bad2
[[1, 2], [3, 6], [6, 10], [2, 1], [6, 3], [10, 6]]

For HP Prime, close equivalent for '+' is extend.
But, instead of extending to end of list, it extend to items inside list.

Cas> bad := [[1,2],[3,6],[6,10]]
Cas> bad2 := map(bad, revserse)       → [[2,1],[6,3],[10,6]]
Cas> extend(bad, bad2)                     → [[1,2,2,1],[3,6,6,3],[6,10,10,6]]

To get what we wanted, we need to use transpose (3 times !)
This work, but ugly. Better solution welcome.

Cas> transpose(extend(transpose(bad), transpose(bad2))
[[1,2],[3,6],[6,10],[2,1],[6,3],[10,6]]


Version 5:

Code:

#cas
mstTree(points, badEdges):=
BEGIN
LOCAL tree, j, k, v, d1, d2, pointsd, n, temp;
LOCAL m, mj, mk;                // mapped indexes
n, pointsd := 1, len(points);   // n = len(vertices)
tree, m := {}, range(1,pointsd+1);
WHILE n < pointsd DO
  d1 := temp := inf;
  FOR j FROM 1 TO n DO
    v := points[mj := m[j]];    // filled vertices
    FOR k FROM n+1 TO pointsd DO
      d2 = v .- points[mk := m[k]];
      IF (d2 *= d2) >= d1 THEN continue END;            
      IF member([mj,mk], badEdges) THEN continue END;
      d1, temp := d2, [k, [mj,mk]];
    END
  END
  k, tree[n] =< temp;
  m[k] =< m[n += 1];            // swap m[k], m[n+1]
  m[n] =< temp[2][2];
END
return tree;
END;
#end

Since it is more natural to read from left to right, I flip the order.
Also, HP Prime always create a new list, even for just updating element.
We might as well update list element by reference, '=<' instead of ':='

OP example:

Cas> pts := [[3,4.5],[1,5],[5,6.5],[2.5,5],[7,3],[1.5,1],[4,2.5],[2,2],[8,4],[3,9],[4.5,0.5],[5,2]]
Cas> mstTree(pts, [])
[[1,4],[4,2],[1,7],[7,12],[12,11],[7,8],[8,6],[12,5],[5,9],[1,3],[3,10]]

It is easier to sketch the tree.

3 10
|
1 4 2
|
7 12 11
| |
| 5 9
|
8 6

Indeed, it is easier to view the entire tree

[Albert Chan]

It seems XCas already has "safe" extend built-in, called semi-augment

semi_augment(A,B) := when(len(A)!=len(B), extend(A,B), extend(append(A,head(B)),tail(B)))

Cas> semi_augment([[1,2],[3,6],[6,10]] , [[2,1],[6,3],[10,6]])

[[1,2], [3,6], [6,10], [2,1], [6,3], [10,6]]       // 3+3 = 6 items

Cas> semi_augment([[1,2]] , [[3,4,5]])

[[1,2], [3,4,5]]       // 1 + 1 = 2 items


Had we use extend instead, we get these:

Cas> extend([[1,2],[3,6],[6,10]] , [[2,1],[6,3],[10,6]])

[[1,2,2,1], [3,6,6,3], [6,10,10,6]]

Cas> extend([[1,2]] , [[3,4,5]])

[[1,2,3,4,5]]