7
$\begingroup$

There are 1000 points in the space. The coordinates are (0,0,0) - (9,9,9). Any point can cover all points on the three straight lines (parallel to the X, Y and Z axes respectively) where this point is located.

Find the minimum number of points n that can cover 1000 points in the space. How can I use Mathematica to quickly find the coordinates of these n(n=50?) points?

Enter image description here

Ben Izd's answer

In the case of 6x6x6(n=18). we can verify that the result is correct. The result of 6x6x6 is obtained by running Ben Izd's algorithm with Intel Core i5-4310M for about 73 minutes. However, the 5x5x5(n=13) case took only 30 seconds.

strip[{x_, y_, z_}, tuples_] := 
 Cases[tuples, Except[{x, y, _} | {x, _, z} | {_, y, z}]]
keyp = {{0, 4, 0}, {4, 0, 1}, {2, 0, 2}, {0, 1, 3}, {0, 5, 4}, {5, 0, 
    5}, {1, 1, 0}, {3, 1, 4}, {2, 2, 1}, {5, 2, 2}, {4, 2, 5}, {5, 3, 
    1}, {4, 3, 2}, {2, 3, 5}, {3, 4, 3}, {1, 4, 4}, {3, 5, 0}, {1, 5, 
    3}};
pts[0] = Tuples[Range[0, 5, 1], 3];
pts[1] = strip[keyp[[1]], pts[0]];
pts[i_] := pts[i] = strip[keyp[[i]], pts[i - 1]]
len = Map[Length, Table[pts[i], {i, 0, 18, 1}]]
Sort[-1*Differences[len], Greater]
(*out1: {216, 200, 184, 170, 156, 144, 132, 119, 108, 95, 83, 72, 62, 52, 42, 30, 20, 10, 0}
out2: {16, 16, 14, 14, 13, 13, 12, 12, 12, 12, 11, 11, 10, 10, 10, 10, 10, 10} *)

In fact, eliminating the most points each time is not the optimal algorithm.

Improve this question
$\endgroup$
10
  • 1
    $\begingroup$ Please explicitly say what problem you encountered in your code. You posted quite a bit of code, probably related to your question, but it is not clear what the code does and where it goes wrong. It would be easier to help if you pinpointed the problem at least. $\endgroup$
    – MarcoB
    Commented May 12, 2022 at 12:46
  • 1
    $\begingroup$ Apparently the OP has no Mathematica related difficulty. Consequently, I am casting a close vote. $\endgroup$
    – Syed
    Commented May 12, 2022 at 13:07
  • $\begingroup$ It's a little hard to understand initially. Maybe: "What is the quickest way to find find the minimum number of points such that their associated orthogonal axes cover all points in the space." Can you first develop a routine that counts all the points on the associated axes given a point (x,y,z) in the space? That would be a start. $\endgroup$
    – josh
    Commented May 12, 2022 at 13:31
  • $\begingroup$ The RegionIntersection function may give coordinates of the relevant points $\endgroup$
    – miss
    Commented May 12, 2022 at 13:50
  • 6
    $\begingroup$ What is the definition of "cover"? Is it like the Rooks Problem, math.stackexchange.com/questions/803276/rooks-in-3d-chess-board ? Seems like a mathematics question, not a Mathematica question, though. If you have trouble implementing an algorithm, describe the algorithm. $\endgroup$
    – Michael E2
    Commented May 12, 2022 at 17:39

3 Answers 3

Reset to default
11
$\begingroup$

Here the problem was solved with a different perspective. You may find other solutions but what made this interesting to me was transforming one problem into another.

Here are the steps:

  1. Create random points in 3d space
  2. For each point, find points it covers
  3. Build a graph out of these connections
  4. Find a minimum set of points which can cover all the vertices by their neighborhood (aka cover all the pairs)

Assume we have a set of pairs:

SeedRandom[75];

pairs = RandomInteger[{0, 3}, {20, 3}];

n = Length[pairs];

pairs that have at least two similar axes will be associated for each pair, we can find them by:

Position[pairs, 
 Alternatives @@ 
  Table[ReplacePart[SOMEPAIR, i -> _], {i, 3}], {1}]

