ILP Part 60 – Fencing Match – Random IT Utensils

This is the sixtieth part of the ILP series. For your convenience you can find other parts in the table of contents in Part 1 – Boolean algebra

Today we are going to solve Fencing Match. We have a board with some numbers on it. Each number indicates how many edges of that field must be included in the road. We need to build the closed road (loop) meeting the requirements.

Code:

var solver = new CplexMilpSolver(15);

var board = "23 1 11\n  3 2  \n  3 2  \n3     2\n  3 1  \n  0 1  \n23 3 20".Split(new []{"\n"}, StringSplitOptions.RemoveEmptyEntries).ToList();

var width = board[0].Count();
var height = board.Count();

Console.WriteLine("Create edges");
var edges = new List<List<System.Dynamic.ExpandoObject>>();
for(var row = 0; row < height;++row){
	var rowVars = new List<System.Dynamic.ExpandoObject>();
	
	for(var column = 0; column < width;++column){
		dynamic expando = new System.Dynamic.ExpandoObject();
		expando.Right = solver.CreateAnonymous(Domain.BinaryInteger);
		expando.Left = solver.CreateAnonymous(Domain.BinaryInteger);
		expando.Bottom = solver.CreateAnonymous(Domain.BinaryInteger);
		expando.Top = solver.CreateAnonymous(Domain.BinaryInteger);
		
		rowVars.Add(expando);
	}
	
	edges.Add(rowVars);
}

Func<int, int, Action<System.Dynamic.ExpandoObject>, System.Dynamic.ExpandoObject> extractField = (row, column, action) => {
	if(row >= 0 && row < height && column >= 0 && column < width){
		dynamic field = edges[row][column];
		action(field);
		return field;
	}
	
	return null;
};
Func<int, int, Action<dynamic>, System.Dynamic.ExpandoObject> upField = (row, column, action) => extractField(row-1, column, action);
Func<int, int, Action<dynamic>, System.Dynamic.ExpandoObject> downField = (row, column, action) => extractField(row+1, column, action);
Func<int, int, Action<dynamic>, System.Dynamic.ExpandoObject> leftField = (row, column, action) => extractField(row, column-1, action);
Func<int, int, Action<dynamic>, System.Dynamic.ExpandoObject> rightField = (row, column, action) => extractField(row, column+1, action);

Console.WriteLine("Make sure fields neighbours match");
for(var row = 0; row < height;++row){
	for(var column = 0; column < width;++column){
		dynamic self = edges[row][column];
		upField(row, column, t => solver.Set(ConstraintType.Equal, self.Top, t.Bottom));
		leftField(row, column, t => solver.Set(ConstraintType.Equal, self.Left, t.Right));
		rightField(row, column, t => solver.Set(ConstraintType.Equal, self.Right, t.Left));
		downField(row, column, t => solver.Set(ConstraintType.Equal, self.Bottom, t.Top));
	}
}

Console.WriteLine("Create crossing points");
var points = new List<List<System.Dynamic.ExpandoObject>>();
for(var row = 0; row <= height;++row){
	var rowVars = new List<System.Dynamic.ExpandoObject>();
	
	for(var column = 0; column <= width;++column){
		dynamic expando = new System.Dynamic.ExpandoObject();
		expando.IdVar = solver.CreateAnonymous(Domain.PositiveOrZeroInteger);
		if(column < width){
			expando.Right = solver.CreateAnonymous(Domain.BinaryInteger);
		}
		if(column > 0){
			expando.Left = solver.CreateAnonymous(Domain.BinaryInteger);
		}
		if(row > 0){
			expando.Up = solver.CreateAnonymous(Domain.BinaryInteger);
		}
		if(row < height){
			expando.Down = solver.CreateAnonymous(Domain.BinaryInteger);
		}
		
		rowVars.Add(expando);
	}
	
	points.Add(rowVars);
}

