discrete-optimization-coursera/puzzle/magicSquare/Program.cs

using System;
using System.Linq;
using Google.OrTools.ConstraintSolver;

namespace magicSquare
{
	class Program
	{
		static void Main(string[] args)
		{
			var n = int.Parse(args[0]);

			var solution = PatternSolve(n);

			PrintSolution(n, solution);
		}

		private static void PrintSolution(int n, int[,] solution)
		{
			Console.WriteLine(n);
			var rowValues =
				from row in Enumerable.Range(0, n)
				let rowVals =
					(from col in Enumerable.Range(0, n)
					 select solution[n - 1 - row, col])
				select rowVals;
			foreach (var rowValue in rowValues)
			{
				Console.WriteLine(string.Join("\t", rowValue.ToArray()));
			}
		}

		public static int[,] LSSolve(SquareMatrix squareMatrix)
		{
			// lock down diagonals

			while (!squareMatrix.IsFeasible())
			{
				// select item
				var maxColumnDeltaPairs = squareMatrix.MaxColumnDeltaPairs().ToList();
				var maxRowDeltaPairs = squareMatrix.MaxRowDeltaPairs().ToList();

				var swaps =
					// build neighborhood
					from cols in maxColumnDeltaPairs
					from rows in maxRowDeltaPairs
					let coord1 = new Tuple<int, int>(cols.Item1, rows.Item1)
					let coord2 = new Tuple<int, int>(cols.Item2, rows.Item2)
					// ignore items on the diagonal
					where !squareMatrix.IsDiagonal(coord1) && !squareMatrix.IsDiagonal(coord2)
					select new {coord1, coord2};
				// order descending by sum of delta difference
				swaps = swaps.OrderByDescending(arg => squareMatrix.AbsDeltaDiff(arg.coord1, arg.coord2));

				// swap
				var swap = swaps.First();
				squareMatrix.Swap(swap.coord1, swap.coord2);
			}
			return null;
		}

		public static int[,] PatternSolve(int n)
		{
			var solution = new int[n, n];

			var nover2Floor = Convert.ToInt32(Math.Floor(n / 2.0));
			for (var row = 1; row <= n; row++)
			{
				for (var col = 1; col <= n; col++)
				{
					var v = n * ((row + col - 1 + nover2Floor) % n) + ((row + 2 * col - 2) % n) + 1;
					solution[row - 1, col - 1] = v;
				}
			}
			return solution;
		}

		public static int[,] CPSolve(int n)
		{

			var n2 = n * n;
			var magicNumber = (n * (n * n + 1)) / 2;

			var solver = new Solver("magicSquares");

			var vars = solver.MakeIntVarMatrix(n, n, 1, n2, "hi");

			var varsFlattened = vars.Flatten();

			var cols =
				from i in Enumerable.Range(0, n)
				let col =
					(from j in Enumerable.Range(0, n)
					 select vars[i, j]).ToArray()
				select col;

			var rows =
				from j in Enumerable.Range(0, n)
				let row =
					(from i in Enumerable.Range(0, n)
					 select vars[i, j]).ToArray()
				select row;

			var diag1 =
				(from i in Enumerable.Range(0, n)
				 select vars[i, i]).ToArray();

			var diag2 =
				(from i in Enumerable.Range(0, n)
				 select vars[i, n - 1 - i]).ToArray();

			foreach (var pieces in
					cols
					.Union(rows)
					.Union(new[] { diag1, diag2 }))
			{
				solver.Add(pieces.Sum() == magicNumber);
			}

			solver.Add(varsFlattened.AllDifferent());

			// TODO: symetry break
			for (var i = 1; i < n; i++)
				solver.Add(diag1[i - 1] == diag1[i] - 1);

			// middle number
			solver.Add(vars[(n - 1) / 2, (n - 1) / 2] == (n * n + 1) / 2);
			solver.Add(vars[(n - 1) / 2, 0] == 1);
			solver.Add(vars[(n - 1) / 2, n - 1] == n * n);

			//var nover2floor = Convert.ToInt32(Math.Floor(n / 2.0));
			//for (var i = 1; i <= n; i++)
			//{
			//	for (var j = 1; j <= n; j++)
			//	{
			//		var v = n * ((i + j - 1 + nover2floor) % n) + ((i + 2 * j - 2) % n) + 1;
			//		solver.Add(vars[i - 1, j - 1] == v);
			//	}
			//}

			var db = solver.MakePhase(varsFlattened, Solver.CHOOSE_FIRST_UNBOUND, Solver.ASSIGN_CENTER_VALUE);

			solver.NewSearch(db);

			var solution = new int[n,n];

			if (solver.NextSolution())
			{
				for (var i = 0; i < n; i++)
				{
					for (var j = 0; j < n; j++)
					{
						solution[i, j] = Convert.ToInt32(vars[i, j].Value());
					}
				}
			}
			solver.EndSearch();

			return solution;
		}
	}
}