If we extend this code further, we could build an AdjacencyMatrix in which we want to select a minimum set of nodes that we can reach all the vertices with at most 1 step(like FindIndependentVertexSet but in minimum terms and little refinement, vertices can be two sides of an edge):

graphRules = 
  DeleteDuplicatesBy[Sort]@
   Catenate@
    Table[Thread[
      index \[UndirectedEdge] 
       Catenate@
        Position[pairs, 
         Alternatives @@ 
          Table[ReplacePart[pairs[[index]], i -> _], {i, 
            3}], {1}]], {index, n}];

adjacency = Unitize@AdjacencyMatrix[graphRules];

Showing the graph:

Graph[graphRules]

enter image description here

Now we find the vertex set:

result = LinearOptimization[
  ConstantArray[1,n],
 {Join[DiagonalMatrix@ConstantArray[1, n], adjacency], 
   Join[ConstantArray[0, n], ConstantArray[-1, n]]}, Integers]

Highlighting the vertex set:

HighlightGraph[graphRules, 
 VertexList[graphRules][[Catenate@Position[result, 1]]]]

enter image description here

The red vertices in the above graph are the pairs you're looking for. Each vertex in the graph is either red or is in the neighborhood of a red vertex (aka is a chosen pair or is covered by a chosen pair). If we get the pairs:

pairs[[VertexList[graphRules][[Catenate@Position[result, 1]]]]]

(* Out: {{2, 3, 0}, {1, 2, 0}, {3, 0, 2}, {2, 1, 1}, {0, 1, 2}, {1, 1, 3}} *)

Visualizing the pairs:

enter image description here

Update 1

Since the above solution for large pairs like (Tuples[Range[0, 6, 1], 3]) takes quite a long time (mainly because of LinearOptimization), it could be faster in a static language. Here we'll discuss an alternative in Julia that is much faster than Mathematica in this particular case.

Requirement:

Now, start a Julia session:

juliaSession = 
 StartExternalSession[<|"System" -> "Julia", 
   "SessionProlog" -> "using JuMP, HiGHS"|>]

Define a function that accepts adjacency as input and returns the result:

