1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * https://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17
18 /*
19 * This is not the original file distributed by the Apache Software Foundation
20 * It has been modified by the Hipparchus project
21 */
22
23 package org.hipparchus.ode.events;
24
25 import java.util.Arrays;
26
27 import org.hipparchus.analysis.UnivariateFunction;
28 import org.hipparchus.analysis.solvers.BracketedUnivariateSolver;
29 import org.hipparchus.ode.ODEState;
30 import org.hipparchus.ode.ODEStateAndDerivative;
31
32 /** Wrapper used to detect only increasing or decreasing events.
33 *
34 * <p>General {@link ODEEventDetector events} are defined implicitly
35 * by a {@link ODEEventDetector#g(ODEStateAndDerivative) g function} crossing
36 * zero. This function needs to be continuous in the event neighborhood,
37 * and its sign must remain consistent between events. This implies that
38 * during an ODE integration, events triggered are alternately events
39 * for which the function increases from negative to positive values,
40 * and events for which the function decreases from positive to
41 * negative values.
42 * </p>
43 *
44 * <p>Sometimes, users are only interested in one type of event (say
45 * increasing events for example) and not in the other type. In these
46 * cases, looking precisely for all events location and triggering
47 * events that will later be ignored is a waste of computing time.</p>
48 *
49 * <p>Users can wrap a regular {@link ODEEventDetector event detector} in
50 * an instance of this class and provide this wrapping instance to
51 * the {@link org.hipparchus.ode.ODEIntegrator ODE solver}
52 * in order to avoid wasting time looking for uninteresting events.
53 * The wrapper will intercept the calls to the {@link
54 * ODEEventDetector#g(ODEStateAndDerivative) g function} and to the {@link
55 * ODEEventHandler#eventOccurred(ODEStateAndDerivative, ODEEventDetector, boolean)
56 * eventOccurred} method in order to ignore uninteresting events. The
57 * wrapped regular {@link ODEEventHandler event handler} will the see only
58 * the interesting events, i.e. either only {@code increasing} events or
59 * {@code decreasing} events. the number of calls to the {@link
60 * ODEEventDetector#g(ODEStateAndDerivative) g function} will also be reduced.</p>
61 * @param <T> type of the event detector
62 * @since 3.0
63 */
64 public class EventSlopeFilter<T extends ODEEventDetector> extends AbstractODEDetector<EventSlopeFilter<T>> {
65
66 /** Number of past transformers updates stored. */
67 private static final int HISTORY_SIZE = 100;
68
69 /** Wrapped event detector.
70 * @since 3.0
71 */
72 private final T rawDetector;
73
74 /** Filter to use. */
75 private final FilterType filter;
76
77 /** Transformers of the g function. */
78 private final Transformer[] transformers;
79
80 /** Update time of the transformers. */
81 private final double[] updates;
82
83 /** Indicator for forward integration. */
84 private boolean forward;
85
86 /** Extreme time encountered so far. */
87 private double extremeT;
88
89 /** Wrap an {@link ODEEventDetector event detector}.
90 * @param rawDetector event detector to wrap
91 * @param filter filter to use
92 * @since 3.0
93 */
94 public EventSlopeFilter(final T rawDetector, final FilterType filter) {
95 this(rawDetector.getMaxCheckInterval(), rawDetector.getMaxIterationCount(),
96 rawDetector.getSolver(), new LocalHandler<>(rawDetector.getHandler()),
97 rawDetector, filter);
98 }
99
100 /** Private constructor with full parameters.
101 * <p>
102 * This constructor is private as users are expected to use the builder
103 * API with the various {@code withXxx()} methods to set up the instance
104 * in a readable manner without using a huge amount of parameters.
105 * </p>
106 * @param maxCheck maximum checking interval (s)
107 * @param maxIter maximum number of iterations in the event time search
108 * @param solver root-finding algorithm to use to detect state events
109 * @param handler event handler to call at event occurrences
110 * @param rawDetector event detector to wrap
111 * @param filter filter to use
112 */
113 private EventSlopeFilter(final AdaptableInterval maxCheck, final int maxIter,
114 final BracketedUnivariateSolver<UnivariateFunction> solver,
115 final ODEEventHandler handler,
116 final T rawDetector, final FilterType filter) {
117 super(maxCheck, maxIter, solver, handler);
118 this.rawDetector = rawDetector;
119 this.filter = filter;
120 this.transformers = new Transformer[HISTORY_SIZE];
121 this.updates = new double[HISTORY_SIZE];
122 }
123
124 /** {@inheritDoc} */
125 @Override
126 protected EventSlopeFilter<T> create(final AdaptableInterval newMaxCheck, final int newMaxIter,
127 final BracketedUnivariateSolver<UnivariateFunction> newSolver,
128 final ODEEventHandler newHandler) {
129 return new EventSlopeFilter<>(newMaxCheck, newMaxIter, newSolver, newHandler,
130 rawDetector, filter);
131 }
132
133 /**
134 * Get the wrapped raw detector.
135 * @return the wrapped raw detector
136 */
137 public T getDetector() {
138 return rawDetector;
139 }
140
141 /** {@inheritDoc} */
142 @Override
143 public void init(final ODEStateAndDerivative initialState, double finalTime) {
144 super.init(initialState, finalTime);
145
146 // delegate to raw handler
147 rawDetector.init(initialState, finalTime);
148
149 // initialize events triggering logic
150 forward = finalTime >= initialState.getTime();
151 extremeT = forward ? Double.NEGATIVE_INFINITY : Double.POSITIVE_INFINITY;
152 Arrays.fill(transformers, Transformer.UNINITIALIZED);
153 Arrays.fill(updates, extremeT);
154
155 }
156
157 /** {@inheritDoc} */
158 @Override
159 public void reset(final ODEStateAndDerivative intermediateState, final double finalTime) {
160 super.reset(intermediateState, finalTime);
161 rawDetector.reset(intermediateState, finalTime);
162 }
163
164 /** {@inheritDoc} */
165 @Override
166 public boolean isForward() {
167 return forward;
168 }
169
170 /** {@inheritDoc} */
171 @Override
172 public double g(final ODEStateAndDerivative state) {
173
174 final double rawG = rawDetector.g(state);
175
176 // search which transformer should be applied to g
177 if (isForward()) {
178 final int last = transformers.length - 1;
179 if (extremeT < state.getTime()) {
180 // we are at the forward end of the history
181
182 // check if a new rough root has been crossed
183 final Transformer previous = transformers[last];
184 final Transformer next = filter.selectTransformer(previous, rawG, forward);
185 if (next != previous) {
186 // there is a root somewhere between extremeT and t.
187 // the new transformer is valid for t (this is how we have just computed
188 // it above), but it is in fact valid on both sides of the root, so
189 // it was already valid before t and even up to previous time. We store
190 // the switch at extremeT for safety, to ensure the previous transformer
191 // is not applied too close of the root
192 System.arraycopy(updates, 1, updates, 0, last);
193 System.arraycopy(transformers, 1, transformers, 0, last);
194 updates[last] = extremeT;
195 transformers[last] = next;
196 }
197
198 extremeT = state.getTime();
199
200 // apply the transform
201 return next.transformed(rawG);
202
203 } else {
204 // we are in the middle of the history
205
206 // select the transformer
207 for (int i = last; i > 0; --i) {
208 if (updates[i] <= state.getTime()) {
209 // apply the transform
210 return transformers[i].transformed(rawG);
211 }
212 }
213
214 return transformers[0].transformed(rawG);
215
216 }
217 } else {
218 if (state.getTime() < extremeT) {
219 // we are at the backward end of the history
220
221 // check if a new rough root has been crossed
222 final Transformer previous = transformers[0];
223 final Transformer next = filter.selectTransformer(previous, rawG, forward);
224 if (next != previous) {
225 // there is a root somewhere between extremeT and t.
226 // the new transformer is valid for t (this is how we have just computed
227 // it above), but it is in fact valid on both sides of the root, so
228 // it was already valid before t and even up to previous time. We store
229 // the switch at extremeT for safety, to ensure the previous transformer
230 // is not applied too close of the root
231 System.arraycopy(updates, 0, updates, 1, updates.length - 1);
232 System.arraycopy(transformers, 0, transformers, 1, transformers.length - 1);
233 updates[0] = extremeT;
234 transformers[0] = next;
235 }
236
237 extremeT = state.getTime();
238
239 // apply the transform
240 return next.transformed(rawG);
241
242 } else {
243 // we are in the middle of the history
244
245 // select the transformer
246 for (int i = 0; i < updates.length - 1; ++i) {
247 if (state.getTime() <= updates[i]) {
248 // apply the transform
249 return transformers[i].transformed(rawG);
250 }
251 }
252
253 return transformers[updates.length - 1].transformed(rawG);
254
255 }
256 }
257
258 }
259
260 /** Local handler.
261 * @param <T> type of the event detector
262 */
263 private static class LocalHandler<T extends ODEEventDetector> implements ODEEventHandler {
264
265 /** Raw handler. */
266 private final ODEEventHandler rawHandler;
267
268 /** Simple constructor.
269 * @param rawHandler raw handler
270 */
271 LocalHandler(final ODEEventHandler rawHandler) {
272 this.rawHandler = rawHandler;
273 }
274
275 /** {@inheritDoc} */
276 @Override
277 public Action eventOccurred(final ODEStateAndDerivative state, final ODEEventDetector detector, final boolean increasing) {
278 // delegate to raw handler, fixing increasing status on the fly
279 @SuppressWarnings("unchecked")
280 final EventSlopeFilter<T> esf = (EventSlopeFilter<T>) detector;
281 return rawHandler.eventOccurred(state, esf, esf.filter.isTriggeredOnIncreasing());
282 }
283
284 /** {@inheritDoc} */
285 @Override
286 public ODEState resetState(final ODEEventDetector detector, final ODEStateAndDerivative state) {
287 // delegate to raw handler
288 @SuppressWarnings("unchecked")
289 final EventSlopeFilter<T> esf = (EventSlopeFilter<T>) detector;
290 return rawHandler.resetState(esf, state);
291 }
292
293 }
294
295 }