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

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 + 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). 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_ensemb... 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@lists.boost.org http://lists.boost.org/mailman/listinfo.cgi/boost-users