Conway's game of life simulates the behavior of a colony of living beings represented by pixels in a plane using a set of simple rules. It is amazing the complex interactions that result from such simple implementations. The Langton Ant is another such creation (for those interested, there is a not too simple Project Euler problem involving this). These examples of what we might generally interpret as "artificial life" are interesting, but can ideas like this be used to solve real problems?
Genetic Algorithms (GA for short henceforth) emulate the propagation of life in an attempt to leverage randomness and simple rules of reproduction and mutation to solve difficult problems more easily in some cases than more traditional techniques. These belong to a class of what are commonly referred to as Randomized Algorithms, which includes mathematical approximations of other physical phenomena as simulated annealing for example ("annealing" is the science of gradually cooling glass from its molten state when it is blown, to room temperature without fracturing it in the process).
Genetic Algorithms have various applications in solving complex dynamic optimization problems, etc, but to illustrate the principles and have some fun in the process, we use these to solve simple Sudoku puzzles in our implementation below. As almost everyone knows, Sudoku is a mathematical puzzle whereby one is presented with a 9x9 grid of cells, some with numbers filled in. This grid is overlaid by a 3x3 grid of cells where each of the 3x3 grids is composed of the smaller 3x3 cells. The goal is to fill in the entire grid with numbers subject to the condition that every 3x3 grid uses the numbers 1 through 9 only once, with the same condition applying to each row and each column of the original 9x9 grid as well. There are thus a total of 9+9+9=27 conditions to the problem, and we can solve this effectively through constrained optimization or through creative use of GAs.
One issue with GAs is that given they are a randomized process, the time to convergence to a complete solution cannot be guaranteed. We see this with our implementation - simple Sudoku puzzles are solved easily by our program, but for the more devious or "evil" level ones, the program thrashes around trying to find a solution, but convergence is very slow. Performance can be improved a lot by building a hybrid method - i.e. employ mathematical or other logical heuristics to fill in some other cells in the harder problems before handing it off to the GAs to solve, to ensure faster convergence, but we do not do that here because our focus is on building a sample implementation of GAs, not boiling the ocean. So, let's stick to the idea at hand, and get right to solving the problem.
How to solve Sudoku with GAs? Here is the approach we will use:
Before going too far though, we first have to explain what mutation and crossover are. Mutation is a small change in the chromosome. Some small subset of bits can be changed arbitrarily - i.e. some digits may be randomly changed to other digits in one or more locations in the chromosome. Crossover is a simulation of sexual reproduction. We take two chromosomes (call them "parents") and generate two off-spring from them by slicing the parents into parts and splicing these parts together keeping in mind the required length constraints of the generated off-spring. Recall from (5) above that we only allow the top 5% of the gene pool to reproduce. What % of the gene pool is newly randomly generated, vs. generated by mutation or by crossover is a design choice, and some tweaking here with a view to optimizing run time, could greatly enhance the performance of the GA implementation.
Code follows along with a sample input and output file. A couple of optimizations can be made to the code below to make it even more efficient. For instance, each new chromosome can be generated by first filling in the empty spaces with the smallest number of options available, then moving on to the harder empty cells, just like humans solve the puzzle. This is left as an exercise to the reader :-)
import os, sys, operator;
from itertools import *;
from random import *;
def getCol(x,n): # get column n from grid x
r=[];
r+=[i[n] for i in x];
return r;
def mkInt(c): # convert all numbers into ints in a list
r=[];
for i in c: r+=[[int(j) if j!='x' else j for j in i]];
return r;
def genConstraints(x):
# collect the constraints now
constraints=[i for i in x]; # horizontal constraints
for i in range(9): constraints+=[getCol(x,i)]; # vertical constraints
# finally, the block constraints
b=[[0,1,2],[3,4,5],[6,7,8]];
for i in product([0,1,2],repeat=2):
p,q=i;
#print "p,q => ",p,q;
c=[];
for i in b[p]:
t=[];
for j in b[q]: t+=[x[i][j]];
#print t;
c+=t;
constraints+=[c];
# now we have all 27 constraints: 9 horizontal, 9 vertical, 9 block
return mkInt(constraints);
def mkGrid(s): # convert a puzzle string to a grid
r=[];
for i in range(9): r+=[s.split(",")[i*9:(i+1)*9]];
return r;
def grabDigit(x): # get a random digit from a list of digits
shuffle(x);
return x[0];
def genChr(x): # generate a chromosome
r=[];
for j in range(9):
nums=[int(k) for k in x[j] if k!='x'];
missing_nums=list(set(range(1,10)).difference(set(nums)));
t=[];
for i in range(9):
if x[j][i]!='x': t+=[int(x[j][i])]
else:
m=list(set(fx[(j,i)]).intersection(set(missing_nums)));
if len(m)==0: m=fx[(j,i)];
a=grabDigit(m);
t+=[a];
if a in missing_nums: missing_nums.remove(a);
r+=[t];
return r;
def validChr(x): # is it a valid chromosome? not used for efficiency
#for i in range(9):
# p=getCol(x,i);
# if sum(p)!=45: return False;
# if list(set(p))!=9: return False;
return True;
def scoreChr(x): # score the chromosome
t=genConstraints(x);
score=0;
for i in t:
if len(list(set(i)))==9: score+=1;
#if sum(i)==45: score+=1;
return score;
def printGrid(x): # print the grid
for i in x:
print " ".join(map(str,i));
return;
def xover(a,b): # cross-over two chromosomes, get two descendants
splice_loc=randint(1,7);
ta=a[:splice_loc]+b[splice_loc:];
tb=b[:splice_loc]+a[splice_loc:];
return (ta,tb);
def xover2(a,b): # another method for cross-over
speed=randint(1,6);
ta,tb=a[:],b[:];
for j in range(speed):
x_loc=randint(1,7);
t=ta[x_loc];
ta[x_loc]=tb[x_loc];
tb[x_loc]=t;
return (ta,tb);
def mutate(a): # mutation
speed=randint(1,6); # how many mutations in one go
for i in range(speed):
gene_seg=randint(0,8); # which gene segment to change
x=genChr(g);
a[gene_seg]=x[gene_seg];
return a;
p=open("s4.txt","r"); # open the puzzle file
pc=p.readlines(); # read it
pc=[i.replace("\n","") for i in pc];
pstr=",".join([i for i in pc]);
g=mkGrid(pstr);
mx=genConstraints(g);
# collect all the constraints by cell...
# find which numbers might fit in each cell. use it for
b=[[0,1,2],[3,4,5],[6,7,8]];
cx,fx={},{};
for i in range(0,9):
for j in range(0,9):
for k in range(0,3):
if i in b[k]: x=k;
if j in b[k]: y=k;
block_id=x*3+y;
cx[(i,j)]=(i,9+j,18+block_id);
fx[(i,j)]=list(set(range(1,10)).difference(set(mx[i]).union(set(mx[9+j])).union(set(mx[18+block_id]))));
# now fx is a dict that contains all possible values for each empty cell
t=fx.values();
for i in range(len(t)):
if len(t[i])==1:
posx,posy=fx.keys()[i];
if g[posx][posy]=='x': g[posx][posy]=t[i][0];
# all constraint strings must have all digits from 1 to 9 just once
# all constraint strings must add up to 9*10/2=45
poolsz=5000; # size of gene pool
gene_pool=[];
while len(gene_pool)!=poolsz: # generate chromosomes to fill pool
seed();
t=genChr(g);
while not validChr(t): t=genChr(g);
if t not in gene_pool: gene_pool+=[t];
seen=gene_pool[:]; # keeps track of all seen chromosomes
# set the parameters for new generations
retain_pct=0.05;
mutate_pct=0.35;
xover_pct=0.55;
new_dna=1-retain_pct-mutate_pct-xover_pct;
generation,solved=0,0;
scores=[]; # store best scores in each generation
# solve the puzzle
while solved==0:
genes=[];
for i in gene_pool: genes+=[(i,scoreChr(i))];
genes=sorted(genes,key=operator.itemgetter(1),reverse=True);
if genes[0][1]==27:
print "puzzle solved";
printGrid(genes[0][0]);
solved=1;
else:
if len(scores)>3 and len(list(set(scores[-3:])))==1:
seed();
jumpahead(randint(1,100000));
poolsz=randint(poolsz,poolsz+500);
retain_pct=0.1;
mutate_pct=round(uniform(0.25,0.3),2);
xover_pct=0.95-mutate_pct;
new_dna=1-retain_pct-mutate_pct-xover_pct;
print "generation: ",generation,"pool size: ",len(genes)," max score: ",genes[0][1];
#printGrid(genes[0][0]);
scores+=[genes[0][1]];
generation+=1;
seed();
chr_retained=int(poolsz*retain_pct);
new_pool=[j[0] for j in genes[:chr_retained]]; # 250 genes to retain
for k in range(int(xover_pct*poolsz/2.)): # generate 450 cross-over genes
seed();
a,b=new_pool[randint(0,len(new_pool)-1)][:],new_pool[randint(0,len(new_pool)-1)][:];
toss=randint(0,1);
if toss==0: ta,tb=xover(a,b);
else: ta,tb=xover2(a,b);
if ta not in seen: new_pool+=[ta];
if tb not in seen: new_pool+=[tb];
for k in range(int(mutate_pct*poolsz)): # generate 150 mutated genes
seed();
a=new_pool[randint(0,len(new_pool)-1)][:];
ta=mutate(a);
if ta not in seen: new_pool+=[ta];
while len(new_pool)<poolsz: # new genetic material
seed();
ta=genChr(g);
if ta not in seen: new_pool+=[ta];
gene_pool=new_pool[:];
seen+=gene_pool[:];
sys.exit(0);
Input File:
7,9,x,x,x,x,3,x,x
x,x,x,x,x,6,9,x,x
8,x,x,x,3,x,x,7,6
x,x,x,x,x,5,x,x,2
x,x,5,4,1,8,7,x,x
4,x,x,7,x,x,x,x,x
6,1,x,x,9,x,x,x,8
x,x,2,3,x,x,x,x,x
x,x,9,x,x,x,x,5,4
Output File:
generation: 0 pool size: 5000 max score: 12
generation: 1 pool size: 5000 max score: 15
generation: 2 pool size: 5000 max score: 17
generation: 3 pool size: 5000 max score: 19
generation: 4 pool size: 5000 max score: 22
generation: 5 pool size: 5000 max score: 25
generation: 6 pool size: 5000 max score: 25
generation: 7 pool size: 5000 max score: 25
puzzle solved
7 9 6 8 5 4 3 2 1
2 4 3 1 7 6 9 8 5
8 5 1 2 3 9 4 7 6
1 3 7 9 6 5 8 4 2
9 2 5 4 1 8 7 6 3
4 6 8 7 2 3 5 1 9
6 1 4 5 9 7 2 3 8
5 8 2 3 4 1 6 9 7
3 7 9 6 8 2 1 5 4
Genetic Algorithms (GA for short henceforth) emulate the propagation of life in an attempt to leverage randomness and simple rules of reproduction and mutation to solve difficult problems more easily in some cases than more traditional techniques. These belong to a class of what are commonly referred to as Randomized Algorithms, which includes mathematical approximations of other physical phenomena as simulated annealing for example ("annealing" is the science of gradually cooling glass from its molten state when it is blown, to room temperature without fracturing it in the process).
Genetic Algorithms have various applications in solving complex dynamic optimization problems, etc, but to illustrate the principles and have some fun in the process, we use these to solve simple Sudoku puzzles in our implementation below. As almost everyone knows, Sudoku is a mathematical puzzle whereby one is presented with a 9x9 grid of cells, some with numbers filled in. This grid is overlaid by a 3x3 grid of cells where each of the 3x3 grids is composed of the smaller 3x3 cells. The goal is to fill in the entire grid with numbers subject to the condition that every 3x3 grid uses the numbers 1 through 9 only once, with the same condition applying to each row and each column of the original 9x9 grid as well. There are thus a total of 9+9+9=27 conditions to the problem, and we can solve this effectively through constrained optimization or through creative use of GAs.
One issue with GAs is that given they are a randomized process, the time to convergence to a complete solution cannot be guaranteed. We see this with our implementation - simple Sudoku puzzles are solved easily by our program, but for the more devious or "evil" level ones, the program thrashes around trying to find a solution, but convergence is very slow. Performance can be improved a lot by building a hybrid method - i.e. employ mathematical or other logical heuristics to fill in some other cells in the harder problems before handing it off to the GAs to solve, to ensure faster convergence, but we do not do that here because our focus is on building a sample implementation of GAs, not boiling the ocean. So, let's stick to the idea at hand, and get right to solving the problem.
How to solve Sudoku with GAs? Here is the approach we will use:
- generate a "board" in string form. this is a string that is 81 characters long, with numbers where there are numbers, and "x" elsewhere to denote empty spaces. a correctly solved board should have all digits from 1 to 9 exactly 9 times, so the sum of the digits in the string in a solved board must always equal (1+2+...+9)*9 =45*9=405.
- "unfold" the string to gather all constraints into a sequential string - this string now has 81+81+81 characters (the first 81 indicates the row constraints, the second 81, the column constraints, and the last 81, the block constraints) for a total of 243 characters.
- generate random "chromosomes", let's say 1000 at a time for our base gene pool. these chromosomes fill in the "x" locations from our base gene in (2) with numbers from 1 through 9, randomly. collect these chromosomes. collect only those chromosomes that meet the requirements of a proper solution - let's say the sum of all digits within should be 405. more stringent conditions here will lead to better chromosomes that take longer to compute. lenient conditions here will generate worse chromosomes that are computed quickly.
- rank the chromosomes in the gene pool based on how well they approximate what might be a correct solution to the Sudoku problem - are all digits used 9 times each, none more, none less? is the grid sum close to 405? what conditions we use here in the "fitness" constraint depend on what conditions we used in (3) above in the gene pool generation process.
- after the gene pool is created and chromosomes ranked for fitness, use Darwin's idea of "survival of the fittest" to keep only the top 25% of the pool, discarding the rest. then, use mutation and crossover to generate the next version of the pool.
- at each iteration of the pool, check to see if any chromosome meets all the criteria of the correct solution. if it does, then that chromosome is the correct solution to the puzzle and we stop there. if none does, keep iterating with more generations of the gene pool till a complete solution is found.
- print the final solution and exit.
Before going too far though, we first have to explain what mutation and crossover are. Mutation is a small change in the chromosome. Some small subset of bits can be changed arbitrarily - i.e. some digits may be randomly changed to other digits in one or more locations in the chromosome. Crossover is a simulation of sexual reproduction. We take two chromosomes (call them "parents") and generate two off-spring from them by slicing the parents into parts and splicing these parts together keeping in mind the required length constraints of the generated off-spring. Recall from (5) above that we only allow the top 5% of the gene pool to reproduce. What % of the gene pool is newly randomly generated, vs. generated by mutation or by crossover is a design choice, and some tweaking here with a view to optimizing run time, could greatly enhance the performance of the GA implementation.
Code follows along with a sample input and output file. A couple of optimizations can be made to the code below to make it even more efficient. For instance, each new chromosome can be generated by first filling in the empty spaces with the smallest number of options available, then moving on to the harder empty cells, just like humans solve the puzzle. This is left as an exercise to the reader :-)
import os, sys, operator;
from itertools import *;
from random import *;
def getCol(x,n): # get column n from grid x
r=[];
r+=[i[n] for i in x];
return r;
def mkInt(c): # convert all numbers into ints in a list
r=[];
for i in c: r+=[[int(j) if j!='x' else j for j in i]];
return r;
def genConstraints(x):
# collect the constraints now
constraints=[i for i in x]; # horizontal constraints
for i in range(9): constraints+=[getCol(x,i)]; # vertical constraints
# finally, the block constraints
b=[[0,1,2],[3,4,5],[6,7,8]];
for i in product([0,1,2],repeat=2):
p,q=i;
#print "p,q => ",p,q;
c=[];
for i in b[p]:
t=[];
for j in b[q]: t+=[x[i][j]];
#print t;
c+=t;
constraints+=[c];
# now we have all 27 constraints: 9 horizontal, 9 vertical, 9 block
return mkInt(constraints);
def mkGrid(s): # convert a puzzle string to a grid
r=[];
for i in range(9): r+=[s.split(",")[i*9:(i+1)*9]];
return r;
def grabDigit(x): # get a random digit from a list of digits
shuffle(x);
return x[0];
def genChr(x): # generate a chromosome
r=[];
for j in range(9):
nums=[int(k) for k in x[j] if k!='x'];
missing_nums=list(set(range(1,10)).difference(set(nums)));
t=[];
for i in range(9):
if x[j][i]!='x': t+=[int(x[j][i])]
else:
m=list(set(fx[(j,i)]).intersection(set(missing_nums)));
if len(m)==0: m=fx[(j,i)];
a=grabDigit(m);
t+=[a];
if a in missing_nums: missing_nums.remove(a);
r+=[t];
return r;
def validChr(x): # is it a valid chromosome? not used for efficiency
#for i in range(9):
# p=getCol(x,i);
# if sum(p)!=45: return False;
# if list(set(p))!=9: return False;
return True;
def scoreChr(x): # score the chromosome
t=genConstraints(x);
score=0;
for i in t:
if len(list(set(i)))==9: score+=1;
#if sum(i)==45: score+=1;
return score;
def printGrid(x): # print the grid
for i in x:
print " ".join(map(str,i));
return;
def xover(a,b): # cross-over two chromosomes, get two descendants
splice_loc=randint(1,7);
ta=a[:splice_loc]+b[splice_loc:];
tb=b[:splice_loc]+a[splice_loc:];
return (ta,tb);
def xover2(a,b): # another method for cross-over
speed=randint(1,6);
ta,tb=a[:],b[:];
for j in range(speed):
x_loc=randint(1,7);
t=ta[x_loc];
ta[x_loc]=tb[x_loc];
tb[x_loc]=t;
return (ta,tb);
def mutate(a): # mutation
speed=randint(1,6); # how many mutations in one go
for i in range(speed):
gene_seg=randint(0,8); # which gene segment to change
x=genChr(g);
a[gene_seg]=x[gene_seg];
return a;
p=open("s4.txt","r"); # open the puzzle file
pc=p.readlines(); # read it
pc=[i.replace("\n","") for i in pc];
pstr=",".join([i for i in pc]);
g=mkGrid(pstr);
mx=genConstraints(g);
# collect all the constraints by cell...
# find which numbers might fit in each cell. use it for
b=[[0,1,2],[3,4,5],[6,7,8]];
cx,fx={},{};
for i in range(0,9):
for j in range(0,9):
for k in range(0,3):
if i in b[k]: x=k;
if j in b[k]: y=k;
block_id=x*3+y;
cx[(i,j)]=(i,9+j,18+block_id);
fx[(i,j)]=list(set(range(1,10)).difference(set(mx[i]).union(set(mx[9+j])).union(set(mx[18+block_id]))));
# now fx is a dict that contains all possible values for each empty cell
t=fx.values();
for i in range(len(t)):
if len(t[i])==1:
posx,posy=fx.keys()[i];
if g[posx][posy]=='x': g[posx][posy]=t[i][0];
# all constraint strings must have all digits from 1 to 9 just once
# all constraint strings must add up to 9*10/2=45
poolsz=5000; # size of gene pool
gene_pool=[];
while len(gene_pool)!=poolsz: # generate chromosomes to fill pool
seed();
t=genChr(g);
while not validChr(t): t=genChr(g);
if t not in gene_pool: gene_pool+=[t];
seen=gene_pool[:]; # keeps track of all seen chromosomes
# set the parameters for new generations
retain_pct=0.05;
mutate_pct=0.35;
xover_pct=0.55;
new_dna=1-retain_pct-mutate_pct-xover_pct;
generation,solved=0,0;
scores=[]; # store best scores in each generation
# solve the puzzle
while solved==0:
genes=[];
for i in gene_pool: genes+=[(i,scoreChr(i))];
genes=sorted(genes,key=operator.itemgetter(1),reverse=True);
if genes[0][1]==27:
print "puzzle solved";
printGrid(genes[0][0]);
solved=1;
else:
if len(scores)>3 and len(list(set(scores[-3:])))==1:
seed();
jumpahead(randint(1,100000));
poolsz=randint(poolsz,poolsz+500);
retain_pct=0.1;
mutate_pct=round(uniform(0.25,0.3),2);
xover_pct=0.95-mutate_pct;
new_dna=1-retain_pct-mutate_pct-xover_pct;
print "generation: ",generation,"pool size: ",len(genes)," max score: ",genes[0][1];
#printGrid(genes[0][0]);
scores+=[genes[0][1]];
generation+=1;
seed();
chr_retained=int(poolsz*retain_pct);
new_pool=[j[0] for j in genes[:chr_retained]]; # 250 genes to retain
for k in range(int(xover_pct*poolsz/2.)): # generate 450 cross-over genes
seed();
a,b=new_pool[randint(0,len(new_pool)-1)][:],new_pool[randint(0,len(new_pool)-1)][:];
toss=randint(0,1);
if toss==0: ta,tb=xover(a,b);
else: ta,tb=xover2(a,b);
if ta not in seen: new_pool+=[ta];
if tb not in seen: new_pool+=[tb];
for k in range(int(mutate_pct*poolsz)): # generate 150 mutated genes
seed();
a=new_pool[randint(0,len(new_pool)-1)][:];
ta=mutate(a);
if ta not in seen: new_pool+=[ta];
while len(new_pool)<poolsz: # new genetic material
seed();
ta=genChr(g);
if ta not in seen: new_pool+=[ta];
gene_pool=new_pool[:];
seen+=gene_pool[:];
sys.exit(0);
Input File:
7,9,x,x,x,x,3,x,x
x,x,x,x,x,6,9,x,x
8,x,x,x,3,x,x,7,6
x,x,x,x,x,5,x,x,2
x,x,5,4,1,8,7,x,x
4,x,x,7,x,x,x,x,x
6,1,x,x,9,x,x,x,8
x,x,2,3,x,x,x,x,x
x,x,9,x,x,x,x,5,4
Output File:
generation: 0 pool size: 5000 max score: 12
generation: 1 pool size: 5000 max score: 15
generation: 2 pool size: 5000 max score: 17
generation: 3 pool size: 5000 max score: 19
generation: 4 pool size: 5000 max score: 22
generation: 5 pool size: 5000 max score: 25
generation: 6 pool size: 5000 max score: 25
generation: 7 pool size: 5000 max score: 25
puzzle solved
7 9 6 8 5 4 3 2 1
2 4 3 1 7 6 9 8 5
8 5 1 2 3 9 4 7 6
1 3 7 9 6 5 8 4 2
9 2 5 4 1 8 7 6 3
4 6 8 7 2 3 5 1 9
6 1 4 5 9 7 2 3 8
5 8 2 3 4 1 6 9 7
3 7 9 6 8 2 1 5 4