Console.WriteLine("Make sure crossing points match");
for(var row = 0; row < height;++row){
	for(var column = 0; column < width;++column){
		dynamic self = points[row][column];
		dynamic field = edges[row][column];
		solver.Set(ConstraintType.Equal, self.Right, field.Top);
	}
}
for(var column = 0; column < width;++column){
	dynamic self = points[height][column];
	dynamic field = edges[height-1][column];
	solver.Set(ConstraintType.Equal, self.Right, field.Bottom);
}
for(var row = 0; row < height;++row){
	for(var column = 1; column <= width;++column){
		dynamic self = points[row][column];
		dynamic field = edges[row][column-1];
		solver.Set(ConstraintType.Equal, self.Left, field.Top);
	}
}
for(var column = 1; column <= width;++column){
	dynamic self = points[height][column];
	dynamic field = edges[height-1][column-1];
	solver.Set(ConstraintType.Equal, self.Left, field.Bottom);
}
for(var column = 0; column < width;++column){
	for(var row = 0; row < height;++row){
		dynamic self = points[row][column];
		dynamic field = edges[row][column];
		solver.Set(ConstraintType.Equal, self.Down, field.Left);
	}
}
for(var row = 0; row < height;++row){
	dynamic self = points[row][width];
	dynamic field = edges[row][width-1];
	solver.Set(ConstraintType.Equal, self.Down, field.Right);
}
for(var column = 0; column < width;++column){
	for(var row = 1; row <= height;++row){
		dynamic self = points[row][column];
		dynamic field = edges[row-1][column];
		solver.Set(ConstraintType.Equal, self.Up, field.Left);
	}
}
for(var row = 1; row <= height;++row){
	dynamic self = points[row][width];
	dynamic field = edges[row-1][width-1];
	solver.Set(ConstraintType.Equal, self.Up, field.Right);
}

Console.WriteLine("Make sure edges match the number");
for(var row = 0; row < height;++row){
	for(var column = 0; column < width;++column){
		if(board[row][column] != ' '){
			var count = board[row][column] - '0';
			dynamic self = edges[row][column];
			solver.Set(ConstraintType.Equal, solver.Operation(OperationType.Addition, self.Top, self.Bottom, self.Left, self.Right), solver.FromConstant(count));
		}
	}
}

Console.WriteLine("Make sure line is continuous");
for(var row = 0; row <= height;++row){
	for(var column = 0; column <= width;++column){
		var neighbours = new List<IVariable>();
		dynamic self = points[row][column];
		
		var selfDictionary = ((IDictionary<String, object>)self);
		
		if(selfDictionary.ContainsKey("Up")) neighbours.Add(self.Up);
		if(selfDictionary.ContainsKey("Down")) neighbours.Add(self.Down);
		if(selfDictionary.ContainsKey("Left")) neighbours.Add(self.Left);
		if(selfDictionary.ContainsKey("Right")) neighbours.Add(self.Right);
		
		var sum = solver.Operation(OperationType.Addition, neighbours.ToArray());
		
		var empty = solver.Operation(OperationType.IsEqual, sum, solver.FromConstant(0));
		var connection = solver.Operation(OperationType.IsEqual, sum, solver.FromConstant(2));
		
		solver.Set(ConstraintType.Equal, solver.Operation(OperationType.Disjunction, empty, connection), solver.FromConstant(1));
	}
}

Console.WriteLine("Make sure there is one line");
var isTopLeftSum = solver.FromConstant(0);
var isBottomRightSum = solver.FromConstant(0);

for(var row = 0; row <= height;++row){
	for(var column = 0; column <= width;++column){
		dynamic self = points[row][column];
		var selfDictionary = ((IDictionary<String, object>)self);
		
		self.IsTopLeft = solver.CreateAnonymous(Domain.BinaryInteger);
		self.IsBottomRight = solver.CreateAnonymous(Domain.BinaryInteger);
		
		isTopLeftSum = solver.Operation(OperationType.Addition, isTopLeftSum, self.IsTopLeft);
		isBottomRightSum = solver.Operation(OperationType.Addition, isBottomRightSum, self.IsBottomRight);
		var isSpecial = solver.Operation(OperationType.Disjunction, self.IsTopLeft, self.IsBottomRight);
		
		solver.Set(ConstraintType.Equal, solver.FromConstant(0), solver.Operation(OperationType.Condition, self.IsBottomRight, self.IdVar, solver.CreateAnonymous(Domain.PositiveOrZeroInteger)));
		
		var neighbours = new List<Tuple<CplexVariable, CplexVariable>>();
		
		if(selfDictionary.ContainsKey("Up")) {
			dynamic up = points[row-1][column];
			neighbours.Add(Tuple.Create(up.IdVar, self.Up));
		}
		if(selfDictionary.ContainsKey("Down")) {
			dynamic down = points[row+1][column];
			neighbours.Add(Tuple.Create(down.IdVar, self.Down));
		}
		if(selfDictionary.ContainsKey("Left")) {
			dynamic left = points[row][column-1];
			neighbours.Add(Tuple.Create(left.IdVar, self.Left));
		}
		if(selfDictionary.ContainsKey("Right")) {
			dynamic right = points[row][column+1];
			neighbours.Add(Tuple.Create(right.IdVar, self.Right));
		}
		
		foreach(dynamic neighbour in neighbours){
			solver.Set(ConstraintType.Equal, solver.FromConstant(1), solver.Operation(OperationType.AbsoluteValue, solver.Operation(OperationType.Subtraction, self.IdVar, solver.Operation(OperationType.Condition, neighbour.Item2, neighbour.Item1, solver.CreateAnonymous(Domain.PositiveOrZeroInteger)))));
		}
		
		for(int i=0;i<neighbours.Count();++i){
			for(int j=i+1;j<neighbours.Count();++j){
				var firstIdVar = solver.Operation(OperationType.Condition, neighbours[i].Item2, neighbours[i].Item1, solver.CreateAnonymous(Domain.PositiveOrZeroInteger));
				var secondIdVar = solver.Operation(OperationType.Condition, neighbours[j].Item2, neighbours[j].Item1, solver.CreateAnonymous(Domain.PositiveOrZeroInteger));
				
				solver.Set(ConstraintType.Equal, solver.Operation(OperationType.AbsoluteValue, solver.Operation(OperationType.Subtraction, firstIdVar, secondIdVar)), solver.Operation(OperationType.Condition, isSpecial, solver.FromConstant(0), solver.FromConstant(2)));
			}
		}
		
		solver.Set(ConstraintType.Equal, solver.FromConstant(1), solver.Operation(OperationType.MaterialImplication, isSpecial, solver.Operation(OperationType.Disjunction, neighbours.Select(n => n.Item2).ToArray())));
	}
}

