Sunday, September 14, 2008

Lights Out Puzzle

Problem Intro:
Related discussion:
topcoder Problem:
ZOJ Problem: TBD

可以把这个问题转化为一个graph 问题,用adjancency matrix(dim = row x col) 可以表示出网格(cell)间的邻接关系。
每个网格的状态是一个finite field。如果只有on/off,则finite field的运算均mod 2.
我们在求解中用到的运算有求倒数和求模,所以这里给出了finite field的 invert(), modulate().

// ----------------------------------------------------------------------------
// ----------------------------------------------------------------------------
#define N 10 // maximum vertex #
#define MAX N*N // maximum edge #

// --- global variables ---
int colcount; // integer, number of columns
int rowcount; // integer, number of rows
int mod; // integer, number of states of a tile (flip: binary, mod = 2)

int mat[MAX][MAX]; // integer[i][j]
int cols[MAX];
int m; // count of rows of the matrix
int n; // count of columns of the matrix
int np; // count of columns of the enlarged matrix
int r; // minimum rank of the matrix
int maxr; // maximum rank of the matrix

// integer[row][col], current states of tiles
int cells[N][N] = {
{1,1,1,1,1,1},
{1,0,2,0,1,0},
{1,0,0,0,1,0},
{1,0,0,0,1,0},
{1,1,1,1,1,0},
{1,2,2,1,1,1}
};

// --- finite field algebra solver
int modulate(int x) {
// returns z such that 0 <= z < x ="=">= 0) return x % mod;
x = (-x) % mod;
if (x == 0) return 0;
return mod - x;
}

int gcd(int x, int y) { // call when: x >= 0 and y >= 0
if (y == 0) return x;
if (x == y) return x;
if (x > y) x = x % y; // x < y
while (x > 0) {
y = y % x; // y < x
if (y == 0) return x;
x = x % y; // x < y
}
return y;
}

int invert(int value) { // call when: 0 <= value < mod
// returns z such that value * z == 1 (mod mod), or 0 if no such z
if (value <= 1) return value;
int seed = gcd(value,mod);
if (seed != 1) return 0;
int a = 1, b = 0, x = value; // invar: a * value + b * mod == x
int c = 0, d = 1, y = mod; // invar: c * value + d * mod == y
while (x > 1) {
int tmp = floor(y / x + 0.0);
y -= x * tmp;
c -= a * tmp;
d -= b * tmp;
tmp = a; a = c; c = tmp;
tmp = b; b = d; d = tmp;
tmp = x; x = y; y = tmp;
}
return a;
}
// ----------------------------------------------------------------------------
// ----------------------------------------------------------------------------


在mathWorld中讲解的很清楚,给出初始状态和目的状态,最终找可行路径归结为求解 linear equation.
求解有三种可能:
  1. 无解
  2. 有唯一解 (矩阵满秩)
  3. 有多个解
由于这是一个finite field 元素的 linear equation,我们用了相应的gaussian elimination method, O(N^4)

几个子函数:
getmat(...), setmat(...) 用来存取adj matrix, 以及作gaussian elimination操作matrix


// ----------------------------------------------------------------------------
// ----------------------------------------------------------------------------
int getmat(int i,int j) { return mat[i][cols[j]]; }
void setmat(int i, int j, int val) { mat[i][cols[j]] = modulate(val); }
// ----------------------------------------------------------------------------
// ----------------------------------------------------------------------------


initMatrix() 根据问题创建adj matrix,不同问题需要改写这一部分
cols[] 用来记录消去过程中,列之间的交换


// ----------------------------------------------------------------------------
// ----------------------------------------------------------------------------
void initMatrix() {
maxr = m < n ? m : n;
for(int i = 0; i < m; ++i)
for(int j = 0; j < n; ++j)
mat[i][j] = 0;
for (int row = 0; row < rowcount; row++)
for (int col = 0; col < colcount; col++){
int i = row * colcount + col;

mat[i][i] = 1;
// 4-connected
if (col > 0) mat[i][i - 1] = 1;
if (row > 0) mat[i][i - colcount] = 1;
if (col < colcount - 1) mat[i][i + 1] = 1;
if (row < rowcount - 1) mat[i][i + colcount] = 1;

// diagonal connection
//if (col > 0 && row > 0) mat[i][i-colcount-1] = 1;
//if (col > 0 && row < rowcount - 1) mat[i][i+colcount-1] = 1;
//if (col < colcount-1 && row > 0) mat[i][i-colcount+1] = 1;
//if (col < colcount-1 && row < rowcount-1) mat[i][i+colcount+1] = 1;
}

for (int j = 0; j < np; j++) cols[j] = j;
return;
}
// ----------------------------------------------------------------------------
// ----------------------------------------------------------------------------


sweep(), sweepStep(), doBasicSweep()构成O(N^4) 高斯消去法
// ----------------------------------------------------------------------------
// ----------------------------------------------------------------------------
void doBasicSweep(int pivoti, int pivotj) {

if(r != pivoti) {
for(int j = 0; j < np; j++) {
int tmp = getmat(r, j);
setmat(r,j, getmat(pivoti,j));
setmat(pivoti,j, tmp);
}
}

if(r != pivotj) {
int tmp = cols[r];
cols[r] = cols[pivotj];
cols[pivotj] = tmp;
}

for (int i = 0; i < m; i++) {
if (i != r) {
int air = getmat(i,r);
if (air != 0)
for (int j = r; j < np; j++)
setmat(i,j, getmat(i,j) - getmat(r,j) * air);
}
}
}

