Coloring edges of a grid

This is a little coloring problem from [1]:

Given an m * n grid, each edge must be colored. However, there are 2 constraints:

1. Entire grid must use only 3 unique colors
2. Each square in the grid must be colored using exactly 2 colors (2 edges per color)

Let's try to build a GAMS model for this. The question itself is not so interesting, but the modeling for setting up these grids is.

Data structures

First, we need to represent the grid. 

 

*------------------------------------------------- * sets to represent grid and edges. * representation of squares is by left-lower point *------------------------------------------------- set    i 'grid x-coordinate'  /i1*i10/    j 'grid y-coordinate'  /j1*j10/    edges(i,j,i,j) 'connects two points'    squares(i,j) 'numbering of squares'    semap(i,j,i,j,i,j) 'map between squares and edges'    c 'colors' /red,green,blue/ ; alias (i,ii),(j,jj);   edges(i,j,i+1,j) = yes; edges(i,j,i,j+1) = yes;   squares(i-1,j-1) = yes; semap(squares(i,j),edges(i,j,ii,jj)) = yes; semap(squares(i,j),edges(ii,jj,i+1,j+1)) = yes;  option edges:0:0:8,semap:0:0:6; display i,j,edges,squares,semap;

There are many ways to do this. Here is one. Each grid point is resented by its \(x\) and \(y\) coordinate, here denoted by \((i,j)\). Edges connect two points, so we get \((i,j,i',j')\) to represent \((i,j) \rightarrow (i',j')\). We only need two edges leaving for each point: one to the right and one to above. This is done by:

edges(i,j,i+1,j) = yes;

edges(i,j,i,j+1) = yes;

When we display this (for a smaller example) we see:

----     31 SET i  grid x-coordinate

i1,    i2,    i3


----     31 SET j  grid y-coordinate

j1,    j2,    j3


----     31 SET edges  connects two points

i1.j1.i1.j2,    i1.j1.i2.j1,    i1.j2.i1.j3,    i1.j2.i2.j2,    i1.j3.i2.j3,    i2.j1.i2.j2,    i2.j1.i3.j1,    i2.j2.i2.j3
i2.j2.i3.j2,    i2.j3.i3.j3,    i3.j1.i3.j2,    i3.j2.i3.j3

I don't use a different numbering scheme for the squares but instead use the coordinates of the left-lower point. We drop the points on the right and upper edge of the grid. This is done by:

squares(i-1,j-1) = yes;

The tricky indexing removes the last row and column. The display confirms this:

----     31 SET squares  numbering of squares

            j1          j2

i1         YES         YES
i2         YES         YES

I also want to map between squares \((i,j)\) and the corresponding four edges \((i',j')\rightarrow (i'',j'')\). So I generate this crazy 6-dimensional structure semap (for squares-edge map). A square \((i,j)\) has two edges leaving from \((i,j\)) and two edges arriving at \((i+1,j+1)\). That is what I use here:

semap(squares(i,j),edges(i,j,ii,jj)) = yes;
semap(squares(i,j),edges(ii,jj,i+1,j+1)) = yes;

I could have dropped the word edges here, as they don't filter anything out. I left it for clarity. 

----     31 SET semap  map between squares and edges

i1.j1.i1.j1.i1.j2,    i1.j1.i1.j1.i2.j1,    i1.j1.i1.j2.i2.j2,    i1.j1.i2.j1.i2.j2,    i1.j2.i1.j2.i1.j3,    i1.j2.i1.j2.i2.j2
i1.j2.i1.j3.i2.j3,    i1.j2.i2.j2.i2.j3,    i2.j1.i2.j1.i2.j2,    i2.j1.i2.j1.i3.j1,    i2.j1.i2.j2.i3.j2,    i2.j1.i3.j1.i3.j2
i2.j2.i2.j2.i2.j3,    i2.j2.i2.j2.i3.j2,    i2.j2.i2.j3.i3.j3,    i2.j2.i3.j2.i3.j3

Of course, this is just one possible approach. Another would be to use a one-dimensional numbering scheme for edges and squares. You still need to specify which edges are used by a square. 

Model

The model itself is not super complicated. 

