1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23 package org.hipparchus.ode.events;
24
25
26 import java.lang.reflect.Field;
27 import org.hipparchus.analysis.UnivariateFunction;
28 import org.hipparchus.analysis.solvers.BracketedUnivariateSolver;
29 import org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver;
30 import org.hipparchus.exception.MathIllegalArgumentException;
31 import org.hipparchus.exception.MathIllegalStateException;
32 import org.hipparchus.ode.EquationsMapper;
33 import org.hipparchus.ode.ExpandableODE;
34 import org.hipparchus.ode.ODEState;
35 import org.hipparchus.ode.ODEStateAndDerivative;
36 import org.hipparchus.ode.OrdinaryDifferentialEquation;
37 import org.hipparchus.ode.SecondaryODE;
38 import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
39 import org.hipparchus.ode.nonstiff.LutherIntegrator;
40 import org.hipparchus.ode.nonstiff.interpolators.ClassicalRungeKuttaStateInterpolator;
41 import org.hipparchus.ode.sampling.DummyStepInterpolator;
42 import org.hipparchus.ode.sampling.ODEStateInterpolator;
43 import org.junit.jupiter.api.Assertions;
44 import org.junit.jupiter.api.Test;
45 import org.junit.jupiter.params.ParameterizedTest;
46 import org.junit.jupiter.params.provider.ValueSource;
47 import org.mockito.Mockito;
48
49 import static org.junit.jupiter.api.Assertions.assertEquals;
50 import static org.junit.jupiter.api.Assertions.assertFalse;
51 import static org.junit.jupiter.api.Assertions.assertTrue;
52
53 class DetectorBasedEventStateTest {
54
55
56 @Test
57 void testEpochComparisonAtLeastSignificantBit() throws NoSuchFieldException, IllegalAccessException {
58
59 final double eventTime = 17016.237999999998;
60
61
62
63 final ODEStateAndDerivative globalCurrent = new ODEStateAndDerivative(17016.238, new double[2], new double[2]);
64 final ODEStateAndDerivative globalPrevious = new ODEStateAndDerivative(17016.237999999998, new double[2], new double[2]);
65 final ClassicalRungeKuttaStateInterpolator interpolator = new ClassicalRungeKuttaStateInterpolator(true, new double[2][2], globalPrevious, globalCurrent,
66 globalPrevious, globalCurrent, null);
67 final ODEStateAndDerivative interpolatedState = interpolator.getInterpolatedState(eventTime);
68 Assertions.assertEquals(interpolatedState.getTime(), globalCurrent.getTime());
69 Assertions.assertNotEquals(interpolatedState.getTime(), globalPrevious.getTime());
70
71
72
73 final DetectorBasedEventState es = new DetectorBasedEventState(new CloseEventsGenerator(16436.238, 17016.238, 720, 1e-10, 100));
74 final Field pendingEvent = DetectorBasedEventState.class.getDeclaredField("pendingEvent");
75 pendingEvent.setAccessible(true);
76 pendingEvent.set(es, true);
77 final Field pendingEventTime = DetectorBasedEventState.class.getDeclaredField("pendingEventTime");
78 pendingEventTime.setAccessible(true);
79 pendingEventTime.set(es, eventTime);
80 final Field afterG = DetectorBasedEventState.class.getDeclaredField("afterG");
81 afterG.setAccessible(true);
82 afterG.set(es, 0.0);
83
84
85 Assertions.assertNotNull(es.doEvent(interpolatedState));
86 }
87
88 @Test
89 void testDoEventThrowsIfTimeMismatch() throws NoSuchFieldException, IllegalAccessException {
90
91 final ODEEventDetector detector = new DummyDetector();
92 final DetectorBasedEventState eventState = new DetectorBasedEventState(detector);
93 final Field pendingEvent = DetectorBasedEventState.class.getDeclaredField("pendingEvent");
94 pendingEvent.setAccessible(true);
95 pendingEvent.set(eventState, true);
96 final Field pendingEventTime = DetectorBasedEventState.class.getDeclaredField("pendingEventTime");
97 pendingEventTime.setAccessible(true);
98 pendingEventTime.set(eventState, 1.0);
99 final ODEStateAndDerivative state = new ODEStateAndDerivative(1.0001, new double[]{0.0}, new double[]{0.0});
100
101 Assertions.assertThrows(org.hipparchus.exception.MathRuntimeException.class, () -> {
102 eventState.doEvent(state);
103 });
104 }
105
106 @Test
107 void testInitSetsInitialValues() {
108
109 final ODEEventDetector detector = new DummyDetector();
110 final DetectorBasedEventState eventState = new DetectorBasedEventState(detector);
111 final ODEStateAndDerivative s0 = new ODEStateAndDerivative(0.0, new double[1], new double[1]);
112
113 eventState.init(s0, 10.0);
114
115 Assertions.assertNotNull(eventState.getEventDetector());
116 Assertions.assertEquals(detector, eventState.getEventDetector());
117 Assertions.assertEquals(Double.NEGATIVE_INFINITY, eventState.getEventTime());
118 }
119
120 @Test
121 void testEvaluateStepNoEvent() {
122
123 final ODEEventDetector detector = new DummyDetector();
124 final DetectorBasedEventState eventState = new DetectorBasedEventState(detector);
125 final ODEStateAndDerivative s0 = new ODEStateAndDerivative(0.0, new double[1], new double[1]);
126 final ODEStateAndDerivative s1 = new ODEStateAndDerivative(1.0, new double[1], new double[1]);
127 final ODEStateInterpolator interpolator = Mockito.mock(ODEStateInterpolator.class);
128 Mockito.when(interpolator.isForward()).thenReturn(true);
129 Mockito.when(interpolator.getPreviousState()).thenReturn(s0);
130 Mockito.when(interpolator.getCurrentState()).thenReturn(s1);
131
132 eventState.init(s0, 1.0);
133 eventState.reinitializeBegin(interpolator);
134 final boolean hasEvent = eventState.evaluateStep(interpolator);
135
136 Assertions.assertFalse(hasEvent);
137 }
138
139
140 @Test
141 void closeEvents() throws MathIllegalArgumentException, MathIllegalStateException {
142
143 final double r1 = 90.0;
144 final double r2 = 135.0;
145 final double gap = r2 - r1;
146
147 final double tolerance = 0.1;
148 DetectorBasedEventState es = new DetectorBasedEventState(new CloseEventsGenerator(r1, r2, 1.5 * gap, tolerance, 100));
149 EquationsMapper mapper = new ExpandableODE(new OrdinaryDifferentialEquation() {
150 @Override
151 public int getDimension() {
152 return 0;
153 }
154 @Override
155 public double[] computeDerivatives(double t, double[] y) {
156 return new double[0];
157 }
158 }).getMapper();
159
160 double[] a = new double[0];
161 ODEStateAndDerivative osdLongBefore = new ODEStateAndDerivative(r1 - 1.5 * gap, a, a);
162 ODEStateAndDerivative osBefore = new ODEStateAndDerivative(r1 - 0.5 * gap, a, a);
163 ODEStateInterpolator interpolatorA = new DummyStepInterpolator(true,
164 osdLongBefore, osBefore,
165 osdLongBefore, osBefore,
166 mapper);
167 es.reinitializeBegin(interpolatorA);
168 assertFalse(es.evaluateStep(interpolatorA));
169
170 ODEStateAndDerivative osdBetween = new ODEStateAndDerivative(0.5 * (r1 + r2), a, a);
171 ODEStateInterpolator interpolatorB = new DummyStepInterpolator(true,
172 osBefore, osdBetween,
173 osBefore, osdBetween,
174 mapper);
175 assertTrue(es.evaluateStep(interpolatorB));
176 assertEquals(r1, es.getEventTime(), tolerance);
177 ODEStateAndDerivative osdAtEvent = new ODEStateAndDerivative(es.getEventTime(), a, a);
178 es.doEvent(osdAtEvent);
179
180 ODEStateAndDerivative osdAfterSecond = new ODEStateAndDerivative(r2 + 0.4 * gap, a, a);
181 ODEStateInterpolator interpolatorC = new DummyStepInterpolator(true,
182 osdAtEvent, osdAfterSecond,
183 osdAtEvent, osdAfterSecond,
184 mapper);
185 assertTrue(es.evaluateStep(interpolatorC));
186 assertEquals(r2, es.getEventTime(), tolerance);
187
188 }
189
190
191 @Test
192 void testIssue695()
193 throws MathIllegalArgumentException, MathIllegalStateException {
194
195 OrdinaryDifferentialEquation equation = new OrdinaryDifferentialEquation() {
196 @Override
197 public int getDimension() {
198 return 1;
199 }
200 @Override
201 public double[] computeDerivatives(double t, double[] y) {
202 return new double[] { 1.0 };
203 }
204 };
205
206 DormandPrince853Integrator integrator = new DormandPrince853Integrator(0.001, 1000, 1.0e-14, 1.0e-14);
207 integrator.addEventDetector(new ResettingEvent(10.99, 0.1, 1.0e-9, 1000));
208 integrator.addEventDetector(new ResettingEvent(11.01, 0.1, 1.0e-9, 1000));
209 integrator.setInitialStepSize(3.0);
210
211 double target = 30.0;
212 ODEStateAndDerivative finalState =
213 integrator.integrate(equation, new ODEState(0.0, new double[1]), target);
214 assertEquals(target, finalState.getTime(), 1.0e-10);
215 assertEquals(32.0, finalState.getPrimaryState()[0], 1.0e-10);
216
217 }
218
219
220 @Test
221 void testIssue965()
222 throws MathIllegalArgumentException, MathIllegalStateException {
223
224 ExpandableODE equation = new ExpandableODE(new OrdinaryDifferentialEquation() {
225 @Override
226 public int getDimension() {
227 return 1;
228 }
229 @Override
230 public double[] computeDerivatives(double t, double[] y) {
231 return new double[] { 2.0 };
232 }
233 });
234 int index = equation.addSecondaryEquations(new SecondaryODE() {
235 @Override
236 public int getDimension() {
237 return 1;
238 }
239 @Override
240 public double[] computeDerivatives(double t, double[] primary,
241 double[] primaryDot, double[] secondary) {
242 return new double[] { -3.0 };
243 }
244 });
245 assertEquals(1, index);
246
247 DormandPrince853Integrator integrator = new DormandPrince853Integrator(0.001, 1000, 1.0e-14, 1.0e-14);
248 integrator.addEventDetector(new SecondaryStateEvent(index, -3.0, 0.1, 1.0e-9, 1000));
249 integrator.setInitialStepSize(3.0);
250
251 ODEState initialState = new ODEState(0.0,
252 new double[] { 0.0 },
253 new double[][] { { 0.0 } });
254 ODEStateAndDerivative finalState = integrator.integrate(equation, initialState, 30.0);
255 assertEquals( 1.0, finalState.getTime(), 1.0e-10);
256 assertEquals( 2.0, finalState.getPrimaryState()[0], 1.0e-10);
257 assertEquals(-3.0, finalState.getSecondaryState(index)[0], 1.0e-10);
258
259 }
260
261 @ParameterizedTest
262 @ValueSource(booleans = {true, false})
263 void testAdaptableInterval(final boolean isForward) {
264
265 final TestDetector detector = new TestDetector();
266 final DetectorBasedEventState eventState = new DetectorBasedEventState(detector);
267 final ODEStateInterpolator mockedInterpolator = Mockito.mock(ODEStateInterpolator.class);
268 final ODEStateAndDerivative stateAndDerivative1 = getStateAndDerivative(1);
269 final ODEStateAndDerivative stateAndDerivative2 = getStateAndDerivative(-1);
270 if (isForward) {
271 Mockito.when(mockedInterpolator.getCurrentState()).thenReturn(stateAndDerivative1);
272 Mockito.when(mockedInterpolator.getPreviousState()).thenReturn(stateAndDerivative2);
273 } else {
274 Mockito.when(mockedInterpolator.getCurrentState()).thenReturn(stateAndDerivative2);
275 Mockito.when(mockedInterpolator.getPreviousState()).thenReturn(stateAndDerivative1);
276 }
277 Mockito.when(mockedInterpolator.isForward()).thenReturn(isForward);
278
279 eventState.evaluateStep(mockedInterpolator);
280
281 if (isForward) {
282 Assertions.assertEquals(1, detector.triggeredForward);
283 Assertions.assertEquals(0, detector.triggeredBackward);
284 } else {
285 Assertions.assertEquals(0, detector.triggeredForward);
286 Assertions.assertEquals(1, detector.triggeredBackward);
287 }
288 }
289
290 @Test
291 void testEventsCloserThanThreshold()
292 throws MathIllegalArgumentException, MathIllegalStateException {
293
294 OrdinaryDifferentialEquation equation = new OrdinaryDifferentialEquation() {
295
296 @Override
297 public int getDimension() {
298 return 1;
299 }
300
301 @Override
302 public double[] computeDerivatives(double t, double[] y) {
303 return new double[] { 1.0 };
304 }
305 };
306
307 LutherIntegrator integrator = new LutherIntegrator(20.0);
308 CloseEventsGenerator eventsGenerator =
309 new CloseEventsGenerator(9.0 - 1.0 / 128, 9.0 + 1.0 / 128, 1.0, 0.02, 1000);
310 integrator.addEventDetector(eventsGenerator);
311 double tEnd = integrator.integrate(equation, new ODEState(0.0, new double[1]), 100.0).getTime();
312 assertEquals( 2, eventsGenerator.getCount());
313 assertEquals( 9.0 + 1.0 / 128, tEnd, 1.0 / 32.0);
314
315 }
316
317 private static class SecondaryStateEvent implements ODEEventDetector {
318
319 private final AdaptableInterval maxCheck;
320 private final int maxIter;
321 private final BracketingNthOrderBrentSolver solver;
322 private final int index;
323 private final double target;
324
325 public SecondaryStateEvent(final int index, final double target,
326 final double maxCheck, final double threshold, final int maxIter) {
327 this.maxCheck = (s, isForward) -> maxCheck;
328 this.maxIter = maxIter;
329 this.solver = new BracketingNthOrderBrentSolver(0, threshold, 0, 5);
330 this.index = index;
331 this.target = target;
332 }
333
334 @Override
335 public AdaptableInterval getMaxCheckInterval() {
336 return maxCheck;
337 }
338
339 @Override
340 public int getMaxIterationCount() {
341 return maxIter;
342 }
343
344 @Override
345 public BracketedUnivariateSolver<UnivariateFunction> getSolver() {
346 return solver;
347 }
348
349
350 @Override
351 public ODEEventHandler getHandler() {
352 return (state, detector, increasing) -> Action.STOP;
353 }
354
355 @Override
356 public double g(ODEStateAndDerivative s) {
357 return s.getSecondaryState(index)[0] - target;
358 }
359
360 }
361
362 private static class DummyDetector implements ODEEventDetector {
363 @Override
364 public AdaptableInterval getMaxCheckInterval() {
365 return (state, isForward) -> 1.0;
366 }
367 @Override
368 public int getMaxIterationCount() {
369 return 10;
370 }
371 @Override
372 public BracketingNthOrderBrentSolver getSolver() {
373 return new BracketingNthOrderBrentSolver();
374 }
375 @Override
376 public ODEEventHandler getHandler() {
377 return (state, detector, increasing) -> Action.CONTINUE;
378 }
379 @Override
380 public double g(ODEStateAndDerivative state) {
381 return 1.0;
382 }
383 }
384
385 private static class ResettingEvent implements ODEEventDetector {
386
387 private static double lastTriggerTime = Double.NEGATIVE_INFINITY;
388 private final AdaptableInterval maxCheck;
389 private final int maxIter;
390 private final BracketingNthOrderBrentSolver solver;
391 private final double tEvent;
392
393 public ResettingEvent(final double tEvent,
394 final double maxCheck, final double threshold, final int maxIter) {
395 this.maxCheck = (s, isForward) -> maxCheck;
396 this.maxIter = maxIter;
397 this.solver = new BracketingNthOrderBrentSolver(0, threshold, 0, 5);
398 this.tEvent = tEvent;
399 }
400
401 @Override
402 public AdaptableInterval getMaxCheckInterval() {
403 return maxCheck;
404 }
405
406 @Override
407 public int getMaxIterationCount() {
408 return maxIter;
409 }
410
411 @Override
412 public BracketedUnivariateSolver<UnivariateFunction> getSolver() {
413 return solver;
414 }
415
416 @Override
417 public double g(ODEStateAndDerivative s) {
418
419
420
421
422 assertTrue(s.getTime() >= lastTriggerTime,
423 "going backard in time! (" + s.getTime() + " < " + lastTriggerTime + ")");
424 return s.getTime() - tEvent;
425 }
426
427 @Override
428 public ODEEventHandler getHandler() {
429 return new ODEEventHandler() {
430 @Override
431 public Action eventOccurred(ODEStateAndDerivative s, ODEEventDetector detector, boolean increasing) {
432
433 lastTriggerTime = s.getTime();
434 return Action.RESET_STATE;
435 }
436
437 @Override
438 public ODEStateAndDerivative resetState(ODEEventDetector detector, ODEStateAndDerivative s) {
439 double[] y = s.getPrimaryState();
440 y[0] += 1.0;
441 return new ODEStateAndDerivative(s.getTime(), y, s.getPrimaryDerivative());
442 }
443 };
444 }
445
446 }
447
448 private static class CloseEventsGenerator implements ODEEventDetector {
449
450 private final AdaptableInterval maxCheck;
451 private final int maxIter;
452 private final BracketingNthOrderBrentSolver solver;
453 final double r1;
454 final double r2;
455 int count;
456
457 public CloseEventsGenerator(final double r1, final double r2,
458 final double maxCheck, final double threshold, final int maxIter) {
459 this.maxCheck = (s, isForward) -> maxCheck;
460 this.maxIter = maxIter;
461 this.solver = new BracketingNthOrderBrentSolver(0, threshold, 0, 5);
462 this.r1 = r1;
463 this.r2 = r2;
464 this.count = 0;
465 }
466
467 @Override
468 public AdaptableInterval getMaxCheckInterval() {
469 return maxCheck;
470 }
471
472 @Override
473 public int getMaxIterationCount() {
474 return maxIter;
475 }
476
477 @Override
478 public BracketedUnivariateSolver<UnivariateFunction> getSolver() {
479 return solver;
480 }
481
482 @Override
483 public double g(ODEStateAndDerivative s) {
484 return (s.getTime() - r1) * (r2 - s.getTime());
485 }
486
487 @Override
488 public ODEEventHandler getHandler() {
489 return (state, detector, increasing) -> ++count < 2 ? Action.CONTINUE : Action.STOP;
490 }
491
492 public int getCount() {
493 return count;
494 }
495
496 }
497
498 private static ODEStateAndDerivative getStateAndDerivative(final double time) {
499 return new ODEStateAndDerivative(time, new double[] {time}, new double[1]);
500 }
501
502 private static class TestDetector implements ODEEventDetector {
503
504 int triggeredForward = 0;
505 int triggeredBackward = 0;
506
507 @Override
508 public AdaptableInterval getMaxCheckInterval() {
509 return (state, isForward) -> {
510 if (isForward) {
511 triggeredForward++;
512 } else {
513 triggeredBackward++;
514 }
515 return 1.;
516 };
517 }
518
519 @Override
520 public int getMaxIterationCount() {
521 return 10;
522 }
523
524 @Override
525 public BracketedUnivariateSolver<UnivariateFunction> getSolver() {
526 return new BracketingNthOrderBrentSolver();
527 }
528
529 @Override
530 public ODEEventHandler getHandler() {
531 return null;
532 }
533
534 @Override
535 public double g(ODEStateAndDerivative state) {
536 return 0.;
537 }
538 }
539
540 }