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
145 // delegate to raw handler
146 rawDetector.init(initialState, finalTime);
147
148 // initialize events triggering logic
149 forward = finalTime >= initialState.getTime();
150 extremeT = forward ? Double.NEGATIVE_INFINITY : Double.POSITIVE_INFINITY;
151 Arrays.fill(transformers, Transformer.UNINITIALIZED);
152 Arrays.fill(updates, extremeT);
153
154 }
155
156 /** {@inheritDoc} */
157 @Override
158 public double g(final ODEStateAndDerivative state) {
159
160 final double rawG = rawDetector.g(state);
161
162 // search which transformer should be applied to g
163 if (forward) {
164 final int last = transformers.length - 1;
165 if (extremeT < state.getTime()) {
166 // we are at the forward end of the history
167
168 // check if a new rough root has been crossed
169 final Transformer previous = transformers[last];
170 final Transformer next = filter.selectTransformer(previous, rawG, forward);
171 if (next != previous) {
172 // there is a root somewhere between extremeT and t.
173 // the new transformer is valid for t (this is how we have just computed
174 // it above), but it is in fact valid on both sides of the root, so
175 // it was already valid before t and even up to previous time. We store
176 // the switch at extremeT for safety, to ensure the previous transformer
177 // is not applied too close of the root
178 System.arraycopy(updates, 1, updates, 0, last);
179 System.arraycopy(transformers, 1, transformers, 0, last);
180 updates[last] = extremeT;
181 transformers[last] = next;
182 }
183
184 extremeT = state.getTime();
185
186 // apply the transform
187 return next.transformed(rawG);
188
189 } else {
190 // we are in the middle of the history
191
192 // select the transformer
193 for (int i = last; i > 0; --i) {
194 if (updates[i] <= state.getTime()) {
195 // apply the transform
196 return transformers[i].transformed(rawG);
197 }
198 }
199
200 return transformers[0].transformed(rawG);
201
202 }
203 } else {
204 if (state.getTime() < extremeT) {
205 // we are at the backward end of the history
206
207 // check if a new rough root has been crossed
208 final Transformer previous = transformers[0];
209 final Transformer next = filter.selectTransformer(previous, rawG, forward);
210 if (next != previous) {
211 // there is a root somewhere between extremeT and t.
212 // the new transformer is valid for t (this is how we have just computed
213 // it above), but it is in fact valid on both sides of the root, so
214 // it was already valid before t and even up to previous time. We store
215 // the switch at extremeT for safety, to ensure the previous transformer
216 // is not applied too close of the root
217 System.arraycopy(updates, 0, updates, 1, updates.length - 1);
218 System.arraycopy(transformers, 0, transformers, 1, transformers.length - 1);
219 updates[0] = extremeT;
220 transformers[0] = next;
221 }
222
223 extremeT = state.getTime();
224
225 // apply the transform
226 return next.transformed(rawG);
227
228 } else {
229 // we are in the middle of the history
230
231 // select the transformer
232 for (int i = 0; i < updates.length - 1; ++i) {
233 if (state.getTime() <= updates[i]) {
234 // apply the transform
235 return transformers[i].transformed(rawG);
236 }
237 }
238
239 return transformers[updates.length - 1].transformed(rawG);
240
241 }
242 }
243
244 }
245
246 /** Local handler.
247 * @param <T> type of the event detector
248 */
249 private static class LocalHandler<T extends ODEEventDetector> implements ODEEventHandler {
250
251 /** Raw handler. */
252 private final ODEEventHandler rawHandler;
253
254 /** Simple constructor.
255 * @param rawHandler raw handler
256 */
257 LocalHandler(final ODEEventHandler rawHandler) {
258 this.rawHandler = rawHandler;
259 }
260
261 /** {@inheritDoc} */
262 @Override
263 public Action eventOccurred(final ODEStateAndDerivative state, final ODEEventDetector detector, final boolean increasing) {
264 // delegate to raw handler, fixing increasing status on the fly
265 @SuppressWarnings("unchecked")
266 final EventSlopeFilter<T> esf = (EventSlopeFilter<T>) detector;
267 return rawHandler.eventOccurred(state, esf, esf.filter.isTriggeredOnIncreasing());
268 }
269
270 /** {@inheritDoc} */
271 @Override
272 public ODEState resetState(final ODEEventDetector detector, final ODEStateAndDerivative state) {
273 // delegate to raw handler
274 @SuppressWarnings("unchecked")
275 final EventSlopeFilter<T> esf = (EventSlopeFilter<T>) detector;
276 return rawHandler.resetState(esf, state);
277 }
278
279 }
280
281 }