How can I use Scipy to fit data generated from a C++ model?

Question

I currently have a functioning and blazing fast model written in C++ and CUDA. However, I'd like to use Scipy.minimize to fit the model to some experimental data. I was hoping it would be easy, but when I try to describe my problem to Google, it keeps telling me that I want to extend Python with C++, the documentation of which makes me want to cry.

I used to use PyCUDA back in the day, but I don't think I am able to wrap the entire model, so I need a little nudge in the right direction.

What I'm trying to accomplish.

from scipy.optimize import minimize
import os

#Initial starting parameters
ics = [ 1, 0, 0 ... 0]

#The objective function is called from my model, which is a built .exe that will accept parameters as maybe a string and somehow output something
model = os.popen("directory/mymodel.exe" + params)

difference = minimize(model, ics, method='Nelder-Mead')

Is something like this possible? I'm having a hard time getting around the interface between C++ and Python and how I can make my .exe generate data interpretable by python so that Scipy can adjust the parameters and iterate.

Answer 1

Other comments have suggested a file-based interface, using an actual C/C++ optimization library, or extending Python with C++. Those are probably better ways to solve your problem, but here's a more narrow minded answer.

To begin, here's a C program called mymodel.c which implements the 2d Rosenbrock function.

#include "stdio.h"
#include "assert.h"
#include "stdlib.h"

/* minimum is at (a, a*a) = (0.123, 0.015129) */
double rosenbrock_2d(double x, double y) {
  const double a = 0.123;
  const double b = 100.0;
  double u, v;
  u = a - x;
  v = y - x*x;
  return u*u + b*v*v;
}

int main(int argc, char **argv)
{
  double x, y;
  assert(argc == 3);
  x = atof(argv[1]);
  y = atof(argv[2]);
  printf("%.16g\n", rosenbrock_2d(x, y));
  return 0;
}

You can build it using something like the following command.

$ gcc -o mymodel.exe mymodel.c

Here's a python script named opt.py which repeatedly calls the executable implementing your model, while attempting to find the minimum of the function.

from scipy.optimize import minimize
import subprocess

#Initial starting parameters
ics = [ 1, 0]

#The objective function is called from my model, which is a built .exe that will accept parameters as maybe a string and somehow output something
def model(x):
    params = ['%.16g' % a for a in x]
    s = subprocess.check_output(['./mymodel.exe'] + params)
    return float(s)

print minimize(model, ics, method='Nelder-Mead')

Then run it:

$ python opt.py                            
  status: 0
    nfev: 86
 success: True
     fun: 2.2359521022451379e-09
       x: array([ 0.12296347,  0.01512302])
 message: 'Optimization terminated successfully.'
     nit: 46

This gives something close to the answer. If you want a closer answer you can use different tolerances or a better method like 'L-BFGS-B'.