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

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)



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.

• Why not just write the data to a file and read it with open()? – Bill Barth Sep 3 '14 at 20:46
• I can do that, but won't that introduce alot of overhead into it? I'm not sure how many times scipy will have to iterate my model to minimize the objective function, so I wanted to achieve the data transfer the fastest possible way on the first try instead of spending a lot of time messing with it. I already have the .exe printing out large matrices for plotting, so it wouldn't be too difficult. – Hair of Slytherin Sep 3 '14 at 20:49
• Might be better to get a C/C++ callable optimization library like TAO. – Bill Barth Sep 3 '14 at 20:59
• You want to extend Python with C++... Wrapping a single C++ function to be callable from Python is not all that complicated, but I agree that it looks scary at first. Thing is, it is well invested time: one you learn how to do it, you will find yourself doing it over and over again... – Jaime Sep 3 '14 at 21:18

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'.

• Alright! That looks like precisely what I was trying to accomplish. If I can bend this to my will, then this will keep my adviser content while I learn how to extend python! I'll give it a shot and accept if it works out. – Hair of Slytherin Sep 4 '14 at 2:06
• One question : could you explain what "'%.16g' % a for a in x" does? To me it looks like you're formatting each element in the list x with maybe 16 digits after the decimal, but I'm not sure how to describe the question to google to get the answer. – Hair of Slytherin Sep 4 '14 at 6:17
• You are correct. See this reference. g is the most compact choice of either floating point or exponential notation with precision of 16 digits. – moyner Sep 4 '14 at 11:30
• @KarstenChu ['%.16g' % a for a in x] is a Python list comprehension that creates a list of strings given a sequence x of floating point numbers. The float->string conversion ('%.16g' % a) converts the floating point number to a string using a standard format specifier from the format specification mini-language. I hope this arms you with some keywords to help with your googling! – k20 Sep 4 '14 at 15:22
• @k20 Yes, it does! Thank you. I've been out of touch with Python for a while so I'd totally forgotten about list comprehension so I was frustrated to find that I couldn't even formulate my question to be answered. – Hair of Slytherin Sep 4 '14 at 22:13