* @see ca.nengo.dynamics.Integrator#integrate(ca.nengo.dynamics.DynamicalSystem, ca.nengo.util.TimeSeries)
*/
public TimeSeries integrate(DynamicalSystem system, TimeSeries input) {
MU.VectorExpander times = new MU.VectorExpander();
MU.MatrixExpander values = new MU.MatrixExpander();
LinearInterpolatorND interpolator = new LinearInterpolatorND(input);
float t0 = input.getTimes()[0];
float tfinal = input.getTimes()[input.getTimes().length - 1];
float hmax = (tfinal - t0) / 2.5f;
float hmin = (tfinal - t0) / 1e9f;
float h = (tfinal - t0) / 100f; //initial guess at step size
float t = t0;
float[] x = system.getState();
times.add(t);
values.add(x);
float[][] k = new float[7][]; //7 stages for each step (although one of them is shared by adjacent steps)
for (int i = 0; i < k.length; i++) {
k[i] = new float[x.length];
}
//Compute the first stage prior to the main loop (subsequently the first stage is assigned from the previous
//step's last stage).
float[] u = interpolator.interpolate(t);
k[0] = system.f(t, u);
while (t < tfinal && h >= hmin) {
if (t + h > tfinal) h = tfinal - t;
for (int j = 0; j < 6; j++) {
float stageTime = t + c[j+1]*h;
u = interpolator.interpolate(stageTime);
float[] ka = new float[x.length];
for (int q = 0; q < x.length; q++) {
for (int r = 0; r <= j; r++) {
ka[q] += k[r][q] * a[j+1][r]; //note: this doesn't look like a matrix-vector product because our k is transposed