*------------------------------------------------- * MIP model *-------------------------------------------------   variables     z                'dummy objective'     x(i,j,i,j,c)     'edge colors'     sqcol(i,j,c)     'square colors' ; binary variable x,sqcol;   Equations      dummyobj      onecolor(i,j,i,j)  'an edge has exactly one color'      implication1(i,j,i,j,i,j,c) 'ecol(edge,c)=1 for any edge => sqcol(square,c)=1'      implication2(i,j,c) 'ecol(edge,c)=0 for all edges => sqcol(square,c)=0'      squarecolors(i,j)   'exactly 2 colors per square' ;   dummyobj.. z =e= 0;   onecolor(edges).. sum(c, x(edges,c)) =e= 1; implication1(squares,edges,c)$semap(squares,edges).. sqcol(squares,c) =g= x(edges,c); implication2(squares,c).. sqcol(squares,c) =l= sum(semap(squares,edges), x(edges,c)); squarecolors(squares).. sum(c, sqcol(squares,c)) =e= 2;   model coloring /all/; solve coloring using mip minimizing z;

We don't have an objective: we only need to find a feasible solution. In GAMS that means that we use a dummy objective.

In the equations, I heavily use the sets squares and edges. This makes the equations simpler, but a little bit misleading. It looks like we have two-dimensional variables and equations, while, in fact, these are multi-dimensional structures. That can be seen in the declarations. Counting the colors in a square is done by two sets of implications:

1. If any of the edges has color \(c\), then the square includes color \(c\)
2. If all of the edges don't have color \(c\), then the square does not include color \(c\) 

I ignore here the requirement "2 edges per color". 

Solution printing

A quick and dirty solution printout shows:

----     77 PARAMETER solution  

i1.j1.i1.j2 GREEN,    i1.j1.i2.j1 GREEN,    i1.j2.i1.j3   RED,    i1.j2.i2.j2   RED,    i1.j3.i2.j3  BLUE,    i2.j1.i2.j2   RED
i2.j1.i3.j1  BLUE,    i2.j2.i2.j3  BLUE,    i2.j2.i3.j2   RED,    i2.j3.i3.j3   RED,    i3.j1.i3.j2   RED,    i3.j2.i3.j3   RED

This is only slightly better than just displaying variable ecol. It may be better to draw a picture. We can write some HTML and SVG to draw the grid and color the edges. For a 25x25 example, this can look like:

A little bit of a low-tech solution: just write out the HTML using the PUT statement. However, it is difficult to beat the simplicity of this approach.

We can solve even a 100x100 grid very quickly (0.89 seconds to find a feasible integer solution). The picture becomes a bit difficult to grasp.

2-edges per color

If you want to use 2 edges per color we need to say actually 

    the number of edges in a square with color \(c\) is 0 or 2.

This can be modeled with a binary variable multiplied by 2. I.e.:

binary variables ZeroOrTwo(i,j,c) 'zero or two edges per color';

equation squareColors2(i,j,c) 'exactly 2 colors per square, and two edges per color';

squareColors2(squares,c).. sum(semap(squares,edges),x(edges,c)) =e= 2*ZeroOrTwo(squares,c);

This constraint implies that only 2 colors are used  for the four edges for each square. So this replaces the three constraints implication1, implication2, SquareColors in the previous model. I.e. we added extra binary variables, but the number of constraints and their complexity decreases. 

This is called Model 2 in the GAMS code.

When I run this, we see for the 25x25 problem:

This model is more difficult. The 25x25 case solves fast, but the 100x100 case takes a while (I never finished it).

A faster 100x100 model

One rule in MIP modeling is: never give up. The previous attempt at solving the 100x100 problem with the 2 edges per color constraint was not very successful. A variant of the model uses the constraints of model 1 and model 2. We just say:

model coloring3 /all/;

solve coloring3 using mip minimizing z;

Somehow this is very fast.

This model produces a rather regular picture:

This picture is for the 25x25 case, but it is similar for the 100x100 case.

It is interesting that combining models 1 and 2 produces faster results than just model 2.

Counting solutions

The original post asked for the total number of solutions. When throwing this model (just 4x4 grid) at the solution pool I see:

--- Cplex Time: 5.94sec (det. 4774.48 ticks)

--- Dumping 373248 solutions from the solution pool...

 

For a 4x5 grid this becomes:

--- Cplex Time: 440.67sec (det. 167328.64 ticks)
--- Dumping 8957952 solutions from the solution pool... 

Obviously, there are way too many solutions to answer this question by a numerical algorithm for a larger grid. This may not be such an interesting question.

 

Conclusion

The model is interesting because of its data structures and indexing techniques. Using sparse multi-dimensional sets, we can populate somewhat complex sets in just a few lines while still being understandable. This would be done very differently in, say, a Python environment. This is a good example where a line-by-line translation does not really work.

