001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017
018package org.apache.commons.math3.ode.events;
019
020import java.util.Arrays;
021
022/** Wrapper used to detect only increasing or decreasing events.
023 *
024 * <p>General {@link EventHandler events} are defined implicitly
025 * by a {@link EventHandler#g(double, double[]) g function} crossing
026 * zero. This function needs to be continuous in the event neighborhood,
027 * and its sign must remain consistent between events. This implies that
028 * during an ODE integration, events triggered are alternately events
029 * for which the function increases from negative to positive values,
030 * and events for which the function decreases from positive to
031 * negative values.
032 * </p>
033 *
034 * <p>Sometimes, users are only interested in one type of event (say
035 * increasing events for example) and not in the other type. In these
036 * cases, looking precisely for all events location and triggering
037 * events that will later be ignored is a waste of computing time.</p>
038 *
039 * <p>Users can wrap a regular {@link EventHandler event handler} in
040 * an instance of this class and provide this wrapping instance to
041 * the {@link org.apache.commons.math3.ode.FirstOrderIntegrator ODE solver}
042 * in order to avoid wasting time looking for uninteresting events.
043 * The wrapper will intercept the calls to the {@link
044 * EventHandler#g(double, double[]) g function} and to the {@link
045 * EventHandler#eventOccurred(double, double[], boolean)
046 * eventOccurred} method in order to ignore uninteresting events. The
047 * wrapped regular {@link EventHandler event handler} will the see only
048 * the interesting events, i.e. either only {@code increasing} events or
049 * {@code decreasing} events. the number of calls to the {@link
050 * EventHandler#g(double, double[]) g function} will also be reduced.</p>
051 *
052 * @since 3.2
053 */
054
055public class EventFilter implements EventHandler {
056
057    /** Number of past transformers updates stored. */
058    private static final int HISTORY_SIZE = 100;
059
060    /** Wrapped event handler. */
061    private final EventHandler rawHandler;
062
063    /** Filter to use. */
064    private final FilterType filter;
065
066    /** Transformers of the g function. */
067    private final Transformer[] transformers;
068
069    /** Update time of the transformers. */
070    private final double[] updates;
071
072    /** Indicator for forward integration. */
073    private boolean forward;
074
075    /** Extreme time encountered so far. */
076    private double extremeT;
077
078    /** Wrap an {@link EventHandler event handler}.
079     * @param rawHandler event handler to wrap
080     * @param filter filter to use
081     */
082    public EventFilter(final EventHandler rawHandler, final FilterType filter) {
083        this.rawHandler   = rawHandler;
084        this.filter       = filter;
085        this.transformers = new Transformer[HISTORY_SIZE];
086        this.updates      = new double[HISTORY_SIZE];
087    }
088
089    /**  {@inheritDoc} */
090    public void init(double t0, double[] y0, double t) {
091
092        // delegate to raw handler
093        rawHandler.init(t0, y0, t);
094
095        // initialize events triggering logic
096        forward  = t >= t0;
097        extremeT = forward ? Double.NEGATIVE_INFINITY : Double.POSITIVE_INFINITY;
098        Arrays.fill(transformers, Transformer.UNINITIALIZED);
099        Arrays.fill(updates, extremeT);
100
101    }
102
103    /**  {@inheritDoc} */
104    public double g(double t, double[] y) {
105
106        final double rawG = rawHandler.g(t, y);
107
108        // search which transformer should be applied to g
109        if (forward) {
110            final int last = transformers.length - 1;
111            if (extremeT < t) {
112                // we are at the forward end of the history
113
114                // check if a new rough root has been crossed
115                final Transformer previous = transformers[last];
116                final Transformer next     = filter.selectTransformer(previous, rawG, forward);
117                if (next != previous) {
118                    // there is a root somewhere between extremeT and t.
119                    // the new transformer is valid for t (this is how we have just computed
120                    // it above), but it is in fact valid on both sides of the root, so
121                    // it was already valid before t and even up to previous time. We store
122                    // the switch at extremeT for safety, to ensure the previous transformer
123                    // is not applied too close of the root
124                    System.arraycopy(updates,      1, updates,      0, last);
125                    System.arraycopy(transformers, 1, transformers, 0, last);
126                    updates[last]      = extremeT;
127                    transformers[last] = next;
128                }
129
130                extremeT = t;
131
132                // apply the transform
133                return next.transformed(rawG);
134
135            } else {
136                // we are in the middle of the history
137
138                // select the transformer
139                for (int i = last; i > 0; --i) {
140                    if (updates[i] <= t) {
141                        // apply the transform
142                        return transformers[i].transformed(rawG);
143                    }
144                }
145
146                return transformers[0].transformed(rawG);
147
148            }
149        } else {
150            if (t < extremeT) {
151                // we are at the backward end of the history
152
153                // check if a new rough root has been crossed
154                final Transformer previous = transformers[0];
155                final Transformer next     = filter.selectTransformer(previous, rawG, forward);
156                if (next != previous) {
157                    // there is a root somewhere between extremeT and t.
158                    // the new transformer is valid for t (this is how we have just computed
159                    // it above), but it is in fact valid on both sides of the root, so
160                    // it was already valid before t and even up to previous time. We store
161                    // the switch at extremeT for safety, to ensure the previous transformer
162                    // is not applied too close of the root
163                    System.arraycopy(updates,      0, updates,      1, updates.length - 1);
164                    System.arraycopy(transformers, 0, transformers, 1, transformers.length - 1);
165                    updates[0]      = extremeT;
166                    transformers[0] = next;
167                }
168
169                extremeT = t;
170
171                // apply the transform
172                return next.transformed(rawG);
173
174            } else {
175                // we are in the middle of the history
176
177                // select the transformer
178                for (int i = 0; i < updates.length - 1; ++i) {
179                    if (t <= updates[i]) {
180                        // apply the transform
181                        return transformers[i].transformed(rawG);
182                    }
183                }
184
185                return transformers[updates.length - 1].transformed(rawG);
186
187            }
188       }
189
190    }
191
192    /**  {@inheritDoc} */
193    public Action eventOccurred(double t, double[] y, boolean increasing) {
194        // delegate to raw handler, fixing increasing status on the fly
195        return rawHandler.eventOccurred(t, y, filter.getTriggeredIncreasing());
196    }
197
198    /**  {@inheritDoc} */
199    public void resetState(double t, double[] y) {
200        // delegate to raw handler
201        rawHandler.resetState(t, y);
202    }
203
204}