Fixed wrong array dimensions in secondary equations handling.

JIRA: MATH-961

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1463684 13f79535-47bb-0310-9956-ffa450edef68

src/changes/changes.xml

@@ -78,6 +78,9 @@ Users are encouraged to upgrade to this version as this release not
   2. A few methods in the FastMath class are in fact slower that their
   counterpart in either Math or StrictMath (cf. MATH-740 and MATH-901).
 ">
+      <action dev="luc" type="fix" issue="MATH-961" >
+        Fixed wrong array dimensions in secondary equations handling in some cases.
+      </action>
       <action dev="luc" type="fix" issue="MATH-960" >
         Fixed missing side effects of secondary equations on main state in
         Ordinary Differential Equations integration.

src/main/java/org/apache/commons/math3/ode/AbstractIntegrator.java

@@ -206,6 +206,22 @@ public abstract class AbstractIntegrator implements FirstOrderIntegrator {
         this.expandable = equations;
     }
 
+    /** Get the differential equations to integrate.
+     * @return differential equations to integrate
+     * @since 3.2
+     */
+    protected ExpandableStatefulODE getExpandable() {
+        return expandable;
+    }
+
+    /** Get the evaluations counter.
+     * @return evaluations counter
+     * @since 3.2
+     */
+    protected Incrementor getEvaluationsCounter() {
+        return evaluations;
+    }
+
     /** {@inheritDoc} */
     public double integrate(final FirstOrderDifferentialEquations equations,
                             final double t0, final double[] y0, final double t, final double[] y)
@@ -340,11 +356,19 @@ public abstract class AbstractIntegrator implements FirstOrderIntegrator {
 
                 // get state at event time
                 interpolator.setInterpolatedTime(eventT);
-                final double[] eventY = interpolator.getInterpolatedState().clone();
+                final double[] eventYPrimary  = interpolator.getInterpolatedState().clone();
+                final double[] eventYComplete = new double[y.length];
+                expandable.getPrimaryMapper().insertEquationData(interpolator.getInterpolatedState(),
+                                                                 eventYComplete);
+                int index = 0;
+                for (EquationsMapper secondary : expandable.getSecondaryMappers()) {
+                    secondary.insertEquationData(interpolator.getInterpolatedSecondaryState(index++),
+                                                 eventYComplete);
+                }
 
                 // advance all event states to current time
                 for (final EventState state : eventsStates) {
-                    state.stepAccepted(eventT, eventY);
+                    state.stepAccepted(eventT, eventYPrimary);
                     isLastStep = isLastStep || state.stop();
                 }
 
@@ -355,18 +379,19 @@ public abstract class AbstractIntegrator implements FirstOrderIntegrator {
 
                 if (isLastStep) {
                     // the event asked to stop integration
-                    System.arraycopy(eventY, 0, y, 0, y.length);
+                    System.arraycopy(eventYComplete, 0, y, 0, y.length);
                     return eventT;
                 }
 
                 boolean needReset = false;
                 for (final EventState state : eventsStates) {
-                    needReset =  needReset || state.reset(eventT, eventY);
+                    needReset =  needReset || state.reset(eventT, eventYComplete);
                 }
                 if (needReset) {
                     // some event handler has triggered changes that
                     // invalidate the derivatives, we need to recompute them
-                    System.arraycopy(eventY, 0, y, 0, y.length);
+                    interpolator.setInterpolatedTime(eventT);
+                    System.arraycopy(eventYComplete, 0, y, 0, y.length);
                     computeDerivatives(eventT, y, yDot);
                     resetOccurred = true;
                     return eventT;

src/main/java/org/apache/commons/math3/ode/ContinuousOutputModel.java

@@ -333,11 +333,27 @@ public class ContinuousOutputModel
    * Get the state vector of the interpolated point.
    * @return state vector at time {@link #getInterpolatedTime}
    * @exception MaxCountExceededException if the number of functions evaluations is exceeded
+   * @see #getInterpolatedSecondaryState(int)
    */
   public double[] getInterpolatedState() throws MaxCountExceededException {
     return steps.get(index).getInterpolatedState();
   }
 
+  /** Get the interpolated secondary state corresponding to the secondary equations.
+   * @param secondaryStateIndex index of the secondary set, as returned by {@link
+   * org.apache.commons.math3.ode.ExpandableStatefulODE#addSecondaryEquations(
+   * org.apache.commons.math3.ode.SecondaryEquations)
+   * ExpandableStatefulODE.addSecondaryEquations(SecondaryEquations)}
+   * @return interpolated secondary state at the current interpolation date
+   * @see #getInterpolatedState()
+   * @since 3.2
+   * @exception MaxCountExceededException if the number of functions evaluations is exceeded
+   */
+  public double[] getInterpolatedSecondaryState(final int secondaryStateIndex)
+    throws MaxCountExceededException {
+    return steps.get(index).getInterpolatedSecondaryState(secondaryStateIndex);
+  }
+
   /** Compare a step interval and a double.
    * @param time point to locate
    * @param interval step interval

src/main/java/org/apache/commons/math3/ode/MultistepIntegrator.java

@@ -228,10 +228,31 @@ public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator {
 
         // start integration, expecting a InitializationCompletedMarkerException
         try {
-            starter.integrate(new CountingDifferentialEquations(y0.length),
+
-                              t0, y0, t, new double[y0.length]);
+            if (starter instanceof AbstractIntegrator) {
+                ((AbstractIntegrator) starter).integrate(getExpandable(), t);
+            } else {
+                starter.integrate(new FirstOrderDifferentialEquations() {
+
+                    /** {@inheritDoc} */
+                    public int getDimension() {
+                        return getExpandable().getTotalDimension();
+                    }
+
+                    /** {@inheritDoc} */
+                    public void computeDerivatives(double t, double[] y, double[] yDot) {
+                        getExpandable().computeDerivatives(t, y, yDot);
+                    }
+
+                }, t0, y0, t, new double[y0.length]);
+            }
+
         } catch (InitializationCompletedMarkerException icme) { // NOPMD
             // this is the expected nominal interruption of the start integrator
+
+            // count the evaluations used by the starter
+            getEvaluationsCounter().incrementCount(starter.getEvaluations());
+
         }
 
         // remove the specific step handler
