View Javadoc
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  package org.hipparchus.ode.events;
23  
24  import org.hipparchus.analysis.UnivariateFunction;
25  import org.hipparchus.analysis.solvers.BracketedUnivariateSolver;
26  import org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver;
27  import org.hipparchus.exception.MathIllegalArgumentException;
28  import org.hipparchus.exception.MathIllegalStateException;
29  import org.hipparchus.ode.ODEIntegrator;
30  import org.hipparchus.ode.ODEState;
31  import org.hipparchus.ode.ODEStateAndDerivative;
32  import org.hipparchus.ode.OrdinaryDifferentialEquation;
33  
34  import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
35  import org.hipparchus.random.RandomGenerator;
36  import org.hipparchus.random.Well19937a;
37  import org.hipparchus.util.FastMath;
38  import org.junit.jupiter.api.Test;
39  import org.junit.jupiter.params.ParameterizedTest;
40  import org.junit.jupiter.params.provider.EnumSource;
41  import org.mockito.Mockito;
42  
43  import static org.junit.jupiter.api.Assertions.assertEquals;
44  import static org.junit.jupiter.api.Assertions.assertFalse;
45  import static org.junit.jupiter.api.Assertions.assertTrue;
46  
47  class EventSlopeFilterTest {
48  
49      @Test
50      void testInit() {
51          // GIVEN
52          final TestDetector detector = new TestDetector();
53          final EventSlopeFilter<?> eventSlopeFilter = new EventSlopeFilter<>(detector,
54                  FilterType.TRIGGER_ONLY_DECREASING_EVENTS);
55          // WHEN
56          eventSlopeFilter.init(Mockito.mock(ODEStateAndDerivative.class), 1);
57          // THEN
58          assertTrue(detector.initialized);
59          assertFalse(detector.resetted);
60      }
61  
62      @ParameterizedTest
63      @EnumSource(FilterType.class)
64      void testReset(final FilterType type) {
65          // GIVEN
66          final TestDetector detector = new TestDetector();
67          final EventSlopeFilter<?> eventSlopeFilter = new EventSlopeFilter<>(detector, type);
68          // WHEN
69          eventSlopeFilter.reset(Mockito.mock(ODEStateAndDerivative.class), 1);
70          // THEN
71          assertTrue(detector.resetted);
72          assertFalse(detector.initialized);
73      }
74  
75      private static class TestDetector extends AbstractODEDetector<TestDetector> {
76          boolean initialized = false;
77          boolean resetted = false;
78  
79          protected TestDetector() {
80              super((state, isForward) -> 0, 1, new BracketingNthOrderBrentSolver(), (s, e, d) -> Action.CONTINUE);
81          }
82  
83          @Override
84          public void init(ODEStateAndDerivative s0, double t) {
85              super.init(s0, t);
86              initialized = true;
87          }
88  
89          @Override
90          public void reset(ODEStateAndDerivative intermediateState, double finalTime) {
91              super.reset(intermediateState, finalTime);
92              resetted = true;
93          }
94  
95          @Override
96          protected TestDetector create(AdaptableInterval newMaxCheck, int newmaxIter, BracketedUnivariateSolver<UnivariateFunction> newSolver, ODEEventHandler newHandler) {
97              return null;
98          }
99  
100         @Override
101         public double g(ODEStateAndDerivative state) {
102             return 1;
103         }
104     }
105 
106     @Test
107     void testHistoryIncreasingForward() {
108 
109         // start point: g > 0
110         testHistory(FilterType.TRIGGER_ONLY_INCREASING_EVENTS,
111                 0.5 * FastMath.PI, 30.5 * FastMath.PI, FastMath.PI, -1);
112 
113         // start point: g = 0
114         testHistory(FilterType.TRIGGER_ONLY_INCREASING_EVENTS,
115                 0, 30.5 * FastMath.PI, FastMath.PI, -1);
116 
117         // start point: g < 0
118         testHistory(FilterType.TRIGGER_ONLY_INCREASING_EVENTS,
119                 1.5 * FastMath.PI, 30.5 * FastMath.PI, FastMath.PI, +1);
120 
121     }
122 
123     @Test
124     void testHistoryIncreasingBackward() {
125 
126         // start point: g > 0
127         testHistory(FilterType.TRIGGER_ONLY_INCREASING_EVENTS,
128                     0.5 * FastMath.PI, -30.5 * FastMath.PI, FastMath.PI, -1);
129 
130         // start point: g = 0
131         testHistory(FilterType.TRIGGER_ONLY_INCREASING_EVENTS,
132                     0, -30.5 * FastMath.PI, FastMath.PI, +1);
133 
134         // start point: g < 0
135         testHistory(FilterType.TRIGGER_ONLY_INCREASING_EVENTS,
136                     1.5 * FastMath.PI, -30.5 * FastMath.PI, FastMath.PI, -1);
137 
138     }
139 
140     @Test
141     void testHistoryDecreasingForward() {
142 
143         // start point: g > 0
144         testHistory(FilterType.TRIGGER_ONLY_DECREASING_EVENTS,
145                     0.5 * FastMath.PI, 30.5 * FastMath.PI, 0, +1);
146 
147         // start point: g = 0
148         testHistory(FilterType.TRIGGER_ONLY_DECREASING_EVENTS,
149                     0, 30.5 * FastMath.PI, 0, +1);
150 
151         // start point: g < 0
152         testHistory(FilterType.TRIGGER_ONLY_DECREASING_EVENTS,
153                     1.5 * FastMath.PI, 30.5 * FastMath.PI, 0, +1);
154 
155     }
156 
157     @Test
158     void testHistoryDecreasingBackward() {
159 
160         // start point: g > 0
161         testHistory(FilterType.TRIGGER_ONLY_DECREASING_EVENTS,
162                     0.5 * FastMath.PI, -30.5 * FastMath.PI, 0, -1);
163 
164         // start point: g = 0
165         testHistory(FilterType.TRIGGER_ONLY_DECREASING_EVENTS,
166                     0, -30.5 * FastMath.PI, 0, -1);
167 
168         // start point: g < 0
169         testHistory(FilterType.TRIGGER_ONLY_DECREASING_EVENTS,
170                     1.5 * FastMath.PI, -30.5 * FastMath.PI, 0, +1);
171 
172     }
173 
174     private void testHistory(FilterType type, double t0, double t1, double refSwitch, double signEven) {
175         Event onlyIncreasing = new Event(Double.POSITIVE_INFINITY, 1.0e-10, 1000, false, true);
176         EventSlopeFilter<Event> eventFilter = new EventSlopeFilter<>(onlyIncreasing, type);
177         eventFilter.init(new ODEStateAndDerivative(t0,
178                                                    new double[] { FastMath.sin(t0),  FastMath.cos(t0) },
179                                                    new double[] { FastMath.cos(t0), -FastMath.sin(t0) }),
180                          t1);
181 
182         // first pass to set up switches history for a long period
183         double h = FastMath.copySign(0.05, t1 - t0);
184         double n = (int) FastMath.floor((t1 - t0) / h);
185         for (int i = 0; i < n; ++i) {
186             double t = t0 + i * h;
187             eventFilter.g(new ODEStateAndDerivative(t,
188                                                     new double[] { FastMath.sin(t),  FastMath.cos(t) },
189                                                     new double[] { FastMath.cos(t), -FastMath.sin(t) }));
190         }
191 
192         // verify old events are preserved, even if randomly accessed
193         RandomGenerator rng = new Well19937a(0xb0e7401265af8cd3l);
194         for (int i = 0; i < 5000; i++) {
195             double t = t0 + (t1 - t0) * rng.nextDouble();
196             double g = eventFilter.g(new ODEStateAndDerivative(t,
197                                                                new double[] { FastMath.sin(t),  FastMath.cos(t) },
198                                                                new double[] { FastMath.cos(t), -FastMath.sin(t) }));
199             int turn = (int) FastMath.floor((t - refSwitch) / (2 * FastMath.PI));
200             if (turn % 2 == 0) {
201                 assertEquals( signEven * FastMath.sin(t), g, 1.0e-10);
202             } else {
203                 assertEquals(-signEven * FastMath.sin(t), g, 1.0e-10);
204             }
205         }
206 
207     }
208 
209     @Test
210     void testIncreasingOnly()
211         throws MathIllegalArgumentException, MathIllegalStateException {
212         ODEIntegrator integrator = new DormandPrince853Integrator(1.0e-3, 100.0, 1e-7, 1e-7);
213         Event allEvents = new Event(0.1, 1.0e-7, 100, true, true);
214         integrator.addEventDetector(allEvents);
215         Event onlyIncreasing = new Event(0.1, 1.0e-7, 100, false, true);
216         integrator.addEventDetector(new EventSlopeFilter<>(onlyIncreasing, FilterType.TRIGGER_ONLY_INCREASING_EVENTS));
217         double t0 = 0.5 * FastMath.PI;
218         double tEnd = 5.5 * FastMath.PI;
219         double[] y = { 0.0, 1.0 };
220         assertEquals(tEnd,
221                             integrator.integrate(new SineCosine(), new ODEState(t0, y), tEnd).getTime(),
222                             1.0e-7);
223 
224         assertEquals(5, allEvents.getEventCount());
225         assertEquals(2, onlyIncreasing.getEventCount());
226 
227     }
228 
229     @Test
230     void testDecreasingOnly()
231         throws MathIllegalArgumentException, MathIllegalStateException {
232         ODEIntegrator integrator = new DormandPrince853Integrator(1.0e-3, 100.0, 1e-7, 1e-7);
233         Event allEvents = new Event(0.1, 1.0e-7, 1000, true, true);
234         integrator.addEventDetector(allEvents);
235         Event onlyDecreasing = new Event(0.1, 1.0e-7, 1000, true, false);
236         integrator.addEventDetector(new EventSlopeFilter<>(onlyDecreasing, FilterType.TRIGGER_ONLY_DECREASING_EVENTS));
237         double t0 = 0.5 * FastMath.PI;
238         double tEnd = 5.5 * FastMath.PI;
239         double[] y = { 0.0, 1.0 };
240         assertEquals(tEnd,
241                             integrator.integrate(new SineCosine(), new ODEState(t0, y), tEnd).getTime(),
242                             1.0e-7);
243 
244         assertEquals(5, allEvents.getEventCount());
245         assertEquals(3, onlyDecreasing.getEventCount());
246 
247     }
248 
249     @Test
250     void testTwoOppositeFilters()
251         throws MathIllegalArgumentException, MathIllegalStateException {
252         ODEIntegrator integrator = new DormandPrince853Integrator(1.0e-3, 100.0, 1e-7, 1e-7);
253         Event allEvents = new Event(0.1, 1.0e-7, 100, true, true);
254         integrator.addEventDetector(allEvents);
255         Event onlyIncreasing = new Event(0.1, 1.0e-7, 100, false, true);
256         integrator.addEventDetector(new EventSlopeFilter<>(onlyIncreasing, FilterType.TRIGGER_ONLY_INCREASING_EVENTS));
257         Event onlyDecreasing = new Event(0.1, 1.0e-7, 100, true, false);
258         integrator.addEventDetector(new EventSlopeFilter<>(onlyDecreasing, FilterType.TRIGGER_ONLY_DECREASING_EVENTS));
259         double t0 = 0.5 * FastMath.PI;
260         double tEnd = 5.5 * FastMath.PI;
261         double[] y = { 0.0, 1.0 };
262         assertEquals(tEnd,
263                             integrator.integrate(new SineCosine(), new ODEState(t0, y), tEnd).getTime(),
264                             1.0e-7);
265 
266         assertEquals(5, allEvents.getEventCount());
267         assertEquals(2, onlyIncreasing.getEventCount());
268         assertEquals(3, onlyDecreasing.getEventCount());
269 
270     }
271 
272     private static class SineCosine implements OrdinaryDifferentialEquation {
273         public int getDimension() {
274             return 2;
275         }
276 
277         public double[] computeDerivatives(double t, double[] y) {
278             return new double[] { y[1], -y[0] };
279         }
280     }
281 
282     /** State events for this unit test. */
283     protected static class Event implements ODEEventDetector {
284 
285         private final AdaptableInterval                             maxCheck;
286         private final int                                           maxIter;
287         private final BracketedUnivariateSolver<UnivariateFunction> solver;
288         private final boolean                                       expectDecreasing;
289         private final boolean                                       expectIncreasing;
290         private int                                                 eventCount;
291 
292         public Event(final double maxCheck, final double threshold, final int maxIter,
293                      boolean expectDecreasing, boolean expectIncreasing) {
294             this.maxCheck         = (s, isForward) -> maxCheck;
295             this.maxIter          = maxIter;
296             this.solver           = new BracketingNthOrderBrentSolver(0, threshold, 0, 5);
297             this.expectDecreasing = expectDecreasing;
298             this.expectIncreasing = expectIncreasing;
299         }
300 
301         public AdaptableInterval getMaxCheckInterval() {
302             return maxCheck;
303         }
304 
305         public int getMaxIterationCount() {
306             return maxIter;
307         }
308 
309         public BracketedUnivariateSolver<UnivariateFunction> getSolver() {
310             return solver;
311         }
312 
313         public int getEventCount() {
314             return eventCount;
315         }
316 
317         @Override
318         public void init(ODEStateAndDerivative s0, double t) {
319             eventCount = 0;
320         }
321 
322         public double g(ODEStateAndDerivative s) {
323             return s.getPrimaryState()[0];
324         }
325 
326         public ODEEventHandler getHandler() {
327             return (ODEStateAndDerivative s, ODEEventDetector detector, boolean increasing) -> {
328                 if (increasing) {
329                     assertTrue(expectIncreasing);
330                 } else {
331                     assertTrue(expectDecreasing);
332                 }
333                 eventCount++;
334                 return Action.RESET_STATE;
335             };
336         }
337 
338     }
339 
340 }