commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject svn commit: r1176734 [1/2] - in /commons/proper/math/trunk/src: main/java/org/apache/commons/math/exception/util/ main/java/org/apache/commons/math/ode/ main/java/org/apache/commons/math/ode/nonstiff/ main/java/org/apache/commons/math/ode/sampling/ mai...
Date Wed, 28 Sep 2011 05:56:43 GMT
Author: luc
Date: Wed Sep 28 05:56:42 2011
New Revision: 1176734

URL: http://svn.apache.org/viewvc?rev=1176734&view=rev
Log:
added support for secondary state in ODE step interpolators

Added:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/EquationsMapper.java   (with props)
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ExpandableStatefulODE.java   (with props)
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/SecondaryEquations.java   (contents, props changed)
      - copied, changed from r1175683, commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AdditionalEquations.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/sampling/DummyStepInterpolator.java   (contents, props changed)
      - copied, changed from r1175683, commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/sampling/DummyStepInterpolator.java
Removed:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AdditionalEquations.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AdditionalStateAndEquations.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ExpandableFirstOrderDifferentialEquations.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ExpandableFirstOrderIntegrator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ExtendedFirstOrderDifferentialEquations.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/MainStateJacobianWrapper.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/sampling/DummyStepInterpolator.java
Modified:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/exception/util/LocalizedFormats.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AbstractIntegrator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/JacobianMatrices.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterJacobianProvider.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterJacobianWrapper.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterizedODE.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsIntegrator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdaptiveStepsizeIntegrator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince54StepInterpolator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince853StepInterpolator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerStepInterpolator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaStepInterpolator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/sampling/AbstractStepInterpolator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/sampling/NordsieckStepInterpolator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/sampling/StepInterpolator.java
    commons/proper/math/trunk/src/main/resources/META-INF/localization/LocalizedFormats_fr.properties
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/JacobianMatricesTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/EulerStepInterpolatorTest.java

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/exception/util/LocalizedFormats.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/exception/util/LocalizedFormats.java?rev=1176734&r1=1176733&r2=1176734&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/exception/util/LocalizedFormats.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/exception/util/LocalizedFormats.java Wed Sep 28 05:56:42 2011
@@ -315,9 +315,8 @@ public enum LocalizedFormats implements 
     UNABLE_TO_SOLVE_SINGULAR_PROBLEM("unable to solve: singular problem"),
     UNBOUNDED_SOLUTION("unbounded solution"),
     UNKNOWN_MODE("unknown mode {0}, known modes: {1} ({2}), {3} ({4}), {5} ({6}), {7} ({8}), {9} ({10}) and {11} ({12})"),
-    UNKNOWN_ADDITIONAL_EQUATION("unknown additional equation"),
     UNKNOWN_PARAMETER("unknown parameter {0}"),
-    UNMATCHED_ODE_IN_EXTENDED_SET("ode does not match the main ode set in the extended set"),
+    UNMATCHED_ODE_IN_EXPANDED_SET("ode does not match the main ode set in the extended set"),
     CANNOT_PARSE_AS_TYPE("string {0} unparseable (from position {1}) as an object of type {2}"), /* keep */
     CANNOT_PARSE("string {0} unparseable (from position {1})"), /* keep */
     UNPARSEABLE_3D_VECTOR("unparseable 3D vector: \"{0}\""),

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AbstractIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AbstractIntegrator.java?rev=1176734&r1=1176733&r2=1176734&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AbstractIntegrator.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AbstractIntegrator.java Wed Sep 28 05:56:42 2011
@@ -47,7 +47,7 @@ import org.apache.commons.math.util.Math
  * @version $Id$
  * @since 2.0
  */