A second feature is that I use multi-dimensional sets to simplify the constraints. Constraints are often the most difficult to debug, so making them as simple as possible is beneficial. 

Finally, I used some HTML to visualize the solution. Clearly, that is useful up to a point.

References

1. Given an m * n grid, each edge must be colored. However, there are constraints, https://stackoverflow.com/questions/76645751/given-an-m-n-grid-each-edge-must-be-colored-however-there-are-constraints

Appendix: complete GAMS model 

$ontext      Color edges of a grid.      From https://stackoverflow.com/questions/76645751/given-an-m-n-grid-each-edge-must-be-colored-however-there-are-constraints          Given an m * n grid, each edge must be colored. However, there are 2 constraints:            1.  Entire grid must use only 3 unique colors          2.  Each square in the grid must be colored using exactly 2 colors (2 edges per color)     Model 1: forget 2 edges per color Model 2: add  2 edges per color Model 3: combine models 1 and 2     $offText   *------------------------------------------------- * sets to represent grid and edges. * representation of squares is by left-lower point *-------------------------------------------------   set    i 'grid x-coordinate'  /i1*i10/    j 'grid y-coordinate'  /j1*j10/    edges(i,j,i,j) 'connects two points'    squares(i,j) 'numbering of squares'    semap(i,j,i,j,i,j) 'map between squares and edges'    c 'colors' /red,green,blue/ ;   alias (i,ii),(j,jj);   edges(i,j,i+1,j) = yes; edges(i,j,i,j+1) = yes;   squares(i-1,j-1) = yes;   semap(squares(i,j),edges(i,j,ii,jj)) = yes; semap(squares(i,j),edges(ii,jj,i+1,j+1)) = yes;   option edges:0:0:8,semap:0:0:6; display i,j,edges,squares,semap;   *------------------------------------------------- * MIP model 1 * ignore two edges per color in a square *-------------------------------------------------   variables   z                'dummy objective'   x(i,j,i,j,c)     'edge colors'   sqCol(i,j,c)     'square colors' ; binary variable x,sqCol;   Equations    dummyObj    oneColor(i,j,i,j)  'an edge has exactly one color'    implication1(i,j,i,j,i,j,c) 'ecol(edge,c)=1 for any edge => sqcol(square,c)=1'    implication2(i,j,c) 'ecol(edge,c)=0 for all edges => sqcol(square,c)=0'    squareColors(i,j)   'exactly 2 colors per square' ;   dummyobj.. z =e= 0;   onecolor(edges).. sum(c, x(edges,c)) =e= 1; implication1(squares,edges,c)$semap(squares,edges).. sqcol(squares,c) =g= x(edges,c); implication2(squares,c).. sqcol(squares,c) =l= sum(semap(squares,edges), x(edges,c)); squarecolors(squares).. sum(c, sqcol(squares,c)) =e= 2;   model coloring1 /all/; solve coloring1 using mip minimizing z;   *------------------------------------------------- * MIP model 2 * add two edges per color in a square *-------------------------------------------------   binary variables ZeroOrTwo(i,j,c) 'zero or two edges per color';   equation squareColors2(i,j,c) 'exactly 2 colors per square, and two edges per color';   squareColors2(squares,c).. sum(semap(squares,edges),x(edges,c)) =e= 2*ZeroOrTwo(squares,c);   model coloring2 /dummyObj,oneColor,squareColors2/; option threads=0; *solve coloring2 using mip minimizing z;   *------------------------------------------------- * MIP model 3 * add two edges per color in a square *-------------------------------------------------   model coloring3 /all/; *solve coloring3 using mip minimizing z;     *------------------------------------------------- * print solution *-------------------------------------------------   acronym RED,GREEN,BLUE; parameter   col(c) / red RED, green GREEN, blue BLUE/   solution(i,j,i,j) ;   loop(c,   solution(edges)$(x.l(edges,c)>0.5) = col(c); );   option solution:0:0:8; display solution;     *------------------------------------------------- * draw solution *-------------------------------------------------     file f /solution.html/; put f;   put '<svg width=800 height=800 viewbox = "0 0 ',(card(i)+1):0:0,' ',(card(j)+1):0:0,'">'/; loop(edges(i,j,ii,jj),    loop(c$(x.l(edges,c)>0.5),        put '<line x1=',ord(i):0:0,' x2=',ord(ii):0:0,' y1=',ord(j):0:0,' y2=',ord(jj):0:0,' stroke="',c.tl:0,'" stroke-width="0.07"/>'/;    ) ) put '</svg>'/;   execute "shellexecute solution.html";