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.CalculusFieldElement;
28 import org.hipparchus.Field;
29 import org.hipparchus.analysis.solvers.BracketedRealFieldUnivariateSolver;
30 import org.hipparchus.ode.FieldODEState;
31 import org.hipparchus.ode.FieldODEStateAndDerivative;
32 import org.hipparchus.util.MathArrays;
33
34 /** Wrapper used to detect only increasing or decreasing events.
35 *
36 * <p>General {@link FieldODEEventDetector events} are defined implicitly
37 * by a {@link FieldODEEventDetector#g(FieldODEStateAndDerivative) g function} crossing
38 * zero. This function needs to be continuous in the event neighborhood,
39 * and its sign must remain consistent between events. This implies that
40 * during an ODE integration, events triggered are alternately events
41 * for which the function increases from negative to positive values,
42 * and events for which the function decreases from positive to
43 * negative values.
44 * </p>
45 *
46 * <p>Sometimes, users are only interested in one type of event (say
47 * increasing events for example) and not in the other type. In these
48 * cases, looking precisely for all events location and triggering
49 * events that will later be ignored is a waste of computing time.</p>
50 *
51 * <p>Users can wrap a regular {@link FieldODEEventDetector event detector} in
52 * an instance of this class and provide this wrapping instance to
53 * the {@link org.hipparchus.ode.FieldODEIntegrator ODE solver}
54 * in order to avoid wasting time looking for uninteresting events.
55 * The wrapper will intercept the calls to the {@link
56 * FieldODEEventDetector#g(FieldODEStateAndDerivative) g function} and to the {@link
57 * FieldODEEventHandler#eventOccurred(FieldODEStateAndDerivative, FieldODEEventDetector, boolean)
58 * eventOccurred} method in order to ignore uninteresting events. The
59 * wrapped regular {@link FieldODEEventDetector event detector} will the see only
60 * the interesting events, i.e. either only {@code increasing} events or
61 * {@code decreasing} events. the number of calls to the {@link
62 * FieldODEEventDetector#g(FieldODEStateAndDerivative) g function} will also be reduced.</p>
63 *
64 * @param <T> type of the event detector
65 * @param <E> the type of the field elements
66 * @since 3.0
67 */
68
69 public class FieldEventSlopeFilter<T extends FieldODEEventDetector<E>, E extends CalculusFieldElement<E>>
70 extends AbstractFieldODEDetector<FieldEventSlopeFilter<T, E>, E> {
71
72 /** Number of past transformers updates stored. */
73 private static final int HISTORY_SIZE = 100;
74
75 /** Wrapped event detector.
76 * @since 3.0
77 */
78 private final T rawDetector;
79
80 /** Filter to use. */
81 private final FilterType filter;
82
83 /** Transformers of the g function. */
84 private final Transformer[] transformers;
85
86 /** Update time of the transformers. */
87 private final E[] updates;
88
89 /** Indicator for forward integration. */
90 private boolean forward;
91
92 /** Extreme time encountered so far. */
93 private E extremeT;
94
95 /** Wrap a {@link FieldODEEventDetector eve,t detector}.
96 * @param field field to which array elements belong
97 * @param rawDetector event detector to wrap
98 * @param filter filter to use
99 */
100 public FieldEventSlopeFilter(final Field<E> field, final T rawDetector, final FilterType filter) {
101 this(field,
102 rawDetector.getMaxCheckInterval(), rawDetector.getMaxIterationCount(),
103 rawDetector.getSolver(), new LocalHandler<>(rawDetector.getHandler()),
104 rawDetector, filter);
105 }
106
107 /** Private constructor with full parameters.
108 * <p>
109 * This constructor is private as users are expected to use the builder
110 * API with the various {@code withXxx()} methods to set up the instance
111 * in a readable manner without using a huge amount of parameters.
112 * </p>
113 * @param field field to which array elements belong
114 * @param maxCheck maximum checking interval (s)
115 * @param maxIter maximum number of iterations in the event time search
116 * @param solver solver to user for locating event
117 * @param handler event handler to call at event occurrences
118 * @param rawDetector event detector to wrap
119 * @param filter filter to use
120 */
121 private FieldEventSlopeFilter(final Field<E> field,
122 final FieldAdaptableInterval<E> maxCheck, final int maxIter,
123 final BracketedRealFieldUnivariateSolver<E> solver,
124 final FieldODEEventHandler<E> handler,
125 final T rawDetector, final FilterType filter) {
126 super(maxCheck, maxIter, solver, handler);
127 this.rawDetector = rawDetector;
128 this.filter = filter;
129 this.transformers = new Transformer[HISTORY_SIZE];
130 this.updates = MathArrays.buildArray(field, HISTORY_SIZE);
131 }
132
133 /** {@inheritDoc} */
134 @Override
135 protected FieldEventSlopeFilter<T, E> create(final FieldAdaptableInterval<E> newMaxCheck, final int newMaxIter,
136 final BracketedRealFieldUnivariateSolver<E> newSolver,
137 final FieldODEEventHandler<E> newHandler) {
138 return new FieldEventSlopeFilter<>(newSolver.getAbsoluteAccuracy().getField(), newMaxCheck, newMaxIter,
139 newSolver, newHandler, rawDetector, filter);
140 }
141
142 /**
143 * Get the wrapped raw detector.
144 * @return the wrapped raw detector
145 */
146 public T getDetector() {
147 return rawDetector;
148 }
149
150 /** {@inheritDoc} */
151 @Override
152 public void init(final FieldODEStateAndDerivative<E> initialState, E finalTime) {
153 super.init(initialState, finalTime);
154
155 // delegate to raw handler
156 rawDetector.init(initialState, finalTime);
157
158 // initialize events triggering logic
159 forward = finalTime.subtract(initialState.getTime()).getReal() >= 0;
160 extremeT = finalTime.getField().getZero().newInstance(forward ? Double.NEGATIVE_INFINITY : Double.POSITIVE_INFINITY);
161 Arrays.fill(transformers, Transformer.UNINITIALIZED);
162 Arrays.fill(updates, extremeT);
163
164 }
165
166 /** {@inheritDoc} */
167 @Override
168 public void reset(final FieldODEStateAndDerivative<E> intermediateState, final E finalTime) {
169 super.reset(intermediateState, finalTime);
170 rawDetector.reset(intermediateState, finalTime);
171 }
172
173 /** {@inheritDoc} */
174 @Override
175 public boolean isForward() {
176 return forward;
177 }
178
179 /** {@inheritDoc} */
180 @Override
181 public E g(final FieldODEStateAndDerivative<E> state) {
182
183 final E rawG = rawDetector.g(state);
184
185 // search which transformer should be applied to g
186 if (isForward()) {
187 final int last = transformers.length - 1;
188 if (extremeT.subtract(state.getTime()).getReal() < 0) {
189 // we are at the forward end of the history
190
191 // check if a new rough root has been crossed
192 final Transformer previous = transformers[last];
193 final Transformer next = filter.selectTransformer(previous, rawG.getReal(), forward);
194 if (next != previous) {
195 // there is a root somewhere between extremeT and t.
196 // the new transformer is valid for t (this is how we have just computed
197 // it above), but it is in fact valid on both sides of the root, so
198 // it was already valid before t and even up to previous time. We store
199 // the switch at extremeT for safety, to ensure the previous transformer
200 // is not applied too close of the root
201 System.arraycopy(updates, 1, updates, 0, last);
202 System.arraycopy(transformers, 1, transformers, 0, last);
203 updates[last] = extremeT;
204 transformers[last] = next;
205 }
206
207 extremeT = state.getTime();
208
209 // apply the transform
210 return next.transformed(rawG);
211
212 } else {
213 // we are in the middle of the history
214
215 // select the transformer
216 for (int i = last; i > 0; --i) {
217 if (updates[i].subtract(state.getTime()).getReal() <= 0) {
218 // apply the transform
219 return transformers[i].transformed(rawG);
220 }
221 }
222
223 return transformers[0].transformed(rawG);
224
225 }
226 } else {
227 if (state.getTime().subtract(extremeT).getReal() < 0) {
228 // we are at the backward end of the history
229
230 // check if a new rough root has been crossed
231 final Transformer previous = transformers[0];
232 final Transformer next = filter.selectTransformer(previous, rawG.getReal(), forward);
233 if (next != previous) {
234 // there is a root somewhere between extremeT and t.
235 // the new transformer is valid for t (this is how we have just computed
236 // it above), but it is in fact valid on both sides of the root, so
237 // it was already valid before t and even up to previous time. We store
238 // the switch at extremeT for safety, to ensure the previous transformer
239 // is not applied too close of the root
240 System.arraycopy(updates, 0, updates, 1, updates.length - 1);
241 System.arraycopy(transformers, 0, transformers, 1, transformers.length - 1);
242 updates[0] = extremeT;
243 transformers[0] = next;
244 }
245
246 extremeT = state.getTime();
247
248 // apply the transform
249 return next.transformed(rawG);
250
251 } else {
252 // we are in the middle of the history
253
254 // select the transformer
255 for (int i = 0; i < updates.length - 1; ++i) {
256 if (state.getTime().subtract(updates[i]).getReal() <= 0) {
257 // apply the transform
258 return transformers[i].transformed(rawG);
259 }
260 }
261
262 return transformers[updates.length - 1].transformed(rawG);
263
264 }
265 }
266
267 }
268
269 /** Local handler.
270 * @param <T> type of the event detector
271 * @param <E> the type of the field elements
272 */
273 private static class LocalHandler<T extends FieldODEEventDetector<E>, E extends CalculusFieldElement<E>>
274 implements FieldODEEventHandler<E> {
275
276 /** Raw handler. */
277 private final FieldODEEventHandler<E> rawHandler;
278
279 /** Simple constructor.
280 * @param rawHandler raw handler
281 */
282 LocalHandler(final FieldODEEventHandler<E> rawHandler) {
283 this.rawHandler = rawHandler;
284 }
285
286 /** {@inheritDoc} */
287 @Override
288 public Action eventOccurred(final FieldODEStateAndDerivative<E> state,
289 final FieldODEEventDetector<E> detector,
290 final boolean increasing) {
291 // delegate to raw handler, fixing increasing status on the fly
292 @SuppressWarnings("unchecked")
293 final FieldEventSlopeFilter<T, E> esf = (FieldEventSlopeFilter<T, E>) detector;
294 return rawHandler.eventOccurred(state, esf, esf.filter.isTriggeredOnIncreasing());
295 }
296
297 /** {@inheritDoc} */
298 @Override
299 public FieldODEState<E> resetState(final FieldODEEventDetector<E> detector,
300 final FieldODEStateAndDerivative<E> state) {
301 // delegate to raw handler
302 @SuppressWarnings("unchecked")
303 final FieldEventSlopeFilter<T, E> esf = (FieldEventSlopeFilter<T, E>) detector;
304 return rawHandler.resetState(esf, state);
305 }
306
307 }
308
309 }