How to create a gsl_odeiv_system structure?

[hakonhagland]

I am following this example: https://www.math.utah.edu/software/gsl/gsl-ref_367.html :

#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>

int
func (double t, const double y[], double f[],
      void *params)
{
  double mu = *(double *)params;
  f[0] = y[1];
  f[1] = -y[0] - mu*y[1]*(y[0]*y[0] - 1);
  return GSL_SUCCESS;
}

int
jac (double t, const double y[], double *dfdy, 
     double dfdt[], void *params)
{
  double mu = *(double *)params;
  gsl_matrix_view dfdy_mat 
    = gsl_matrix_view_array (dfdy, 2, 2);
  gsl_matrix * m = &dfdy_mat.matrix; 
  gsl_matrix_set (m, 0, 0, 0.0);
  gsl_matrix_set (m, 0, 1, 1.0);
  gsl_matrix_set (m, 1, 0, -2.0*mu*y[0]*y[1] - 1.0);
  gsl_matrix_set (m, 1, 1, -mu*(y[0]*y[0] - 1.0));
  dfdt[0] = 0.0;
  dfdt[1] = 0.0;
  return GSL_SUCCESS;
}

int
main (void)
{
  const gsl_odeiv_step_type * T 
    = gsl_odeiv_step_rk8pd;

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

  double mu = 10;
  gsl_odeiv_system sys = {func, jac, 2, &mu};

  double t = 0.0, t1 = 100.0;
  double h = 1e-6;
  double y[2] = { 1.0, 0.0 };

  gsl_ieee_env_setup();

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

      if (status != GSL_SUCCESS)
          break;

      printf("%.5e %.5e %.5e\n", t, y[0], y[1]);
    }

  gsl_odeiv_evolve_free(e);
  gsl_odeiv_control_free(c);
  gsl_odeiv_step_free(s);
  return 0;
}

I would like to write this in Perl using Math::GSL. I tried:

use strict;
use warnings;

# Trying to implement a Perl Math::GSL::ODEIV version of the
# Gnu GSL C code given here:
#
#   https://www.math.utah.edu/software/gsl/gsl-ref_367.html
#
use Math::GSL::Errno qw($GSL_SUCCESS);
use Math::GSL::ODEIV qw/ :all /;
use Math::GSL::Matrix qw/:all/;
use Math::GSL::IEEEUtils qw/ :all /;

sub func {
    my ($t, $y, $f, $params) = @_;
    my $mu = $params->{mu};

    $f->[0] = $y->[1];
    $f->[1] = -$y->[0] - $mu*$y->[1]*(($y->[0])**2 - 1);
    return $GSL_SUCCESS;
}

sub jac {
    my ($t, $y, $dfdy, $dfdt, $params) = @_;

    my $mu = $params->{mu};
    my $m = gsl_matrix_view_array($dfdy, 2, 2);
    gsl_matrix_set( $m, 0, 0, 0.0 );
    gsl_matrix_set( $m, 0, 1, 1.0 );
    gsl_matrix_set( $m, 1, 0, -2.0 * $mu * $y->[0] * $y->[1] - 1.0 );
    gsl_matrix_set( $m, 1, 1, -$mu * ($y->[0]**2 - 1.0) );
    $dfdt->[0] = 0.0;
    $dfdt->[1] = 0.0;
    return $GSL_SUCCESS;
}

my $T = $gsl_odeiv_step_rk8pd;
my $s = gsl_odeiv_step_alloc($T, 2);
my $c = gsl_odeiv_control_y_new(1e-6, 0.0);
my $e = gsl_odeiv_evolve_alloc(2);
my $mu = 10;
# Now I would like to do this, but it does not work:
#
#    my $sys = Math::GSL::ODEIV::gsl_odeiv_system->new( \&func, \&jac, 2, \$mu );
#
#
# but this however works:
my $sys = Math::GSL::ODEIV::gsl_odeiv_system->new( );
# but then how can I set the function, jacobian, dimension and parameters of
# the $sys object (struct)?
# 

[leto]

@hakonhagland Thanks for this detailed bug report. Can you show exactly what error you get?

Also, which version of GSL and Math::GSL ?

[hakonhagland]

This was run on my Ubuntu laptop:

