Author: luc
Date: Mon Apr 8 14:41:32 2013
New Revision: 1465654
URL: http://svn.apache.org/r1465654
Log:
Fixed inconsistency preventing use of secondary states in ODE events.
JIRA: MATH-965
Modified:
commons/proper/math/trunk/src/changes/changes.xml
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/AbstractIntegrator.java
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/events/EventState.java
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/ode/events/EventStateTest.java
Modified: commons/proper/math/trunk/src/changes/changes.xml
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/changes/changes.xml?rev=1465654&r1=1465653&r2=1465654&view=diff
==============================================================================
--- commons/proper/math/trunk/src/changes/changes.xml (original)
+++ commons/proper/math/trunk/src/changes/changes.xml Mon Apr 8 14:41:32 2013
@@ -51,6 +51,10 @@ If the output is not quite correct, chec
</properties>
<body>
<release version="x.y" date="TBD" description="TBD">
+ <action dev="luc" type="fix" issue="MATH-965" >
+ Fixed inconsistent dimensions preventing use of secondary states
+ in ODE events.
+ </action>
</release>
<release version="3.2" date="2013-04-06" description="
This is a minor release: It combines bug fixes and new features.
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/AbstractIntegrator.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/AbstractIntegrator.java?rev=1465654&r1=1465653&r2=1465654&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/AbstractIntegrator.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/AbstractIntegrator.java
Mon Apr 8 14:41:32 2013
@@ -188,6 +188,7 @@ public abstract class AbstractIntegrator
evaluations.resetCount();
for (final EventState state : eventsStates) {
+ state.setExpandable(expandable);
state.getEventHandler().init(t0, y0, t);
}
@@ -356,7 +357,6 @@ public abstract class AbstractIntegrator
// get state at event time
interpolator.setInterpolatedTime(eventT);
- final double[] eventYPrimary =
interpolator.getInterpolatedState().clone();
final double[] eventYComplete = new double[y.length];
expandable.getPrimaryMapper().insertEquationData(interpolator.getInterpolatedState(),
eventYComplete);
@@ -368,7 +368,7 @@ public abstract class AbstractIntegrator
// advance all event states to current time
for (final EventState state : eventsStates) {
- state.stepAccepted(eventT, eventYPrimary);
+ state.stepAccepted(eventT, eventYComplete);
isLastStep = isLastStep || state.stop();
}
@@ -412,7 +412,14 @@ public abstract class AbstractIntegrator
// last part of the step, after the last event
interpolator.setInterpolatedTime(currentT);
- final double[] currentY = interpolator.getInterpolatedState();
+ final double[] currentY = new double[y.length];
+
expandable.getPrimaryMapper().insertEquationData(interpolator.getInterpolatedState(),
+ currentY);
+ int index = 0;
+ for (EquationsMapper secondary : expandable.getSecondaryMappers())
{
+
secondary.insertEquationData(interpolator.getInterpolatedSecondaryState(index++),
+ currentY);
+ }
for (final EventState state : eventsStates) {
state.stepAccepted(currentT, currentY);
isLastStep = isLastStep || state.stop();
Modified:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/events/EventState.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/events/EventState.java?rev=1465654&r1=1465653&r2=1465654&view=diff
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/events/EventState.java
(original)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/events/EventState.java
Mon Apr 8 14:41:32 2013
@@ -25,6 +25,8 @@ import org.apache.commons.math3.analysis
import org.apache.commons.math3.analysis.solvers.UnivariateSolverUtils;
import org.apache.commons.math3.exception.MaxCountExceededException;
import org.apache.commons.math3.exception.NoBracketingException;
+import org.apache.commons.math3.ode.EquationsMapper;
+import org.apache.commons.math3.ode.ExpandableStatefulODE;
import org.apache.commons.math3.ode.sampling.StepInterpolator;
import org.apache.commons.math3.util.FastMath;
@@ -55,6 +57,9 @@ public class EventState {
/** Upper limit in the iteration count for event localization. */
private final int maxIterationCount;
+ /** Equation being integrated. */
+ private ExpandableStatefulODE expandable;
+
/** Time at the beginning of the step. */
private double t0;
@@ -107,6 +112,7 @@ public class EventState {
this.solver = solver;
// some dummy values ...
+ expandable = null;
t0 = Double.NaN;
g0 = Double.NaN;
g0Positive = true;
@@ -125,6 +131,13 @@ public class EventState {
return handler;
}
+ /** Set the equation.
+ * @param expandable equation being integrated
+ */
+ public void setExpandable(final ExpandableStatefulODE expandable) {
+ this.expandable = expandable;
+ }
+
/** Get the maximal time interval between events handler checks.
* @return maximal time interval between events handler checks
*/
@@ -156,7 +169,7 @@ public class EventState {
t0 = interpolator.getPreviousTime();
interpolator.setInterpolatedTime(t0);
- g0 = handler.g(t0, interpolator.getInterpolatedState());
+ g0 = handler.g(t0, getCompleteState(interpolator));
if (g0 == 0) {
// excerpt from MATH-421 issue:
// If an ODE solver is setup with an EventHandler that return STOP
@@ -175,12 +188,32 @@ public class EventState {
FastMath.abs(solver.getRelativeAccuracy() * t0));
final double tStart = t0 + 0.5 * epsilon;
interpolator.setInterpolatedTime(tStart);
- g0 = handler.g(tStart, interpolator.getInterpolatedState());
+ g0 = handler.g(tStart, getCompleteState(interpolator));
}
g0Positive = g0 >= 0;
}
+ /** Get the complete state (primary and secondary).
+ * @param interpolator interpolator to use
+ * @return complete state
+ */
+ private double[] getCompleteState(final StepInterpolator interpolator) {
+
+ final double[] complete = new double[expandable.getTotalDimension()];
+
+
expandable.getPrimaryMapper().insertEquationData(interpolator.getInterpolatedState(),
+ complete);
+ int index = 0;
+ for (EquationsMapper secondary : expandable.getSecondaryMappers()) {
+
secondary.insertEquationData(interpolator.getInterpolatedSecondaryState(index++),
+ complete);
+ }
+
+ return complete;
+
+ }
+
/** Evaluate the impact of the proposed step on the event handler.
* @param interpolator step interpolator for the proposed step
* @return true if the event handler triggers an event before
@@ -207,7 +240,7 @@ public class EventState {
public double value(final double t) throws
LocalMaxCountExceededException {
try {
interpolator.setInterpolatedTime(t);
- return handler.g(t,
interpolator.getInterpolatedState());
+ return handler.g(t, getCompleteState(interpolator));
} catch (MaxCountExceededException mcee) {
throw new LocalMaxCountExceededException(mcee);
}
@@ -221,7 +254,7 @@ public class EventState {
// evaluate handler value at the end of the substep
final double tb = t0 + (i + 1) * h;
interpolator.setInterpolatedTime(tb);
- final double gb = handler.g(tb,
interpolator.getInterpolatedState());
+ final double gb = handler.g(tb,
getCompleteState(interpolator));
// check events occurrence
if (g0Positive ^ (gb >= 0)) {
Modified:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/ode/events/EventStateTest.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/ode/events/EventStateTest.java?rev=1465654&r1=1465653&r2=1465654&view=diff
==============================================================================
---
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/ode/events/EventStateTest.java
(original)
+++
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/ode/events/EventStateTest.java
Mon Apr 8 14:41:32 2013
@@ -23,7 +23,9 @@ import org.apache.commons.math3.exceptio
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.ode.ExpandableStatefulODE;
import org.apache.commons.math3.ode.FirstOrderDifferentialEquations;
+import org.apache.commons.math3.ode.SecondaryEquations;
import org.apache.commons.math3.ode.nonstiff.DormandPrince853Integrator;
import org.apache.commons.math3.ode.sampling.AbstractStepInterpolator;
import org.apache.commons.math3.ode.sampling.DummyStepInterpolator;
@@ -56,6 +58,13 @@ public class EventStateTest {
EventState es = new EventState(closeEventsGenerator, 1.5 * gap,
tolerance, 100,
new BrentSolver(tolerance));
+ es.setExpandable(new ExpandableStatefulODE(new
FirstOrderDifferentialEquations() {
+ public int getDimension() {
+ return 0;
+ }
+ public void computeDerivatives(double t, double[] y, double[]
yDot) {
+ }
+ }));
AbstractStepInterpolator interpolator =
new DummyStepInterpolator(new double[0], new double[0], true);
@@ -145,5 +154,74 @@ public class EventStateTest {
}
+ // Jira: MATH-965
+ @Test
+ public void testIssue965()
+ throws DimensionMismatchException, NumberIsTooSmallException,
+ MaxCountExceededException, NoBracketingException {
+
+ ExpandableStatefulODE equation =
+ new ExpandableStatefulODE(new
FirstOrderDifferentialEquations() {
+
+ public int getDimension() {
+ return 1;
+ }
+
+ public void computeDerivatives(double t, double[] y, double[]
yDot) {
+ yDot[0] = 2.0;
+ }
+ });
+ equation.setTime(0.0);
+ equation.setPrimaryState(new double[1]);
+ equation.addSecondaryEquations(new SecondaryEquations() {
+
+ public int getDimension() {
+ return 1;
+ }
+
+ public void computeDerivatives(double t, double[] primary,
+ double[] primaryDot, double[]
secondary,
+ double[] secondaryDot) {
+ secondaryDot[0] = -3.0;
+ }
+ });
+ int index = equation.getSecondaryMappers()[0].getFirstIndex();
+
+ DormandPrince853Integrator integrator = new
DormandPrince853Integrator(0.001, 1000, 1.0e-14, 1.0e-14);
+ integrator.addEventHandler(new SecondaryStateEvent(index, -3.0), 0.1,
1.0e-9, 1000);
+ integrator.setInitialStepSize(3.0);
+
+ integrator.integrate(equation, 30.0);
+ Assert.assertEquals( 1.0, equation.getTime(), 1.0e-10);
+ Assert.assertEquals( 2.0, equation.getPrimaryState()[0], 1.0e-10);
+ Assert.assertEquals(-3.0, equation.getSecondaryState(0)[0], 1.0e-10);
+
+ }
+
+ private static class SecondaryStateEvent implements EventHandler {
+
+ private int index;
+ private final double target;
+
+ public SecondaryStateEvent(final int index, final double target) {
+ this.index = index;
+ this.target = target;
+ }
+
+ public void init(double t0, double[] y0, double t) {
+ }
+
+ public double g(double t, double[] y) {
+ return y[index] - target;
+ }
+
+ public Action eventOccurred(double t, double[] y, boolean increasing) {
+ return Action.STOP;
+ }
+
+ public void resetState(double t, double[] y) {
+ }
+
+ }
}
