Python code to run a C program that calculates whether points are in the mandelbrot set

Question

This is part of a project to explore the Mandelbrot set. The code to graphically navigate the complex plane and allow the user to adjust the colormap and number of iterations is all written in Python, because that sort of user interface work is easiest in a scripting language. However, Python is not very good at inner loops that run hundreds of times for each pixel, so I wrote a C program that does the actual computation. The python function runs the C executable in a separate process that passes the result back from C to python as a human-readable string. Parsing the string back to integers is slow. What is a more pythonic way of creating the interface between C and python? Also, any general tips for writing faster, better C are appreciated, I'm new to the language.

Python code:

import numpy as np
from subprocess import Popen
import subprocess

def getcounts(xmin, xmax, ymin, ymax, dim, depth):
    """render a view of the mandelbrot set on the given rectangle of the
    complex plane, by calling a C script and interpreting stdout as a
    series of integers
    output is a numpy array reshaped into a dim*dim square"""

    p = Popen(["./a.out", repr(xmin), repr(xmax), repr(ymin),
    repr(ymax), str(dim), str(depth)], stdout=subprocess.PIPE)

    counts = np.fromstring(p.stdout.readline()[:-1], dtype = np.int, sep = ' ')

    return counts.reshape([dim, dim]).transpose()

C code:

#include <stdio.h>
#include <stdlib.h> 
int main(int argc, char *argv[]){
  int dim = 100;        /*defaults produce a nice overview of the mandelbrot set*/
  int depth = 300;
  double minx = -2;
  double maxx = 2;
  double miny = -1.5;
  double maxy = 1.5;     

  if(argc > 6){
      minx = atof(argv[1]);
      maxx = atof(argv[2]);
      miny = atof(argv[3]);
      maxy = atof(argv[4]);
      dim = atoi(argv[5]);
      depth = atoi(argv[6]);
  }

  double xstep = (maxx - minx)/dim;
  double ystep = (maxy - miny)/dim;

  int i, j;
  for(i = 0; i < dim; i++){
    double c_real = minx + i * xstep;
    for(j = 0; j < dim; j++){
      double c_imag = miny + j * ystep;
      int count = 0;
      double z_real = 0;
      double z_imag = 0;
      double zreal_temp;
      while((z_real * z_real + z_imag * z_imag < 4) && count < depth){
          count ++;
          zreal_temp = z_real * z_real - z_imag * z_imag + c_real;
          z_imag = z_real * z_imag * 2 + c_imag;
          z_real = zreal_temp;
      }
      printf("%d ", count);
    }   
  }
}

Finally, here's a python script to test the above function: I've cut out all the bits for moving around and adjusting settings since this question is about the C - Python interface

import numpy as np
import matplotlib.pyplot as plt
xmin = -2
xmax = 2
ymin = -2
ymax = 2
depth = 60   
dim = 600
counts = getcounts(xmin, xmax, ymin, ymax, dim, depth)
plt.imshow(counts, extent = [xmin, xmax, ymax, ymin], interpolation = 'nearest')
plt.show()

Answer 1

In my opinion, Popen is simply the wrong approach to this problem. Everything else follows from that. Python provides an API to directly expose C functions to Python. Boost.Python is a fantastic library to make it easy to expose C++ classes without having to deal with the messy direct Python API. Boost.Numpy makes it easy to return numpy arrays directly.

Regardless of which one you choose, you should choose one and use it to create something to the effect of mandelbrot.so which will have one function that returns a list of lists. You can then just use that function directly. Rather than write getcounts() to do some work doing string processing, just call mandelbrot.getcounts() with all the same args. Something like:

import numpy as np
import matplotlib.pyplot as plt
import mandelbrot

xmin = -2
xmax = 2
ymin = -2
ymax = 2
depth = 60   
dim = 600
counts = mandelbrot.getcounts(xmin, xmax, ymin, ymax, dim, depth)
plt.imshow(counts, extent = [xmin, xmax, ymax, ymin], interpolation = 'nearest')
plt.show()

There are several huge advantages to doing it this way:

• The ability to actually do any sort of error checking.
• Not having to do any output processing converting between one program's output and another program's desired input.
• Not having to start a new process and do input handling on that one will make the code perform dramatically better.
• It'll read more Pythonic on one end (just call a function that returns a value) and more C/C++-ic on the other end (have a function that returns a value there too, instead of printing output). This has the added benefit of being able to actually write tests for your C/C++ code.

Answer 2

This is nitpicky, but your Popen constructor is odd. Why are you calling a mix of repr and str functions for your number arguments? For integers there's not much of a difference except if you get large, where repr is less desireable:

>>> repr(10000000000000000000000000000)
'10000000000000000000000000000L'
>>> str(10000000000000000000000000000)
'10000000000000000000000000000'

I don't think you want that L. For floats, repr offers more precision but you don't use any floats in your sample. If it's an intentional consideration, note it in a comment so that it's clear for another user.

Even aside from that, you have a long Popen constructor. You could shorten it if you only needed one function, because you could call map.

["./a.out", str(xmin), str(xmax), str(ymin), str(ymax), str(dim), str(depth)]

map(str, ["./a.out", xmin, xmax, ymin, ymax, dim, depth])

Note that there's no adverse result of calling str("./a.out") so this will just give you a list of strings the way you want them.