$ uname -a
Linux hakon-Vostro-5568 4.15.0-36-generic #39-Ubuntu SMP Mon Sep 24 16:19:09 UTC 2018 x86_64 x86_64 x86_64 GNU/Linux
$ lsb_release -a
No LSB modules are available.
Distributor ID:	Ubuntu
Description:	Ubuntu 18.04.1 LTS
Release:	18.04
Codename:	bionic
$ perl --version | head -2
This is perl 5, version 26, subversion 1 (v5.26.1) built for x86_64-linux
$ swig -version
SWIG Version 3.0.12
Compiled with g++ [x86_64-pc-linux-gnu]
Configured options: +pcre
$ gsl-config --version
2.5
$  perl -MMath::GSL -E 'say $Math::GSL::VERSION'
0.39

The error I get when including the commented out line:

my $sys = Math::GSL::ODEIV::gsl_odeiv_system->new( \&func, \&jac, 2, \$mu );

is

RuntimeError Usage: new_gsl_odeiv_system(); at /home/hakon/perlbrew/perls/perl-5.26.1/lib/site_perl/5.26.1/x86_64-linux/Math/GSL/ODEIV.pm line 96.

[leto]

@hakonhagland have you tried creating the object, and then trying to set the jacobian? There should be some way to set it after creation time

[leto]

@hakonhagland

This function can set the jacobian after creation time:

Math-GSL-0.39/lib/Math/GSL/ODEIV.pm:89:*swig_jacobian_set = *Math::GSL::ODEIVc::gsl_odeiv_system_jacobian_set;

[hakonhagland]

@leto I tried both:

Math::GSL::ODEIV::gsl_odeiv_system::swig_jacobian_set( $sys, \&jac );

and

$sys->swig_jacobian_set( \&jac );

both of these give error:

TypeError in method 'gsl_odeiv_system_jacobian_set', argument 2 of type 'int (*)(double,double const [],double *,double [],void *)'

[leto]

We may need more swig/xs glue, similar to the fdf solver stuff that was just merged

[morgantomj]

leto, hakonhagland,

My thanks to both of you for looking into gsl_odeiv_system. It may be asking a bit much, but would it be possible to add the 'driver' GSL functions to any new distro of Math::GSL::ODEIV? Note that they also require a useful gsl_odeiv_system method.

One other question: I have seen reference to 'odeiv2' in the GSL lib and in the Math::GSL::ODEIV package. How does it relate to the functions in the perl libs?

Thanx,
Tom Morgan

[leto]

@morgantomj it looks like odeiv2 was added in 1.15 and Math::GSL has never had any glue for it. Just created an issue for that: #157

I haven't used these 'driver' functions before, can you give an example of what you want to do with them?

[morgantomj]

Duke, I'm interested in using the Math::GSL::ODEIV module functions to replace some code that now uses Numerical Recipes C-code (which is copyrighted). There is an example of using the gsl driver functions in the gsl example for the ODEIV C-code. The example shows that the driver functions simplify solving the ODE system. It still needs the gsl_odeiv_system data type, tho. I don't quite follow how SWIG is used to glue the gsl C-code together with perl, but I would guess that the SWIG code for the driver funstions is very similar to the code for the control or evolution SWIG code. So my naive view is that extending the code to the driver functions should be fairly easy and quick...Famous Last Words and Easy for me to say. Thanx again for you work on this. Tom Morgan

…

On Sat, Oct 20, 2018 at 11:31 AM Duke Leto ***@***.***> wrote: @morgantomj <https://github.com/morgantomj> it looks like odeiv2 was added in 1.15 and Math::GSL has never had any glue for it. Just created an issue for that: #157 <#157> I haven't used these 'driver' functions before, can you give an example of what you want to do with them? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#156 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AqQEXDlOBJISNCJkbAfC5BThOUv9sUKUks5um0HCgaJpZM4XZ_ja> .

[morgantomj]

Duke,
I believe that the problem with defining a system structure in the Math::GSL::ODEIV is that there is no C function that returns the 'gsl_odeiv2_system sys' structure - in the GSL C-code example you find this C syntax to define the sys:

gsl_odeiv2_system sys = {func, jac, 2, &mu};

