starting point for a solution to this question: http://stats.stackexchange.com/questions/59674/parameters-estimation-of-ode-system

gistfile1.py

import numpy as np
from PyGMO import *
from PyGMO.problem import base
import scipy
import scipy.weave as weave
import random
import operator
import sys

c_function = r'''
int
func (double t, const double y[], double f[],
      void *params)
{

  double a1 = ((double*)params)[0];
  double a2 = ((double*)params)[1];
  double a3 = ((double*)params)[2];
  double a4 = ((double*)params)[3];
  double a5 = ((double*)params)[4];
  double a6 = ((double*)params)[5];
  double a7 = ((double*)params)[6];
  double a8 = ((double*)params)[7];
  double a9 = ((double*)params)[8];

  // S = y[0], D = y[1], F = y[2]
  f[0] = a1*y[0]+a2*y[1]+a3*y[2];
  f[1] = a4*y[0]+a5*y[1]+a6*y[2];
  f[2] = a7*y[0]+a8*y[1]+a9*y[2];

  return GSL_SUCCESS;
}

'''

c_main = r'''
 const gsl_odeiv_step_type * T = gsl_odeiv_step_rkf45;

  gsl_odeiv_step * s 
    = gsl_odeiv_step_alloc (T, 3);
  gsl_odeiv_control * c 
    = gsl_odeiv_control_y_new (1e-6, 0.0);
  gsl_odeiv_evolve * e 
    = gsl_odeiv_evolve_alloc (3);

  gsl_odeiv_system sys = {func, NULL, 3, params};

  double t = 0.0, t1 = 11.;
  double h = 1e-6;
  double y[3] = { 17110., 1157.33, 104181. };


  while (t < t1)
    {
      int status = gsl_odeiv_evolve_apply (e, c, s,
                                           &sys, 
                                           &t, t1,
                                           &h, y);

      if (status != GSL_SUCCESS)
          break;

    }

  results[0] = y[0];
  results[1] = y[1];
  results[2] = y[2];

  gsl_odeiv_evolve_free(e);
  gsl_odeiv_control_free(c);
  gsl_odeiv_step_free(s);

'''


class lpa_problem(base):
    def __init__(self):
        super(lpa_problem, self).__init__(9, 0, 1)

        self.set_bounds(0.,10.)  # set bounds of parameters here
        
        self.c_function = c_function
        self.c_main = c_main

    def _objfun_impl(self, variable_parameters):
        results = scipy.zeros(3)
        params = scipy.array(variable_parameters)

        compiler = 'gcc'
        variables = ['results', 'params']
        libs = ['gsl','gslcblas','m']
        headers = ['<stdio.h>','<gsl/gsl_errno.h>', '<gsl/gsl_matrix.h>', '<gsl/gsl_odeiv.h>']
        libpath = ['/usr/lib64/']
        headerpath = ['PATH/TO/PYTHON.H']

        res = weave.inline(c_main, variables, compiler=compiler,
                           libraries=libs, headers=headers,
                           include_dirs = headerpath, library_dirs = libpath,
                           support_code = c_function)

        S_final = results[0]
        D_final = results[1]
        F_final = results[2]

        return ( (21162. - S_final)**2 + (1541.52-D_final)**2 + (242864.-F_final)**2, )


def main():
    algo = algorithm.de(gen=10)
    prob = lpa_problem()
    print 'fetched problem'

    
    archi = archipelago(algo,prob, 100, 20, topology = topology.ageing_clustered_ba(a=100))

    archi.evolve(1)
    archi.join()
    
    for i in range(2000):
        archi.evolve(1); archi.join()
        if i % 10 == 0:
            print '--------------'
            print 'round',i
            for island in random.sample(archi, 10):
                print island.population.champion.x, island.population.champion.f, island.population.champion.c
            print '--------------'
            print min([island.population.champion.f[0] for island in archi])
            print max([island.population.champion.f[0] for island in archi])
            print '--------------'

    for island_i, island in enumerate(archi):
        print 'island',island_i
        print island.population.champion.x, island.population.champion.f, island.population.champion.c


if __name__ == '__main__':
    main()