solver.Set(ConstraintType.Equal, isTopLeftSum, solver.FromConstant(1));
solver.Set(ConstraintType.Equal, isBottomRightSum, solver.FromConstant(1));

solver.AddGoal("Goal", solver.FromConstant(0));
solver.Solve();

for(var row = 0; row < height;++row){
	for(var column = 0; column < width; ++ column){
		dynamic self = edges[row][column];
		Console.Write(solver.GetValue(self.Top) > 0.5 ? " -" : "  ");
	}
	Console.WriteLine();
	
	for(var column = 0; column < width; ++ column){
		dynamic self = edges[row][column];
		Console.Write(solver.GetValue(self.Left) > 0.5 ? "|" : " ");
		Console.Write(board[row][column]);
		
		if(column == width - 1){
			Console.Write(solver.GetValue(self.Right) > 0.5 ? "|" : " ");
		}
	}
	Console.WriteLine();
	
	if(row == height - 1){
		for(var column = 0; column < width; ++ column){
			dynamic self = edges[row][column];
			Console.Write(solver.GetValue(self.Bottom) > 0.5 ? " -" : "  ");
		}
		Console.WriteLine();
	}
}

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

var solver = new CplexMilpSolver(15);

var board = "23 1 11\n 3 2 \n 3 2 \n3 2\n 3 1 \n 0 1 \n23 3 20".Split(new []{"\n"}, StringSplitOptions.RemoveEmptyEntries).ToList();

var width = board[0].Count();

var height = board.Count();

Console.WriteLine("Create edges");

var edges = new List<List<System.Dynamic.ExpandoObject>>();

for(var row = 0; row < height;++row){

var rowVars = new List<System.Dynamic.ExpandoObject>();

for(var column = 0; column < width;++column){

dynamic expando = new System.Dynamic.ExpandoObject();

expando.Right = solver.CreateAnonymous(Domain.BinaryInteger);

expando.Left = solver.CreateAnonymous(Domain.BinaryInteger);

expando.Bottom = solver.CreateAnonymous(Domain.BinaryInteger);

expando.Top = solver.CreateAnonymous(Domain.BinaryInteger);

rowVars.Add(expando);

}

edges.Add(rowVars);

}

Func<int, int, Action<System.Dynamic.ExpandoObject>, System.Dynamic.ExpandoObject> extractField = (row, column, action) => {

if(row >= 0 && row < height && column >= 0 && column < width){

dynamic field = edges[row][column];

action(field);

return field;

}

return null;

};

Func<int, int, Action<dynamic>, System.Dynamic.ExpandoObject> upField = (row, column, action) => extractField(row-1, column, action);

Func<int, int, Action<dynamic>, System.Dynamic.ExpandoObject> downField = (row, column, action) => extractField(row+1, column, action);

Func<int, int, Action<dynamic>, System.Dynamic.ExpandoObject> leftField = (row, column, action) => extractField(row, column-1, action);

Func<int, int, Action<dynamic>, System.Dynamic.ExpandoObject> rightField = (row, column, action) => extractField(row, column+1, action);