If there were a C-function like ( but properly written):
gsl_odeiv2_system *sys = define_gsl_system (
int func (double t, const double y[], double f[],
void *params),
int jac (double t, const double y[], double *dfdy,
double dfdt[], void *params) func (double t, const double y[], double f[],
void *params),
const gsl_odeiv2_step_type * T,
void *params};

My C-code is poor but he idea is to NOT use the construct 'sys = {...}' but to have a function that returns *sys, given the func, jac, etc. as arguments, then there could be a corresponding swig function to implement it in perl Math::GSL::ODEIV.

Ideally this new function would be defined in the gsl C source itself. I think then it would just fall out of the code that builds Math::GSL::ODEIV. I don't know if it's possible to add such a function to the Math::GSL side, so that swig can add it to the pile, assuming the gsl library has already been compiled.

Thanks for your work on this ...

Tom Morgan

[leto]

@morgantomj you might want to look at and test out the latest code from @hakonhagland . it's likely to be merged

[morgantomj]

Thu Jan 30 15:55:39 EST 2020

I have been thinking about how to get around the problem of using the Math::GSL::ODEIV module -
that is, to use the ODEIV (and ODEIV2) Math::GSL module. It seems to me it is necessary to pass
a C-pointer to a C-function from the perl script to the Math::GSL::ODEIV module interface methods.
But the user provided function
(that represents the Ordinary Differential Equation (ODE) system) has not yet been
compiled. So the problem is how to pass a C pointer to a C-function the perl module method(s).

In the C-code for the gsl library for the ODEIV module, there is a C-code example:

/////////////// beginning of C code ///////////////
/* Example for using gsl ODEIV library
https://www.gnu.org/software/gsl/doc/html/ode-initval.html
*/

#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>

int
func (double t, const double y[], double f[],
void params)
{
(void)(t); / avoid unused parameter warning */
double mu = (double )params;
f[0] = y[1];
f[1] = -y[0] - muy[1](y[0]*y[0] - 1);
return GSL_SUCCESS;
}

int
jac (double t, const double y[], double *dfdy,
double dfdt[], void params)
{
(void)(t); / avoid unused parameter warning */
double mu = (double )params;
gsl_matrix_view dfdy_mat
= gsl_matrix_view_array (dfdy, 2, 2);
gsl_matrix * m = &dfdy_mat.matrix;
gsl_matrix_set (m, 0, 0, 0.0);
gsl_matrix_set (m, 0, 1, 1.0);
gsl_matrix_set (m, 1, 0, -2.0muy[0]y[1] - 1.0);
gsl_matrix_set (m, 1, 1, -mu(y[0]*y[0] - 1.0));
dfdt[0] = 0.0;
dfdt[1] = 0.0;
return GSL_SUCCESS;
}

int
main (void)
{
double mu = 10;
gsl_odeiv2_system sys = {func, jac, 2, &mu};

gsl_odeiv2_driver * d =
gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rk8pd,
1e-6, 1e-6, 0.0);
int i;
double t = 0.0, t1 = 100.0;
double y[2] = { 1.0, 0.0 };

for (i = 1; i <= 100; i++)
{
double ti = i * t1 / 100.0;
int status = gsl_odeiv2_driver_apply (d, &t, ti, y);

  if (status != GSL_SUCCESS)
    {
      printf ("error, return value=%d\n", status);
      break;
    }

  printf ("%.5e %.5e %.5e\n", t, y[0], y[1]);
}

gsl_odeiv2_driver_free (d);
return 0;
}

/////////////// End of C-code //////////////

A few things to notice:

1. There are 2 functions defined here 'func' and 'jac' (the Jacobian)

2. The stucture sys is declared as: gsl_odeiv2_system. This is a structure
   which has 4 fields - two of which are C pointers to user provided functions.

3. The code uses 'sys = {func, jac, 2, &mu}' to assign to the 4 fields.

4. 'func' and 'jac' are how C function pointers are defined. Here is the declaration of the
   function 'func':
   int func (double t, const double y[], double f[], void *params);
   There is a similar declaration for 'jac'.

5. the first field MUST have the signature:
   int (*) (double t, const double y[], double f[], void *params);

6. The structure 'sys' is passed to 'driver' function.

