Author: luc
Date: Tue Mar 19 14:09:58 2013
New Revision: 1458298
URL: http://svn.apache.org/r1458298
Log:
Added a way to trigger only increasing or decreasing events in ODE.
The method used is completely different from what was suggested in the
request report. It is fully compatible with previous versions and most
importantly has no side effects on users that do not need this feature,
as it is based on an upstream filtering of the event definition
function. The events detection logic is not modified and it fact it does
not even know events are filtered before they are detected.
JIRA: MATH-811
Added:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/events/EventFilter.java
(with props)
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/events/FilterType.java
(with props)
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/events/Transformer.java
(with props)
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/ode/events/EventFilterTest.java
(with props)
Modified:
commons/proper/math/trunk/src/changes/changes.xml
Modified: commons/proper/math/trunk/src/changes/changes.xml
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/changes/changes.xml?rev=1458298&r1=1458297&r2=1458298&view=diff
==============================================================================
--- commons/proper/math/trunk/src/changes/changes.xml (original)
+++ commons/proper/math/trunk/src/changes/changes.xml Tue Mar 19 14:09:58 2013
@@ -55,6 +55,9 @@ This is a minor release: It combines bug
Changes to existing features were made in a backwards-compatible
way such as to allow drop-in replacement of the v3.1[.1] JAR file.
">
+ <action dev="luc" type="add" issue="MATH-811" >
+ Added a way to trigger only increasing or decreasing events in ODE
integration.
+ </action>
<action dev="luc" type="fix" issue="MATH-950" >
Fixed missing update in ODE event handlers, when a RESET_STATE is
triggered.
</action>
Added:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/events/EventFilter.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/events/EventFilter.java?rev=1458298&view=auto
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/events/EventFilter.java
(added)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/events/EventFilter.java
Tue Mar 19 14:09:58 2013
@@ -0,0 +1,207 @@
+/*
+ * 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.math3.ode.events;
+
+import java.util.Arrays;
+
+import org.apache.commons.math3.ode.FirstOrderIntegrator;
+
+
+/** Wrapper used to detect only increasing or decreasing events.
+ *
+ * <p>General {@link EventHandler events} are defined implicitely
+ * by a {@link EventHandler#g(double, double[]) g function} crossing
+ * zero. This function needs to be continuous in the event neighborhood,
+ * and its sign must remain consistent between events. This implies that
+ * during an ODE integration, events triggered are alternately events
+ * for which the function increases from negative to positive values,
+ * and events for which the function decreases from positive to
+ * negative values.
+ * </p>
+ *
+ * <p>Sometimes, users are only interested in one type of event (say
+ * increasing events for example) and not in the other type. In these
+ * cases, looking precisely for all events location and triggering
+ * events that will later be ignored is a waste of computing time.</p>
+ *
+ * <p>Users can wrap a regular {@link EventHandler event handler} in
+ * an instance of this class and provide this wrapping instance to
+ * the {@link FirstOrderIntegrator ODE solver} in order to avoid wasting
+ * time looking for uninteresting events. The wrapper will intercept
+ * the calls to the {@link EventHandler#g(double, double[]) g function}
+ * and to the {@link EventHandler#eventOccurred(double, double[], boolean)
+ * eventOccurred} method in order to ignore uninteresting events. The
+ * wrapped regular {@link EventHandler event handler} will the see only
+ * the interesting events, i.e. either only {@code increasing} events or
+ * {@code decreasing} events. the number of calls to the {@link
+ * EventHandler#g(double, double[]) g function} will also be reduced.</p>
+ *
+ * @version $Id$
+ * @since 3.2
+ */
+
+public class EventFilter implements EventHandler {
+
+ /** Number of past transformers updates stored. */
+ private static final int HISTORY_SIZE = 100;
+
+ /** Wrapped event handler. */
+ private final EventHandler rawHandler;
+
+ /** Filter to use. */
+ private final FilterType filter;
+
+ /** Transformers of the g function. */
+ private final Transformer[] transformers;
+
+ /** Update time of the transformers. */
+ private final double[] updates;
+
+ /** Indicator for forward integration. */
+ private boolean forward;
+
+ /** Extreme time encountered so far. */
+ private double extremeT;
+
+ /** Wrap an {@link EventHandler event handler}.
+ * @param rawHandler event handler to wrap
+ * @param filter filter to use
+ */
+ public EventFilter(final EventHandler rawHandler, final FilterType filter)
{
+ this.rawHandler = rawHandler;
+ this.filter = filter;
+ this.transformers = new Transformer[HISTORY_SIZE];
+ this.updates = new double[HISTORY_SIZE];
+ }
+
+ /** {@inheritDoc} */
+ public void init(double t0, double[] y0, double t) {
+
+ // delegate to raw handler
+ rawHandler.init(t0, y0, t);
+
+ // initialize events triggering logic
+ forward = t >= t0;
+ extremeT = forward ? Double.NEGATIVE_INFINITY :
Double.POSITIVE_INFINITY;
+ Arrays.fill(transformers, Transformer.UNINITIALIZED);
+ Arrays.fill(updates, extremeT);
+
+ }
+
+ /** {@inheritDoc} */
+ public double g(double t, double[] y) {
+
+ final double rawG = rawHandler.g(t, y);
+
+ // search which transformer should be applied to g
+ if (forward) {
+ final int last = transformers.length - 1;
+ if (extremeT < t) {
+ // we are at the forward end of the history
+
+ // check if a new rough root has been crossed
+ final Transformer previous = transformers[last];
+ final Transformer next =
filter.selectTransformer(previous, rawG, forward);
+ if (next != previous) {
+ // there is a root somewhere between extremeT end t
+ // the new transformer, which is valid on both sides of
the root,
+ // so it is valid for t (this is how we have just computed
it above),
+ // but it was already valid before, so we store the switch
at extremeT
+ // for safety, to ensure the previous transformer is not
applied too
+ // close of the root
+ System.arraycopy(updates, 1, updates, 0, last);
+ System.arraycopy(transformers, 1, transformers, 0, last);
+ updates[last] = extremeT;
+ transformers[last] = next;
+ }
+
+ extremeT = t;
+
+ // apply the transform
+ return next.transformed(rawG);
+
+ } else {
+ // we are in the middle of the history
+
+ // select the transformer
+ for (int i = last; i > 0; --i) {
+ if (updates[i] <= t) {
+ // apply the transform
+ return transformers[i].transformed(rawG);
+ }
+ }
+
+ return transformers[0].transformed(rawG);
+
+ }
+ } else {
+ if (t < extremeT) {
+ // we are at the backward end of the history
+
+ // check if a new rough root has been crossed
+ final Transformer previous = transformers[0];
+ final Transformer next =
filter.selectTransformer(previous, rawG, forward);
+ if (next != previous) {
+ // there is a root somewhere between extremeT end t
+ // the new transformer, which is valid on both sides of
the root,
+ // so it is valid for t (this is how we have just computed
it above),
+ // but it was already valid before, so we store the switch
at extremeT
+ // for safety, to ensure the previous transformer is not
applied too
+ // close of the root
+ System.arraycopy(updates, 0, updates, 1,
updates.length - 1);
+ System.arraycopy(transformers, 0, transformers, 1,
transformers.length - 1);
+ updates[0] = extremeT;
+ transformers[0] = next;
+ }
+
+ extremeT = t;
+
+ // apply the transform
+ return next.transformed(rawG);
+
+ } else {
+ // we are in the middle of the history
+
+ // select the transformer
+ for (int i = 0; i < updates.length - 1; ++i) {
+ if (t <= updates[i]) {
+ // apply the transform
+ return transformers[i].transformed(rawG);
+ }
+ }
+
+ return transformers[updates.length - 1].transformed(rawG);
+
+ }
+ }
+
+ }
+
+ /** {@inheritDoc} */
+ public Action eventOccurred(double t, double[] y, boolean increasing) {
+ // delegate to raw handler, fixing increasing status on the fly
+ return rawHandler.eventOccurred(t, y, filter.getTriggeredIncreasing());
+ }
+
+ /** {@inheritDoc} */
+ public void resetState(double t, double[] y) {
+ // delegate to raw handler
+ rawHandler.resetState(t, y);
+ }
+
+}
Propchange:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/events/EventFilter.java
------------------------------------------------------------------------------
svn:eol-style = native
Propchange:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/events/EventFilter.java
------------------------------------------------------------------------------
svn:keywords = "Author Date Id Revision"
Added:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/events/FilterType.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/events/FilterType.java?rev=1458298&view=auto
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/events/FilterType.java
(added)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/events/FilterType.java
Tue Mar 19 14:09:58 2013
@@ -0,0 +1,399 @@
+/*
+ * 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.math3.ode.events;
+
+import org.apache.commons.math3.exception.MathInternalError;
+import org.apache.commons.math3.util.Precision;
+
+
+/** Enumerate for {@link EventFilter filtering events}.
+ *
+ * @version $Id$
+ * @since 3.2
+ */
+
+public enum FilterType {
+
+ /** Constant for triggering only decreasing events.
+ * <p>When this filter is used, the wrapped {@link EventHandler
+ * event handler} {@link EventHandler#eventOccurred(double, double[],
+ * boolean) eventOccurred} method will be called <em>only</em> with
+ * its {@code increasing} argument set to false.</p>
+ */
+ TRIGGER_ONLY_DECREASING_EVENTS {
+
+ /** {@inheritDoc} */
+ protected boolean getTriggeredIncreasing() {
+ return false;
+ }
+
+ /** {@inheritDoc}
+ * <p>
+ * states scheduling for computing h(t,y) as an altered version of
g(t, y)
+ * <ul>
+ * <li>0 are triggered events for which a zero is produced (here
decreasing events)</li>
+ * <li>X are ignored events for which zero is masked (here increasing
events)</li>
+ * </ul>
+ * </p>
+ * <pre>
+ * g(t)
+ * ___ ___ ___
+ * / \ / \ / \
+ * / \ / \ / \
+ * / g>0 \ / g>0 \ / g>0 \
+ * / \ / \ / \
+ * ----- X --------- 0 --------- X --------- 0 --------- X ---------
0 ---
+ * / \ / \ /
\
+ * / \ g<0 / \ g<0 /
\ g<0
+ * / \ / \ /
\ /
+ * ___/ \___/ \___/
\___/
+ * </pre>
+ * <pre>
+ * h(t,y)) as an alteration of g(t,y)
+ * ___ ___ ___
+ * \ / \ / \ / \
+ * \ / \ h=+g / \ / \
+ * \ / \ h=min(-s,-g,+g) / \ / \
+ * \_/ \ / \_/ \
+ * ------ ---------- 0 ----------_---------- 0 ---------------------
0 ---
+ * \ / \ /
\
+ * h=max(+s,-g,+g) \ / \ / h=max(+s,-g,+g)
\
+ * \ / \ / h=-g
\ /
+ * \___/ \___/
\___/
+ * </pre>
+ * <p>
+ * As shown by the figure above, several expressions are used to
compute h,
+ * depending on the current state:
+ * <ul>
+ * <li>h = max(+s,-g,+g)</li>
+ * <li>h = +g</li>
+ * <li>h = min(-s,-g,+g)</li>
+ * <li>h = -g</li>
+ * </ul>
+ * where s is a tiny positive value: {@link Precision#SAFE_MIN}.
+ * </p>
+ */
+ protected Transformer selectTransformer(final Transformer previous,
+ final double g, final boolean
forward) {
+ if (forward) {
+ switch (previous) {
+ case UNINITIALIZED :
+ // we are initializing the first point
+ if (g > 0) {
+ // initialize as if previous root (i.e. backward
one) was an ignored increasing event
+ return Transformer.MAX;
+ } else if (g < 0) {
+ // initialize as if previous root (i.e. backward
one) was a triggered decreasing event
+ return Transformer.PLUS;
+ } else {
+ // we are exactly at a root, we don't know if it
is an increasing
+ // or a decreasing event, we remain in
uninitialized state
+ return Transformer.UNINITIALIZED;
+ }
+ case PLUS :
+ if (g >= 0) {
+ // we have crossed the zero line on an ignored
increasing event,
+ // we must change the transformer
+ return Transformer.MIN;
+ } else {
+ // we are still in the same status
+ return previous;
+ }
+ case MINUS :
+ if (g >= 0) {
+ // we have crossed the zero line on an ignored
increasing event,
+ // we must change the transformer
+ return Transformer.MAX;
+ } else {
+ // we are still in the same status
+ return previous;
+ }
+ case MIN :
+ if (g <= 0) {
+ // we have crossed the zero line on a triggered
decreasing event,
+ // we must change the transformer
+ return Transformer.MINUS;
+ } else {
+ // we are still in the same status
+ return previous;
+ }
+ case MAX :
+ if (g <= 0) {
+ // we have crossed the zero line on a triggered
decreasing event,
+ // we must change the transformer
+ return Transformer.PLUS;
+ } else {
+ // we are still in the same status
+ return previous;
+ }
+ default :
+ // this should never happen
+ throw new MathInternalError();
+ }
+ } else {
+ switch (previous) {
+ case UNINITIALIZED :
+ // we are initializing the first point
+ if (g > 0) {
+ // initialize as if previous root (i.e. forward
one) was a triggered decreasing event
+ return Transformer.MINUS;
+ } else if (g < 0) {
+ // initialize as if previous root (i.e. forward
one) was an ignored increasing event
+ return Transformer.MIN;
+ } else {
+ // we are exactly at a root, we don't know if it
is an increasing
+ // or a decreasing event, we remain in
uninitialized state
+ return Transformer.UNINITIALIZED;
+ }
+ case PLUS :
+ if (g <= 0) {
+ // we have crossed the zero line on an ignored
increasing event,
+ // we must change the transformer
+ return Transformer.MAX;
+ } else {
+ // we are still in the same status
+ return previous;
+ }
+ case MINUS :
+ if (g <= 0) {
+ // we have crossed the zero line on an ignored
increasing event,
+ // we must change the transformer
+ return Transformer.MIN;
+ } else {
+ // we are still in the same status
+ return previous;
+ }
+ case MIN :
+ if (g >= 0) {
+ // we have crossed the zero line on a triggered
decreasing event,
+ // we must change the transformer
+ return Transformer.PLUS;
+ } else {
+ // we are still in the same status
+ return previous;
+ }
+ case MAX :
+ if (g >= 0) {
+ // we have crossed the zero line on a triggered
decreasing event,
+ // we must change the transformer
+ return Transformer.MINUS;
+ } else {
+ // we are still in the same status
+ return previous;
+ }
+ default :
+ // this should never happen
+ throw new MathInternalError();
+ }
+ }
+ }
+
+ },
+
+ /** Constant for triggering only increasing events.
+ * <p>When this filter is used, the wrapped {@link EventHandler
+ * event handler} {@link EventHandler#eventOccurred(double, double[],
+ * boolean) eventOccurred} method will be called <em>only</em> with
+ * its {@code increasing} argument set to true.</p>
+ */
+ TRIGGER_ONLY_INCREASING_EVENTS {
+
+ /** {@inheritDoc} */
+ protected boolean getTriggeredIncreasing() {
+ return true;
+ }
+
+ /** {@inheritDoc}
+ * <p>
+ * states scheduling for computing h(t,y) as an altered version of
g(t, y)
+ * <ul>
+ * <li>0 are triggered events for which a zero is produced (here
increasing events)</li>
+ * <li>X are ignored events for which zero is masked (here decreasing
events)</li>
+ * </ul>
+ * </p>
+ * <pre>
+ * g(t)
+ * ___ ___ ___
+ * / \ / \ / \
+ * / \ / \ / \
+ * / g>0 \ / g>0 \ / g>0 \
+ * / \ / \ / \
+ * ----- 0 --------- X --------- 0 --------- X --------- 0 ---------
X ---
+ * / \ / \ /
\
+ * / \ g<0 / \ g<0 /
\ g<0
+ * / \ / \ /
\ /
+ * ___/ \___/ \___/
\___/
+ * </pre>
+ * <pre>
+ * h(t,y)) as an alteration of g(t,y)
+ * ___ ___
+ * \ / \ / \
+ * \ h=-g / \ / \ h=-g
+ * \ h=min(-s,-g,+g) / \ / \
h=min(-s,-g,+g)
+ * \ / \_/ \
+ * ------0 ----------_---------- 0 --------------------- 0 ---------
_ ---
+ * \ / \ / \ /
\
+ * \ / \ / h=max(+s,-g,+g) \ /
\
+ * \ / \ / h=+g \ /
\ /
+ * \___/ \___/ \___/
\___/
+ * </pre>
+ * <p>
+ * As shown by the figure above, several expressions are used to
compute h,
+ * depending on the current state:
+ * <ul>
+ * <li>h = max(+s,-g,+g)</li>
+ * <li>h = +g</li>
+ * <li>h = min(-s,-g,+g)</li>
+ * <li>h = -g</li>
+ * </ul>
+ * where s is a tiny positive value: {@link Precision#SAFE_MIN}.
+ * </p>
+ */
+ protected Transformer selectTransformer(final Transformer previous,
+ final double g, final boolean
forward) {
+ if (forward) {
+ switch (previous) {
+ case UNINITIALIZED :
+ // we are initializing the first point
+ if (g > 0) {
+ // initialize as if previous root (i.e. backward
one) was a triggered increasing event
+ return Transformer.PLUS;
+ } else if (g < 0) {
+ // initialize as if previous root (i.e. backward
one) was an ignored decreasing event
+ return Transformer.MIN;
+ } else {
+ // we are exactly at a root, we don't know if it
is an increasing
+ // or a decreasing event, we remain in
uninitialized state
+ return Transformer.UNINITIALIZED;
+ }
+ case PLUS :
+ if (g <= 0) {
+ // we have crossed the zero line on an ignored
decreasing event,
+ // we must change the transformer
+ return Transformer.MAX;
+ } else {
+ // we are still in the same status
+ return previous;
+ }
+ case MINUS :
+ if (g <= 0) {
+ // we have crossed the zero line on an ignored
decreasing event,
+ // we must change the transformer
+ return Transformer.MIN;
+ } else {
+ // we are still in the same status
+ return previous;
+ }
+ case MIN :
+ if (g >= 0) {
+ // we have crossed the zero line on a triggered
increasing event,
+ // we must change the transformer
+ return Transformer.PLUS;
+ } else {
+ // we are still in the same status
+ return previous;
+ }
+ case MAX :
+ if (g >= 0) {
+ // we have crossed the zero line on a triggered
increasing event,
+ // we must change the transformer
+ return Transformer.MINUS;
+ } else {
+ // we are still in the same status
+ return previous;
+ }
+ default :
+ // this should never happen
+ throw new MathInternalError();
+ }
+ } else {
+ switch (previous) {
+ case UNINITIALIZED :
+ // we are initializing the first point
+ if (g > 0) {
+ // initialize as if previous root (i.e. forward
one) was an ignored decreasing event
+ return Transformer.MAX;
+ } else if (g < 0) {
+ // initialize as if previous root (i.e. forward
one) was a triggered increasing event
+ return Transformer.MINUS;
+ } else {
+ // we are exactly at a root, we don't know if it
is an increasing
+ // or a decreasing event, we remain in
uninitialized state
+ return Transformer.UNINITIALIZED;
+ }
+ case PLUS :
+ if (g >= 0) {
+ // we have crossed the zero line on an ignored
decreasing event,
+ // we must change the transformer
+ return Transformer.MIN;
+ } else {
+ // we are still in the same status
+ return previous;
+ }
+ case MINUS :
+ if (g >= 0) {
+ // we have crossed the zero line on an ignored
decreasing event,
+ // we must change the transformer
+ return Transformer.MAX;
+ } else {
+ // we are still in the same status
+ return previous;
+ }
+ case MIN :
+ if (g <= 0) {
+ // we have crossed the zero line on a triggered
increasing event,
+ // we must change the transformer
+ return Transformer.MINUS;
+ } else {
+ // we are still in the same status
+ return previous;
+ }
+ case MAX :
+ if (g <= 0) {
+ // we have crossed the zero line on a triggered
increasing event,
+ // we must change the transformer
+ return Transformer.PLUS;
+ } else {
+ // we are still in the same status
+ return previous;
+ }
+ default :
+ // this should never happen
+ throw new MathInternalError();
+ }
+ }
+ }
+
+ };
+
+ /** Get the increasing status of triggered events.
+ * @return true if triggered events are increasing events
+ */
+ protected abstract boolean getTriggeredIncreasing();
+
+ /** Get next function transformer in the specified direction.
+ * @param previous transformer active on the previous point with respect
+ * to integration direction (may be null if no previous point is known)
+ * @param g current value of the g function
+ * @param forward true if integration goes forward
+ * @return next transformer transformer
+ */
+ protected abstract Transformer selectTransformer(Transformer previous,
+ double g, boolean
forward);
+
+}
Propchange:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/events/FilterType.java
------------------------------------------------------------------------------
svn:eol-style = native
Propchange:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/events/FilterType.java
------------------------------------------------------------------------------
svn:keywords = "Author Date Id Revision"
Added:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/events/Transformer.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/events/Transformer.java?rev=1458298&view=auto
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/events/Transformer.java
(added)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/events/Transformer.java
Tue Mar 19 14:09:58 2013
@@ -0,0 +1,103 @@
+/*
+ * 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.math3.ode.events;
+
+import org.apache.commons.math3.util.FastMath;
+import org.apache.commons.math3.util.Precision;
+
+
+/** Transformer for {@link EventHandler#g(double, double[]) g functions}.
+ * @see EventFilter
+ * @see FilterType
+ * @version $Id$
+ * @since 3.2
+ */
+enum Transformer {
+
+ /** Transformer computing transformed = 0.
+ * <p>
+ * This transformer is used when we initialize the filter, until we get at
+ * least one non-zero value to select the proper transformer.
+ * </p>
+ */
+ UNINITIALIZED {
+ /** {@inheritDoc} */
+ protected double transformed(final double g) {
+ return 0;
+ }
+ },
+
+ /** Transformer computing transformed = g.
+ * <p>
+ * When this transformer is applied, the roots of the original function
+ * are preserved, with the same {@code increasing/decreasing} status.
+ * </p>
+ */
+ PLUS {
+ /** {@inheritDoc} */
+ protected double transformed(final double g) {
+ return g;
+ }
+ },
+
+ /** Transformer computing transformed = -g.
+ * <p>
+ * When this transformer is applied, the roots of the original function
+ * are preserved, with reversed {@code increasing/decreasing} status.
+ * </p>
+ */
+ MINUS {
+ /** {@inheritDoc} */
+ protected double transformed(final double g) {
+ return -g;
+ }
+ },
+
+ /** Transformer computing transformed = min(-{@link Precision#SAFE_MIN},
-g, +g).
+ * <p>
+ * When this transformer is applied, the transformed function is
+ * guaranteed to be always strictly negative (i.e. there are no roots).
+ * </p>
+ */
+ MIN {
+ /** {@inheritDoc} */
+ protected double transformed(final double g) {
+ return FastMath.min(-Precision.SAFE_MIN, FastMath.min(-g, +g));
+ }
+ },
+
+ /** Transformer computing transformed = max(+{@link Precision#SAFE_MIN},
-g, +g).
+ * <p>
+ * When this transformer is applied, the transformed function is
+ * guaranteed to be always strictly positive (i.e. there are no roots).
+ * </p>
+ */
+ MAX {
+ /** {@inheritDoc} */
+ protected double transformed(final double g) {
+ return FastMath.max(+Precision.SAFE_MIN, FastMath.max(-g, +g));
+ }
+ };
+
+ /** Transform value of function g.
+ * @param g raw value of function g
+ * @return transformed value of function g
+ */
+ protected abstract double transformed(double g);
+
+}
Propchange:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/events/Transformer.java
------------------------------------------------------------------------------
svn:eol-style = native
Propchange:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/ode/events/Transformer.java
------------------------------------------------------------------------------
svn:keywords = "Author Date Id Revision"
Added:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/ode/events/EventFilterTest.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/ode/events/EventFilterTest.java?rev=1458298&view=auto
==============================================================================
---
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/ode/events/EventFilterTest.java
(added)
+++
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/ode/events/EventFilterTest.java
Tue Mar 19 14:09:58 2013
@@ -0,0 +1,279 @@
+/*
+ * 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.math3.ode.events;
+
+import java.io.IOException;
+import java.io.PrintStream;
+
+import org.apache.commons.math3.analysis.solvers.BracketingNthOrderBrentSolver;
+import org.apache.commons.math3.exception.DimensionMismatchException;
+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.FirstOrderDifferentialEquations;
+import org.apache.commons.math3.ode.FirstOrderIntegrator;
+import org.apache.commons.math3.ode.nonstiff.DormandPrince853Integrator;
+import org.apache.commons.math3.random.RandomGenerator;
+import org.apache.commons.math3.random.Well19937a;
+import org.apache.commons.math3.util.FastMath;
+import org.junit.Assert;
+import org.junit.Test;
+
+public class EventFilterTest {
+
+ @Test
+ public void testHistoryIncreasingForward() {
+
+ // start point: g > 0
+ testHistory(FilterType.TRIGGER_ONLY_INCREASING_EVENTS,
+ 0.5 * FastMath.PI, 30.5 * FastMath.PI, FastMath.PI, -1);
+
+ // start point: g = 0
+ testHistory(FilterType.TRIGGER_ONLY_INCREASING_EVENTS,
+ 0, 30.5 * FastMath.PI, FastMath.PI, -1);
+
+ // start point: g < 0
+ testHistory(FilterType.TRIGGER_ONLY_INCREASING_EVENTS,
+ 1.5 * FastMath.PI, 30.5 * FastMath.PI, FastMath.PI, +1);
+
+ }
+
+ @Test
+ public void testHistoryIncreasingBackward() {
+
+ // start point: g > 0
+ testHistory(FilterType.TRIGGER_ONLY_INCREASING_EVENTS,
+ 0.5 * FastMath.PI, -30.5 * FastMath.PI, FastMath.PI, -1);
+
+ // start point: g = 0
+ testHistory(FilterType.TRIGGER_ONLY_INCREASING_EVENTS,
+ 0, -30.5 * FastMath.PI, FastMath.PI, +1);
+
+ // start point: g < 0
+ testHistory(FilterType.TRIGGER_ONLY_INCREASING_EVENTS,
+ 1.5 * FastMath.PI, -30.5 * FastMath.PI, FastMath.PI, -1);
+
+ }
+
+ @Test
+ public void testHistoryDecreasingForward() {
+
+ // start point: g > 0
+ testHistory(FilterType.TRIGGER_ONLY_DECREASING_EVENTS,
+ 0.5 * FastMath.PI, 30.5 * FastMath.PI, 0, +1);
+
+ // start point: g = 0
+ testHistory(FilterType.TRIGGER_ONLY_DECREASING_EVENTS,
+ 0, 30.5 * FastMath.PI, 0, +1);
+
+ // start point: g < 0
+ testHistory(FilterType.TRIGGER_ONLY_DECREASING_EVENTS,
+ 1.5 * FastMath.PI, 30.5 * FastMath.PI, 0, +1);
+
+ }
+
+ @Test
+ public void testHistoryDecreasingBackward() {
+
+ // start point: g > 0
+ testHistory(FilterType.TRIGGER_ONLY_DECREASING_EVENTS,
+ 0.5 * FastMath.PI, -30.5 * FastMath.PI, 0, -1);
+
+ // start point: g = 0
+ testHistory(FilterType.TRIGGER_ONLY_DECREASING_EVENTS,
+ 0, -30.5 * FastMath.PI, 0, -1);
+
+ // start point: g < 0
+ testHistory(FilterType.TRIGGER_ONLY_DECREASING_EVENTS,
+ 1.5 * FastMath.PI, -30.5 * FastMath.PI, 0, +1);
+
+ }
+
+ public void testHistory(FilterType type, double t0, double t1, double
refSwitch, double signEven) {
+ Event onlyIncreasing = new Event(false, true);
+ EventFilter eventFilter =
+ new EventFilter(onlyIncreasing, type);
+ eventFilter.init(t0, new double[] {1.0, 0.0}, t1);
+
+ // first pass to set up switches history for a long period
+ double h = FastMath.copySign(0.05, t1 - t0);
+ double n = (int) FastMath.floor((t1 - t0) / h);
+ for (int i = 0; i < n; ++i) {
+ double t = t0 + i * h;
+ eventFilter.g(t, new double[] { FastMath.sin(t), FastMath.cos(t)
});
+ }
+
+ // verify old events are preserved, even if randomly accessed
+ RandomGenerator rng = new Well19937a(0xb0e7401265af8cd3l);
+ try {
+ PrintStream out = new PrintStream("/home/luc/x.dat");
+ for (int i = 0; i < 5000; i++) {
+ double t = t0 + (t1 - t0) * rng.nextDouble();
+ double g = eventFilter.g(t, new double[] { FastMath.sin(t),
FastMath.cos(t) });
+ int turn = (int) FastMath.floor((t - refSwitch) / (2 *
FastMath.PI));
+ out.println(t + " " + g);
+ if (turn % 2 == 0) {
+ Assert.assertEquals( signEven * FastMath.sin(t), g, 1.0e-10);
+ } else {
+ Assert.assertEquals(-signEven * FastMath.sin(t), g, 1.0e-10);
+ }
+ }
+ out.close();
+ } catch (IOException ioe) {
+ Assert.fail(ioe.getLocalizedMessage());
+ }
+
+ }
+
+ @Test
+ public void testIncreasingOnly()
+ throws DimensionMismatchException, NumberIsTooSmallException,
+ MaxCountExceededException, NoBracketingException {
+ double e = 1e-15;
+ FirstOrderIntegrator integrator;
+ integrator = new DormandPrince853Integrator(1.0e-3, 100.0, 1e-7, 1e-7);
+ Event allEvents = new Event(true, true);
+ integrator.addEventHandler(allEvents, 0.1, e, 1000,
+ new BracketingNthOrderBrentSolver(1.0e-7,
5));
+ Event onlyIncreasing = new Event(false, true);
+ integrator.addEventHandler(new EventFilter(onlyIncreasing,
+
FilterType.TRIGGER_ONLY_INCREASING_EVENTS),
+ 0.1, e, 100,
+ new BracketingNthOrderBrentSolver(1.0e-7,
5));
+ double t0 = 0.5 * FastMath.PI;
+ double tEnd = 5.5 * FastMath.PI;
+ double[] y = { 0.0, 1.0 };
+ Assert.assertEquals(tEnd,
+ integrator.integrate(new SineCosine(), t0, y,
tEnd, y),
+ 1.0e-7);
+
+ Assert.assertEquals(5, allEvents.getEventCount());
+ Assert.assertEquals(2, onlyIncreasing.getEventCount());
+
+ }
+
+ @Test
+ public void testDecreasingOnly()
+ throws DimensionMismatchException, NumberIsTooSmallException,
+ MaxCountExceededException, NoBracketingException {
+ double e = 1e-15;
+ FirstOrderIntegrator integrator;
+ integrator = new DormandPrince853Integrator(1.0e-3, 100.0, 1e-7, 1e-7);
+ Event allEvents = new Event(true, true);
+ integrator.addEventHandler(allEvents, 0.1, e, 1000,
+ new BracketingNthOrderBrentSolver(1.0e-7,
5));
+ Event onlyDecreasing = new Event(true, false);
+ integrator.addEventHandler(new EventFilter(onlyDecreasing,
+
FilterType.TRIGGER_ONLY_DECREASING_EVENTS),
+ 0.1, e, 1000,
+ new BracketingNthOrderBrentSolver(1.0e-7,
5));
+ double t0 = 0.5 * FastMath.PI;
+ double tEnd = 5.5 * FastMath.PI;
+ double[] y = { 0.0, 1.0 };
+ Assert.assertEquals(tEnd,
+ integrator.integrate(new SineCosine(), t0, y,
tEnd, y),
+ 1.0e-7);
+
+ Assert.assertEquals(5, allEvents.getEventCount());
+ Assert.assertEquals(3, onlyDecreasing.getEventCount());
+
+ }
+
+ @Test
+ public void testTwoOppositeFilters()
+ throws DimensionMismatchException, NumberIsTooSmallException,
+ MaxCountExceededException, NoBracketingException {
+ double e = 1e-15;
+ FirstOrderIntegrator integrator;
+ integrator = new DormandPrince853Integrator(1.0e-3, 100.0, 1e-7, 1e-7);
+ Event allEvents = new Event(true, true);
+ integrator.addEventHandler(allEvents, 0.1, e, 1000,
+ new BracketingNthOrderBrentSolver(1.0e-7,
5));
+ Event onlyIncreasing = new Event(false, true);
+ integrator.addEventHandler(new EventFilter(onlyIncreasing,
+
FilterType.TRIGGER_ONLY_INCREASING_EVENTS),
+ 0.1, e, 1000,
+ new BracketingNthOrderBrentSolver(1.0e-7,
5));
+ Event onlyDecreasing = new Event(true, false);
+ integrator.addEventHandler(new EventFilter(onlyDecreasing,
+
FilterType.TRIGGER_ONLY_DECREASING_EVENTS),
+ 0.1, e, 1000,
+ new BracketingNthOrderBrentSolver(1.0e-7,
5));
+ double t0 = 0.5 * FastMath.PI;
+ double tEnd = 5.5 * FastMath.PI;
+ double[] y = { 0.0, 1.0 };
+ Assert.assertEquals(tEnd,
+ integrator.integrate(new SineCosine(), t0, y,
tEnd, y),
+ 1.0e-7);
+
+ Assert.assertEquals(5, allEvents.getEventCount());
+ Assert.assertEquals(2, onlyIncreasing.getEventCount());
+ Assert.assertEquals(3, onlyDecreasing.getEventCount());
+
+ }
+
+ private static class SineCosine implements FirstOrderDifferentialEquations
{
+ public int getDimension() {
+ return 2;
+ }
+
+ public void computeDerivatives(double t, double[] y, double[] yDot) {
+ yDot[0] = y[1];
+ yDot[1] = -y[0];
+ }
+ }
+
+ /** State events for this unit test. */
+ protected static class Event implements EventHandler {
+
+ private final boolean expectDecreasing;
+ private final boolean expectIncreasing;
+ private int eventCount;
+
+ public Event(boolean expectDecreasing, boolean expectIncreasing) {
+ this.expectDecreasing = expectDecreasing;
+ this.expectIncreasing = expectIncreasing;
+ }
+
+ public int getEventCount() {
+ return eventCount;
+ }
+
+ public void init(double t0, double[] y0, double t) {
+ eventCount = 0;
+ }
+
+ public double g(double t, double[] y) {
+ return y[0];
+ }
+
+ public Action eventOccurred(double t, double[] y, boolean increasing) {
+ if (increasing) {
+ Assert.assertTrue(expectIncreasing);
+ } else {
+ Assert.assertTrue(expectDecreasing);
+ }
+ eventCount++;
+ return Action.RESET_STATE;
+ }
+
+ public void resetState(double t, double[] y) {
+ // in fact, we don't really reset anything for this test
+ }
+
+ }
+}
Propchange:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/ode/events/EventFilterTest.java
------------------------------------------------------------------------------
svn:eol-style = native
Propchange:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/ode/events/EventFilterTest.java
------------------------------------------------------------------------------
svn:keywords = "Author Date Id Revision"
