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