7. The gsl library is already compiled, but the function 'func' is user provided,
   and needs to be compiled, and the gsl library linked in to make an executable.
   Also for 'jac'.

The above all works to solve the user provided ODE system difined by 'func' and 'jac'.

The problem, as I see it, occurs when trying to use the perl module Math::GSL::ODEIV
to solve the same problem.

Somehow the 4 'sys' structure fields must be correctly assigned to. In particular,
the 2 C-function pointersare pointers to C functions. As far as I can tell the user
MUST create a pointer to a C-function from within the perl script that is using the
module Math::GSL::ODEIV.

So it appears that in order to use the Math::GSL::ODEIV module, the user needs to write
the C-code for the function 'func', for example, and SOMEHOW compile it and link that
code to the shared object ODEIV.so that the SWIG process has produced and which is used
in Math::GSL::ODEIV.

By the way, there are other parts of the gls library that require the user to define a C function,
and pass a pointer to that function to some other library function. For example, finding roots,
finding numerical derivatives, numerical integration. I suspect all of these areas will have
the same type of problem as I have described above for Math::GSL::ODEIV

I am at the limit of my knowledge of both C, perl, and Dynaloader. Is there some way to use
swig to compile the user provided functions (probably with some wrapper functions, too)
and then somehow link in the ODEIV.so library so that in a perl script, the user could
assign to the 'sys' structure, and then pass that to the driver function, and thereby solve the
system of ODE's.

If not SWIG, would it be feasible to use something like Inlne::C to do the same thing.

By the way I wrote a C-function which can take the place of the 'sys = {func, jac, 2, &mu}'.
I vainly hoped that it would allow the assignment of 'func' to the first field, etc.
That was the point when I realized that the real problem is how to assign a C pointer to a C function
to a field in the 'sys' structure from within a perl script.

///////////// Beginning of C code system.c //////////////
/* ode-initval2/system.c
*
*

• This program is free software; you can redistribute it and/or modify
• it under the terms of the GNU General Public License as published by
• the Free Software Foundation; either version 3 of the License, or (at
• your option) any later version.
• This program is distributed in the hope that it will be useful, but
• WITHOUT ANY WARRANTY; without even the implied warranty of
• MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
• General Public License for more details.
• You should have received a copy of the GNU General Public License
• along with this program; if not, write to the Free Software
• Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  */

/* System routine for odeiv2.
TJM 2019-12-6: I want to use SWIG to allow access to perl scripts of ODEIV2
In the example C-code examples, this construct is used to define/declare the
system:
gsl_odeiv2_system sys = {func, jac, 2, &mu};
This construct directly sets (initializes) the 4 fields of the gsl_odeiv2_system struct.
I do not know how to perform the equivalent in the perl/SWIG equivalent code.
So I am defining a function to perform the assignment, which SWIG should be able to
incorporate into the perl, allowing the system to be defined via the new function
'gsl_odeiv2_system_define'.

*/

#include <config.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv2.h>
#include <gsl/gsl_machine.h>

static gsl_odeiv2_system *
gsl_odeiv2_system_alloc_new ( int (*function) (double t, const double y[], double dydt[], void *params),
int (*jacobian) (double t, const double y[], double *dfdy, double dfdt[],
void *params),
size_t dimension,
void *params
)
{
gsl_odeiv2_system *sys;

sys = (gsl_odeiv2_system *) calloc (1, sizeof (gsl_odeiv2_system)); // this should allocate the space for struct 'sys'

if (sys == NULL)
{
GSL_ERROR_NULL ("failed to allocate space for 'system' struct",
GSL_ENOMEM);
}

sys->function = function;
sys->jacobian = jacobian;
sys->dimension = dimension;
sys->params = params;
return sys;
}

void
gsl_odeiv2_system_free (gsl_odeiv2_system * sys)
{
free (sys);
}

/////////////// End of system.c //////////

So I'm stumped as to how to get around this problem.

Any thoughts?

Tom Morgan [email protected]

[hakonhagland]

Hi Tom,

nice to see that you are still working on this issue! Have you tried the code in PR #172 ? I think it solves the problem you are referring to. Can you try the following:

git clone https://github.com/leto/math--gsl.git
cd math--gsl
git fetch origin pull/172/head
git checkout -b odeiv FETCH_HEAD
perl Build.PL
./Build installdeps --cpan_client 'cpanm'
./Build install  # NOTE: I am using perlbrew, so I do not need to run this as sudo

Then create a this Perl test script:

use strict;
use warnings;
use Math::GSL::Errno qw($GSL_SUCCESS);
use Math::GSL::ODEIV qw/ :all /;
use Math::GSL::Matrix qw/:all/;
use Math::GSL::IEEEUtils qw/ :all /;

sub func {
    my ($t, $y, $dydt, $params) = @_;
    my $mu = $params->{mu};
    $dydt->[0] = $y->[1];
    $dydt->[1] = -$y->[0] - $mu*$y->[1]*(($y->[0])**2 - 1);
    return $GSL_SUCCESS;
}

sub jac {
    my ($t, $y, $dfdy, $dfdt, $params) = @_;

    my $mu = $params->{mu};
    my $m = gsl_matrix_view_array($dfdy, 2, 2);
    gsl_matrix_set( $m, 0, 0, 0.0 );
    gsl_matrix_set( $m, 0, 1, 1.0 );
    gsl_matrix_set( $m, 1, 0, (-2.0 * $mu * $y->[0] * $y->[1]) - 1.0 );
    gsl_matrix_set( $m, 1, 1, -$mu * (($y->[0])**2 - 1.0) );
    $dfdt->[0] = 0.0;
    $dfdt->[1] = 0.0;
    return $GSL_SUCCESS;
}

my $T = $gsl_odeiv_step_rk8pd;
my $s = gsl_odeiv_step_alloc($T, 2);
my $c = gsl_odeiv_control_y_new(1e-6, 0.0);
my $e = gsl_odeiv_evolve_alloc(2);
my $params = { mu => 10 };
my $sys = Math::GSL::ODEIV::gsl_odeiv_system->new(\&func, \&jac, 2, $params );
my $t = 0.0;
my $t1 = 100.0;
my $h = 1e-6;
my $y = [ 1.0, 0.0 ];
gsl_ieee_env_setup;
while ($t < $t1) {
    my $status = gsl_odeiv_evolve_apply ($e, $c, $s, $sys, \$t, $t1, \$h, $y);
    last if $status != $GSL_SUCCESS;
    printf "%.5e %.5e %.5e\n", $t, $y->[0], $y->[1];
}
gsl_odeiv_evolve_free($e);
gsl_odeiv_control_free($c);
gsl_odeiv_step_free($s);

When I run this test script on my Ubuntu laptop I get the expected output.

[leto]

@hakonhagland I would like to get some of your stuff merged, but we ran into test failures on GSL 2.6 stuff. What is the current state of things? I haven't looked at this code in a while.

[morgantomj]

Håkon Hægland, duke, Good job ... I was able to run your code successfully after following your instructions. Then I modified my perl code to mimic yours, and was able to get mine to run as well. I had to hide my gsl-2.5 library, and my Math::GSL, but things seem to be working as advertised. My sincere thanx to the both of you for getting this to work. I'm not sure how you did it, but I'm very happy you solved the problem I was having passing ref's to functions, and fill out the '$sys' structure. Thanx again to both of you.. Tom Morgan

…

On Sun, Feb 2, 2020 at 10:53 PM Duke Leto ***@***.***> wrote: @hakonhagland <https://github.com/hakonhagland> I would like to get some of your stuff merged, but we ran into test failures on GSL 2.6 stuff. What is the current state of things? I haven't looked at this code in a while. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#156?email_source=notifications&email_token=AKSAIXABWZOKN4VMYYQHPJTRA6IMHA5CNFSM4F3H7DNKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEKSM24I#issuecomment-581225841>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AKSAIXH3WCPEFSYTYXMGJQDRA6IMHANCNFSM4F3H7DNA> .

[leto]

@morgantomj all the thanks go to @hakonhagland for this fix, but I am glad that Math::GSL exists and is useful for you! Hopefully we can get this cleaned up and merged and working with recent versions of GSL

[hakonhagland]

@leto Thanks! I have also been occupied by other things the last months so I have not had time to look at the GSL 2.6 issue we discussed earlier. I will see if I can try revivify issue #173 soon, as it would be nice if we could merge this..