juliaLinearOptimizer = ExternalFunction[juliaSession,
"function temp(data)
    
    n=size(data)[1];

    model = Model();
    set_optimizer(model,HiGHS.Optimizer);

    @variable(model,v[1:n],Bin)

    @objective(model,Min,sum(v[1:n]))

    for row in eachrow(data)
        @constraint(model,sum(v.*row)>=1)
    end

    optimize!(model)

    return round.(value.(v))
end"]

It solved the 6x6x6 in 40 seconds!

Floor @ Flatten @ juliaLinearOptimizer[Normal[adjacency]]

Note that the first call may take a little more time. Also, don't forget to delete the object after you're done:

DeleteObject[juliaSession]

Comparison for 5x5x5

Software/Language Time (second)
Mathematica 13.0.1 (LinearOptimization) 21
Matlab 2021a (intlinprog) 11
Julia 1.7.2 (HiGHS 1.1.3 - 2nd time) 1.5

Update 2

For checking the correctness you can use:

DeleteCases[pairs, 
 Alternatives @@ (({#1, #2, _} | {#1, _, #3} | {_, #2, #3}) & @@@ chosen)]

Where pairs are all the points and chosen are the picked ones.

Improve this answer
$\endgroup$
7
  • 1
    $\begingroup$ @fcwyq don't be so bombastic. Ben Izd isn't claiming to be a graph theory expert. No one is stupid for not understanding a single piece of code. Perhaps Ben hasn't solved your problem. Perhaps his answer is still useful. $\endgroup$
    – Adam
    Commented May 12, 2022 at 18:17
  • 2
    $\begingroup$ @fcwyq: Given the points Ben creates above in pairs, and the red vertices he computes, those points: {{2, 3, 0}, {1, 2, 0}, {3, 0, 2}, {2, 1, 1}, {0, 1, 2}, {1, 1, 3}}, when extended orthogonally, do indeed cover all points in pairs. Try it. Maybe we can then extend Ben's code to create the points in your space to subsequently compute the necessary vertices. $\endgroup$
    – josh
    Commented May 12, 2022 at 18:30
  • 1
    $\begingroup$ @josh: Thank you. I understand. Great answer ! $\endgroup$
    – miss
    Commented May 12, 2022 at 18:32
  • 1
    $\begingroup$ @Ben: Found typo: need Graph[graphRules] $\endgroup$
    – josh
    Commented May 12, 2022 at 18:37
  • $\begingroup$ Hi to everyone, First of all, I'm not an expert in any field let alone graph theory. As @Adam said no one is stupid for not understanding a code, especially here (with this belief code gold community would be filled with desperate people ;), I'm responsible for not conveying the main point clearly, so I'll edit the answer. $\endgroup$
    – Ben Izd
    Commented May 12, 2022 at 18:38
7
$\begingroup$

A single point covers at most 28 points. So $\frac{1000}{28}\approx35.7<n$.

Now here's some code that selects random points from the remaining ones until it covers all 1000:

NestWhile[{Append[#, #3], 
DeleteCases[#2,Alternatives @@ Table[ReplacePart[#3, i -> _], {i, 3}]]} & @@ 
Append[#, RandomChoice@#[[2]]] &, {{}, 
Flatten[Outer[List, #, #, #] &@Range[0, 9], 2]},Length@#[[2]] > 0 &][[1]]

I ran that a few times and found a list of 78 points:

{{1,5,3},{9,5,5},{5,6,9},{4,4,5},{4,2,1},{1,1,1},
 {4,0,2},{2,9,7},{1,8,7},{9,9,9},{4,1,7},{6,5,7},
 {2,3,9},{3,1,3},{2,6,0},{3,7,2},{8,8,0},{5,2,5},
 {4,9,8},{6,7,9},{8,1,4},{6,6,1},{7,6,8},{7,7,5},
 {2,4,1},{2,7,6},{1,3,4},{9,4,2},{0,0,5},{5,7,8},
 {0,9,3},{8,5,1},{7,3,1},{1,6,6},{9,0,1},{5,9,1},
 {0,2,4},{0,4,9},{9,3,7},{8,3,6},{8,4,3},{6,4,4},
 {5,5,0},{3,5,6},{5,0,3},{0,7,0},{8,6,7},{3,6,4},
 {2,2,2},{7,1,2},{7,0,0},{9,8,6},{1,2,0},{6,3,0},
 {4,6,3},{6,1,6},{2,8,4},{3,9,0},{8,0,9},{3,8,9},
 {3,0,8},{4,5,9},{6,2,8},{5,4,6},{2,1,8},{7,2,3},
 {7,9,4},{1,9,5},{8,9,2},{3,4,7},{6,8,5},{0,3,2},
 {0,8,8},{9,7,4},{9,1,0},{5,8,2},{1,4,8},{3,3,5}}

So we can say for sure that $36\leq n\leq78$.

I ran the code for smaller grids. The pattern seems to be you can cover a $\{0\ldots m\}^3$ grid with just under $m^2$ points, like 1 or 2 or 3 fewer.

blegh

Those are coverings of $\{0,1\}^3$, $\{0,1,2\}^3$,..., $\{0,...8\}^3$. They have 2,5,8,15,25,35,47,62 points respectively. Here's the code to make them, so you can check for yourself.

MapIndexed[({l,c}\[Function]Graphics3D[{Point@#,#2},Boxed->False]&[#,Table[Line@{ReplacePart[#,i->0],ReplacePart[#,i->c]},{i,3}]&/@l])[#,First@#2]&,{{{0,0,0},{1,1,1}},{{0,0,0},{0,2,2},{2,0,2},{2,2,0},{1,1,1}},{{0,0,0},{0,3,3},{3,0,3},{3,3,0},{1,2,2},{2,1,2},{1,1,1},{2,2,1}},{{0,0,0},{0,4,4},{4,0,4},{4,4,0},{2,3,1},{3,1,1},{3,3,2},{2,2,2},{1,2,1},{1,0,3},{3,2,3},{1,1,2},{1,3,4},{2,1,3},{4,3,3}},{{0,0,0},{0,5,5},{5,0,5},{5,5,0},{2,3,3},{2,4,0},{3,3,2},{0,4,4},{2,5,1},{4,0,2},{1,5,2},{5,4,1},{5,2,4},{1,2,3},{4,3,4},{3,2,5},{3,1,3},{0,2,1},{3,0,4},{2,1,2},{4,2,0},{4,5,3},{4,4,5},{1,1,1},{1,3,5}},{{0,0,0},{0,6,6},{6,0,6},{6,6,0},{2,4,2},{2,2,5},{5,5,5},{4,4,0},{6,1,3},{0,4,1},{3,6,4},{4,5,2},{6,3,4},{3,5,1},{4,0,4},{5,3,2},{3,4,5},{3,2,0},{3,0,3},{0,5,3},{4,2,3},{1,1,0},{1,0,2},{6,2,1},{2,5,4},{1,4,6},{2,0,1},{0,1,2},{5,2,4},{2,3,6},{1,3,5},{5,6,3},{4,1,6},{4,6,1},{5,1,1}},{{0,0,0},{0,7,7},{7,0,7},{7,7,0},{3,3,1},{7,3,4},{1,1,5},{1,0,3},{2,4,4},{6,3,2},{5,6,7},{6,1,7},{4,7,2},{6,2,0},{5,2,2},{2,7,6},{7,5,6},{1,5,1},{2,6,1},{4,3,6},{4,4,1},{7,6,2},{5,1,0},{3,5,4},{5,0,5},{0,4,5},{2,3,7},{3,2,5},{6,6,3},{3,0,2},{0,1,3},{0,2,4},{2,5,2},{6,7,5},{4,5,7},{1,6,6},{1,4,0},{5,7,3},{3,6,0},{2,2,3},{7,2,1},{4,6,4},{3,4,6},{6,0,4},{7,4,3},{1,2,7},{1,7,4}},{{0,0,0},{0,8,8},{8,0,8},{8,8,0},{4,4,8},{6,3,4},{6,7,5},{3,8,3},{0,3,5},{6,1,7},{7,3,7},{6,0,3},{3,2,5},{5,1,6},{1,7,6},{3,1,1},{4,8,4},{5,5,3},{7,1,2},{8,3,3},{4,0,2},{6,2,6},{8,7,2},{2,7,7},{0,6,3},{4,6,0},{2,3,6},{2,2,1},{1,6,1},{1,4,7},{0,5,1},{8,2,7},{7,0,4},{1,1,8},{6,6,2},{7,5,6},{7,4,5},{4,7,1},{7,2,0},{1,2,4},{1,5,0},{5,7,4},{5,6,8},{4,5,5},{3,5,8},{2,4,4},{5,8,5},{3,6,6},{4,1,3},{5,3,2},{5,0,7},{6,8,1},{0,4,2},{8,1,4},{5,4,0},{8,4,6},{2,0,5},{7,7,3},{2,8,2},{8,6,5},{2,1,0},{3,7,0}}}]

This pattern supports the idea that a GLB for $n$ is something like 76.

Improve this answer
$\endgroup$
5
$\begingroup$

An alternate approach to structuring the problem. Since this is Mathematica, we're required to use patterns. I think it is possible to make it more "pattern-ish"

Build up the list of points...

tuples = Tuples[Range[10], 3];

Define a function that strips out the points covered given {x,y,z}

strip[{x_, y_, z_}, tuples_] := Cases[tuples, Except[{x, y, _} | {x, _, z} | {_, y, z}]];

Do a random reduction

reduceRandom[tuples_] := Module[{x, y, z},
                        {x, y, z} = RandomChoice[tuples];
                        strip[{x, y, z}, tuples]
                        ]

Or a continual local optimization (greedy algorithm) maximizing the number of points removed at each step.

reduceMaximal[tuples_] := Module[{x, y, z},
                  {x, y, z} = First@MinimalBy[tuples, Length@strip[#, tuples] &];
                  strip[{x, y, z}, tuples]
                  ]

Try them. Greedy algorithm first

res = Cases[FixedPointList[reduceMaximal, tuples], Except[{}]];
Length[res]
(* 76 *)

Pretty sure this is not optimal. Random efforts...

Sort@Table[Length@Cases[FixedPointList[reduceRandom, tuples], Except[{}]], 300]

(* {80, 82, 82, 82, 82, 82, 83, 83, 83, 83, 83, 83,... *)
Improve this answer
$\endgroup$

Not the answer you're looking for? Browse other questions tagged or ask your own question.