int sweepStep() {
int i;
int j;
bool finished = true;
for (j = r; j < n; j++) {
for (i = r; i < m; i++) {
int aij = getmat(i,j);
if (aij != 0) finished = false;
int inv = invert(aij);
if (inv != 0) {
for (int jj = r; jj < np; jj++)
setmat(i,jj, getmat(i,jj) * inv);
doBasicSweep(i,j);
return 1;
}
}
}
if (finished) { // we have: 0x = b (every matrix element is 0)
maxr = r; // rank(A) == maxr
for (j = n; j < np; j++)
for (i = r; i < m; i++)
if (getmat(i,j) != 0) {
//std::cout << "no solution\n";
return -1; // no solution since b != 0
}
return 1; // 0x = 0 has solutions including x = 0
}
return -1; // failed in finding a solution
}

int sweep() {
for (r = 0; r < maxr; r++) {
int status = sweepStep();
if(status < 0) return -1;
if(status == 0) return 0;
if (r == maxr) break;
}
return find_min_flip();
}// ----------------------------------------------------------------------------
// ----------------------------------------------------------------------------


solve 是整个问题入口,初始化各变量,调用initMatrix生成adj matrix,并填入目标状态
mat[i][n] (i = 1...n) --- mat矩阵最后一列

// ----------------------------------------------------------------------------
// ----------------------------------------------------------------------------
int solve(int goal) {

int size = colcount * rowcount;
m = size;
n = size;
np = n + 1;
initMatrix();
for (int row = 0; row < rowcount; row++)
for (int col = 0; col < colcount; col++)
mat[row * colcount + col][n] = modulate(goal - cells[row][col]);
return sweep();
}// ----------------------------------------------------------------------------
// ----------------------------------------------------------------------------

find_min_flip根据矩阵的秩r,
如果满秩,直接求出解,否则
可知有n-r个自由变量,构成mod^(n-r)种组合
e.g., mod = 2, n-r = 4, 则有2^4 = 16种组合。
对每一种组合代入原方程组,求解其它变量值,对于所有可能,再统计其中步数最少的输出
这里用了iterative backtracking, 两重循环,
外重是当前处理第k个自由变量,
内重是当前处理的自由变量的第(0...mod-1)种取值可能

// ----------------------------------------------------------------------------
// ----------------------------------------------------------------------------
int find_min_flip() {

vector sol(n,-1);
vector best(n, -1);

if(r == n) {
int flips = 0;
for(int j = 0; j < n; ++j) {
if(j > 0 && j % colcount == 0) cout << endl;
if(getmat(j,n) > 0) flips += getmat(j,n);
cout << getmat(j,n) << " ";
}
cout << endl;
return flips;
}
else if(r < n) {
int minflips = n;
int k = r;
while (k >= r) {
while(sol[k] < mod) {
sol[k] = sol[k] + 1;
if(k == n-1) {
for(int j = r-1; j >= 0; --j) {
int sum = 0;
for(int jj = j+1; jj < n; ++jj)
sum += sol[jj] * getmat(j,jj);
sol[j] = modulate(getmat(j,n) - sum);
}

int fl = 0;
for(int j = 0; j < n; ++j) {
if(sol[j] == 1) fl++;
//cout << sol[j] << " ";
}
//cout << endl;
if(fl < minflips) {
minflips = fl;
best = sol;
}
}
else k = k + 1;
}
sol[k] = -1;
k = k - 1;
}
for(int j = 0; j < n; ++j) {
if(j > 0 && j % colcount == 0) cout << endl;
cout << best[j] << " ";
}
cout << endl;
return minflips;
}
return -1;
}// ----------------------------------------------------------------------------
// ----------------------------------------------------------------------------

test in main(), 6x6 grid, on/off two states, goal state 0 or 1, initial state in cells. 输出的两个矩阵分别对应目标状态为全关和全开,矩阵中为“1”的是应该按的开关。
// ----------------------------------------------------------------------------
// ----------------------------------------------------------------------------
int main() {
mod = 2;
rowcount = colcount = 6;
for(int i = 0; i < mod ; ++i){
solve(i); cout << endl;
}
}// ----------------------------------------------------------------------------
// ----------------------------------------------------------------------------

outputs:
1 1 1 1 0 0
1 0 0 1 0 1
1 0 0 0 1 1
1 1 0 0 1 1
0 0 1 1 0 1
0 1 1 1 0 0

0 1 0 0 0 1
1 1 1 0 1 1
0 1 1 1 0 0
0 0 1 1 0 0
0 1 0 0 1 1
1 1 0 0 0 1


For the topCoder problem: LightedPanels. We only need to modify the initMatrix() to support 6-connected adj matrix, then ...

int LightedPanels::minTouch(vector  board) {

colcount = board[0].length();
rowcount = board.size();
imgcount = 2;

for (int col = 0; col < colcount; col++)
for (int row = 0; row < rowcount; row++) {
cells[row][col] = board[row][col] == '.' ? 0 : 1;
}

return solve(1);
}