Re: [Boost-users] A couple of odeint questions

Subject: Re: [Boost-users] A couple of odeint questions
From: Karsten Ahnert (karsten.ahnert_at_[hidden])
Date: 2012-10-02 15:02:17

Hi,

On 10/02/2012 08:01 PM, Tiago Gehring wrote:
> Hello,
> 
> I've been using the odeint library for some time now  - and I really
> like it, thanks for the great work! I know that the odeint library is
> currently being reviewed but I'm actually writing here to ask a couple
> of questions about my usage scenario, hope that I'm in the right place.
> 
> My usage scenario involves neural networks where I have a set of n
> equations for the N neurons themselves and another m equations for all
> the M synapses (connections) where usually M >> N and the two sets of
> equations depends on each other. What I'm trying to do is to put all
> these N + M equations into one large system and to integrate them at
> once while also keeping the state vectors separated. This is because of
> the different dimensionality of the two sets of equations and because my
> intention is to later use OpenCL to do the calculations (which as I
> understand works better if you use something like a
> vexcl::multivector<double, M> + another vexcl::multivector<double N> and
> I don't see how to pack all of this in one large state vector);

I think vex::multivector< double , M > is something different then a
vector of M elements. I think M play more or less the role of a
dimension. You have to resize such a multivector with length L. a
multivector encapsulates then M vectors with length L.

> So my first question is: can I somehow "combine" two sets of state
> vectors into one thing or integrate two systems in parallel
> (synchronously)? I tried the former (one combined state object) and even
> defined my own combined algebra (a std::pair of std::vectors - later
> vexcl::multivectors) but then stumbled into problems with my value type
> (which I defined as a std::array<double, 2>). It seems that I have to
> define all possible mathematical and comparison operations for my pairs
> (arrays).

This is in principle possible. You can use a fusion sequence, for exmaple

typedef fusion::vector< state_type1, state_type2 > state_type;

You then need a resize_impl specialization which resizes state_type1 and
state_type2. And I think you do not need to write a resize_impl for your
overall state_type. Odeint will dispatch it to the resizers of the
underlying state_types.

For the algebra part it depends on your state_type's. If they are
ordinary range you need a algebra which calls the range_algebra on all
element of the fusion sequence. If your type supports expression
template (like the types from boost.ublas, mtl) you can simply use the
odeint::fusion_algebra. We have done something similar with viennacl,
see
https://github.com/ddemidov/vexcl_odeint_paper/blob/master/src/lorenz_ensemble/viennacl_lorenz_ensemble.cpp

Your value_type should not be a pair. It should be double. Nevertheless,
you can have a pair as state_type or vector of pairs as state_type. In
the tutorial we show how one can define a basic type Point with all
needed operators which can be used in vector< Point > for the state_type.

> 
> And my second question: my problem involves integrating the system until
> one of the neurons fires - at which point it's state variable (membrane
> potential) jumps to a fixed value, i.e., I have a sharp discontinuity in
> the state variable. My way of dealing with this right now was to
> implement a custom integration function (based on the adaptive
> integration one) which has a built in function that checks if the
> threshold was reached at every integration step. When this happens I
> break the integration, reset the state of the neuron in question to the
> fixed value and start integrating again. I wonder if it's possible to
> solve this in an easier way... Is it possible for example to change the
> state on the fly and can the system handle discontinuities? (i.e. can I
> integrate without stopping each time a neuron fires).

A custom integration routine is maybe not a bad idea. Some people
already suggested such a feature. I think it is time for us to seriously
think about implementing this kind of controller.

But odeint supports iterators, have a look at the iterators section in
odeint docs and at the iterators examples. In principle you can do
something like

boost::find_if( make_const_step_range(
    stepper , sys , state , t_start , t_end , dt ) ,
    []( const state_type &x ) { // your break condition } );

> 
> Thanks a lot,
> Tiago
> 
> 
> _______________________________________________
> Boost-users mailing list
> Boost-users_at_[hidden]
> http://lists.boost.org/mailman/listinfo.cgi/boost-users
>