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_fft, len_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");',

          headers=['<stdio.h>'],

          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'],

headers=['<fftw3.h>'],

libraries=['fftw3'],

        include_dirs=fftw_info['include_dirs'],

library_dirs=fftw_info['library_dirs'],

compiler='gcc' )

print comp

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)
    #vpython_adapter = DBAdapters.DBVPythonGraph(identify="queens", frequency=1)
    #ga.setDBAdapter(vpython_adapter)
    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()