@@ -353,20 +374,33 @@ public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator {
                 // first step, we need to store also the beginning of the step
                 interpolator.setInterpolatedTime(prev);
                 t[0] = prev;
-                System.arraycopy(interpolator.getInterpolatedState(), 0,
+                final ExpandableStatefulODE expandable = getExpandable();
-                                 y[0],    0, y[0].length);
+                final EquationsMapper primary = expandable.getPrimaryMapper();
-                System.arraycopy(interpolator.getInterpolatedDerivatives(), 0,
+                primary.insertEquationData(interpolator.getInterpolatedState(), y[count]);
-                                 yDot[0], 0, yDot[0].length);
+                primary.insertEquationData(interpolator.getInterpolatedDerivatives(), yDot[count]);
+                int index = 0;
+                for (final EquationsMapper secondary : expandable.getSecondaryMappers()) {
+                    secondary.insertEquationData(interpolator.getInterpolatedSecondaryState(index), y[count]);
+                    secondary.insertEquationData(interpolator.getInterpolatedSecondaryDerivatives(index), yDot[count]);
+                    ++index;
+                }
             }
 
             // store the end of the step
             ++count;
             interpolator.setInterpolatedTime(curr);
             t[count] = curr;
-            System.arraycopy(interpolator.getInterpolatedState(), 0,
+
-                             y[count],    0, y[count].length);
+            final ExpandableStatefulODE expandable = getExpandable();
-            System.arraycopy(interpolator.getInterpolatedDerivatives(), 0,
+            final EquationsMapper primary = expandable.getPrimaryMapper();
-                             yDot[count], 0, yDot[count].length);
+            primary.insertEquationData(interpolator.getInterpolatedState(), y[count]);
+            primary.insertEquationData(interpolator.getInterpolatedDerivatives(), yDot[count]);
+            int index = 0;
+            for (final EquationsMapper secondary : expandable.getSecondaryMappers()) {
+                secondary.insertEquationData(interpolator.getInterpolatedSecondaryState(index), y[count]);
+                secondary.insertEquationData(interpolator.getInterpolatedSecondaryDerivatives(index), yDot[count]);
+                ++index;
+            }
 
             if (count == t.length - 1) {
 
@@ -411,30 +445,4 @@ public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator {
 
     }
 
-    /** Wrapper for differential equations, ensuring start evaluations are counted. */
-    private class CountingDifferentialEquations implements FirstOrderDifferentialEquations {
-
-        /** Dimension of the problem. */
-        private final int dimension;
-
-        /** Simple constructor.
-         * @param dimension dimension of the problem
-         */
-        public CountingDifferentialEquations(final int dimension) {
-            this.dimension = dimension;
-        }
-
-        /** {@inheritDoc} */
-        public void computeDerivatives(double t, double[] y, double[] dot)
-            throws MaxCountExceededException, DimensionMismatchException {
-            MultistepIntegrator.this.computeDerivatives(t, y, dot);
-        }
-
-        /** {@inheritDoc} */
-        public int getDimension() {
-            return dimension;
-        }
-
-    }
-
 }

