1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22 package org.hipparchus.ode.events;
23
24 import org.hipparchus.CalculusFieldElement;
25 import org.hipparchus.Field;
26 import org.hipparchus.analysis.UnivariateFunction;
27 import org.hipparchus.analysis.solvers.BracketedRealFieldUnivariateSolver;
28 import org.hipparchus.analysis.solvers.BracketedUnivariateSolver;
29 import org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver;
30 import org.hipparchus.analysis.solvers.FieldBracketingNthOrderBrentSolver;
31 import org.hipparchus.exception.MathIllegalArgumentException;
32 import org.hipparchus.exception.MathIllegalStateException;
33 import org.hipparchus.ode.FieldExpandableODE;
34 import org.hipparchus.ode.FieldODEIntegrator;
35 import org.hipparchus.ode.FieldODEState;
36 import org.hipparchus.ode.FieldODEStateAndDerivative;
37 import org.hipparchus.ode.FieldOrdinaryDifferentialEquation;
38 import org.hipparchus.ode.ODEIntegrator;
39 import org.hipparchus.ode.ODEState;
40 import org.hipparchus.ode.ODEStateAndDerivative;
41 import org.hipparchus.ode.OrdinaryDifferentialEquation;
42 import org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator;
43 import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
44 import org.hipparchus.random.RandomGenerator;
45 import org.hipparchus.random.Well19937a;
46 import org.hipparchus.util.Binary64Field;
47 import org.hipparchus.util.FastMath;
48 import org.hipparchus.util.MathArrays;
49 import org.junit.Assert;
50 import org.junit.Test;
51
52 public class EventSlopeFilterTest {
53
54 @Test
55 public void testHistoryIncreasingForward() {
56
57
58 testHistory(FilterType.TRIGGER_ONLY_INCREASING_EVENTS,
59 0.5 * FastMath.PI, 30.5 * FastMath.PI, FastMath.PI, -1);
60
61
62 testHistory(FilterType.TRIGGER_ONLY_INCREASING_EVENTS,
63 0, 30.5 * FastMath.PI, FastMath.PI, -1);
64
65
66 testHistory(FilterType.TRIGGER_ONLY_INCREASING_EVENTS,
67 1.5 * FastMath.PI, 30.5 * FastMath.PI, FastMath.PI, +1);
68
69 }
70
71 @Test
72 public void testHistoryIncreasingForwardField() {
73
74
75 testHistoryField(Binary64Field.getInstance(), FilterType.TRIGGER_ONLY_INCREASING_EVENTS,
76 0.5 * FastMath.PI, 30.5 * FastMath.PI, FastMath.PI, -1);
77
78
79 testHistoryField(Binary64Field.getInstance(), FilterType.TRIGGER_ONLY_INCREASING_EVENTS,
80 0, 30.5 * FastMath.PI, FastMath.PI, -1);
81
82
83 testHistoryField(Binary64Field.getInstance(), FilterType.TRIGGER_ONLY_INCREASING_EVENTS,
84 1.5 * FastMath.PI, 30.5 * FastMath.PI, FastMath.PI, +1);
85
86 }
87
88 @Test
89 public void testHistoryIncreasingBackward() {
90
91
92 testHistory(FilterType.TRIGGER_ONLY_INCREASING_EVENTS,
93 0.5 * FastMath.PI, -30.5 * FastMath.PI, FastMath.PI, -1);
94
95
96 testHistory(FilterType.TRIGGER_ONLY_INCREASING_EVENTS,
97 0, -30.5 * FastMath.PI, FastMath.PI, +1);
98
99
100 testHistory(FilterType.TRIGGER_ONLY_INCREASING_EVENTS,
101 1.5 * FastMath.PI, -30.5 * FastMath.PI, FastMath.PI, -1);
102
103 }
104
105 @Test
106 public void testHistoryIncreasingBackwardField() {
107
108
109 testHistoryField(Binary64Field.getInstance(), FilterType.TRIGGER_ONLY_INCREASING_EVENTS,
110 0.5 * FastMath.PI, -30.5 * FastMath.PI, FastMath.PI, -1);
111
112
113 testHistoryField(Binary64Field.getInstance(), FilterType.TRIGGER_ONLY_INCREASING_EVENTS,
114 0, -30.5 * FastMath.PI, FastMath.PI, +1);
115
116
117 testHistoryField(Binary64Field.getInstance(), FilterType.TRIGGER_ONLY_INCREASING_EVENTS,
118 1.5 * FastMath.PI, -30.5 * FastMath.PI, FastMath.PI, -1);
119
120 }
121
122 @Test
123 public void testHistoryDecreasingForward() {
124
125
126 testHistory(FilterType.TRIGGER_ONLY_DECREASING_EVENTS,
127 0.5 * FastMath.PI, 30.5 * FastMath.PI, 0, +1);
128
129
130 testHistory(FilterType.TRIGGER_ONLY_DECREASING_EVENTS,
131 0, 30.5 * FastMath.PI, 0, +1);
132
133
134 testHistory(FilterType.TRIGGER_ONLY_DECREASING_EVENTS,
135 1.5 * FastMath.PI, 30.5 * FastMath.PI, 0, +1);
136
137 }
138
139 @Test
140 public void testHistoryDecreasingForwardField() {
141
142
143 testHistoryField(Binary64Field.getInstance(), FilterType.TRIGGER_ONLY_DECREASING_EVENTS,
144 0.5 * FastMath.PI, 30.5 * FastMath.PI, 0, +1);
145
146
147 testHistoryField(Binary64Field.getInstance(), FilterType.TRIGGER_ONLY_DECREASING_EVENTS,
148 0, 30.5 * FastMath.PI, 0, +1);
149
150
151 testHistoryField(Binary64Field.getInstance(), FilterType.TRIGGER_ONLY_DECREASING_EVENTS,
152 1.5 * FastMath.PI, 30.5 * FastMath.PI, 0, +1);
153
154 }
155
156 @Test
157 public void testHistoryDecreasingBackward() {
158
159
160 testHistory(FilterType.TRIGGER_ONLY_DECREASING_EVENTS,
161 0.5 * FastMath.PI, -30.5 * FastMath.PI, 0, -1);
162
163
164 testHistory(FilterType.TRIGGER_ONLY_DECREASING_EVENTS,
165 0, -30.5 * FastMath.PI, 0, -1);
166
167
168 testHistory(FilterType.TRIGGER_ONLY_DECREASING_EVENTS,
169 1.5 * FastMath.PI, -30.5 * FastMath.PI, 0, +1);
170
171 }
172
173 @Test
174 public void testHistoryDecreasingBackwardField() {
175
176
177 testHistoryField(Binary64Field.getInstance(), FilterType.TRIGGER_ONLY_DECREASING_EVENTS,
178 0.5 * FastMath.PI, -30.5 * FastMath.PI, 0, -1);
179
180
181 testHistoryField(Binary64Field.getInstance(), FilterType.TRIGGER_ONLY_DECREASING_EVENTS,
182 0, -30.5 * FastMath.PI, 0, -1);
183
184
185 testHistoryField(Binary64Field.getInstance(), FilterType.TRIGGER_ONLY_DECREASING_EVENTS,
186 1.5 * FastMath.PI, -30.5 * FastMath.PI, 0, +1);
187
188 }
189
190 private void testHistory(FilterType type, double t0, double t1, double refSwitch, double signEven) {
191 Event onlyIncreasing = new Event(Double.POSITIVE_INFINITY, 1.0e-10, 1000, false, true);
192 EventSlopeFilter<Event> eventFilter = new EventSlopeFilter<>(onlyIncreasing, type);
193 eventFilter.init(new ODEStateAndDerivative(t0,
194 new double[] { FastMath.sin(t0), FastMath.cos(t0) },
195 new double[] { FastMath.cos(t0), -FastMath.sin(t0) }),
196 t1);
197
198
199 double h = FastMath.copySign(0.05, t1 - t0);
200 double n = (int) FastMath.floor((t1 - t0) / h);
201 for (int i = 0; i < n; ++i) {
202 double t = t0 + i * h;
203 eventFilter.g(new ODEStateAndDerivative(t,
204 new double[] { FastMath.sin(t), FastMath.cos(t) },
205 new double[] { FastMath.cos(t), -FastMath.sin(t) }));
206 }
207
208
209 RandomGenerator rng = new Well19937a(0xb0e7401265af8cd3l);
210 for (int i = 0; i < 5000; i++) {
211 double t = t0 + (t1 - t0) * rng.nextDouble();
212 double g = eventFilter.g(new ODEStateAndDerivative(t,
213 new double[] { FastMath.sin(t), FastMath.cos(t) },
214 new double[] { FastMath.cos(t), -FastMath.sin(t) }));
215 int turn = (int) FastMath.floor((t - refSwitch) / (2 * FastMath.PI));
216 if (turn % 2 == 0) {
217 Assert.assertEquals( signEven * FastMath.sin(t), g, 1.0e-10);
218 } else {
219 Assert.assertEquals(-signEven * FastMath.sin(t), g, 1.0e-10);
220 }
221 }
222
223 }
224
225 private <T extends CalculusFieldElement<T>> void testHistoryField(Field<T> field, FilterType type, double t0, double t1, double refSwitch, double signEven) {
226 FieldEvent<T> onlyIncreasing = new FieldEvent<>(Double.POSITIVE_INFINITY,
227 field.getZero().newInstance(1.0e-10), 1000, false, true);
228 FieldEventSlopeFilter<FieldEvent<T>, T> eventFilter = new FieldEventSlopeFilter<>(field, onlyIncreasing, type);
229 eventFilter.init(buildStateAndDerivative(field, t0), field.getZero().add(t1));
230
231
232 double h = FastMath.copySign(0.05, t1 - t0);
233 double n = (int) FastMath.floor((t1 - t0) / h);
234 for (int i = 0; i < n; ++i) {
235 double t = t0 + i * h;
236 eventFilter.g(buildStateAndDerivative(field, t));
237 }
238
239
240 RandomGenerator rng = new Well19937a(0xb0e7401265af8cd3l);
241 for (int i = 0; i < 5000; i++) {
242 double t = t0 + (t1 - t0) * rng.nextDouble();
243 double g = eventFilter.g(buildStateAndDerivative(field, t)).getReal();
244 int turn = (int) FastMath.floor((t - refSwitch) / (2 * FastMath.PI));
245 if (turn % 2 == 0) {
246 Assert.assertEquals( signEven * FastMath.sin(t), g, 1.0e-10);
247 } else {
248 Assert.assertEquals(-signEven * FastMath.sin(t), g, 1.0e-10);
249 }
250 }
251
252 }
253
254 private <T extends CalculusFieldElement<T>> FieldODEStateAndDerivative<T> buildStateAndDerivative(Field<T> field, double t0) {
255 T t0F = field.getZero().add(t0);
256 T[] y0F = MathArrays.buildArray(field, 2);
257 y0F[0] = FastMath.sin(t0F);
258 y0F[1] = FastMath.cos(t0F);
259 T[] y0DotF = MathArrays.buildArray(field, 2);
260 y0DotF[0] = FastMath.cos(t0F);
261 y0DotF[1] = FastMath.sin(t0F).negate();
262 return new FieldODEStateAndDerivative<>(t0F, y0F, y0DotF);
263 }
264
265 @Test
266 public void testIncreasingOnly()
267 throws MathIllegalArgumentException, MathIllegalStateException {
268 ODEIntegrator integrator = new DormandPrince853Integrator(1.0e-3, 100.0, 1e-7, 1e-7);
269 Event allEvents = new Event(0.1, 1.0e-7, 100, true, true);
270 integrator.addEventDetector(allEvents);
271 Event onlyIncreasing = new Event(0.1, 1.0e-7, 100, false, true);
272 integrator.addEventDetector(new EventSlopeFilter<>(onlyIncreasing, FilterType.TRIGGER_ONLY_INCREASING_EVENTS));
273 double t0 = 0.5 * FastMath.PI;
274 double tEnd = 5.5 * FastMath.PI;
275 double[] y = { 0.0, 1.0 };
276 Assert.assertEquals(tEnd,
277 integrator.integrate(new SineCosine(), new ODEState(t0, y), tEnd).getTime(),
278 1.0e-7);
279
280 Assert.assertEquals(5, allEvents.getEventCount());
281 Assert.assertEquals(2, onlyIncreasing.getEventCount());
282
283 }
284
285 @Test
286 public void testIncreasingOnlyField() {
287 doTestIncreasingOnlyField(Binary64Field.getInstance());
288 }
289
290 private <T extends CalculusFieldElement<T>> void doTestIncreasingOnlyField(Field<T> field) {
291 FieldODEIntegrator<T> integrator = new DormandPrince853FieldIntegrator<>(field, 1.0e-3, 100.0, 1e-7, 1e-7);
292 FieldEvent<T> allEvents = new FieldEvent<>(0.1, field.getZero().newInstance(1.0e-7), 100, true, true);
293 integrator.addEventDetector(allEvents);
294 FieldEvent<T> onlyIncreasing = new FieldEvent<>(0.1, field.getZero().newInstance(1.0e-7), 100, false, true);
295 integrator.addEventDetector(new FieldEventSlopeFilter<>(field, onlyIncreasing,
296 FilterType.TRIGGER_ONLY_INCREASING_EVENTS));
297 T t0 = field.getZero().add(0.5 * FastMath.PI);
298 T tEnd = field.getZero().add(5.5 * FastMath.PI);
299 T[] y = MathArrays.buildArray(field, 2);
300 y[0] = field.getZero();
301 y[1] = field.getOne();
302 Assert.assertEquals(tEnd.getReal(),
303 integrator.integrate(new FieldExpandableODE<>(new FieldSineCosine<T>()),
304 new FieldODEState<>(t0, y), tEnd).getTime().getReal(),
305 1.0e-7);
306
307 Assert.assertEquals(5, allEvents.getEventCount());
308 Assert.assertEquals(2, onlyIncreasing.getEventCount());
309
310 }
311
312 @Test
313 public void testDecreasingOnly()
314 throws MathIllegalArgumentException, MathIllegalStateException {
315 ODEIntegrator integrator = new DormandPrince853Integrator(1.0e-3, 100.0, 1e-7, 1e-7);
316 Event allEvents = new Event(0.1, 1.0e-7, 1000, true, true);
317 integrator.addEventDetector(allEvents);
318 Event onlyDecreasing = new Event(0.1, 1.0e-7, 1000, true, false);
319 integrator.addEventDetector(new EventSlopeFilter<>(onlyDecreasing, FilterType.TRIGGER_ONLY_DECREASING_EVENTS));
320 double t0 = 0.5 * FastMath.PI;
321 double tEnd = 5.5 * FastMath.PI;
322 double[] y = { 0.0, 1.0 };
323 Assert.assertEquals(tEnd,
324 integrator.integrate(new SineCosine(), new ODEState(t0, y), tEnd).getTime(),
325 1.0e-7);
326
327 Assert.assertEquals(5, allEvents.getEventCount());
328 Assert.assertEquals(3, onlyDecreasing.getEventCount());
329
330 }
331
332 @Test
333 public void testDecreasingOnlyField() {
334 doTestDecreasingOnlyField(Binary64Field.getInstance());
335 }
336
337 private <T extends CalculusFieldElement<T>> void doTestDecreasingOnlyField(Field<T> field) {
338 FieldODEIntegrator<T> integrator = new DormandPrince853FieldIntegrator<>(field, 1.0e-3, 100.0, 1e-7, 1e-7);
339 FieldEvent<T> allEvents = new FieldEvent<>(0.1, field.getZero().newInstance(1.0e-7), 100, true, true);
340 integrator.addEventDetector(allEvents);
341 FieldEvent<T> onlyDecreasing = new FieldEvent<>(0.1, field.getZero().newInstance(1.0e-7), 100, true, false);
342 integrator.addEventDetector(new FieldEventSlopeFilter<>(field, onlyDecreasing,
343 FilterType.TRIGGER_ONLY_DECREASING_EVENTS));
344 T t0 = field.getZero().add(0.5 * FastMath.PI);
345 T tEnd = field.getZero().add(5.5 * FastMath.PI);
346 T[] y = MathArrays.buildArray(field, 2);
347 y[0] = field.getZero();
348 y[1] = field.getOne();
349 Assert.assertEquals(tEnd.getReal(),
350 integrator.integrate(new FieldExpandableODE<>(new FieldSineCosine<T>()),
351 new FieldODEState<>(t0, y), tEnd).getTime().getReal(),
352 1.0e-7);
353
354 Assert.assertEquals(5, allEvents.getEventCount());
355 Assert.assertEquals(3, onlyDecreasing.getEventCount());
356
357 }
358
359 @Test
360 public void testTwoOppositeFilters()
361 throws MathIllegalArgumentException, MathIllegalStateException {
362 ODEIntegrator integrator = new DormandPrince853Integrator(1.0e-3, 100.0, 1e-7, 1e-7);
363 Event allEvents = new Event(0.1, 1.0e-7, 100, true, true);
364 integrator.addEventDetector(allEvents);
365 Event onlyIncreasing = new Event(0.1, 1.0e-7, 100, false, true);
366 integrator.addEventDetector(new EventSlopeFilter<>(onlyIncreasing, FilterType.TRIGGER_ONLY_INCREASING_EVENTS));
367 Event onlyDecreasing = new Event(0.1, 1.0e-7, 100, true, false);
368 integrator.addEventDetector(new EventSlopeFilter<>(onlyDecreasing, FilterType.TRIGGER_ONLY_DECREASING_EVENTS));
369 double t0 = 0.5 * FastMath.PI;
370 double tEnd = 5.5 * FastMath.PI;
371 double[] y = { 0.0, 1.0 };
372 Assert.assertEquals(tEnd,
373 integrator.integrate(new SineCosine(), new ODEState(t0, y), tEnd).getTime(),
374 1.0e-7);
375
376 Assert.assertEquals(5, allEvents.getEventCount());
377 Assert.assertEquals(2, onlyIncreasing.getEventCount());
378 Assert.assertEquals(3, onlyDecreasing.getEventCount());
379
380 }
381
382 @Test
383 public void testTwoOppositeFiltersField() {
384 doestTwoOppositeFiltersField(Binary64Field.getInstance());
385 }
386
387 private <T extends CalculusFieldElement<T>> void doestTwoOppositeFiltersField(Field<T> field)
388 throws MathIllegalArgumentException, MathIllegalStateException {
389 FieldODEIntegrator<T> integrator = new DormandPrince853FieldIntegrator<>(field, 1.0e-3, 100.0, 1e-7, 1e-7);
390 FieldEvent<T> allEvents = new FieldEvent<>(0.1, field.getZero().newInstance(1.0e-7), 100, true, true);
391 integrator.addEventDetector(allEvents);
392 FieldEvent<T> onlyIncreasing = new FieldEvent<>(0.1, field.getZero().newInstance(1.0e-7), 100, false, true);
393 integrator.addEventDetector(new FieldEventSlopeFilter<>(field, onlyIncreasing,
394 FilterType.TRIGGER_ONLY_INCREASING_EVENTS));
395 FieldEvent<T> onlyDecreasing = new FieldEvent<>(0.1, field.getZero().newInstance(1.0e-7), 100, true, false);
396 integrator.addEventDetector(new FieldEventSlopeFilter<>(field, onlyDecreasing,
397 FilterType.TRIGGER_ONLY_DECREASING_EVENTS));
398 T t0 = field.getZero().add(0.5 * FastMath.PI);
399 T tEnd = field.getZero().add(5.5 * FastMath.PI);
400 T[] y = MathArrays.buildArray(field, 2);
401 y[0] = field.getZero();
402 y[1] = field.getOne();
403 Assert.assertEquals(tEnd.getReal(),
404 integrator.integrate(new FieldExpandableODE<>(new FieldSineCosine<T>()),
405 new FieldODEState<>(t0, y), tEnd).getTime().getReal(),
406 1.0e-7);
407
408 Assert.assertEquals(5, allEvents.getEventCount());
409 Assert.assertEquals(2, onlyIncreasing.getEventCount());
410 Assert.assertEquals(3, onlyDecreasing.getEventCount());
411
412 }
413
414 private static class SineCosine implements OrdinaryDifferentialEquation {
415 public int getDimension() {
416 return 2;
417 }
418
419 public double[] computeDerivatives(double t, double[] y) {
420 return new double[] { y[1], -y[0] };
421 }
422 }
423
424 private static class FieldSineCosine<T extends CalculusFieldElement<T>> implements FieldOrdinaryDifferentialEquation<T> {
425 public int getDimension() {
426 return 2;
427 }
428
429 public T[] computeDerivatives(T t, T[] y) {
430 final T[] yDot = MathArrays.buildArray(t.getField(), getDimension());
431 yDot[0] = y[1];
432 yDot[1] = y[0].negate();
433 return yDot;
434 }
435 }
436
437
438 protected static class Event implements ODEEventDetector {
439
440 private final AdaptableInterval maxCheck;
441 private final int maxIter;
442 private final BracketedUnivariateSolver<UnivariateFunction> solver;
443 private final boolean expectDecreasing;
444 private final boolean expectIncreasing;
445 private int eventCount;
446
447 public Event(final double maxCheck, final double threshold, final int maxIter,
448 boolean expectDecreasing, boolean expectIncreasing) {
449 this.maxCheck = s -> maxCheck;
450 this.maxIter = maxIter;
451 this.solver = new BracketingNthOrderBrentSolver(0, threshold, 0, 5);
452 this.expectDecreasing = expectDecreasing;
453 this.expectIncreasing = expectIncreasing;
454 }
455
456 public AdaptableInterval getMaxCheckInterval() {
457 return maxCheck;
458 }
459
460 public int getMaxIterationCount() {
461 return maxIter;
462 }
463
464 public BracketedUnivariateSolver<UnivariateFunction> getSolver() {
465 return solver;
466 }
467
468 public int getEventCount() {
469 return eventCount;
470 }
471
472 public void init(ODEStateAndDerivative s0, double t) {
473 eventCount = 0;
474 }
475
476 public double g(ODEStateAndDerivative s) {
477 return s.getPrimaryState()[0];
478 }
479
480 public ODEEventHandler getHandler() {
481 return (ODEStateAndDerivative s, ODEEventDetector detector, boolean increasing) -> {
482 if (increasing) {
483 Assert.assertTrue(expectIncreasing);
484 } else {
485 Assert.assertTrue(expectDecreasing);
486 }
487 eventCount++;
488 return Action.RESET_STATE;
489 };
490 }
491
492 }
493
494
495 protected static class FieldEvent<T extends CalculusFieldElement<T>> implements FieldODEEventDetector<T> {
496
497 private final FieldAdaptableInterval<T> maxCheck;
498 private final int maxIter;
499 private final BracketedRealFieldUnivariateSolver<T> solver;
500 private final boolean expectDecreasing;
501 private final boolean expectIncreasing;
502 private int eventCount;
503
504 public FieldEvent(final double maxCheck, final T threshold, final int maxIter,
505 boolean expectDecreasing, boolean expectIncreasing) {
506 this.maxCheck = s -> maxCheck;
507 this.maxIter = maxIter;
508 this.solver = new FieldBracketingNthOrderBrentSolver<>(threshold.getField().getZero(),
509 threshold,
510 threshold.getField().getZero(),
511 5);
512 this.expectDecreasing = expectDecreasing;
513 this.expectIncreasing = expectIncreasing;
514 }
515
516 public FieldAdaptableInterval<T> getMaxCheckInterval() {
517 return maxCheck;
518 }
519
520 public int getMaxIterationCount() {
521 return maxIter;
522 }
523
524 public BracketedRealFieldUnivariateSolver<T> getSolver() {
525 return solver;
526 }
527
528 public int getEventCount() {
529 return eventCount;
530 }
531
532 public void init(FieldODEStateAndDerivative<T> s0, T t) {
533 eventCount = 0;
534 }
535
536 public T g(FieldODEStateAndDerivative<T> s) {
537 return s.getPrimaryState()[0];
538 }
539
540 public FieldODEEventHandler<T> getHandler() {
541 return (FieldODEStateAndDerivative<T> s, FieldODEEventDetector<T> detector, boolean increasing) -> {
542 if (increasing) {
543 Assert.assertTrue(expectIncreasing);
544 } else {
545 Assert.assertTrue(expectDecreasing);
546 }
547 eventCount++;
548 return Action.RESET_STATE;
549 };
550 }
551
552 }
553 }