Console.WriteLine("Make sure fields neighbours match");

for(var row = 0; row < height;++row){

for(var column = 0; column < width;++column){

dynamic self = edges[row][column];

upField(row, column, t => solver.Set(ConstraintType.Equal, self.Top, t.Bottom));

leftField(row, column, t => solver.Set(ConstraintType.Equal, self.Left, t.Right));

rightField(row, column, t => solver.Set(ConstraintType.Equal, self.Right, t.Left));

downField(row, column, t => solver.Set(ConstraintType.Equal, self.Bottom, t.Top));

}

Console.WriteLine("Create crossing points");

var points = new List<List<System.Dynamic.ExpandoObject>>();

for(var row = 0; row <= height;++row){

var rowVars = new List<System.Dynamic.ExpandoObject>();

for(var column = 0; column <= width;++column){

dynamic expando = new System.Dynamic.ExpandoObject();

expando.IdVar = solver.CreateAnonymous(Domain.PositiveOrZeroInteger);

if(column < width){

expando.Right = solver.CreateAnonymous(Domain.BinaryInteger);

}

if(column > 0){

expando.Left = solver.CreateAnonymous(Domain.BinaryInteger);

}

if(row > 0){

expando.Up = solver.CreateAnonymous(Domain.BinaryInteger);

}

if(row < height){

expando.Down = solver.CreateAnonymous(Domain.BinaryInteger);

}

rowVars.Add(expando);

}

points.Add(rowVars);

}

Console.WriteLine("Make sure crossing points match");

for(var row = 0; row < height;++row){

for(var column = 0; column < width;++column){

dynamic self = points[row][column];

dynamic field = edges[row][column];

solver.Set(ConstraintType.Equal, self.Right, field.Top);

}

for(var column = 0; column < width;++column){

dynamic self = points[height][column];

dynamic field = edges[height-1][column];

solver.Set(ConstraintType.Equal, self.Right, field.Bottom);

}

for(var row = 0; row < height;++row){

for(var column = 1; column <= width;++column){

dynamic self = points[row][column];

dynamic field = edges[row][column-1];

solver.Set(ConstraintType.Equal, self.Left, field.Top);

}

for(var column = 1; column <= width;++column){

dynamic self = points[height][column];

dynamic field = edges[height-1][column-1];

solver.Set(ConstraintType.Equal, self.Left, field.Bottom);

}

for(var column = 0; column < width;++column){

for(var row = 0; row < height;++row){

dynamic self = points[row][column];

dynamic field = edges[row][column];

solver.Set(ConstraintType.Equal, self.Down, field.Left);

}

for(var row = 0; row < height;++row){

dynamic self = points[row][width];

dynamic field = edges[row][width-1];

solver.Set(ConstraintType.Equal, self.Down, field.Right);

}

for(var column = 0; column < width;++column){

for(var row = 1; row <= height;++row){

dynamic self = points[row][column];

dynamic field = edges[row-1][column];

solver.Set(ConstraintType.Equal, self.Up, field.Left);

}

for(var row = 1; row <= height;++row){

dynamic self = points[row][width];

dynamic field = edges[row-1][width-1];

solver.Set(ConstraintType.Equal, self.Up, field.Right);

}

Console.WriteLine("Make sure edges match the number");

for(var row = 0; row < height;++row){

for(var column = 0; column < width;++column){

if(board[row][column] != ' '){

var count = board[row][column] - '0';

dynamic self = edges[row][column];

solver.Set(ConstraintType.Equal, solver.Operation(OperationType.Addition, self.Top, self.Bottom, self.Left, self.Right), solver.FromConstant(count));

}

Console.WriteLine("Make sure line is continuous");

for(var row = 0; row <= height;++row){

for(var column = 0; column <= width;++column){

var neighbours = new List<IVariable>();

dynamic self = points[row][column];

var selfDictionary = ((IDictionary<String, object>)self);

if(selfDictionary.ContainsKey("Up")) neighbours.Add(self.Up);

if(selfDictionary.ContainsKey("Down")) neighbours.Add(self.Down);

if(selfDictionary.ContainsKey("Left")) neighbours.Add(self.Left);

if(selfDictionary.ContainsKey("Right")) neighbours.Add(self.Right);

var sum = solver.Operation(OperationType.Addition, neighbours.ToArray());

var empty = solver.Operation(OperationType.IsEqual, sum, solver.FromConstant(0));

var connection = solver.Operation(OperationType.IsEqual, sum, solver.FromConstant(2));

solver.Set(ConstraintType.Equal, solver.Operation(OperationType.Disjunction, empty, connection), solver.FromConstant(1));

}

Console.WriteLine("Make sure there is one line");

var isTopLeftSum = solver.FromConstant(0);

var isBottomRightSum = solver.FromConstant(0);

for(var row = 0; row <= height;++row){

for(var column = 0; column <= width;++column){

dynamic self = points[row][column];

var selfDictionary = ((IDictionary<String, object>)self);

self.IsTopLeft = solver.CreateAnonymous(Domain.BinaryInteger);

self.IsBottomRight = solver.CreateAnonymous(Domain.BinaryInteger);

isTopLeftSum = solver.Operation(OperationType.Addition, isTopLeftSum, self.IsTopLeft);

isBottomRightSum = solver.Operation(OperationType.Addition, isBottomRightSum, self.IsBottomRight);

var isSpecial = solver.Operation(OperationType.Disjunction, self.IsTopLeft, self.IsBottomRight);

solver.Set(ConstraintType.Equal, solver.FromConstant(0), solver.Operation(OperationType.Condition, self.IsBottomRight, self.IdVar, solver.CreateAnonymous(Domain.PositiveOrZeroInteger)));

var neighbours = new List<Tuple<CplexVariable, CplexVariable>>();

if(selfDictionary.ContainsKey("Up")) {

dynamic up = points[row-1][column];

neighbours.Add(Tuple.Create(up.IdVar, self.Up));

}

if(selfDictionary.ContainsKey("Down")) {

dynamic down = points[row+1][column];

neighbours.Add(Tuple.Create(down.IdVar, self.Down));

}

if(selfDictionary.ContainsKey("Left")) {

dynamic left = points[row][column-1];

neighbours.Add(Tuple.Create(left.IdVar, self.Left));

}

if(selfDictionary.ContainsKey("Right")) {

dynamic right = points[row][column+1];

neighbours.Add(Tuple.Create(right.IdVar, self.Right));

}

foreach(dynamic neighbour in neighbours){

solver.Set(ConstraintType.Equal, solver.FromConstant(1), solver.Operation(OperationType.AbsoluteValue, solver.Operation(OperationType.Subtraction, self.IdVar, solver.Operation(OperationType.Condition, neighbour.Item2, neighbour.Item1, solver.CreateAnonymous(Domain.PositiveOrZeroInteger)))));

}

for(int i=0;i<neighbours.Count();++i){

for(int j=i+1;j<neighbours.Count();++j){

var firstIdVar = solver.Operation(OperationType.Condition, neighbours[i].Item2, neighbours[i].Item1, solver.CreateAnonymous(Domain.PositiveOrZeroInteger));

var secondIdVar = solver.Operation(OperationType.Condition, neighbours[j].Item2, neighbours[j].Item1, solver.CreateAnonymous(Domain.PositiveOrZeroInteger));

solver.Set(ConstraintType.Equal, solver.Operation(OperationType.AbsoluteValue, solver.Operation(OperationType.Subtraction, firstIdVar, secondIdVar)), solver.Operation(OperationType.Condition, isSpecial, solver.FromConstant(0), solver.FromConstant(2)));

}

solver.Set(ConstraintType.Equal, solver.FromConstant(1), solver.Operation(OperationType.MaterialImplication, isSpecial, solver.Operation(OperationType.Disjunction, neighbours.Select(n => n.Item2).ToArray())));

}

solver.Set(ConstraintType.Equal, isTopLeftSum, solver.FromConstant(1));

solver.Set(ConstraintType.Equal, isBottomRightSum, solver.FromConstant(1));

solver.AddGoal("Goal", solver.FromConstant(0));

solver.Solve();

for(var row = 0; row < height;++row){

for(var column = 0; column < width; ++ column){

dynamic self = edges[row][column];

Console.Write(solver.GetValue(self.Top) > 0.5 ? " -" : " ");

}

Console.WriteLine();

for(var column = 0; column < width; ++ column){

dynamic self = edges[row][column];

Console.Write(solver.GetValue(self.Left) > 0.5 ? "|" : " ");

Console.Write(board[row][column]);

if(column == width - 1){

Console.Write(solver.GetValue(self.Right) > 0.5 ? "|" : " ");

}

Console.WriteLine();

if(row == height - 1){

for(var column = 0; column < width; ++ column){

dynamic self = edges[row][column];

Console.Write(solver.GetValue(self.Bottom) > 0.5 ? " -" : " ");

}

Console.WriteLine();

}

In line 3 we encode the board as a string.

We are going to model two things: fields and crossing points. Each field has four edges (top, bottom, left, right) represented by binary variables. You can see that in lines 8-24.

In lines 26-38 we define some helpers, and then in lines 40-49 we make sure that edges of adjacent fields match (so right edge of some field is equal to the left edge of the field on the right).

Second thing we model is crossing points. For each field there is one crossing point in the top left corner, also, fields on the right and bottom boundary have additional crossing points in the bottom right corner. This is in lines 51-76. We also define another variable IdVar for numbering the points in a BFS-like style later.

Next, we make sure crossing points match. Each crossing point includes up to four edges (left, right, up, down) which need to match with field edges. A lot of juggling in lines 78-137.

That was “simple” but cumbersome. Once we have the correct model we can easily model our requirements.

First, in lines 128-137 we make sure that edges of fields with numbers are selected properly.

In lines 139-159 we make sure the road is continuous. The same trick as always — we take edges for a given crossing point and make sure that either none of them or exactly two of them are selected for the road.

Finally, in lines 161-2016 we make sure there is exactly one road. We ignored this requirement in previous parts so it’s high time to takle it.

We base our idea on BFS-like approach. We want to cut the road in two equal halves, identify one “end” of it as a starting point and then number all points on the road based on the Manhattan distance.

We designate two points which we call “top left” and “bottom right”. They don’t need to be actual top-left-most and bottom-right-most points on the road. We need to make sure there is exactly one top left and one bottom right. This will guarantee there is exactly one road. So we start with lines 162 and 163 which we use to count the points and then in lines 215-216 we set the constraint.

Starting in line 165 we iterate through all the crossing points. We define binaries for top left and bottom right points in lines 170-171. We then add them in lines 173 and 174. Effectively only one road point will be a “top left” and some other will be “bottom right”.

We also calculate if given point is special – either top left or bottom right – in line 175.

Now we need to implement the BFS-like logic. In line 177 we use conditions to specify that if this point is a bottom-right one then its IdVar must be equal to zero. This effectively starts the numbering. Notice that IdVar is a natural number. If this point is not a bottom-right corner then line 177 set its IdVar to a free variable effectively neutralizing the constraint. This is actually how you make “if” condition in ILP (instead of “if-else” expression).

Next, in lines 181-196 we find neighbors. We group tuples of neighbours’ IdVars and our edges. Then in lines 198-200 we add a constraint that if there is a road (edge is equal to one) then neighbour’s IdVar must differ by one from our value. Since we operate on natural numbers and we start numbering from zero, neighbours will have increasing values.

In lines 202-209 we make some magic. If you take a continuous road and number points starting from the corner, you can see that for each “middle” point neighbours’ IdVars differ by two. However, for two special points (top left and bottom right corners) the neighbours are equal. This is what we set here — we conditionally take IdVars for neighbours (if there are edges) and then set the difference to zero or two, depending on whether given point is a special corner or not.

Finally, in line 211 we make sure that the special point can be on the road only — if is special is one (so given point is a corner) then we use material implication to indicate that at we must have some edges set for this field.

Output:

Tried aggregator 5 times.
MIP Presolve eliminated 15467 rows and 2645 columns.
MIP Presolve modified 36628 coefficients.
Aggregator did 16957 substitutions.
Reduced MIP has 28421 rows, 17654 columns, and 84682 nonzeros.
Reduced MIP has 13739 binaries, 3915 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.59 sec. (334.22 ticks)
Probing fixed 4272 vars, tightened 1764 bounds.
Probing changed sense of 3 constraints.
Probing time = 1.98 sec. (381.40 ticks)
Cover probing fixed 0 vars, tightened 8 bounds.
Tried aggregator 3 times.
MIP Presolve eliminated 10986 rows and 6579 columns.
MIP Presolve modified 5755 coefficients.
Aggregator did 2635 substitutions.
Reduced MIP has 14800 rows, 8440 columns, and 42921 nonzeros.
Reduced MIP has 5666 binaries, 2774 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.13 sec. (129.00 ticks)
Probing time = 0.08 sec. (12.19 ticks)
Tried aggregator 1 time.
MIP Presolve eliminated 1573 rows and 944 columns.
Reduced MIP has 13227 rows, 7496 columns, and 38740 nonzeros.
Reduced MIP has 5028 binaries, 2468 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.03 sec. (18.68 ticks)
Probing time = 0.05 sec. (7.86 ticks)
Clique table members: 33932.
MIP emphasis: balance optimality and feasibility.
MIP search method: dynamic search.
Parallel mode: deterministic, using up to 4 threads.
Root relaxation solution time = 0.22 sec. (274.85 ticks)

        Nodes                                         Cuts/
   Node  Left     Objective  IInf  Best Integer    Best Bound    ItCnt     Gap

      0     0        0.0000  1857                      0.0000     3767
      0     0        0.0000  3096                  Cuts: 1249     5833
      0     0        0.0000  2658                   Cuts: 533     7307
      0     0        0.0000  2459                   Cuts: 351     7947
      0     0        0.0000  2214                   Cuts: 197     8257
      0     0        0.0000  2047                    Cuts: 42     8477
      0     0        0.0000  2033                    Cuts: 58     8531
      0     0        0.0000  2020                    Cuts: 58     8592
      0     0        0.0000  1991                    Cuts: 24     8629
      0     0        0.0000  1984                     Cuts: 7     8662
      0     2        0.0000  1984                      0.0000     8662
Elapsed time = 11.20 sec. (4434.32 ticks, tree = 0.01 MB, solutions = 0)
      5     7        0.0000  1582                      0.0000    10209
     11    13        0.0000  2542                      0.0000    14980
     43    45        0.0000  2576                      0.0000    25907
    113    93        0.0000  2801                      0.0000    35995
    259   183        0.0000   607                      0.0000    39929
    335   213        0.0000  1703                      0.0000    43877
    547   235        0.0000  1255                      0.0000    48099
    700   236    infeasible                            0.0000    53520
    716   230    infeasible                            0.0000    54089
   1060   263        0.0000   250                      0.0000    67348
Elapsed time = 17.25 sec. (7827.67 ticks, tree = 7.18 MB, solutions = 0)
   1409   266    infeasible                            0.0000    81492
   1776   260    infeasible                            0.0000    97601
   2184   267    infeasible                            0.0000   115041
   2573   274    infeasible                            0.0000   130798
   3018   277    infeasible                            0.0000   145637
   3263   245    infeasible                            0.0000   156266
   3685   283    infeasible                            0.0000   171075
   4175   263    infeasible                            0.0000   186666
   4727   273        0.0000   320                      0.0000   201106
   4777   273        0.0000  2125                      0.0000   211968
Elapsed time = 49.45 sec. (22799.02 ticks, tree = 24.51 MB, solutions = 0)
   4848   289        0.0000  1670                      0.0000   226062
   5209   212        0.0000   821                      0.0000   248000
   5264   157    infeasible                            0.0000   265899
   5339    82    infeasible                            0.0000   284511
   5374    49        0.0000   524                     -0.0000   306065
   5391    39        0.0000  1265                     -0.0000   320383
   5411    27    infeasible                           -0.0000   332759
   5437    23    infeasible                           -0.0000   344390
*  5457+   17                            0.0000       -0.0000   346785    0.00%

GUB cover cuts applied:  28
Clique cuts applied:  25
Cover cuts applied:  758
Implied bound cuts applied:  7
Flow cuts applied:  384
Mixed integer rounding cuts applied:  84
Zero-half cuts applied:  174
Gomory fractional cuts applied:  5

Root node processing (before b&c):
  Real time             =   11.16 sec. (4393.59 ticks)
Parallel b&c, 4 threads:
  Real time             =   57.55 sec. (26592.93 ticks)
  Sync time (average)   =    5.86 sec.
  Wait time (average)   =    0.00 sec.
                          ------------
Total (root+branch&cut) =   68.70 sec. (30986.52 ticks)
   -     -
 2|3|  1| |1 1
 -   -       -
|    3| |2| | |
     -
|   |3  |2| | |
 -   - -   -
 3|          2|
 -   -     - -
|   |3|  1|
   -   -   -
| |  0  |1  |
   -   -   -
|2 3| |3  |2 0
 - -   - -

100

101

102

103

104

105

106

107

108

109

110

Tried aggregator 5 times.

MIP Presolve eliminated 15467 rows and 2645 columns.

MIP Presolve modified 36628 coefficients.

Aggregator did 16957 substitutions.

Reduced MIP has 28421 rows, 17654 columns, and 84682 nonzeros.

Reduced MIP has 13739 binaries, 3915 generals, 0 SOSs, and 0 indicators.

Presolve time = 0.59 sec. (334.22 ticks)

Probing fixed 4272 vars, tightened 1764 bounds.

Probing changed sense of 3 constraints.

Probing time = 1.98 sec. (381.40 ticks)

Cover probing fixed 0 vars, tightened 8 bounds.

Tried aggregator 3 times.

MIP Presolve eliminated 10986 rows and 6579 columns.

MIP Presolve modified 5755 coefficients.

Aggregator did 2635 substitutions.

Reduced MIP has 14800 rows, 8440 columns, and 42921 nonzeros.

Reduced MIP has 5666 binaries, 2774 generals, 0 SOSs, and 0 indicators.

Presolve time = 0.13 sec. (129.00 ticks)

Probing time = 0.08 sec. (12.19 ticks)

Tried aggregator 1 time.

MIP Presolve eliminated 1573 rows and 944 columns.

Reduced MIP has 13227 rows, 7496 columns, and 38740 nonzeros.

Reduced MIP has 5028 binaries, 2468 generals, 0 SOSs, and 0 indicators.

Presolve time = 0.03 sec. (18.68 ticks)

Probing time = 0.05 sec. (7.86 ticks)

Clique table members: 33932.

MIP emphasis: balance optimality and feasibility.

MIP search method: dynamic search.

Parallel mode: deterministic, using up to 4 threads.

Root relaxation solution time = 0.22 sec. (274.85 ticks)

Nodes Cuts/

Node Left Objective IInf Best Integer Best Bound ItCnt Gap

0 0 0.0000 1857 0.0000 3767

0 0 0.0000 3096 Cuts: 1249 5833

0 0 0.0000 2658 Cuts: 533 7307

0 0 0.0000 2459 Cuts: 351 7947

0 0 0.0000 2214 Cuts: 197 8257

0 0 0.0000 2047 Cuts: 42 8477

0 0 0.0000 2033 Cuts: 58 8531

0 0 0.0000 2020 Cuts: 58 8592

0 0 0.0000 1991 Cuts: 24 8629

0 0 0.0000 1984 Cuts: 7 8662

0 2 0.0000 1984 0.0000 8662

Elapsed time = 11.20 sec. (4434.32 ticks, tree = 0.01 MB, solutions = 0)

5 7 0.0000 1582 0.0000 10209

11 13 0.0000 2542 0.0000 14980

43 45 0.0000 2576 0.0000 25907

113 93 0.0000 2801 0.0000 35995

259 183 0.0000 607 0.0000 39929

335 213 0.0000 1703 0.0000 43877

547 235 0.0000 1255 0.0000 48099

700 236 infeasible 0.0000 53520

716 230 infeasible 0.0000 54089

1060 263 0.0000 250 0.0000 67348

Elapsed time = 17.25 sec. (7827.67 ticks, tree = 7.18 MB, solutions = 0)

1409 266 infeasible 0.0000 81492

1776 260 infeasible 0.0000 97601

2184 267 infeasible 0.0000 115041

2573 274 infeasible 0.0000 130798

3018 277 infeasible 0.0000 145637

3263 245 infeasible 0.0000 156266

3685 283 infeasible 0.0000 171075

4175 263 infeasible 0.0000 186666

4727 273 0.0000 320 0.0000 201106

4777 273 0.0000 2125 0.0000 211968

Elapsed time = 49.45 sec. (22799.02 ticks, tree = 24.51 MB, solutions = 0)

4848 289 0.0000 1670 0.0000 226062

5209 212 0.0000 821 0.0000 248000

5264 157 infeasible 0.0000 265899

5339 82 infeasible 0.0000 284511

5374 49 0.0000 524 -0.0000 306065

5391 39 0.0000 1265 -0.0000 320383

5411 27 infeasible -0.0000 332759

5437 23 infeasible -0.0000 344390

* 5457+ 17 0.0000 -0.0000 346785 0.00%

GUB cover cuts applied: 28

Clique cuts applied: 25

Cover cuts applied: 758

Implied bound cuts applied: 7

Flow cuts applied: 384

Mixed integer rounding cuts applied: 84

Zero-half cuts applied: 174

Gomory fractional cuts applied: 5

Root node processing (before b&c):

Real time = 11.16 sec. (4393.59 ticks)

Parallel b&c, 4 threads:

Real time = 57.55 sec. (26592.93 ticks)

Sync time (average) = 5.86 sec.

Wait time (average) = 0.00 sec.

------------

Total (root+branch&cut) = 68.70 sec. (30986.52 ticks)

- -

2|3| 1| |1 1

- - -

| 3| |2| | |

| |3 |2| | |

- - - -

3| 2|

- - - -

| |3| 1|

- - -

| | 0 |1 |

- - -

|2 3| |3 |2 0

- - - -

Printing may be a little misleading at first sight but it shows the correct answer.

Also, this solution is prone to numerical instability.