src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsBashforthIntegrator.java

@@ -22,6 +22,7 @@ import org.apache.commons.math3.exception.MaxCountExceededException;
 import org.apache.commons.math3.exception.NoBracketingException;
 import org.apache.commons.math3.exception.NumberIsTooSmallException;
 import org.apache.commons.math3.linear.Array2DRowRealMatrix;
+import org.apache.commons.math3.ode.EquationsMapper;
 import org.apache.commons.math3.ode.ExpandableStatefulODE;
 import org.apache.commons.math3.ode.sampling.NordsieckStepInterpolator;
 import org.apache.commons.math3.util.FastMath;
@@ -255,7 +256,14 @@ public class AdamsBashforthIntegrator extends AdamsIntegrator {
             final double stepEnd = stepStart + stepSize;
             interpolator.shift();
             interpolator.setInterpolatedTime(stepEnd);
-            System.arraycopy(interpolator.getInterpolatedState(), 0, y, 0, y0.length);
+            final ExpandableStatefulODE expandable = getExpandable();
+            final EquationsMapper primary = expandable.getPrimaryMapper();
+            primary.insertEquationData(interpolator.getInterpolatedState(), y);
+            int index = 0;
+            for (final EquationsMapper secondary : expandable.getSecondaryMappers()) {
+                secondary.insertEquationData(interpolator.getInterpolatedSecondaryState(index), y);
+                ++index;
+            }
 
             // evaluate the derivative
             computeDerivatives(stepEnd, y, yDot);

src/main/java/org/apache/commons/math3/ode/nonstiff/AdamsMoultonIntegrator.java

@@ -25,6 +25,7 @@ import org.apache.commons.math3.exception.NoBracketingException;
 import org.apache.commons.math3.exception.NumberIsTooSmallException;
 import org.apache.commons.math3.linear.Array2DRowRealMatrix;
 import org.apache.commons.math3.linear.RealMatrixPreservingVisitor;
+import org.apache.commons.math3.ode.EquationsMapper;
 import org.apache.commons.math3.ode.ExpandableStatefulODE;
 import org.apache.commons.math3.ode.sampling.NordsieckStepInterpolator;
 import org.apache.commons.math3.util.FastMath;
@@ -250,7 +251,14 @@ public class AdamsMoultonIntegrator extends AdamsIntegrator {
                 // predict a first estimate of the state at step end (P in the PECE sequence)
                 final double stepEnd = stepStart + stepSize;
                 interpolator.setInterpolatedTime(stepEnd);
-                System.arraycopy(interpolator.getInterpolatedState(), 0, yTmp, 0, y0.length);
+                final ExpandableStatefulODE expandable = getExpandable();
+                final EquationsMapper primary = expandable.getPrimaryMapper();
+                primary.insertEquationData(interpolator.getInterpolatedState(), yTmp);
+                int index = 0;
+                for (final EquationsMapper secondary : expandable.getSecondaryMappers()) {
+                    secondary.insertEquationData(interpolator.getInterpolatedSecondaryState(index), yTmp);
+                    ++index;
+                }
 
                 // evaluate a first estimate of the derivative (first E in the PECE sequence)
                 computeDerivatives(stepEnd, yTmp, yDot);