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.CalculusFieldElement;
25  import org.hipparchus.Field;
26  import org.hipparchus.analysis.solvers.BracketedRealFieldUnivariateSolver;
27  import org.hipparchus.analysis.solvers.FieldBracketingNthOrderBrentSolver;
28  import org.hipparchus.exception.MathIllegalArgumentException;
29  import org.hipparchus.exception.MathIllegalStateException;
30  import org.hipparchus.ode.FieldExpandableODE;
31  import org.hipparchus.ode.FieldODEIntegrator;
32  import org.hipparchus.ode.FieldODEState;
33  import org.hipparchus.ode.FieldODEStateAndDerivative;
34  import org.hipparchus.ode.FieldOrdinaryDifferentialEquation;
35  import org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator;
36  import org.hipparchus.random.RandomGenerator;
37  import org.hipparchus.random.Well19937a;
38  import org.hipparchus.util.Binary64Field;
39  import org.hipparchus.util.FastMath;
40  import org.hipparchus.util.MathArrays;
41  import org.junit.jupiter.api.Test;
42  
43  import static org.junit.jupiter.api.Assertions.assertEquals;
44  import static org.junit.jupiter.api.Assertions.assertTrue;
45  
46  class FieldEventSlopeFilterTest {
47  
48      @Test
49      void testHistoryIncreasingForwardField() {
50  
51          // start point: g > 0
52          testHistoryField(Binary64Field.getInstance(), FilterType.TRIGGER_ONLY_INCREASING_EVENTS,
53                           0.5 * FastMath.PI, 30.5 * FastMath.PI, FastMath.PI, -1);
54  
55          // start point: g = 0
56          testHistoryField(Binary64Field.getInstance(), FilterType.TRIGGER_ONLY_INCREASING_EVENTS,
57                           0, 30.5 * FastMath.PI, FastMath.PI, -1);
58  
59          // start point: g < 0
60          testHistoryField(Binary64Field.getInstance(), FilterType.TRIGGER_ONLY_INCREASING_EVENTS,
61                           1.5 * FastMath.PI, 30.5 * FastMath.PI, FastMath.PI, +1);
62  
63      }
64  
65      @Test
66      void testHistoryIncreasingBackwardField() {
67  
68          // start point: g > 0
69          testHistoryField(Binary64Field.getInstance(), FilterType.TRIGGER_ONLY_INCREASING_EVENTS,
70                           0.5 * FastMath.PI, -30.5 * FastMath.PI, FastMath.PI, -1);
71  
72          // start point: g = 0
73          testHistoryField(Binary64Field.getInstance(), FilterType.TRIGGER_ONLY_INCREASING_EVENTS,
74                           0, -30.5 * FastMath.PI, FastMath.PI, +1);
75  
76          // start point: g < 0
77          testHistoryField(Binary64Field.getInstance(), FilterType.TRIGGER_ONLY_INCREASING_EVENTS,
78                           1.5 * FastMath.PI, -30.5 * FastMath.PI, FastMath.PI, -1);
79  
80      }
81  
82      @Test
83      void testHistoryDecreasingForwardField() {
84  
85          // start point: g > 0
86          testHistoryField(Binary64Field.getInstance(), FilterType.TRIGGER_ONLY_DECREASING_EVENTS,
87                           0.5 * FastMath.PI, 30.5 * FastMath.PI, 0, +1);
88  
89          // start point: g = 0
90          testHistoryField(Binary64Field.getInstance(), FilterType.TRIGGER_ONLY_DECREASING_EVENTS,
91                           0, 30.5 * FastMath.PI, 0, +1);
92  
93          // start point: g < 0
94          testHistoryField(Binary64Field.getInstance(), FilterType.TRIGGER_ONLY_DECREASING_EVENTS,
95                           1.5 * FastMath.PI, 30.5 * FastMath.PI, 0, +1);
96  
97      }
98  
99      @Test
100     void testHistoryDecreasingBackwardField() {
101 
102         // start point: g > 0
103         testHistoryField(Binary64Field.getInstance(), FilterType.TRIGGER_ONLY_DECREASING_EVENTS,
104                          0.5 * FastMath.PI, -30.5 * FastMath.PI, 0, -1);
105 
106         // start point: g = 0
107         testHistoryField(Binary64Field.getInstance(), FilterType.TRIGGER_ONLY_DECREASING_EVENTS,
108                          0, -30.5 * FastMath.PI, 0, -1);
109 
110         // start point: g < 0
111         testHistoryField(Binary64Field.getInstance(), FilterType.TRIGGER_ONLY_DECREASING_EVENTS,
112                          1.5 * FastMath.PI, -30.5 * FastMath.PI, 0, +1);
113 
114     }
115 
116     private <T extends CalculusFieldElement<T>> void testHistoryField(Field<T> field, FilterType type, double t0, double t1, double refSwitch, double signEven) {
117         FieldEvent<T> onlyIncreasing = new FieldEvent<>(Double.POSITIVE_INFINITY,
118                                                         field.getZero().newInstance(1.0e-10), 1000, false, true);
119         FieldEventSlopeFilter<FieldEvent<T>, T> eventFilter = new FieldEventSlopeFilter<>(field, onlyIncreasing, type);
120         eventFilter.init(buildStateAndDerivative(field, t0), field.getZero().add(t1));
121 
122         // first pass to set up switches history for a long period
123         double h = FastMath.copySign(0.05, t1 - t0);
124         double n = (int) FastMath.floor((t1 - t0) / h);
125         for (int i = 0; i < n; ++i) {
126             double t = t0 + i * h;
127             eventFilter.g(buildStateAndDerivative(field, t));
128         }
129 
130         // verify old events are preserved, even if randomly accessed
131         RandomGenerator rng = new Well19937a(0xb0e7401265af8cd3l);
132         for (int i = 0; i < 5000; i++) {
133             double t = t0 + (t1 - t0) * rng.nextDouble();
134             double g = eventFilter.g(buildStateAndDerivative(field, t)).getReal();
135             int turn = (int) FastMath.floor((t - refSwitch) / (2 * FastMath.PI));
136             if (turn % 2 == 0) {
137                 assertEquals( signEven * FastMath.sin(t), g, 1.0e-10);
138             } else {
139                 assertEquals(-signEven * FastMath.sin(t), g, 1.0e-10);
140             }
141         }
142 
143     }
144 
145     private <T extends CalculusFieldElement<T>> FieldODEStateAndDerivative<T> buildStateAndDerivative(Field<T> field, double t0) {
146         T t0F      = field.getZero().add(t0);
147         T[] y0F    = MathArrays.buildArray(field, 2);
148         y0F[0]     = FastMath.sin(t0F);
149         y0F[1]     = FastMath.cos(t0F);
150         T[] y0DotF = MathArrays.buildArray(field, 2);
151         y0DotF[0]  = FastMath.cos(t0F);
152         y0DotF[1]  = FastMath.sin(t0F).negate();
153         return new FieldODEStateAndDerivative<>(t0F, y0F, y0DotF);
154     }
155 
156     @Test
157     void testIncreasingOnlyField() {
158         doTestIncreasingOnlyField(Binary64Field.getInstance());
159     }
160 
161     private <T extends CalculusFieldElement<T>> void doTestIncreasingOnlyField(Field<T> field) {
162         FieldODEIntegrator<T> integrator = new DormandPrince853FieldIntegrator<>(field, 1.0e-3, 100.0, 1e-7, 1e-7);
163         FieldEvent<T> allEvents = new FieldEvent<>(0.1, field.getZero().newInstance(1.0e-7), 100, true, true);
164         integrator.addEventDetector(allEvents);
165         FieldEvent<T> onlyIncreasing = new FieldEvent<>(0.1, field.getZero().newInstance(1.0e-7), 100, false, true);
166         integrator.addEventDetector(new FieldEventSlopeFilter<>(field, onlyIncreasing,
167                                                                 FilterType.TRIGGER_ONLY_INCREASING_EVENTS));
168         T t0   = field.getZero().add(0.5 * FastMath.PI);
169         T tEnd = field.getZero().add(5.5 * FastMath.PI);
170         T[] y  = MathArrays.buildArray(field, 2);
171         y[0]   = field.getZero();
172         y[1]   = field.getOne();
173         assertEquals(tEnd.getReal(),
174                             integrator.integrate(new FieldExpandableODE<>(new FieldSineCosine<T>()),
175                                                  new FieldODEState<>(t0, y), tEnd).getTime().getReal(),
176                             1.0e-7);
177 
178         assertEquals(5, allEvents.getEventCount());
179         assertEquals(2, onlyIncreasing.getEventCount());
180 
181     }
182 
183     @Test
184     void testDecreasingOnlyField() {
185         doTestDecreasingOnlyField(Binary64Field.getInstance());
186     }
187 
188     private <T extends CalculusFieldElement<T>> void doTestDecreasingOnlyField(Field<T> field) {
189         FieldODEIntegrator<T> integrator = new DormandPrince853FieldIntegrator<>(field, 1.0e-3, 100.0, 1e-7, 1e-7);
190         FieldEvent<T> allEvents = new FieldEvent<>(0.1, field.getZero().newInstance(1.0e-7), 100, true, true);
191         integrator.addEventDetector(allEvents);
192         FieldEvent<T> onlyDecreasing = new FieldEvent<>(0.1, field.getZero().newInstance(1.0e-7), 100, true, false);
193         integrator.addEventDetector(new FieldEventSlopeFilter<>(field, onlyDecreasing,
194                                                                 FilterType.TRIGGER_ONLY_DECREASING_EVENTS));
195         T t0   = field.getZero().add(0.5 * FastMath.PI);
196         T tEnd = field.getZero().add(5.5 * FastMath.PI);
197         T[] y  = MathArrays.buildArray(field, 2);
198         y[0]   = field.getZero();
199         y[1]   = field.getOne();
200         assertEquals(tEnd.getReal(),
201                             integrator.integrate(new FieldExpandableODE<>(new FieldSineCosine<T>()),
202                                                  new FieldODEState<>(t0, y), tEnd).getTime().getReal(),
203                             1.0e-7);
204 
205         assertEquals(5, allEvents.getEventCount());
206         assertEquals(3, onlyDecreasing.getEventCount());
207 
208     }
209 
210     @Test
211     void testTwoOppositeFiltersField() {
212         doestTwoOppositeFiltersField(Binary64Field.getInstance());
213     }
214 
215     private <T extends CalculusFieldElement<T>> void doestTwoOppositeFiltersField(Field<T> field)
216         throws MathIllegalArgumentException, MathIllegalStateException {
217         FieldODEIntegrator<T> integrator = new DormandPrince853FieldIntegrator<>(field, 1.0e-3, 100.0, 1e-7, 1e-7);
218         FieldEvent<T> allEvents = new FieldEvent<>(0.1, field.getZero().newInstance(1.0e-7), 100, true, true);
219         integrator.addEventDetector(allEvents);
220         FieldEvent<T> onlyIncreasing = new FieldEvent<>(0.1, field.getZero().newInstance(1.0e-7), 100, false, true);
221         integrator.addEventDetector(new FieldEventSlopeFilter<>(field, onlyIncreasing,
222                                                                 FilterType.TRIGGER_ONLY_INCREASING_EVENTS));
223         FieldEvent<T> onlyDecreasing = new FieldEvent<>(0.1, field.getZero().newInstance(1.0e-7), 100, true, false);
224         integrator.addEventDetector(new FieldEventSlopeFilter<>(field, onlyDecreasing,
225                                                                 FilterType.TRIGGER_ONLY_DECREASING_EVENTS));
226         T t0   = field.getZero().add(0.5 * FastMath.PI);
227         T tEnd = field.getZero().add(5.5 * FastMath.PI);
228         T[] y  = MathArrays.buildArray(field, 2);
229         y[0]   = field.getZero();
230         y[1]   = field.getOne();
231         assertEquals(tEnd.getReal(),
232                             integrator.integrate(new FieldExpandableODE<>(new FieldSineCosine<T>()),
233                                                  new FieldODEState<>(t0, y), tEnd).getTime().getReal(),
234                             1.0e-7);
235 
236         assertEquals(5, allEvents.getEventCount());
237         assertEquals(2, onlyIncreasing.getEventCount());
238         assertEquals(3, onlyDecreasing.getEventCount());
239 
240     }
241 
242     private static class FieldSineCosine<T extends CalculusFieldElement<T>> implements FieldOrdinaryDifferentialEquation<T> {
243         public int getDimension() {
244             return 2;
245         }
246 
247         public T[] computeDerivatives(T t, T[] y) {
248             final T[] yDot = MathArrays.buildArray(t.getField(), getDimension());
249             yDot[0] = y[1];
250             yDot[1] = y[0].negate();
251             return yDot;
252         }
253     }
254 
255     /** State events for this unit test. */
256     protected static class FieldEvent<T extends CalculusFieldElement<T>> implements FieldODEEventDetector<T> {
257 
258         private final FieldAdaptableInterval<T>             maxCheck;
259         private final int                                   maxIter;
260         private final BracketedRealFieldUnivariateSolver<T> solver;
261         private final boolean                               expectDecreasing;
262         private final boolean                               expectIncreasing;
263         private int                                         eventCount;
264 
265         public FieldEvent(final double maxCheck, final T threshold, final int maxIter,
266                           boolean expectDecreasing, boolean expectIncreasing) {
267             this.maxCheck         = (s, isForward) -> maxCheck;
268             this.maxIter          = maxIter;
269             this.solver           = new FieldBracketingNthOrderBrentSolver<>(threshold.getField().getZero(),
270                                                                             threshold,
271                                                                             threshold.getField().getZero(),
272                                                                             5);
273             this.expectDecreasing = expectDecreasing;
274             this.expectIncreasing = expectIncreasing;
275         }
276 
277         public FieldAdaptableInterval<T> getMaxCheckInterval() {
278             return maxCheck;
279         }
280 
281         public int getMaxIterationCount() {
282             return maxIter;
283         }
284 
285         public BracketedRealFieldUnivariateSolver<T> getSolver() {
286             return solver;
287         }
288 
289         public int getEventCount() {
290             return eventCount;
291         }
292 
293         @Override
294         public void init(FieldODEStateAndDerivative<T> s0, T t) {
295             eventCount = 0;
296         }
297 
298         public T g(FieldODEStateAndDerivative<T> s) {
299             return s.getPrimaryState()[0];
300         }
301 
302         public FieldODEEventHandler<T> getHandler() {
303             return (FieldODEStateAndDerivative<T> s, FieldODEEventDetector<T> detector, boolean increasing) -> {
304                 if (increasing) {
305                     assertTrue(expectIncreasing);
306                 } else {
307                     assertTrue(expectDecreasing);
308                 }
309                 eventCount++;
310                 return Action.RESET_STATE;
311             };
312         }
313 
314     }
315 }