scipy.optimize.least_squares not giving correct result when using C++ program from Python

I am trying to use the scipy.optimize.least_squares function in order to minimize a model generated from a C++ code.

Here is a minimal example to demonstrate the problem (taken from the Scipy docs).

I use two functions, one that is pure Python and another one that calls a C++ code to generate the model. Both functions return the same results when tested but not when used in scipy.optimize.least_squares.

import numpy as np
from scipy.optimize import least_squares
import csv
import os
import time

def gen_data(t, a, b, c, noise=0, n_outliers=0, random_state=0):
    y = a + b * np.exp(t * c)
    rnd = np.random.RandomState(random_state)
    error = noise * rnd.randn(t.size)
    outliers = rnd.randint(0, t.size, n_outliers)
    error[outliers] *= 10
    return y + error

a = 0.5
b = 2.0
c = -1
t_min = 0
t_max = 10
n_points = 15
t_train = np.linspace(t_min, t_max, n_points)
y_train = gen_data(t_train, a, b, c, noise=0.1, n_outliers=3)

def fun_cpp(x, t, y):
    a = os.popen("./a.out "+str(x[0])+" "+str(x[1])+" "+str(x[2])).readlines()
    model = [float(i) for i in a[0].split(' ')[:-1]]

    return model - y

def fun_python(x, t, y):
    return (x[0] + x[1] * np.exp(x[2] * t)) - y

x0 = np.array([1.0, 1.0, 0.0])

res_lsq = least_squares(fun_python, x0, args=(t_train, y_train), verbose = 2)
res_lsq2 = least_squares(fun_cpp, x0, args=(t_train, y_train), verbose = 2)

t_test = np.linspace(t_min, t_max, n_points * 10)
y_true = gen_data(t_test, a, b, c)
y_lsq = gen_data(t_test, *res_lsq.x)
y_lsq2 = gen_data(t_test, *res_lsq2.x)

import matplotlib.pyplot as plt
plt.ion() 
plt.plot(t_train, y_train, 'o')
plt.plot(t_test, y_true, 'k', linewidth=2, label='true')
plt.plot(t_test, y_lsq, label='pure Python')
plt.plot(t_test, y_lsq2, label='C++/Python')

plt.legend()
plt.show()

And here is the C++ code (to run :./a.out 0.5 2 -1):

#include<iostream> 
#include <stdlib.h>
#include <math.h>
#include <vector>
using std::vector;

vector<double> linspace(double a, double b, int n) {
    vector<double> array;
    double step = (b-a) / (n-1);

    while(a <= b) {
        array.push_back(a);
        a += step;  
    }
    return array;
}

int main(int argc, char** argv) {  // argv values must be ints or floats. For example 0.5, 2.0, -1
  using namespace std;

  vector<double> x = linspace(0, 10, 15);
  vector<double> model;
  vector<double> parameters (3);

  if(argc<4) {
    parameters[0] = 0.5;
    parameters[1] = 2.0;
    parameters[2] = -1;
  }
  else{
    parameters[0] = atof(argv[1]);
    parameters[1] = atof(argv[2]);
    parameters[2] = atof(argv[3]);
  }

  for (int i=0; i<x.size(); i++){
    model.push_back(parameters[0]+ parameters[1] * exp(parameters[2] * x[i]));
    cout << model[i]<< " ";
  }
}

This code results in the following figure:

Any advice on how to proceed ?

Thanks

Solution

Ok so I think I've identified my problem.

The input function need to be derivable which fun_python is whereas fun_cpp is not as it call upon an external program.

To solve my problem a black box optimization algorithm can be used, for example: https://github.com/paulknysh/blackbox

The example code above can be modified to:

import numpy as np
from scipy.optimize import least_squares
import csv
import os
import time

def gen_data(t, a, b, c, noise=0, n_outliers=0, random_state=0):
    y = a + b * np.exp(t * c)
    rnd = np.random.RandomState(random_state)
    error = noise * rnd.randn(t.size)
    outliers = rnd.randint(0, t.size, n_outliers)
    error[outliers] *= 10
    return y + error

a = 0.5
b = 2.0
c = -1
t_min = 0
t_max = 10
n_points = 15
t_train = np.linspace(t_min, t_max, n_points)
y_train = gen_data(t_train, a, b, c, noise=0.1, n_outliers=3)

def fun_cpp(x):
    global y_train
    a = os.popen("./a.out "+str(x[0])+" "+str(x[1])+" "+str(x[2])).readlines()
    model = [float(i) for i in a[0].split(' ')[:-1]]
    chisquared = np.sum([(y1-y2)**2. for y1,y2 in zip(y_train, model)])
    return chisquared

def fun_python(x, t, y):
    return (x[0] + x[1] * np.exp(x[2] * t)) - y

x0 = np.array([1.0, 1.0, 0.0])

res_lsq = least_squares(fun_python, x0, args=(t_train, y_train))

import blackbox as bb

bb.search(f=fun_cpp,  # given function
              box=[[0., 1.], [0., 3], [-2., 1.]],  # range of values for each parameter
              n=40,  # number of function calls on initial stage (global search)
              m=40,  # number of function calls on subsequent stage (local search)
              batch=4,  # number of calls that will be evaluated in parallel
              resfile='output.csv')  # text file where results will be saved

rowSelect = []
with open('output.csv') as csvfile:
    readCSV = csv.reader(csvfile, delimiter=',')
    for row in readCSV:
        rowSelect.append(row)

bestFitParam = [float(x) for x in rowSelect[1][0:-1]]

t_test = np.linspace(t_min, t_max, n_points * 10)
y_true = gen_data(t_test, a, b, c)
y_lsq = gen_data(t_test, *res_lsq.x)
y_lsq2 = gen_data(t_test, *bestFitParam)

import matplotlib.pyplot as plt
plt.ion() 
plt.plot(t_train, y_train, 'o')
plt.plot(t_test, y_true, 'k', linewidth=2, label='true')
plt.plot(t_test, y_lsq, label='pure Python')
plt.plot(t_test, y_lsq2, label='Black-box C++/Python')

plt.legend()
plt.show()

Which results in the following figure:

The Black box optimization results in better results but takes longer to compute.