1   /* Copyright 2002-2025 CS GROUP
2    * Licensed to CS GROUP (CS) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * CS 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    *   http://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  package org.orekit.estimation.measurements.gnss;
18  
19  import org.hipparchus.util.FastMath;
20  
21  /** Sampler for generating long integers between two limits in an alternating pattern.
22   * <p>
23   * Given a center a and a radius r, this class will generate integers kᵢ such
24   * that a - r ≤ kᵢ ≤ a + r. The generation order will start from the middle
25   * (i.e. k₀ is the long integer closest to a) and go towards the boundaries,
26   * alternating between values lesser than a and values greater than a.
27   * For example, with a = 17.3 and r = 5.2, it will generate: k₀ = 17, k₁ = 18,
28   * k₂ = 16, k₃ = 19, k₄ = 15, k₅ = 20, k₆ = 14, k₇ = 21, k₈ = 13, k₉ = 22.
29   * </p>
30   * <p>
31   * There are no hard limits to the generation, i.e. in the example above, the
32   * generator will happily generate k₁₀ = 12, k₁₁ = 23, k₁₂ = 11... In fact, if
33   * there are no integers at all between {@code a - r} and {@code a + r}, even
34   * the initial k₀ that is implicitly generated at construction will be out of
35   * range. The {@link #inRange()} method can be used to check if the last generator
36   * is still producing numbers within the initial range or if it has already
37   * started generating out of range numbers.
38   * </p>
39   * <p>
40   * If there are integers between {@code a - r} and {@code a + r}, it is guaranteed
41   * that they will all be generated once before {@link #inRange()} starts returning
42   * {@code false}.
43   * </p>
44   * <p>
45   * This allows to explore the range for one integer ambiguity starting
46   * with the most probable values (closest to a) and continuing with
47   * values less probable.
48   * </p>
49   * @see <a href="https://www.researchgate.net/publication/2790708_The_LAMBDA_method_for_integer_ambiguity_estimation_implementation_aspects">
50   * The LAMBDA method for integer ambiguity estimation: implementation aspects</a>
51   * @see <a href="https://oeis.org/A001057">
52   * A001057: Canonical enumeration of integers: interleaved positive and negative integers with zero prepended.</a>
53   * @author Luc Maisonobe
54   * @since 10.0
55   */
56  class AlternatingSampler {
57  
58      /** Range midpoint. */
59      private final double a;
60  
61      /** Offset with respect to A001057. */
62      private final long offset;
63  
64      /** Sign with respect to A001057. */
65      private final long sign;
66  
67      /** Minimum number to generate. */
68      private long min;
69  
70      /** Maximum number to generate. */
71      private long max;
72  
73      /** Previous generated number in A001057. */
74      private long k1;
75  
76      /** Current generated number in A001057. */
77      private long k0;
78  
79      /** Current generated number. */
80      private long current;
81  
82      /** Simple constructor.
83       * <p>
84       * A first initial integer is already generated as a side effect of
85       * construction, so {@link #getCurrent()} can be called even before
86       * calling {@link #generateNext()}. If there are no integers at
87       * all between {@code a - r} and {@code a + r}, then this initial
88       * integer will already be out of range.
89       * </p>
90       * @param a range midpoint
91       * @param r range radius
92       */
93      AlternatingSampler(final double a, final double r) {
94  
95          this.a      = a;
96          this.offset = (long) FastMath.rint(a);
97          this.sign   = offset <= a ? +1 : -1;
98          setRadius(r);
99  
100         this.k1      = 0;
101         this.k0      = 0;
102         this.current = offset;
103     }
104 
105     /** Reset the range radius.
106      * <p>
107      * Resetting radius is allowed during sampling, it simply changes
108      * the boundaries used when calling {@link #inRange()}. Resetting
109      * the radius does not change the sampling itself, neither the
110      * {@link #getCurrent() current} value nor the {@link #generateNext()
111      * next generated} ones.
112      * </p>
113      * <p>
114      * A typical use case for calling {@link #setRadius(double)} during
115      * sampling is to reduce sampling interval. It is used to shrink
116      * the search ellipsoid on the fly in LAMBDA-based methods in order
117      * to speed-up search.
118      * </p>
119      * @param r range radius
120      */
121     public void setRadius(final double r) {
122         this.min = (long) FastMath.ceil(a - r);
123         this.max = (long) FastMath.floor(a + r);
124     }
125 
126     /** Get the range midpoint.
127      * @return range midpoint
128      */
129     public double getMidPoint() {
130         return a;
131     }
132 
133     /** Get current value.
134      * @return current value
135      */
136     public long getCurrent() {
137         return current;
138     }
139 
140     /** Check if the current value is within range.
141      * @return true if current value is within range
142      */
143     public boolean inRange() {
144         return min <= current && current <= max;
145     }
146 
147     /** Generate next value.
148      */
149     public void generateNext() {
150 
151         // apply A001057 recursion
152         final long k2 = k1;
153         k1 = k0;
154         k0 = 1 - (k1 << 1) - k2;
155 
156         // take offset and sign into account
157         current = offset + sign * k0;
158 
159     }
160 
161 }