## Wednesday, October 26, 2011

### Using Scipy To Do Inline FFTs

INSTALLING
You need to do a sudo apt-get install python2.x-dev to get inline to work (I did 2.7). When scipy.weave.inline() runs, there is a call to the C header 'Python.h' which is contained in the dev package.

If there are issues with building the C code, i.e. cryptic errors like "Error Code -1 returned", adding verbose=2 to the arguments of the inline() call will show every line output as it's being built - this is very useful for figuring out why the inline build is breaking!

BACKGROUND
This code pretty much builds on my previous blog post about FFTW. I am a big fan of the speed and simplicity of Python, especially for interacting with the web. However, after my data is fetched from the web, I also want to perform some pretty intense analytics, and even NumPy left something to be desired. After reading this page on heavy computation in Python (http://www.scipy.org/PerformancePython), I found out that inline C can get pretty close to pure C functions, and by interfacing the PyObjects directly, disk I/O can be avoided, which may also save some compute time. And it is also much, much simpler (and faster to code) in Python instead of pure C.

CODES OF NOTE
There are a couple of nice things about inlining these functions. Notice the ['comp', 'num_fft', 'fft_len', 'val'] inside the scipy.weave.inline() call. This exposes those Python objects to the C code, allowing them to be used. I use num_fftlen_fft, and val to set loop sizes and initial values, then after performing the FFTs, the values in out are written to comp using C array indexing. Because the numpy array is declared as complex, the real and imaginary components are combined and written using (re,imag) form.

FINAL THOUGHTS
Here are the statistics I get running two chunks of code against each other, each using dft_plan_many to do 100,000 length 32 FFTs, with FFTW_PLAN set to FFTW_ESTIMATE.

Python
$time ./fftw_plan_many.py real 0m0.738s user 0m0.416s sys 0m0.192s C$ time ./many_dft.o
real 0m0.188s
user 0m0.064s
sys 0m0.088s

Not bad to do 100,000 FFTs, albeit small ones (32 bin) in < 1 sec. However C (like Sebulba) always wins in numerics.

I am planning to continue this investigation thread at a later date. Once my nVidia development computer gets up and running, I will likely attempt to do CUDAFFT with this same interface, and compare that to the current findings.

SIMPLE CODE
#!/usr/bin/python
from scipy.weave import inline

inline( 'printf("Hello World");',
compiler='gcc')
print ""

PYTHON FFTW PLAN MANY CODE
#!/usr/bin/python
import numpy as np
from scipy.weave import inline
from numpy.distutils.system_info import get_info
fftw_plan_many_dft = \
"""
fftw_complex* in;
fftw_complex* out;
fftw_plan plan;

in = (fftw_complex *) fftw_malloc(sizeof(fftw_complex)*fft_len*num_fft);
out = (fftw_complex *) fftw_malloc(sizeof(fftw_complex)*fft_len*num_fft);

plan = fftw_plan_many_dft( 1,
&fft_len,
num_fft,
in,
NULL,
num_fft,
1,
out,
NULL,
num_fft,
1,
FFTW_FORWARD,
FFTW_ESTIMATE );

for (int i = 0; i < num_fft; i++){
val+=1;
for (int j = 0; j < fft_len; j++) {
in[j+i*fft_len][0] = val;
in[j+i*fft_len][1] = val ;
}
}
fftw_execute(plan);
for (int i =0; i < num_fft; i++){
for (int j = 0; j < fft_len; j++){
comp[j+i*fft_len] = std::complex<double>(out[j+i*fft_len][0]/fft_len, out[j+i*fft_len][1]/fft_len);
}
}
fftw_destroy_plan(plan);
"""

fftw_info = get_info( 'fftw3' )
num_fft = 100000
fft_len = 32
val = 0.
comp = np.array(np.zeros(num_fft*fft_len), dtype=complex)
inline( fftw_plan_many_dft,
['comp', 'num_fft', 'fft_len', 'val'],
libraries=['fftw3'],
include_dirs=fftw_info['include_dirs'],
library_dirs=fftw_info['library_dirs'],
compiler='gcc' )
print comp

## Saturday, October 8, 2011

### A Genetic Solution to the N-queens Problem

GETTING IN THE QUEEN'S GENES
Lately, I have become interested in the techniques of genetic algorithms, shortened as GA. After speaking with a friend about an assignment in his algorithms course, I decided that my first GA project was going to be an expandable genetic n-queens solver. I have included the code below.

The N-queens problem is one of the classic problems of computer science. The premise is simple - for any chess board NxN, where N>3, there exists at least one arrangement such that N queens can be placed without conflicting in either a row or a column. More information and a better explanation can be found here - http://en.wikipedia.org/wiki/Nqueens . The reason this problem is popular in algorithms classes around the globe (besides the association with Djikstra) is because for N > 8 or so, the problem has a solution space large enough that brute forcing becomes impractical. This requires the student to think about the parameters of the puzzle and find ways to eliminate "stupid" possibilities,  paring down the number of possibilities to a reasonable volume.

SIMPLE EXAMPLE
Suppose we choose to solve the problem for N =4. This means that the board size is 4^2=16, and the number of queens we can fit inside the board is 4. Picturing an empty chess board, it would look something like this:
[0 0 0 0]
[0 0 0 0]
[0 0 0 0]
[0 0 0 0]

If we were to place the queens randomly around the board:
[1 0 1 0]
[0 0 0 1]
[0 0 0 0]
[0 0 0 1]

Here we can clearly see that the purple queens conflict by row, while the yellow conflict by column. There is also a diagonal conflict between the purple and yellow in the upper right hand corner. Shifting things some more, we could eventually see:
[0 0 1 0]
[1 0 0 0]
[0 0 0 1]
[0 1 0 0]

There are no row, column, or diagonal conflicts. This is one valid solution to the N-queens problem.

CODE PROCESS
Upon searching for existing examples, I found samples of code from http://pyevolve.sourceforge.net/0_6rc1/examples.html , which actually have an example of an N-queens solver (second from the bottom of the page). I didn't like the logic in the sample code, so I decided to create my own flavor for the answer, borrowing from this existing source. I also figured it could be a good way to learn more about genetic programming.

I am a big fan of functional programming, so I have a strong tendency to use (and abuse) the pseudo-functional abilities available in python. Luckily, using map() and filter() tends to speed up a lot of python code, at least in my experience.

The example from sourceforge uses simpler python code, but the logic behind it seems unintuitive to me. When I first read about the problem, my first thought was that minimizing conflicts would eventually lead to the right answer. I also thought of the board exactly as it appears in the above simple example, whereas the sample code used the range function to create a list of numbers (0-board_size) to look for conflicts. I found it simpler to just use ones and zeroes. The sourceforge code also counts collisions, but then inverts the scoring for some reason, requiring the use of "maximize" int the .setMinimax method. Assuredly, I did do some wonky Python tricks to make my code work, but I will explain myself below.

PYTHON EXPLANATION
The very first step that I took was downloading the PyEvolve package. This can be done either by doing sudo apt-get install python-pyevolve in a Debian-based Linux install, or by downloading the source from http://pyevolve.sourceforge.net/ . Once PyEvolve is installed, everything should be ready to go.

Line-by-Line of queens_func
if sum(genome) < BORDER: return BOARD_SIZE
This line is required to keep the GA from ranking bad solutions with high scores. In this case, the GA will create some solutions with no pieces in a row (or no pieces in any rows!). These solutions should be regarded as genetic failures, and given the worst possible score. Otherwise, a solution of all zeroes could be considered to meet the exit criteria, as there would be no conflicts.

array = map( lambda x: list(genome[x::BORDER]), range(BORDER))
The first abuse of map occurs here. map(a,b) applies a given function a to every element in list b. I can can also do something like map(lambda x: x+1, b). By using lambda x:, I am defining an inline function that will apply x+1 to any input x. Basically x++. This is identical to writing def function(x): return x+1 elsewhere in my file, and doing map(function, b). By doing list(genome[x::BORDER], I am indexing through the list genome in steps of size BORDER, starting at x (where x is the number we are currently on in range(BORDER), and creating a list for each grouping. For example:
genome = [1, 2, 3, 4, 5, 6, 7, 8, 9]
BORDER = 3
range(BORDER) -> [0, 1, 2]
map( lambda x: genome[x::BORDER], range(BORDER)) -> [[1, 4, 7], [2, 4, 8], [3, 6, 9]]
This can also interpreted as a list of columns, where array[num] returns a list of the values in column num. Effectively, this creates a 2D list of lists out of a 1D list, without involving NumPy. I didn't want to include NumPy because I may have to use a similar program on work computer that doesn't have access to NumPy.

cols = map(sum, array)
Nothing too special here - we simply want to sum up the number of pieces in each column of the array. Ideally, this result will be all ones, as that means there is only 1 piece in every column.

rows = map(sum, map(list, zip(*array)))

This is a fun one. zip(*array) uses the splat operator to "unlist" every column in array, then regroup based on position in each sublist. Applying this to the earlier example:
array = [[1, 4, 7], [2, 4, 8], [3, 6,9]]
zip(*array) -> [(1, 2, 3), (4, 5, 6), (7, 8, 9)].
This can be generalized to [(all pos(0)), (all pos(1)), (all pos(2)), ...]. Doing map(list, ) turns each tuple into a list, creating a 2D transpose of the earlier defined array. Summing these will give the total number of queens in each row. Once again the ideal result is all ones.

rows_ints = [] for row in rows: rows_ints.append(int(''.join(map(str, row)), 2))
rows_ints = map( lambda x: int(math.log(x,2)) if x !=0 else x, rows_ints)
For this part, I wanted to turn each row into it's binary equivalent, so [1, 0, 0, 0] would become 1000. I did this because I thought an easy way to test for diagonals would be to see if row neighbors are also power of 2 neighbors, by checking the power of 2 exponent. If they are off by one, then a diagonal conflict is present. A brief example:
[1,0,0] -> 100 = 8 -> 2^3 -> 3
[0,1,0] -> 010 = 4 -> 2^2 -> 2
[0,0,1] -> 001 = 2 -> 2^1 -> 1
The lambda expression actually contains a pseudo if-else statement. int(math.log(x,2) is returned if x isn't equal to zero - otherwise x (which must be zero) is returned unchanged. This is because log(x) will not work if x = 0. Gotta love log rules.

diags = filter( lambda x: x==1, map( lambda i,j: abs(i-j), rows_ints[1:], rows_ints[:-1]))
Working from the inside out, map( lambda i,j: abs(i-j), rows_ints[1:], rows_ints[:-1]) takes a shifted slice from rows_ints and point-by-point subracts the values from the default slice. A brief example:
a = [1,2,3]
a[1:] -> [2,3]
a[:-1] -> [1,2]
map( lambda i,j: abs(i-j), a[1:], a[:-1] ) -> [abs(2-1), abs(3-2)] ->[1,1]
Once this operation is complete, the next step is to filter out all the differences of 1 between rows. As discussed previously, a difference of 1 between neighbor rows indicates that there is a diagonal conflict between them.

return max(map(sum, rows)) + max(map(sum, cols)) - 2 + sum(diags)
Our return score will be the worst count per row, the worst count per column, and the total number of diagonal conflicts. This result is scaled down by -2 because the best result possible for max(map(sum))) for rows is 1, and the same for columns. This would indicate that the worst values for queens per column is 1, and the worst value for queens per row is 1, meaning there are no row/col conflicts. If sum(diags) was then 0, this would be a solution that meets the bestrawscore criteria. If genome.setParams was changed to  genome.setParams(bestrawscore=2, rounddecimal=2), then the -2 could be eliminated from the return statement.

In the queens_init function
[1]*N1 + [0]*N2 will actually create one expanded list that is of length N1+N2, where the first N1 values are 1 and the rest are 0. Pretty self explanatory, but I thought it was neat that python understood using the * operator for expansion on a list, rather than multiplying each value in the list by N1 or N2. Of course, that is the same "feature" that annoyed me in the past after coming to python from MATLAB, but I am beginning to understand why this is the default behavior, and scientific pythonistas need to use NumPy to get closer to MATLAB functionality.

CHANGES IN MAIN
I only made a few other changes to the sample code. I changed the genome.setParams(bestrawscore=0), because my methodology was to minimize, rather than maximize - and once I found a score of 0, that was a solution to the problem. This also meant adjusting to ga.setMinimax(Consts.minimaxType["minimize"]), rather than maximize.

One of the more superflous changes was to ga.setPopulationSize(sum(range(NQUEENS))). I felt that increasing or decreasing the size per generation based on the board size might help avoid stagnation of the gene pool - I never ran into this issue but reading various sources about genetic algorithms I became worried that it might become an issue for large values of BORDER. I have no evidence to back this thought, but it didn't appear to hurt anything, either.

I also increased ga.setGenerations(1E7), so that it would effectively run until a solution was found. As far as I know, there is no way to set a true "infinite" generations setting in PyEvolve, but this is pretty close.

FINAL NOTES
Feel free to use, abuse, modify, or otherwise put to use the code below. Please let me know if you know about the best way to set population size, or if any other problems are found.

CODE
#!/usr/bin/python
from pyevolve import *
from random import shuffle
import math

BORDER = 8
NQUEENS = BORDER
BOARD_SIZE = pow(BORDER, 2)

def queens_func(genome):
if sum(genome) < BORDER:
return BOARD_SIZE
array = map( lambda x: list(genome[x::BORDER]), range(BORDER))
cols = array
rows = map(list, zip(*array))
rows_ints = []

for row in rows:
rows_ints.append(int(''.join(map(str, row)), 2))
rows_ints = map( lambda x: int(math.log(x,2)) if x !=0 else x, rows_ints)
diags = filter( lambda x: x==1, map( lambda i,j: abs(i-j), rows_ints[1:], rows_ints[:-1]))
return max(map(sum, rows)) + max(map(sum, cols)) - 2 + sum(diags)

def queens_init(genome, **args):
genome.genomeList = [1]*NQUEENS + [0]*(BOARD_SIZE-NQUEENS)
shuffle(genome.genomeList)

def main_func():
genome = G1DList.G1DList(BOARD_SIZE)
genome.setParams(bestrawscore=0, rounddecimal=2)
genome.initializator.set(queens_init)
genome.mutator.set(Mutators.G1DListMutatorSwap)
genome.crossover.set(Crossovers.G1DListCrossoverCutCrossfill)
genome.evaluator.set(queens_func)
ga = GSimpleGA.GSimpleGA(genome)
ga.terminationCriteria.set(GSimpleGA.RawScoreCriteria)
ga.setMinimax(Consts.minimaxType["minimize"])
ga.setPopulationSize(100)
ga.setGenerations(1E7)
ga.setMutationRate(0.02)
ga.setCrossoverRate(1.0)
ga.evolve(freq_stats=50)
best = ga.bestIndividual()
print best
for i in map(lambda x: list(best[x::BORDER]), range(BORDER)):
print i

print "Best individual score : %.2f\n" %(best.getRawScore())

if __name__ == "__main__":
main_func()

## Friday, September 16, 2011

### Using FFTW plan_many_dft

CHASING FFTW3
The fftw3 library can be installed by doing sudo apt-get install libfftw3-3, and also
sudo ln -s /usr/lib/libfftw3.so.3.2.4 /usr/lib/libfftw3.so.

Apparently gcc can’t figure out that you are linking in a library unless it ends explicitly in .so, so by creating a symlink to the actual library a quick workaround is achieved.

The final gcc line to compile these programs is gcc <file> -o <objectname> -lfftw3 -lm

If libfftw3 gets broken during installation, bad things can happen. For me, this meant hours of trying to hack together a solution.

I finally ended up throwing in the towel, and doing sudo apt-get remove libfftw3-3, sudo apt-get purge libfftw3-3, sudo apt-get install libfftw3-3, sudo apt-get install libfftw3-3-dev, sudo apt-get install libfftw3-3-doc. Now everything works it great.

The single code is pretty simple, and is explained well by the FFTW documentation. READ IT! I mainly wrote the single code to have a reference template for doing other fftw stuff. The plan_many interface was confusing for me, so I wrote some code explaining it and how it is used.

A SIMPLE SCRIPT
Write a (python!... or something less awesome) script to print system time, call the executable, then print system time. EDIT: Or... just do time <your executable> in Linux.
In the dftmany C code -
#include <stdlib.h>
Change num_fft to a large number - say 400, 1000, etc.
Change fft_len to something more “usable”, say 8192
Fill the in array like this
in[j+i*fft_len][0] = (double)rand();
in[j+i*fft_len][1] = (double)rand();

Now laugh, maniacally.

Even with basic FFTW flags like FFTW_ESTIMATE, 400 8192 length FFTs can be  computed in .17s on a 3.2 GHz Core i7 - other processors may be slightly slower, but this application is single threaded so any computer with at least one idle 3+ GHz core could likely get similar results.

To go faster, try these gcc flags on for size - and remember “oh3” NOT “zero3”
CFLAGS = -O3 -ffast-math -mfpmath=sse -funroll-loops

And to go even EVEN faster, one could compile fftw from source using different build flags.
Have fun!

SIMPLE CODE
#include <stdio.h>
#include <fftw3.h>

int main (void) {
int N = 64;
double val = 1.;

fftw_complex* in;
fftw_complex* out;
fftw_plan plan;

in = (fftw_complex *) fftw_malloc(sizeof(fftw_complex) * N);
out = (fftw_complex *) fftw_malloc(sizeof(fftw_complex) * N);

plan = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);

for (int i = 0; i < N; i++) {
in[i][0] = val;
in[i][1] = 0.;
}

fftw_execute(plan);

#Results are non-normalized, need to scale by N. Should see DC=1, then all 0
for (int j = 0; j < N; j++) {
printf("%f+%fj", out[j][0]/N, out[j][1]/N);
}

fftw_destroy_plan(plan);

return 0;
}

PLAN_MANY CODE
#include <stdio.h>
#include <fftw3.h>

int main (void) {
int num_fft = 3;
int fft_len = 64;
double val = 0.;

fftw_complex* in;
fftw_complex* out;
fftw_plan plan;

//Alloc empty arrays the size of the data to be transformed
in = (fftw_complex *) fftw_malloc(sizeof(fftw_complex)*fft_len*num_fft);
out = (fftw_complex *) fftw_malloc(sizeof(fftw_complex)*fft_len*num_fft);

//Rank is 1, because the operation is for 1d arrays.
//
//Pointer to fft_len denotes that each fft is fft_len wide.
//
//Stride parameters are num_fft, because there are num_fft arrays contiguous in memory.
//
//Embed parameters are null because these arrays are not subsets of a larger n-dim
//array.
plan = fftw_plan_many_dft( 1,
&fft_len,
num_fft,
in,
NULL,
num_fft,
1,
out,
NULL,
num_fft,
1,
FFTW_FORWARD,
FFTW_ESTIMATE );

//Create input array of the form [1----][2----][3----][4-----]... etc. Use increasing val to show
//that the fft boundaries are correct. If two different values showed up in the same
//fft chunk, frequency content would be non-zero.
for (int i = 0; i < num_fft; i++){
val+=1;
for (int j = 0; j < fft_len; j++) {
in[j+i*fft_len][0] = val;
in[j+i*fft_len][1] = 0.;
}
}

fftw_execute(plan);

//Read back the output. Should be something like [1+0j 0+0j........... 1+0j 0+0j........] i.e.
//only DC terms show up in the spectrum output.
for (int q = 0; q < num_fft; q++) {
for (int r = 0; r < fft_len; r++) {
printf("%f+%fj\n", out[r+q*fft_len][0]/fft_len, out[r+q*fft_len][1]/fft_len);
}
}

fftw_destroy_plan(plan);

return 0;
}