-public abstract class AbstractIntegrator implements ExpandableFirstOrderIntegrator {
+public abstract class AbstractIntegrator implements FirstOrderIntegrator {
 
     /** Step handler. */
     protected Collection<StepHandler> stepHandlers;
@@ -77,7 +77,7 @@ public abstract class AbstractIntegrator
     private Incrementor evaluations;
 
     /** Differential equations to integrate. */
-    private transient ExpandableFirstOrderDifferentialEquations equations;
+    private transient ExpandableStatefulODE equations;
 
     /** Build an instance.
      * @param name name of the method
@@ -185,21 +185,59 @@ public abstract class AbstractIntegrator
         evaluations.resetCount();
     }
 
-    /** Set the differential equations.
-     * @param equations differential equations to integrate
-     * @see #computeDerivatives(double, double[], double[])
+    /** Set the equations.
+     * @param equations equations to set
      */
-    protected void setEquations(final ExpandableFirstOrderDifferentialEquations equations) {
+    protected void setEquations(final ExpandableStatefulODE equations) {
         this.equations = equations;
     }
 
     /** {@inheritDoc} */
-    public double integrate(FirstOrderDifferentialEquations equations,
-                            double t0, double[] y0, double t, double[] y)
+    public double integrate(final FirstOrderDifferentialEquations equations,
+                            final double t0, final double[] y0, final double t, final double[] y)
         throws MathIllegalStateException, MathIllegalArgumentException {
-        return integrate(new ExpandableFirstOrderDifferentialEquations(equations), t0, y0, t, y);
+
+        if (y0.length != equations.getDimension()) {
+            throw new DimensionMismatchException(y0.length, equations.getDimension());
+        }
+        if (y.length != equations.getDimension()) {
+            throw new DimensionMismatchException(y.length, equations.getDimension());
+        }
+
+        // prepare expandable stateful equations
+        final ExpandableStatefulODE expandable = new ExpandableStatefulODE(equations);
+        expandable.setTime(t0);
+        expandable.setPrimaryState(y0);
+
+        // perform integration
+        integrate(expandable, t);
+
+        // extract results back from the stateful equations
+        System.arraycopy(expandable.getPrimaryState(), 0, y, 0, y.length);
+        return expandable.getTime();
+
     }
 
+    /** Integrate a set of differential equations up to the given time.
+     * <p>This method solves an Initial Value Problem (IVP).</p>
+     * <p>The set of differential equations is composed of a main set, which
+     * can be extended by some sets of secondary equations. The set of
+     * equations must be already set up with initial time and partial states.
+     * At integration completion, the final time and partial states will be
+     * available in the same object.</p>
+     * <p>Since this method stores some internal state variables made
+     * available in its public interface during integration ({@link
+     * #getCurrentSignedStepsize()}), it is <em>not</em> thread-safe.</p>
+     * @param equations complete set of differential equations to integrate
+     * @param t target time for the integration
+     * (can be set to a value smaller than <code>t0</code> for backward integration)
+     * @throws MathIllegalStateException if the integrator cannot perform integration
+     * @throws MathIllegalArgumentException if integration parameters are wrong (typically
+     * too small integration span)
+     */
+    public abstract void integrate(ExpandableStatefulODE equations, double t)
+        throws MathIllegalStateException, MathIllegalArgumentException;
+
     /** Compute the derivatives and check the number of evaluations.
      * @param t current value of the independent <I>time</I> variable
      * @param y array containing the current value of the state vector
@@ -335,33 +373,19 @@ public abstract class AbstractIntegrator
 
     }
 
-    /** Perform some sanity checks on the integration parameters.
-     * @param ode differential equations set
-     * @param t0 start time
-     * @param y0 state vector at t0
+    /** Check the integration span.
      * @param t target time for the integration
-     * @param y placeholder where to put the state vector
-     * @exception DimensionMismatchException if some inconsistency is detected
      * @exception NumberIsTooSmallException if integration span is too small
      */
-    protected void sanityChecks(final ExpandableFirstOrderDifferentialEquations ode,
-                                final double t0, final double[] y0,
-                                final double t, final double[] y)
-        throws DimensionMismatchException, NumberIsTooSmallException {
-
-        if (ode.getMainSetDimension() != y0.length) {
-            throw new DimensionMismatchException(ode.getDimension(), y0.length);
-        }
-
-        if (ode.getMainSetDimension() != y.length) {
-            throw new DimensionMismatchException(ode.getDimension(), y.length);
-        }
+    protected void sanityChecks(final ExpandableStatefulODE equations, final double t)
+        throws NumberIsTooSmallException {
 
-        if (FastMath.abs(t - t0) <= 1.0e-12 * FastMath.max(FastMath.abs(t0), FastMath.abs(t))) {
+        final double threshold = 1000 * FastMath.ulp(FastMath.max(FastMath.abs(equations.getTime()),
+                                                                  FastMath.abs(t)));
+        final double dt = FastMath.abs(equations.getTime() - t);
+        if (dt <= threshold) {
             throw new NumberIsTooSmallException(LocalizedFormats.TOO_SMALL_INTEGRATION_INTERVAL,
-                                                FastMath.abs(t - t0),
-                                                1.0e-12 * FastMath.max(FastMath.abs(t0), FastMath.abs(t)),
-                                                false);
+                                                dt, threshold, false);
         }
 
     }

Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/EquationsMapper.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/EquationsMapper.java?rev=1176734&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/EquationsMapper.java (added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/EquationsMapper.java Wed Sep 28 05:56:42 2011
@@ -0,0 +1,98 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.math.ode;
+
+import java.io.Serializable;
+
+import org.apache.commons.math.exception.DimensionMismatchException;
+
+/**
+ * Class mapping the part of a complete state or derivative that pertains
+ * to a specific differential equation.
+ * <p>
+ * Instances of this class are guaranteed to be immutable.
+ * </p>
+ * @see SecondaryEquations
+ * @version $Id$
+ * @since 3.0
+ */
+public class EquationsMapper implements Serializable {
+
+    /** Serializable UID. */
+    private static final long serialVersionUID = 20110925L;
+
+    /** Index of the first equation element in complete state arrays. */
+    final int firstIndex;
+
+    /** Dimension of the secondary state parameters. */
+    final int dimension;
+
+    /** simple constructor.
+     * @param firstIndex index of the first equation element in complete state arrays
+     * @param dimension dimension of the secondary state parameters
+     */
+    public EquationsMapper(final int firstIndex, final int dimension) {
+        this.firstIndex = firstIndex;
+        this.dimension  = dimension;
+    }
+
+    /** Get the index of the first equation element in complete state arrays.
+     * @return index of the first equation element in complete state arrays
+     */
+    public int getFirstIndex() {
+        return firstIndex;
+    }
+
+    /** Get the dimension of the secondary state parameters.
+     * @return dimension of the secondary state parameters
+     */
+    public int getDimension() {
+        return dimension;
+    }
+
+    /** Extract equation data from a complete state or derivative array.
+     * @param complete complete state or derivative array from which
+     * equation data should be retrieved
+     * @param equationData placeholder where to put equation data
+     * @throws DimensionMismatchException if the dimension of the equation data does not
+     * match the mapper dimension
+     */
+    public void extractEquationData(double[] complete, double[] equationData)
+        throws DimensionMismatchException {
+        if (equationData.length != dimension) {
+            throw new DimensionMismatchException(equationData.length, dimension);
+        }
+        System.arraycopy(complete, firstIndex, equationData, 0, dimension);
+    }
+
+    /** Insert equation data into a complete state or derivative array.
+     * @param equationData equation data to be inserted into the complete array
+     * @param complete placeholder where to put equation data (only the
+     * part corresponding to the equation will be overwritten)
+     * @throws DimensionMismatchException if the dimension of the equation data does not
+     * match the mapper dimension
+     */
+    public void insertEquationData(double[] equationData, double[] complete)
+        throws DimensionMismatchException {
+        if (equationData.length != dimension) {
+            throw new DimensionMismatchException(equationData.length, dimension);
+        }
+        System.arraycopy(equationData, 0, complete, firstIndex, dimension);
+    }
+
+}

Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/EquationsMapper.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/EquationsMapper.java
------------------------------------------------------------------------------
    svn:keywords = "Author Date Id Revision"

Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ExpandableStatefulODE.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ExpandableStatefulODE.java?rev=1176734&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ExpandableStatefulODE.java (added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ExpandableStatefulODE.java Wed Sep 28 05:56:42 2011
@@ -0,0 +1,327 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.math.ode;
+
+import java.util.ArrayList;
+import java.util.List;
+
+import org.apache.commons.math.exception.DimensionMismatchException;
+
+
+/**
+ * This class represents a combined set of first order differential equations,
+ * with at least a primary set of equations expandable by some sets of secondary
+ * equations.
+ * <p>
+ * One typical use case is the computation of the Jacobian matrix for some ODE.
+ * In this case, the primary set of equations corresponds to the raw ODE, and we
+ * add to this set another bunch of secondary equations which represent the Jacobian
+ * matrix of the primary set.
+ * </p>
+ * <p>
+ * We want the integrator to use <em>only</em> the primary set to estimate the
+ * errors and hence the step sizes. It should <em>not</em> use the secondary
+ * equations in this computation. The {@link AbstractIntegrator integrator} will
+ * be able to know where the primary set ends and so where the secondary sets begin.
+ * </p>
+ *
+ * @see FirstOrderDifferentialEquations
+ * @see JacobianMatrices
+ *
+ * @version $Id$
+ * @since 3.0
+ */
+
+public class ExpandableStatefulODE {
+
+    /** Primary differential equation. */
+    private final FirstOrderDifferentialEquations primary;
+
+    /** Mapper for primary equation. */
+    private final EquationsMapper primaryMapper;
+
+    /** Time. */
+    private double time;
+
+    /** State. */
+    private final double[] primaryState;
+
+    /** State derivative. */
+    private final double[] primaryStateDot;
+
+    /** Components of the expandable ODE. */
+    private List<SecondaryComponent> components;
+
+    /** Build an expandable set from its primary ODE set.
+     * @param ode the primary set of differential equations to be integrated.
+     */
+    public ExpandableStatefulODE(final FirstOrderDifferentialEquations primary) {
+        final int n          = primary.getDimension();
+        this.primary         = primary;
+        this.primaryMapper   = new EquationsMapper(0, n);
+        this.time            = Double.NaN;
+        this.primaryState    = new double[n];
+        this.primaryStateDot = new double[n];
+        this.components      = new ArrayList<ExpandableStatefulODE.SecondaryComponent>();
+    }
+
+    /** Get the primary set of differential equations.
+     * @return primary set of differential equations
+     */
+    public FirstOrderDifferentialEquations getPrimary() {
+        return primary;
+    }
+
+    /** Return the dimension of the complete set of equations.
+     * <p>
+     * The complete set of equations correspond to the primary set plus all secondary sets.
+     * </p>
+     * @return dimension of the complete set of equations
+     */
+    public int getTotalDimension() {
+        if (components.isEmpty()) {
+            // there are no secondary equations, the complete set is limited to the primary set
+            return primaryMapper.getDimension();
+        } else {
+            // there are secondary equations, the complete set ends after the last set
+            final EquationsMapper lastMapper = components.get(components.size() - 1).mapper;
+            return lastMapper.getFirstIndex() + lastMapper.getDimension();
+        }
+    }
+
+    /** Get the current time derivative of the complete state vector.
+     * @param t current value of the independent <I>time</I> variable
+     * @param y array containing the current value of the complete state vector
+     * @param yDot placeholder array where to put the time derivative of the complete state vector
+     */
+    public void computeDerivatives(final double t, final double[] y, final double[] yDot) {
+
+        // compute derivatives of the primary equations
+        primaryMapper.extractEquationData(y, primaryState);
+        primary.computeDerivatives(t, primaryState, primaryStateDot);
+        primaryMapper.insertEquationData(primaryStateDot, yDot);
+
+        // Add contribution for secondary equations
+        for (final SecondaryComponent component : components) {
+            component.mapper.extractEquationData(y, component.state);
+            component.equation.computeDerivatives(t, primaryState, primaryStateDot,
+                                                  component.state, component.stateDot);
+            component.mapper.insertEquationData(component.stateDot, yDot);
+        }
+
+    }
+
+    /** Add a set of secondary equations to be integrated along with the primary set.
+     * @param secondary secondary equations set
+     * @return index of the secondary equation in the expanded state
+     */
+    public int addSecondaryEquations(final SecondaryEquations secondary) {
+
+        final int firstIndex;
+        if (components.isEmpty()) {
+            // lazy creation of the components list
+            components = new ArrayList<ExpandableStatefulODE.SecondaryComponent>();
+            firstIndex = primary.getDimension();
+        } else {
+            final SecondaryComponent last = components.get(components.size() - 1);
+            firstIndex = last.mapper.getFirstIndex() + last.mapper.getDimension();
+        }
+
+        components.add(new SecondaryComponent(secondary, firstIndex));
+
+        return components.size() - 1;
+
+    }
+
+    /** Get an equations mapper for the primary equations set.
+     * @return mapper for the primary set
+     * @see #getSecondaryMappers()
+     */
+    public EquationsMapper getPrimaryMapper() {
+        return primaryMapper;
+    }
+
+    /** Get the equations mappers for the secondary equations sets.
+     * @return equations mappers for the secondary equations sets
+     * @see #getPrimaryMapper()
+     */
+    public EquationsMapper[] getSecondaryMappers() {
+        final EquationsMapper[] mappers = new EquationsMapper[components.size()];
+        for (int i = 0; i < mappers.length; ++i) {
+            mappers[i] = components.get(i).mapper;
+        }
+        return mappers;
+    }
+
+    /** Set current time.
+     * @param time current time
+     */
+    public void setTime(final double time) {
+        this.time = time;
+    }
+
+    /** Get current time.
+     * @return current time
+     */
+    public double getTime() {
+        return time;
+    }
+
+    /** Set primary part of the current state.
+     * @param primaryState primary part of the current state
+     * @throws DimensionMismatchException if the dimension of the array does not
+     * match the primary set
+     */
+    public void setPrimaryState(final double[] primaryState) throws DimensionMismatchException {
+
+        // safety checks
+        if (primaryState.length != this.primaryState.length) {
+            throw new DimensionMismatchException(primaryState.length, this.primaryState.length);
+        }
+
+        // set the data
+        System.arraycopy(primaryState, 0, this.primaryState, 0, primaryState.length);
+
+    }
+
+    /** Get primary part of the current state.
+     * @return primary part of the current state
+     */
+    public double[] getPrimaryState() {
+        return primaryState.clone();
+    }
+
+    /** Get primary part of the current state derivative.
+     * @return primary part of the current state derivative
+     */
+    public double[] getPrimaryStateDot() {
+        return primaryStateDot.clone();
+    }
+
+    /** Set secondary part of the current state.
+     * @param index index of the part to set as returned by {@link
+     * #addSecondaryEquations(SecondaryEquations)}
+     * @param secondaryState secondary part of the current state
+     * @throws DimensionMismatchException if the dimension of the partial state does not
+     * match the selected equations set dimension
+     */
+    public void setSecondaryState(final int index, final double[] secondaryState)
+        throws DimensionMismatchException {
+
+        // get either the secondary state
+        double[] localArray = components.get(index).state;
+
+        // safety checks
+        if (secondaryState.length != localArray.length) {
+            throw new DimensionMismatchException(secondaryState.length, localArray.length);
+        }
+
+        // set the data
+        System.arraycopy(secondaryState, 0, localArray, 0, secondaryState.length);
+
+    }
+
+    /** Get secondary part of the current state.
+     * @param index index of the part to set as returned by {@link
+     * #addSecondaryEquations(SecondaryEquations)}
+     * @return secondary part of the current state
+     */
+    public double[] getSecondaryState(final int index) {
+        return components.get(index).state.clone();
+    }
+
+    /** Get secondary part of the current state derivative.
+     * @param index index of the part to set as returned by {@link
+     * #addSecondaryEquations(SecondaryEquations)}
+     * @return secondary part of the current state derivative
+     */
+    public double[] getSecondaryStateDot(final int index) {
+        return components.get(index).stateDot.clone();
+    }
+
+    /** Set the complete current state.
+     * @param completeState complete current state to copy data from
+     * @throws DimensionMismatchException if the dimension of the complete state does not
+     * match the complete equations sets dimension
+     */
+    public void setCompleteState(final double[] completeState)
+        throws DimensionMismatchException {
+
+        // safety checks
+        if (completeState.length != getTotalDimension()) {
+            throw new DimensionMismatchException(completeState.length, getTotalDimension());
+        }
+
+        // set the data
+        primaryMapper.extractEquationData(completeState, primaryState);
+        for (final SecondaryComponent component : components) {
+            component.mapper.extractEquationData(completeState, component.state);
+        }
+
+    }
+
+    /** Get the complete current state.
+     * @return complete current state
+     * @throws DimensionMismatchException if the dimension of the complete state does not
+     * match the complete equations sets dimension
+     */
+    public double[] getCompleteState() {
+
+        // allocate complete array
+        double[] completeState = new double[getTotalDimension()];
+
+        // set the data
+        primaryMapper.insertEquationData(primaryState, completeState);
+        for (final SecondaryComponent component : components) {
+            component.mapper.insertEquationData(component.state, completeState);
+        }
+
+        return completeState;
+
+    }
+
+    /** Components of the compound stateful ODE. */
+    private static class SecondaryComponent {
+
+        
+        /** Secondary differential equation. */
+        private final SecondaryEquations equation;
+
+        /** Mapper between local and complete arrays. */
+        private final EquationsMapper mapper;
+
+        /** State. */
+        private final double[] state;
+
+        /** State derivative. */
+        private final double[] stateDot;
+
+        /** Simple constructor.
+         * @param equation secondary differential equation
+         * @param first index index to use for the first element in the complete arrays
+         */
+        public SecondaryComponent(final SecondaryEquations equation, final int firstIndex) {
+            final int n   = equation.getDimension();
+            this.equation = equation;
+            mapper        = new EquationsMapper(firstIndex, n);
+            state         = new double[n];
+            stateDot      = new double[n];
+        }
+
+    }
+
+}

Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ExpandableStatefulODE.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ExpandableStatefulODE.java
------------------------------------------------------------------------------
    svn:keywords = "Author Date Id Revision"

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/JacobianMatrices.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/JacobianMatrices.java?rev=1176734&r1=1176733&r2=1176734&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/JacobianMatrices.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/JacobianMatrices.java Wed Sep 28 05:56:42 2011
@@ -25,25 +25,25 @@ import org.apache.commons.math.exception
 import org.apache.commons.math.exception.util.LocalizedFormats;
 
 /**
- * This class defines a set of {@link AdditionalEquations additional equations} to
- * compute the jacobian matrices with respect to the initial state vector and, if
- * any, to some parameters of the main ODE set.
+ * This class defines a set of {@link SecondaryEquations secondary equations} to
+ * compute the Jacobian matrices with respect to the initial state vector and, if
+ * any, to some parameters of the primary ODE set.
  * <p>
- * It is intended to be packed into an {@link ExpandableFirstOrderDifferentialEquations}
- * in conjunction with a main set of ODE, which may be:
+ * It is intended to be packed into an {@link ExpandableStatefulODE}
+ * in conjunction with a primary set of ODE, which may be:
  * <ul>
  * <li>a {@link FirstOrderDifferentialEquations}</li>
  * <li>a {@link MainStateJacobianProvider}</li>
  * </ul>
- * In order to compute jacobian matrices with respect to some parameters of the
- * main ODE set, the following parameter jacobian providers may be set:
+ * In order to compute Jacobian matrices with respect to some parameters of the
+ * primary ODE set, the following parameter Jacobian providers may be set:
  * <ul>
  * <li>a {@link ParameterJacobianProvider}</li>
  * <li>a {@link ParameterizedODE}</li>
  * </ul>
  * </p>
  *
- * @see ExpandableFirstOrderDifferentialEquations
+ * @see ExpandableStatefulODE
  * @see FirstOrderDifferentialEquations
  * @see MainStateJacobianProvider
  * @see ParameterJacobianProvider
@@ -52,288 +52,242 @@ import org.apache.commons.math.exception
  * @version $Id$
  * @since 3.0
  */
-public class JacobianMatrices implements AdditionalEquations {
+public class JacobianMatrices {
 
     /** Expandable first order differential equation. */
-    private ExpandableFirstOrderDifferentialEquations efode;
+    private ExpandableStatefulODE efode;
 
-    /** FODE without exact main jacobian computation skill. */
-    private FirstOrderDifferentialEquations fode = null;
+    /** Index of the instance in the expandable set. */
+    private int index;
 
-    /** FODE with exact main jacobian computation skill. */
-    private MainStateJacobianProvider jode = null;
+    /** FODE with exact primary Jacobian computation skill. */
+    private MainStateJacobianProvider jode;
 
-    /** FODE without exact parameter jacobian computation skill. */
-    private ParameterizedODE pode = null;
-
-    /** FODE with exact parameter jacobian computation skill. */
-    private List<ParameterJacobianProvider> pjp = new ArrayList<ParameterJacobianProvider>();;
-
-    /** List of parameters selected for parameter jacobian computation. */
-    private List<ParameterConfiguration> selectedParameters = null;
+    /** FODE without exact parameter Jacobian computation skill. */
+    private ParameterizedODE pode;
 
     /** Main state vector dimension. */
     private int stateDim;
 
+    /** Selected parameters for parameter Jacobian computation. */
+    private ParameterConfiguration[] selectedParameters;
+
+    /** FODE with exact parameter Jacobian computation skill. */
+    private List<ParameterJacobianProvider> jacobianProviders;
+
     /** Parameters dimension. */
-    private int paramDim = 0;
+    private int paramDim;
 
-    /** Current main state jacobian matrix in a row. */
-    private double[] mainJacobianInARow;
+    /** Boolean for selected parameters consistency. */
+    private boolean dirtyParameter;
 
-    /** Current parameters jacobian matrices in a row. */
-    private double[] parameterJacobiansInARow = null;
+    /** State and parameters Jacobian matrices in a row. */
+    private double[] matricesData;
 
-    /** Step used for finite difference computation of jacobian matrix
-     *  w.r.t. main state vector. */
-    private double[] hY = null;
+    /** Simple constructor for a secondary equations set computing Jacobian matrices.
+     * <p>
+     * Parameters must belong to the supported ones given by {@link
+     * Parameterizable#getParametersNames()}, so the primary set of differential
+     * equations must be {@link Parameterizable}.
+     * </p>
+     * <p>Note that each selection clears the previous selected parameters.</p>
+     *
+     * @param fode the primary first order differential equations set to extend
+     * @param hY step used for finite difference computation with respect to state vector
+     * @param parameters parameters to consider for Jacobian matrices processing
+     * (may be null if parameters Jacobians is not desired)
+     * @exception MathIllegalArgumentException if one parameter is not supported
+     * or there is a dimension mismatch with {@code hY}
+     */
+    public JacobianMatrices(final FirstOrderDifferentialEquations fode, final double[] hY,
+                            final String... parameters)
+        throws MathIllegalArgumentException {
+        this(new MainStateJacobianWrapper(fode, hY), parameters);
+    }
 
-    /** Boolean for fode consistency. */
-    private boolean dirtyMainState = false;
+    /** Simple constructor for a secondary equations set computing Jacobian matrices.
+     * <p>
+     * Parameters must belong to the supported ones given by {@link
+     * Parameterizable#getParametersNames()}, so the primary set of differential
+     * equations must be {@link Parameterizable}.
+     * </p>
+     * <p>Note that each selection clears the previous selected parameters.</p>
+     *
+     * @param jode the primary first order differential equations set to extend
+     * @param parameters parameters to consider for Jacobian matrices processing
+     * (may be null if parameters Jacobians is not desired)
+     * @exception MathIllegalArgumentException if one parameter is not supported
+     */
+    public JacobianMatrices(final MainStateJacobianProvider jode,
+                            final String... parameters)
+        throws MathIllegalArgumentException {
 
-    /** Boolean for selected parameters consistency. */
-    private boolean dirtyParameter = false;
+        this.efode = null;
+        this.index = -1;
 
-    /** Simple constructor for an additional equations set computing jacobian matrices.
-     * <p>This additional equations set is added internally to the expandable
-     * first order differential equations set thanks to the
-     * {@link ExpandableFirstOrderDifferentialEquations#addAdditionalEquations(AdditionalEquations)}
-     * method.
-     * @param extended the expandable first order differential equations set
-     * @param jode the main first order differential equations set to extend
-     * @exception IllegalArgumentException if jode does not match the main set to be extended given by
-     *            {@link ExpandableFirstOrderDifferentialEquations#getMainSet() extended.getMainSet()}
-     */
-    public JacobianMatrices(final ExpandableFirstOrderDifferentialEquations extended,
-                            final MainStateJacobianProvider jode)
-        throws IllegalArgumentException {
-
-        checkCompatibility(extended, jode);
-
-        efode = extended;
-        stateDim = efode.getMainSetDimension();
-        mainJacobianInARow = new double[stateDim * stateDim];
         this.jode = jode;
-        efode.addAdditionalEquations(this);
-        setInitialMainStateJacobian();
+        this.pode = null;
+
+        this.stateDim = jode.getDimension();
+
+        if (parameters == null) {
+            selectedParameters = null;
+            paramDim = 0;
+        } else {
+            this.selectedParameters = new ParameterConfiguration[parameters.length];
+            for (int i = 0; i < parameters.length; ++i) {
+                selectedParameters[i] = new ParameterConfiguration(parameters[i], Double.NaN);
+            }
+            paramDim = parameters.length;
+        }
+        this.dirtyParameter = false;
+
+        this.jacobianProviders = new ArrayList<ParameterJacobianProvider>();
+
+        // set the default initial state Jacobian to the identity
+        // and the default initial parameters Jacobian to the null matrix
+        matricesData = new double[(stateDim + paramDim) * stateDim];
+        for (int i = 0; i < stateDim; ++i) {
+            matricesData[i * (stateDim + 1)] = 1.0;
+        }
+
     }
 
-    /** Simple constructor for an additional equations set computing jacobian matrices.
-     * <p>This additional equations set is added internally to the expandable
-     * first order differential equations set thanks to the
-     * {@link ExpandableFirstOrderDifferentialEquations#addAdditionalEquations(AdditionalEquations)}
-     * method.
-     * @param extended the expandable first order differential equations set
-     * @param fode the main first order differential equations set to extend
-     * @exception IllegalArgumentException if fode does not match the main set to be extended given by
-     *            {@link ExpandableFirstOrderDifferentialEquations#getMainSet() extended.getMainSet()}
+    /** Register the variational equations for the Jacobians matrices to the expandable set.
+     * @exception MathIllegalArgumentException if the primary set of the expandable set does
+     * not match the one used to build the instance
+     * @see ExpandableStatefulODE#addSecondaryEquations(SecondaryEquations)
      */
-    public JacobianMatrices(final ExpandableFirstOrderDifferentialEquations extended,
-                            final FirstOrderDifferentialEquations fode)
-        throws IllegalArgumentException {
+    public void registerVariationalEquations(final ExpandableStatefulODE expandable)
+        throws MathIllegalArgumentException {
 
-        checkCompatibility(extended, fode);
+        // safety checks
+        final FirstOrderDifferentialEquations ode = (jode instanceof MainStateJacobianWrapper) ?
+                                                    ((MainStateJacobianWrapper) jode).ode :
+                                                    jode;
+        if (expandable.getPrimary() != ode) {
+            throw new MathIllegalArgumentException(LocalizedFormats.UNMATCHED_ODE_IN_EXPANDED_SET);
+        }
+
+        efode = expandable;
+        index = efode.addSecondaryEquations(new JacobiansSecondaryEquations());
+        efode.setSecondaryState(index, matricesData);
 
-        efode = extended;
-        stateDim = efode.getMainSetDimension();
-        mainJacobianInARow = new double[stateDim * stateDim];
-        this.fode = fode;
-        dirtyMainState = true;
-        efode.addAdditionalEquations(this);
-        setInitialMainStateJacobian();
     }
 
-    /** Add a parameter jacobian provider.
-     * @param pjp the parameter jacobian provider to compute exactly the parameter jacobian matrix
+    /** Add a parameter Jacobian provider.
+     * @param provider the parameter Jacobian provider to compute exactly the parameter Jacobian matrix
      */
-    public void setParameterJacobianProvider(final ParameterJacobianProvider pjp) {
-        this.pjp.add(pjp);
+    public void addParameterJacobianProvider(final ParameterJacobianProvider provider) {
+        jacobianProviders.add(provider);
     }
 
-    /** Add a parameter jacobian provider.
-     * @param pjp the parameterized ODE to compute by finite difference the parameter jacobian matrix
+    /** Add a parameter Jacobian provider.
+     * @param pode the parameterized ODE to compute the parameter Jacobian matrix using finite differences 
      */
     public void setParameterizedODE(final ParameterizedODE pode) {
         this.pode = pode;
         dirtyParameter = true;
     }
 
-    /** Select the parameters to consider for jacobian matrices processing.
-     * <p>
-     * Parameters must belong to the supported ones given by {@link
-     * Parameterizable#getParametersNames()}, so the main set of differential
-     * equations must be {@link Parameterizable}.
-     * </p>
-     * <p>Note that each selection clears the previous selected parameters.</p>
-     *
-     * @param parameters parameters to consider for jacobian matrices processing
-     * @exception IllegalArgumentException if one parameter is not supported
-     */
-    public void selectParameters(final String... parameters) throws IllegalArgumentException {
-        
-        selectedParameters = new ArrayList<ParameterConfiguration>();
-        for (String param : parameters) {
-            selectedParameters.add(new ParameterConfiguration(param, Double.NaN));
-        }
-        paramDim = parameters.length;
-        parameterJacobiansInARow = new double[paramDim * stateDim];
-        setInitialParameterJacobians();
-    }
-
     /** Set the step associated to a parameter in order to compute by finite
-     *  difference the jacobian matrix.
+     *  difference the Jacobian matrix.
      * <p>
-     * Needed if and only if the main ODE set is a {@link ParameterizedODE}
-     * and the parameter has been {@link #selectParameters(String ...) selected}
+     * Needed if and only if the primary ODE set is a {@link ParameterizedODE}.
      * </p>
      * <p>
-     * For pval, a non zero value of the parameter, pval * Math.sqrt(MathUtils.EPSILON)
-     * is a reasonable value for such a step.
+     * Given a non zero parameter value pval for the parameter, a reasonable value
+     * for such a step is {@code pval * FastMath.sqrt(MathUtils.EPSILON)}.
      * </p>
      * <p>
-     * A zero value for such a step doesn't enable to compute the parameter jacobian matrix.
+     * A zero value for such a step doesn't enable to compute the parameter Jacobian matrix.
      * </p>
-     * @param parameter parameter to consider for jacobian processing
-     * @param hP step for jacobian finite difference computation w.r.t. the specified parameter
+     * @param parameter parameter to consider for Jacobian processing
+     * @param hP step for Jacobian finite difference computation w.r.t. the specified parameter
      * @see ParameterizedODE
      * @exception IllegalArgumentException if the parameter is not supported
      */
     public void setParameterStep(final String parameter, final double hP) {
 
-        boolean found = false;
         for (ParameterConfiguration param: selectedParameters) {
             if (parameter.equals(param.getParameterName())) {
                 param.setHP(hP);
-                found = true;
                 dirtyParameter = true;
-                break;
+                return;
             }
         }
-        if (!found) {
-            throw new MathIllegalArgumentException(LocalizedFormats.UNKNOWN_PARAMETER,
-                                                   parameter);
-        }
-    }
 
-    /** Set the steps in order to compute by finite difference the jacobian
-     *  matrix with respect to main state.
-     * <p>
-     * Needed if and only if the main set is a {@link FirstOrderDifferentialEquations}.
-     * </p>
-     * <p>
-     * Zero values for those steps don't enable to compute the main state jacobian matrix.
-     * </p>
-     * @param hY step used for finite difference computation with respect to state vector
-     * @exception IllegalArgumentException if the hY has not the dimension of the main state
-     * given by {@link ExpandableFirstOrderDifferentialEquations#getMainSetDimension()}
-     */
-    public void setMainStateSteps(final double[] hY) {
+        throw new MathIllegalArgumentException(LocalizedFormats.UNKNOWN_PARAMETER, parameter);
 
-        if (fode != null) {
-            // Check dimension
-            checkDimension(stateDim, hY);
-            this.hY = hY.clone();
-            dirtyMainState = true;           
-        }
     }
 
-    /** Set the initial value of the jacobian matrix with respect to state.
-     * @param dYdY0 initial jacobian matrix w.r.t. state
+    /** Set the initial value of the Jacobian matrix with respect to state.
+     * <p>
+     * If this method is not called, the initial value of the Jacobian
+     * matrix with respect to state is set to identity.
+     * </p>
+     * @param dYdY0 initial Jacobian matrix w.r.t. state
      * @exception IllegalArgumentException if matrix dimensions are incorrect
      */
-    public void setInitialMainStateJacobian(final double[][] dYdY0) {
+    public void setInitialMainStateJacobian(final double[][] dYdY0)
+        throws MathIllegalArgumentException {
 
         // Check dimensions
         checkDimension(stateDim, dYdY0);
         checkDimension(stateDim, dYdY0[0]);
 
         // store the matrix in row major order as a single dimension array
-        int index = 0;
+        int i = 0;
         for (final double[] row : dYdY0) {
-            System.arraycopy(row, 0, mainJacobianInARow, index, stateDim);
-            index += stateDim;
+            System.arraycopy(row, 0, matricesData, i, stateDim);
+            i += stateDim;
+        }
+
+        if (efode != null) {
+            efode.setSecondaryState(index, matricesData);
         }
-        // set initial additional state value in expandable fode
-        efode.setInitialAdditionalState(mainJacobianInARow, this);
+
     }
 
-    /** Set the initial value of the jacobian matrix with respect to one parameter.
-     * <p>The parameter must be {@link #selectParameters(String...) selected}.</p>
+    /** Set the initial value of a column of the Jacobian matrix with respect to one parameter.
+     * <p>
+     * If this method is not called for some parameter, the initial value of
+     * the column of the Jacobian matrix with respect to this parameter is set to zero.
+     * </p>
      * @param pName parameter name
-     * @param dYdP initial jacobian matrix w.r.t. the parameter
-     * @exception IllegalArgumentException if matrix dimensions are incorrect
+     * @param dYdP initial Jacobian column vector with respect to the parameter
+     * @exception MathIllegalArgumentException if a parameter is not supported
      */
-    public void setInitialParameterJacobian(final String pName, final double[] dYdP) {
+    public void setInitialParameterJacobian(final String pName, final double[] dYdP)
+        throws MathIllegalArgumentException {
 
         // Check dimensions
         checkDimension(stateDim, dYdP);
 
-        // store the matrix in a global single dimension array
-        boolean found = false;
-        int index = 0;
+        // store the column in a global single dimension array
+        int i = stateDim * stateDim;
         for (ParameterConfiguration param: selectedParameters) {
             if (pName.equals(param.getParameterName())) {
-                System.arraycopy(dYdP, 0, parameterJacobiansInARow, index, stateDim);
-                double[] p = new double[this.getDimension()];
-                index = stateDim * stateDim;
-                System.arraycopy(mainJacobianInARow, 0, p, 0, index);
-                System.arraycopy(parameterJacobiansInARow, 0, p, index, stateDim * paramDim);
-                // set initial additional state value in expandable fode
-                efode.setInitialAdditionalState(p, this);
-                found = true;
-                break;
+                System.arraycopy(dYdP, 0, matricesData, i, stateDim);
+                if (efode != null) {
+                    efode.setSecondaryState(index, matricesData);
+                }
+                return;
             }
-            index += stateDim;
-        }
-        if (! found) {
-            throw new MathIllegalArgumentException(LocalizedFormats.UNKNOWN_PARAMETER,
-                                                   pName);
+            i += stateDim;
         }
-    }
 
-    /** Set the default initial value of the jacobian matrix with respect to state.
-     * <p>dYdY0 is set to the identity matrix.</p>
-     */
-    public void setInitialMainStateJacobian() {
-        final double[][] dYdY0 = new double[stateDim][stateDim];
-        for (int i = 0; i < stateDim; ++i) {
-            dYdY0[i][i] = 1.0;
-        }
-        setInitialMainStateJacobian(dYdY0);
-    }
-
-    /** Set the default initial value of the jacobian matrix with respect to one parameter.
-     * <p>The parameter must be {@link #selectParameters(String...) selected}.</p>
-     * <p>dYdP is set to the null matrix.</p>
-     * @param pName parameter name
-     */
-    public void setInitialParameterJacobian(final String pName) {
-        setInitialParameterJacobian(pName, new double[stateDim]);
-    }
-
-    /** Set the default initial values of jacobian matrices with respect to all parameters.
-     */
-    public void setInitialParameterJacobians() {
-        for (ParameterConfiguration param: selectedParameters) {
-            setInitialParameterJacobian(param.getParameterName());
-        }
-    }
+        throw new MathIllegalArgumentException(LocalizedFormats.UNKNOWN_PARAMETER, pName);
 
-    /** Set default initial values for jacobian matrices.
-     * <p>dYdY0 is set to the identity matrix and all dYdP are set to zero.</p>
-     */
-    public void setInitialJacobians() {
-        setInitialMainStateJacobian();
-        setInitialParameterJacobians();
     }
 
-    /** Get the current value of the jacobian matrix with respect to state.
-     * @param dYdY0 current jacobian matrix with respect to state.
+    /** Get the current value of the Jacobian matrix with respect to state.
+     * @param dYdY0 current Jacobian matrix with respect to state.
      */
     public void getCurrentMainSetJacobian(final double[][] dYdY0) {
 
         // get current state for this set of equations from the expandable fode
-        double[] p = efode.getCurrentAdditionalState(this);
+        double[] p = efode.getSecondaryState(index);
 
         int index = 0;
         for (int i = 0; i < stateDim; i++) {
@@ -343,14 +297,14 @@ public class JacobianMatrices implements
 
     }
 
-    /** Get the current value of the jacobian matrix with respect to one parameter.
-     * @param pName name of the parameter for the computed jacobian matrix 
-     * @param dYdP current jacobian matrix with respect to the named parameter
+    /** Get the current value of the Jacobian matrix with respect to one parameter.
+     * @param pName name of the parameter for the computed Jacobian matrix 
+     * @param dYdP current Jacobian matrix with respect to the named parameter
      */
     public void getCurrentParameterJacobian(String pName, final double[] dYdP) {
 
         // get current state for this set of equations from the expandable fode
-        double[] p = efode.getCurrentAdditionalState(this);
+        double[] p = efode.getSecondaryState(index);
 
         int index = stateDim * stateDim;
         for (ParameterConfiguration param: selectedParameters) {
@@ -363,108 +317,151 @@ public class JacobianMatrices implements
 
     }
 
-    /** {@inheritDoc} */
-    public int getDimension() {
-        return stateDim * (stateDim + paramDim);
+    /** Check array dimensions.
+     * @param expected expected dimension
+     * @param array (may be null if expected is 0)
+     * @throws DimensionMismatchException if the array dimension does not match the expected one
+     */
+    private void checkDimension(final int expected, final Object array)
+        throws DimensionMismatchException {
+        int arrayDimension = (array == null) ? 0 : Array.getLength(array);
+        if (arrayDimension != expected) {
+            throw new DimensionMismatchException(arrayDimension, expected);
+        }
     }
 
-    /** {@inheritDoc} */
-    public void computeDerivatives(final double t, final double[] y, final double[] yDot,
-                                   final double[] z, final double[] zDot) {
+    /** Local implementation of secondary equations.
+     * <p>
+     * This class is an inner class to ensure proper scheduling of calls
+     * by forcing the use of {@link JacobianMatrices#registerVariationalEquations(ExpandableStatefulODE)}.
+     * </p>
+     */
+    private class JacobiansSecondaryEquations implements SecondaryEquations {
 
-        // Lazy initialization
-        if (dirtyMainState) {
-            jode = new MainStateJacobianWrapper(fode, hY);
-            dirtyMainState = false;
-        }
-        if (dirtyParameter && (paramDim != 0)) {
-            pjp.add(new ParameterJacobianWrapper(jode, pode, selectedParameters));
-            dirtyParameter = false;
+        /** {@inheritDoc} */
+        public int getDimension() {
+            return stateDim * (stateDim + paramDim);
         }
 
-        // variational equations:
-        // from d[dy/dt]/dy0 and d[dy/dt]/dp to d[dy/dy0]/dt and d[dy/dp]/dt
+        /** {@inheritDoc} */
+        public void computeDerivatives(final double t, final double[] y, final double[] yDot,
+                                       final double[] z, final double[] zDot) {
 
-        // compute jacobian matrix with respect to main state
-        double[][] dFdY = new double[stateDim][stateDim];
-        jode.computeMainStateJacobian(t, y, yDot, dFdY);
+            // Lazy initialization
+            if (dirtyParameter && (paramDim != 0)) {
+                jacobianProviders.add(new ParameterJacobianWrapper(jode, pode, selectedParameters));
+                dirtyParameter = false;
+            }
 
-        // Dispatch jacobian matrix in the compound additional state vector
-        for (int i = 0; i < stateDim; ++i) {
-            final double[] dFdYi = dFdY[i];
-            for (int j = 0; j < stateDim; ++j) {
-                double s = 0;
-                final int startIndex = j;
-                int zIndex = startIndex;
-                for (int l = 0; l < stateDim; ++l) {
-                    s += dFdYi[l] * z[zIndex];
-                    zIndex += stateDim;
+            // variational equations:
+            // from d[dy/dt]/dy0 and d[dy/dt]/dp to d[dy/dy0]/dt and d[dy/dp]/dt
+
+            // compute Jacobian matrix with respect to primary state
+            double[][] dFdY = new double[stateDim][stateDim];
+            jode.computeMainStateJacobian(t, y, yDot, dFdY);
+
+            // Dispatch Jacobian matrix in the compound secondary state vector
+            for (int i = 0; i < stateDim; ++i) {
+                final double[] dFdYi = dFdY[i];
+                for (int j = 0; j < stateDim; ++j) {
+                    double s = 0;
+                    final int startIndex = j;
+                    int zIndex = startIndex;
+                    for (int l = 0; l < stateDim; ++l) {
+                        s += dFdYi[l] * z[zIndex];
+                        zIndex += stateDim;
+                    }
+                    zDot[startIndex + i * stateDim] = s;
                 }
-                zDot[startIndex + i * stateDim] = s;
             }
-        }
 
-        if (paramDim != 0) {
-            // compute jacobian matrices with respect to parameters
-            double[] dFdP = new double[stateDim];
-            int startIndex = stateDim * stateDim;
-            for (ParameterConfiguration param: selectedParameters) {
-                boolean found = false;
-                for (ParameterJacobianProvider provider: pjp) {
-                    if (provider.isSupported(param.getParameterName())) {
-                        try {
-                            provider.computeParameterJacobian(t, y, yDot, param.getParameterName(), dFdP);
-                            for (int i = 0; i < stateDim; ++i) {
-                                final double[] dFdYi = dFdY[i];
-                                int zIndex = startIndex;
-                                double s = dFdP[i];
-                                for (int l = 0; l < stateDim; ++l) {
-                                    s += dFdYi[l] * z[zIndex];
-                                    zIndex++;
+            if (paramDim != 0) {
+                // compute Jacobian matrices with respect to parameters
+                double[] dFdP = new double[stateDim];
+                int startIndex = stateDim * stateDim;
+                for (ParameterConfiguration param: selectedParameters) {
+                    boolean found = false;
+                    for (ParameterJacobianProvider provider: jacobianProviders) {
+                        if (provider.isSupported(param.getParameterName())) {
+                            try {
+                                provider.computeParameterJacobian(t, y, yDot, param.getParameterName(), dFdP);
+                                for (int i = 0; i < stateDim; ++i) {
+                                    final double[] dFdYi = dFdY[i];
+                                    int zIndex = startIndex;
+                                    double s = dFdP[i];
+                                    for (int l = 0; l < stateDim; ++l) {
+                                        s += dFdYi[l] * z[zIndex];
+                                        zIndex++;
+                                    }
+                                    zDot[startIndex + i] = s;
                                 }
-                                zDot[startIndex + i] = s;
+                                startIndex += stateDim;
+                                found = true;
+                                break;
+                            } catch (IllegalArgumentException iae) {
                             }
-                            startIndex += stateDim;
-                            found = true;
-                            break;
-                        } catch (IllegalArgumentException iae) {
                         }
                     }
-                }
-                if (! found) {
-                    throw new MathIllegalArgumentException(LocalizedFormats.UNKNOWN_PARAMETER,
-                                                           param);
+                    if (! found) {
+                        throw new MathIllegalArgumentException(LocalizedFormats.UNKNOWN_PARAMETER,
+                                                               param);
+                    }
                 }
             }
-        }
 
+        }
     }
 
-    /** Check compatibility between the main set in the expandable ode and an ordinary ode.
-     * @param expended expandable ode containing a main set
-     * @param ode single ode to check 
-     * @throws MathIllegalArgumentException if single ode doesn't match the main ode set in the extended ode
+    /** Wrapper class to compute jacobian matrices by finite differences for ODE
+     *  which do not compute them by themselves.
      */
-    private void checkCompatibility(final ExpandableFirstOrderDifferentialEquations extended,
-                                    final FirstOrderDifferentialEquations ode)
-        throws MathIllegalArgumentException {
+    private static class MainStateJacobianWrapper implements MainStateJacobianProvider {
 
-        if (!(ode == extended.getMainSet())) {
-            throw new MathIllegalArgumentException(LocalizedFormats.UNMATCHED_ODE_IN_EXTENDED_SET);
+        /** Raw ODE without jacobians computation skill to be wrapped into a MainStateJacobianProvider. */
+        private final FirstOrderDifferentialEquations ode;
+
+        /** Steps for finite difference computation of the jacobian df/dy w.r.t. state. */
+        private final double[] hY;
+
+        /** Wrap a {@link FirstOrderDifferentialEquations} into a {@link MainStateJacobianProvider}.
+         * @param ode original ODE problem, without jacobians computation skill
+         * @param hY step sizes to compute the jacobian df/dy
+         * @see JacobianMatrices#setMainStateSteps(double[])
+         */
+        public MainStateJacobianWrapper(final FirstOrderDifferentialEquations ode,
+                                        final double[] hY) {
+            this.ode = ode;
+            this.hY = hY.clone();
         }
-    }
 
-    /** Check array dimensions.
-     * @param expected expected dimension
-     * @param array (may be null if expected is 0)
-     * @throws DimensionMismatchException if the array dimension does not match the expected one
-     */
-    private void checkDimension(final int expected, final Object array)
-        throws DimensionMismatchException {
-        int arrayDimension = (array == null) ? 0 : Array.getLength(array);
-        if (arrayDimension != expected) {
-            throw new DimensionMismatchException(arrayDimension, expected);
+        /** {@inheritDoc} */
+        public int getDimension() {
+            return ode.getDimension();
+        }
+
+        /** {@inheritDoc} */
+        public void computeDerivatives(double t, double[] y, double[] yDot) {
+            ode.computeDerivatives(t, y, yDot);
+        }
+
+        /** {@inheritDoc} */
+        public void computeMainStateJacobian(double t, double[] y, double[] yDot,
+                                             double[][] dFdY) {
+
+            final int n = ode.getDimension();
+            final double[] tmpDot = new double[n];
+
+            for (int j = 0; j < n; ++j) {
+                final double savedYj = y[j];
+                y[j] += hY[j];
+                ode.computeDerivatives(t, y, tmpDot);
+                for (int i = 0; i < n; ++i) {
+                    dFdY[i][j] = (tmpDot[i] - yDot[i]) / hY[j];
+                }
+                y[j] = savedYj;
+            }
         }
+
     }
 
 }

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterJacobianProvider.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterJacobianProvider.java?rev=1176734&r1=1176733&r2=1176734&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterJacobianProvider.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterJacobianProvider.java Wed Sep 28 05:56:42 2011
@@ -18,7 +18,7 @@ package org.apache.commons.math.ode;
 
 import org.apache.commons.math.exception.MathIllegalArgumentException;
 
-/** Interface to compute exactly jacobian matrix for some parameter
+/** Interface to compute exactly Jacobian matrix for some parameter
  *  when computing {@link JacobianMatrices partial derivatives equations}.
  * 
  * @version $Id$
@@ -26,18 +26,19 @@ import org.apache.commons.math.exception
  */
 public interface ParameterJacobianProvider extends Parameterizable {
 
-    /** Compute the jacobian matrix of ODE with respect to one parameter.
+    /** Compute the Jacobian matrix of ODE with respect to one parameter.
      * <p>The parameter must be one given by {@link #getParametersNames()}.</p>
      * @param t current value of the independent <I>time</I> variable
      * @param y array containing the current value of the main state vector
      * @param yDot array containing the current value of the time derivative
      * of the main state vector
      * @param paramName name of the parameter to consider
-     * @param dFdP placeholder array where to put the jacobian matrix of the
+     * @param dFdP placeholder array where to put the Jacobian matrix of the
      * ODE with respect to the parameter
      * @throws MathIllegalArgumentException if the parameter is not supported
      */
     void computeParameterJacobian(double t, double[] y, double[] yDot,
                                   String paramName, double[] dFdP)
         throws MathIllegalArgumentException;
+
 }

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterJacobianWrapper.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterJacobianWrapper.java?rev=1176734&r1=1176733&r2=1176734&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterJacobianWrapper.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterJacobianWrapper.java Wed Sep 28 05:56:42 2011
@@ -20,7 +20,7 @@ import java.util.Collection;
 import java.util.HashMap;
 import java.util.Map;
 
-/** Wrapper class to compute jacobian matrices by finite differences for ODE
+/** Wrapper class to compute Jacobian matrices by finite differences for ODE
  *  which do not compute them by themselves.
  *  
  * @version $Id$
@@ -31,21 +31,21 @@ class ParameterJacobianWrapper implement
     /** Main ODE set. */
     private final FirstOrderDifferentialEquations fode;
 
-    /** Raw ODE without jacobian computation skill to be wrapped into a ParameterJacobianProvider. */
+    /** Raw ODE without Jacobian computation skill to be wrapped into a ParameterJacobianProvider. */
     private final ParameterizedODE pode;
 
-    /** Steps for finite difference computation of the jacobian df/dp w.r.t. parameters. */
+    /** Steps for finite difference computation of the Jacobian df/dp w.r.t. parameters. */
     private final Map<String, Double> hParam;
 
     /** Wrap a {@link ParameterizedODE} into a {@link ParameterJacobianProvider}.
      * @param fode main first order differential equations set
-     * @param pode additional problem, without parametre jacobian computation skill
-     * @param paramsAndSteps parameters and steps to compute the jacobians df/dp
+     * @param pode secondary problem, without parameter Jacobian computation skill
+     * @param paramsAndSteps parameters and steps to compute the Jacobians df/dp
      * @see JacobianMatrices#setParameterStep(String, double)
      */
     public ParameterJacobianWrapper(final FirstOrderDifferentialEquations fode,
                                     final ParameterizedODE pode,
-                                    final Collection<ParameterConfiguration> paramsAndSteps) {
+                                    final ParameterConfiguration[] paramsAndSteps) {
         this.fode = fode;
         this.pode = pode;
         this.hParam = new HashMap<String, Double>();

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterizedODE.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterizedODE.java?rev=1176734&r1=1176733&r2=1176734&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterizedODE.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ParameterizedODE.java Wed Sep 28 05:56:42 2011
@@ -16,7 +16,7 @@
  */
 package org.apache.commons.math.ode;
 
-/** Interface to compute by finite difference jacobian matrix for some parameter
+/** Interface to compute by finite difference Jacobian matrix for some parameter
  *  when computing {@link JacobianMatrices partial derivatives equations}.
  * 
  * @version $Id$

Copied: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/SecondaryEquations.java (from r1175683, commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AdditionalEquations.java)
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/SecondaryEquations.java?p2=commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/SecondaryEquations.java&p1=commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AdditionalEquations.java&r1=1175683&r2=1176734&rev=1176734&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AdditionalEquations.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/SecondaryEquations.java Wed Sep 28 05:56:42 2011
@@ -18,38 +18,39 @@
 package org.apache.commons.math.ode;
 
 /**
- * This interface allows users to add their own differential equations to a main
+ * This interface allows users to add secondary differential equations to a primary
  * set of differential equations.
  * <p>
  * In some cases users may need to integrate some problem-specific equations along
- * with a main set of differential equations. One example is optimal control where
+ * with a primary set of differential equations. One example is optimal control where
  * adjoined parameters linked to the minimized hamiltonian must be integrated.
  * </p>
  * <p>
- * This interface allows users to add such equations to a main set of {@link
+ * This interface allows users to add such equations to a primary set of {@link
  * FirstOrderDifferentialEquations first order differential equations}
  * thanks to the {@link
- * ExpandableFirstOrderDifferentialEquations#addAdditionalEquations(AdditionalEquations)}
+ * ExpandableStatefulODE#addSecondaryEquations(SecondaryEquations)}
  * method.
  * </p>
- * @see ExpandableFirstOrderDifferentialEquations
+ * @see ExpandableStatefulODE
  * @version $Id$
  * @since 3.0
  */
-public interface AdditionalEquations {
+public interface SecondaryEquations {
 
-    /** Get the dimension of the additional state parameters.
-     * @return dimension of the additional state parameters
+    /** Get the dimension of the secondary state parameters.
+     * @return dimension of the secondary state parameters
      */
     int getDimension();
 
-    /** Compute the derivatives related to the additional state parameters.
+    /** Compute the derivatives related to the secondary state parameters.
      * @param t current value of the independent <I>time</I> variable
-     * @param y array containing the current value of the main state vector
-     * @param yDot array containing the derivative of the main state vector
-     * @param z array containing the current value of the additional state vector
-     * @param zDot placeholder array where to put the derivative of the additional state vector
+     * @param primary array containing the current value of the primary state vector
+     * @param primaryDot array containing the derivative of the primary state vector
+     * @param secondary array containing the current value of the secondary state vector
+     * @param secondaryDot placeholder array where to put the derivative of the secondary state vector
      */
-    void computeDerivatives(double t, double[] y, double[] yDot, double[] z, double[] zDot);
+    void computeDerivatives(double t, double[] primary, double[] primaryDot,
+                            double[] secondary, double[] secondaryDot);
 
 }

Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/SecondaryEquations.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/SecondaryEquations.java
------------------------------------------------------------------------------
    svn:keywords = "Author Date Id Revision"

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java?rev=1176734&r1=1176733&r2=1176734&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java Wed Sep 28 05:56:42 2011
@@ -20,7 +20,7 @@ package org.apache.commons.math.ode.nons
 import org.apache.commons.math.exception.MathIllegalArgumentException;
 import org.apache.commons.math.exception.MathIllegalStateException;
 import org.apache.commons.math.linear.Array2DRowRealMatrix;
-import org.apache.commons.math.ode.ExpandableFirstOrderDifferentialEquations;
+import org.apache.commons.math.ode.ExpandableStatefulODE;
 import org.apache.commons.math.ode.sampling.NordsieckStepInterpolator;
 import org.apache.commons.math.ode.sampling.StepHandler;
 import org.apache.commons.math.util.FastMath;
@@ -189,31 +189,23 @@ public class AdamsBashforthIntegrator ex
 
     /** {@inheritDoc} */
     @Override
-    public double integrate(final ExpandableFirstOrderDifferentialEquations equations,
-                            final double t0, final double[] z0,
-                            final double t, final double[] z)
+    public void integrate(final ExpandableStatefulODE equations, final double t)
         throws MathIllegalStateException, MathIllegalArgumentException {
 
-        sanityChecks(equations, t0, z0, t, z);
+        sanityChecks(equations, t);
         setEquations(equations);
         resetEvaluations();
-        final boolean forward = t > t0;
+        final boolean forward = t > equations.getTime();
 
         // initialize working arrays
-        final int totalDim = equations.getDimension();
-        final int mainDim  = equations.getMainSetDimension();
-        final double[] y0  = new double[totalDim];
-        final double[] y   = new double[totalDim];
-        System.arraycopy(z0, 0, y0, 0, mainDim);
-        System.arraycopy(equations.getCurrentAdditionalStates(), 0, y0, mainDim, totalDim - mainDim);
-        if (y != y0) {
-            System.arraycopy(y0, 0, y, 0, totalDim);
-        }
-        final double[] yDot = new double[totalDim];
+        final double[] y0   = equations.getCompleteState();
+        final double[] y    = y0.clone();
+        final double[] yDot = new double[y.length];
 
         // set up an interpolator sharing the integrator arrays
         final NordsieckStepInterpolator interpolator = new NordsieckStepInterpolator();
-        interpolator.reinitialize(y, forward);
+        interpolator.reinitialize(y, forward,
+                                  equations.getPrimaryMapper(), equations.getSecondaryMappers());
 
         // set up integration control objects
         for (StepHandler handler : stepHandlers) {
@@ -222,7 +214,7 @@ public class AdamsBashforthIntegrator ex
         setStateInitialized(false);
 
         // compute the initial Nordsieck vector using the configured starter integrator
-        start(t0, y, t);
+        start(equations.getTime(), y, t);
         interpolator.reinitialize(stepStart, stepSize, scaled, nordsieck);
         interpolator.storeTime(stepStart);
         final int lastRow = nordsieck.getRowDimension() - 1;
@@ -317,13 +309,11 @@ public class AdamsBashforthIntegrator ex
 
         } while (!isLastStep);
 
-        // dispatch result between main and additional states
-        System.arraycopy(y, 0, z, 0, z.length);
-        equations.setCurrentAdditionalState(y);
+        // dispatch results
+        equations.setTime(stepStart);
+        equations.setCompleteState(y);
 
-        final double stopTime = stepStart;
         resetInternalState();
-        return stopTime;
 
     }
 

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsIntegrator.java?rev=1176734&r1=1176733&r2=1176734&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsIntegrator.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsIntegrator.java Wed Sep 28 05:56:42 2011
@@ -20,7 +20,7 @@ package org.apache.commons.math.ode.nons
 import org.apache.commons.math.exception.MathIllegalArgumentException;
 import org.apache.commons.math.exception.MathIllegalStateException;
 import org.apache.commons.math.linear.Array2DRowRealMatrix;
-import org.apache.commons.math.ode.ExpandableFirstOrderDifferentialEquations;
+import org.apache.commons.math.ode.ExpandableStatefulODE;
 import org.apache.commons.math.ode.MultistepIntegrator;
 
 
@@ -86,9 +86,7 @@ public abstract class AdamsIntegrator ex
 
     /** {@inheritDoc} */
     @Override
-    public abstract double integrate(final ExpandableFirstOrderDifferentialEquations equations,
-                                     final double t0, final double[] y0,
-                                     final double t, final double[] y)
+    public abstract void integrate(final ExpandableStatefulODE equations, final double t)
         throws MathIllegalStateException, MathIllegalArgumentException;
 
     /** {@inheritDoc} */

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java?rev=1176734&r1=1176733&r2=1176734&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java Wed Sep 28 05:56:42 2011
@@ -23,7 +23,7 @@ import org.apache.commons.math.exception
 import org.apache.commons.math.exception.MathIllegalStateException;
 import org.apache.commons.math.linear.Array2DRowRealMatrix;
 import org.apache.commons.math.linear.RealMatrixPreservingVisitor;
-import org.apache.commons.math.ode.ExpandableFirstOrderDifferentialEquations;
+import org.apache.commons.math.ode.ExpandableStatefulODE;
 import org.apache.commons.math.ode.sampling.NordsieckStepInterpolator;
 import org.apache.commons.math.ode.sampling.StepHandler;
 import org.apache.commons.math.util.FastMath;
@@ -206,34 +206,26 @@ public class AdamsMoultonIntegrator exte
 
     /** {@inheritDoc} */
     @Override
-    public double integrate(final ExpandableFirstOrderDifferentialEquations equations,
-                            final double t0, final double[] z0,
-                            final double t, final double[] z)
+    public void integrate(final ExpandableStatefulODE equations,final double t)
         throws MathIllegalStateException, MathIllegalArgumentException {
 
-        sanityChecks(equations, t0, z0, t, z);
+        sanityChecks(equations, t);
         setEquations(equations);
         resetEvaluations();
-        final boolean forward = t > t0;
+        final boolean forward = t > equations.getTime();
 
         // initialize working arrays
-        final int totalDim = equations.getDimension();
-        final int mainDim  = equations.getMainSetDimension();
-        final double[] y0  = new double[totalDim];
-        final double[] y   = new double[totalDim];
-        System.arraycopy(z0, 0, y0, 0, mainDim);
-        System.arraycopy(equations.getCurrentAdditionalStates(), 0, y0, mainDim, totalDim - mainDim);
-        if (y != y0) {
-            System.arraycopy(y0, 0, y, 0, totalDim);
-        }
-        final double[] yDot = new double[totalDim];
-        final double[] yTmp = new double[totalDim];
-        final double[] predictedScaled = new double[totalDim];
+        final double[] y0   = equations.getCompleteState();
+        final double[] y    = y0.clone();
+        final double[] yDot = new double[y.length];
+        final double[] yTmp = new double[y.length];
+        final double[] predictedScaled = new double[y.length];
         Array2DRowRealMatrix nordsieckTmp = null;
 
         // set up two interpolators sharing the integrator arrays
         final NordsieckStepInterpolator interpolator = new NordsieckStepInterpolator();
-        interpolator.reinitialize(y, forward);
+        interpolator.reinitialize(y, forward,
+                                  equations.getPrimaryMapper(), equations.getSecondaryMappers());
 
         // set up integration control objects
         for (StepHandler handler : stepHandlers) {
@@ -242,7 +234,7 @@ public class AdamsMoultonIntegrator exte
         setStateInitialized(false);
 
         // compute the initial Nordsieck vector using the configured starter integrator
-        start(t0, y, t);
+        start(equations.getTime(), y, t);
         interpolator.reinitialize(stepStart, stepSize, scaled, nordsieck);
         interpolator.storeTime(stepStart);
 
@@ -295,7 +287,7 @@ public class AdamsMoultonIntegrator exte
             updateHighOrderDerivativesPhase2(predictedScaled, correctedScaled, nordsieckTmp);
 
             // discrete events handling
-            System.arraycopy(yTmp, 0, y, 0, totalDim);
+            System.arraycopy(yTmp, 0, y, 0, y.length);
             interpolator.reinitialize(stepEnd, stepSize, correctedScaled, nordsieckTmp);
             interpolator.storeTime(stepStart);
             interpolator.shift();
@@ -335,14 +327,11 @@ public class AdamsMoultonIntegrator exte
 
         } while (!isLastStep);
 
-        // dispatch result between main and additional states
-        System.arraycopy(y, 0, z, 0, z.length);
-        equations.setCurrentAdditionalState(y);
-
-        final double stopTime  = stepStart;
-        stepStart = Double.NaN;
-        stepSize  = Double.NaN;
-        return stopTime;
+        // dispatch results
+        equations.setTime(stepStart);
+        equations.setCompleteState(y);
+
+        resetInternalState();
 
     }
 

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdaptiveStepsizeIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdaptiveStepsizeIntegrator.java?rev=1176734&r1=1176733&r2=1176734&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdaptiveStepsizeIntegrator.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdaptiveStepsizeIntegrator.java Wed Sep 28 05:56:42 2011
@@ -23,7 +23,7 @@ import org.apache.commons.math.exception
 import org.apache.commons.math.exception.NumberIsTooSmallException;
 import org.apache.commons.math.exception.util.LocalizedFormats;
 import org.apache.commons.math.ode.AbstractIntegrator;
-import org.apache.commons.math.ode.ExpandableFirstOrderDifferentialEquations;
+import org.apache.commons.math.ode.ExpandableStatefulODE;
 import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
 import org.apache.commons.math.util.FastMath;
 
@@ -43,9 +43,9 @@ import org.apache.commons.math.util.Fast
  * relTol which will be used for all components.
  * </p>
  * <p>
- * If the Ordinary Differential Equations is an {@link ExpandableFirstOrderDifferentialEquations
+ * If the Ordinary Differential Equations is an {@link ExpandableStatefulODE
  * extended ODE} rather than a {@link FirstOrderDifferentialEquations basic ODE}, then
- * <em>only</em> the {@link ExpandableFirstOrderDifferentialEquations#getMainSet() main part}
+ * <em>only</em> the {@link ExpandableStatefulODE#getMainSet() main part}
  * of the state vector is used for stepsize control, not the complete state vector.
  * </p>
  *
@@ -213,24 +213,14 @@ public abstract class AdaptiveStepsizeIn
     }
   }
 
-  /** Perform some sanity checks on the integration parameters.
-   * @param equations differential equations set
-   * @param t0 start time
-   * @param y0 state vector at t0
-   * @param t target time for the integration
-   * @param y placeholder where to put the state vector
-   * @exception DimensionMismatchException if some inconsistency is detected
-   * @exception NumberIsTooSmallException if integration span is too small
-   */
+  /** {@inheritDoc} */
   @Override
-  protected void sanityChecks(final ExpandableFirstOrderDifferentialEquations equations,
-                              final double t0, final double[] y0,
-                              final double t, final double[] y)
+  protected void sanityChecks(final ExpandableStatefulODE equations, final double t)
       throws DimensionMismatchException, NumberIsTooSmallException {
 
-      super.sanityChecks(equations, t0, y0, t, y);
+      super.sanityChecks(equations, t);
 
-      mainSetDimension = equations.getMainSetDimension();
+      mainSetDimension = equations.getPrimaryMapper().getDimension();
 
       if ((vecAbsoluteTolerance != null) && (vecAbsoluteTolerance.length != mainSetDimension)) {
           throw new DimensionMismatchException(mainSetDimension, vecAbsoluteTolerance.length);
@@ -349,9 +339,7 @@ public abstract class AdaptiveStepsizeIn
   }
 
   /** {@inheritDoc} */
-  public abstract double integrate (ExpandableFirstOrderDifferentialEquations equations,
-                                    double t0, double[] y0,
-                                    double t, double[] y)
+  public abstract void integrate (ExpandableStatefulODE equations, double t)
     throws MathIllegalStateException, MathIllegalArgumentException;
 
   /** {@inheritDoc} */

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince54StepInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince54StepInterpolator.java?rev=1176734&r1=1176733&r2=1176734&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince54StepInterpolator.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince54StepInterpolator.java Wed Sep 28 05:56:42 2011
@@ -18,6 +18,7 @@
 package org.apache.commons.math.ode.nonstiff;
 
 import org.apache.commons.math.ode.AbstractIntegrator;
+import org.apache.commons.math.ode.EquationsMapper;
 import org.apache.commons.math.ode.sampling.StepInterpolator;
 
 /**
@@ -145,8 +146,10 @@ class DormandPrince54StepInterpolator
   /** {@inheritDoc} */
   @Override
   public void reinitialize(final AbstractIntegrator integrator,
-                           final double[] y, final double[][] yDotK, final boolean forward) {
-    super.reinitialize(integrator, y, yDotK, forward);
+                           final double[] y, final double[][] yDotK, final boolean forward,
+                           final EquationsMapper primaryMapper,
+                           final EquationsMapper[] secondaryMappers) {
+    super.reinitialize(integrator, y, yDotK, forward, primaryMapper, secondaryMappers);
     v1 = null;
     v2 = null;
     v3 = null;



Mime
View raw message