1   /* Copyright 2002-2026 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;
18  
19  import java.util.ArrayList;
20  import java.util.Arrays;
21  import java.util.IdentityHashMap;
22  import java.util.List;
23  import java.util.Map;
24  import java.util.function.Function;
25  
26  import org.hipparchus.linear.MatrixUtils;
27  import org.hipparchus.linear.RealMatrix;
28  import org.orekit.propagation.SpacecraftState;
29  import org.orekit.utils.drivers.ParameterDriver;
30  import org.orekit.utils.drivers.ParameterDriversList;
31  import org.orekit.utils.TimeStampedPVCoordinates;
32  
33  /** Class multiplexing several measurements as one.
34   * <p>
35   * Date comes from the first measurement, observed and estimated
36   * values result from gathering all underlying measurements values.
37   *
38   * @author Luc Maisonobe
39   * @since 10.1
40   */
41  public class MultiplexedMeasurement extends AbstractMeasurement<MultiplexedMeasurement> {
42  
43      /** Type of the measurement. */
44      public static final String MEASUREMENT_TYPE = "MultiplexedMeasurement";
45  
46      /** Multiplexed measurements. */
47      private final List<ObservedMeasurement<?>> observedMeasurements;
48  
49      /** Multiplexed measurements without derivatives.
50       */
51      private final List<EstimatedMeasurementBase<?>> estimatedMeasurementsWithoutDerivatives;
52  
53      /** Multiplexed measurements. */
54      private final List<EstimatedMeasurement<?>> estimatedMeasurements;
55  
56      /** Multiplexed parameters drivers. */
57      private final ParameterDriversList parametersDrivers;
58  
59      /** Total dimension. */
60      private final int dimension;
61  
62      /** Total number of satellites involved. */
63      private final int nbSat;
64  
65      /** States mapping. */
66      private final int[][] multiplexedToUnderlying;
67  
68      /** States mapping. */
69      private final int[][] underlyingToMultiplexed;
70  
71      /** Simple constructor.
72       * @param measurements measurements to multiplex
73       * @since 10.1
74       */
75      public MultiplexedMeasurement(final List<ObservedMeasurement<?>> measurements) {
76          super(measurements.getFirst().getDate(),
77                multiplex(measurements, ComparableMeasurement::getObservedValue),
78                multiplexMeasurementQuality(measurements), multiplex(measurements));
79  
80          this.observedMeasurements                    = measurements;
81          this.estimatedMeasurementsWithoutDerivatives = new ArrayList<>();
82          this.estimatedMeasurements                   = new ArrayList<>();
83          this.parametersDrivers                       = new ParameterDriversList();
84  
85          // gather parameters drivers
86          int dim = 0;
87          for (final ObservedMeasurement<?> m : measurements) {
88              for (final ParameterDriver driver : m.getParametersDrivers()) {
89                  parametersDrivers.add(driver);
90              }
91              dim += m.getDimension();
92          }
93          parametersDrivers.sort();
94          for (final ParameterDriver driver : parametersDrivers.getDrivers()) {
95              addParameterDriver(driver);
96          }
97          this.dimension = dim;
98  
99          // set up states mappings for observed satellites
100         final List<ObservableSatellite> deduplicated = getSatellites();
101         this.nbSat   = deduplicated.size();
102         this.multiplexedToUnderlying = new int[measurements.size()][];
103         this.underlyingToMultiplexed = new int[measurements.size()][deduplicated.size()];
104         for (int i = 0; i < multiplexedToUnderlying.length; ++i) {
105             final List<ObservableSatellite> satellites = measurements.get(i).getSatellites();
106             multiplexedToUnderlying[i] = new int[satellites.size()];
107             for (int j = 0; j < multiplexedToUnderlying[i].length; ++j) {
108                 final int index = satellites.get(j).getPropagatorIndex();
109                 for (int k = 0; k < nbSat; ++k) {
110                     if (deduplicated.get(k).getPropagatorIndex() == index) {
111                         multiplexedToUnderlying[i][j] = k;
112                         underlyingToMultiplexed[i][k] = j;
113                         break;
114                     }
115                 }
116             }
117         }
118 
119     }
120 
121     /** Get the underlying measurements.
122      * @return underlying measurements
123      */
124     public List<ObservedMeasurement<?>> getMeasurements() {
125         return observedMeasurements;
126     }
127 
128     /** Get the underlying estimated measurements without derivatives.
129      * @return underlying estimated measurements without derivatives
130      * @since 12.0
131      */
132     public List<EstimatedMeasurementBase<?>> getEstimatedMeasurementsWithoutDerivatives() {
133         return estimatedMeasurementsWithoutDerivatives;
134     }
135 
136     /** Get the underlying estimated measurements.
137      * @return underlying estimated measurements
138      */
139     public List<EstimatedMeasurement<?>> getEstimatedMeasurements() {
140         return estimatedMeasurements;
141     }
142 
143     /** Get the spacecraft state index in the underlying measurement.
144      * @param measurementIndex index of the underlying measurement
145      * @param multiplexedStateIndex index of the spacecraft state in the multiplexed array
146      * @return spacecraft state index in the underlying measurement
147      * @since 13.0
148      */
149     public int getUnderlyingStateIndex(final int measurementIndex, final int multiplexedStateIndex) {
150         return multiplexedToUnderlying[measurementIndex][multiplexedStateIndex];
151     }
152 
153     /** Get the spacecraft state index in the multiplexed measurement.
154      * @param measurementIndex index of the underlying measurement
155      * @param underlyingStateIndex index of the spacecraft state in the underlying array
156      * @return spacecraft state index in the multiplexed measurement
157      * @since 13.0
158      */
159     public int getMultiplexedStateIndex(final int measurementIndex, final int underlyingStateIndex) {
160         return underlyingToMultiplexed[measurementIndex][underlyingStateIndex];
161     }
162 
163     /** {@inheritDoc} */
164     @Override
165     protected EstimatedMeasurementBase<MultiplexedMeasurement> theoreticalEvaluationWithoutDerivatives(final int iteration,
166                                                                                                        final int evaluation,
167                                                                                                        final SpacecraftState[] states,
168                                                                                                        final boolean fillParticipants) {
169 
170         final SpacecraftState[]              evaluationStates = new SpacecraftState[nbSat];
171         final double[]                       value            = new double[dimension];
172 
173         // loop over all multiplexed measurements
174         estimatedMeasurementsWithoutDerivatives.clear();
175         int index = 0;
176         for (int i = 0; i < observedMeasurements.size(); ++i) {
177 
178             // filter states involved in the current measurement
179             final SpacecraftState[] filteredStates = new SpacecraftState[multiplexedToUnderlying[i].length];
180             for (int j = 0; j < multiplexedToUnderlying[i].length; ++j) {
181                 filteredStates[j] = states[getUnderlyingStateIndex(i, j)];
182             }
183 
184             // perform evaluation
185             final EstimatedMeasurementBase<?> eI = observedMeasurements.get(i).estimateWithoutDerivatives(iteration,
186                     evaluation, filteredStates);
187             estimatedMeasurementsWithoutDerivatives.add(eI);
188 
189             // extract results
190             final double[] valueI = eI.getEstimatedValue();
191             System.arraycopy(valueI, 0, value, index, valueI.length);
192             index += valueI.length;
193 
194             // extract states
195             final SpacecraftState[] statesI = eI.getStates();
196             for (int j = 0; j < multiplexedToUnderlying[i].length; ++j) {
197                 evaluationStates[multiplexedToUnderlying[i][j]] = statesI[j];
198             }
199 
200         }
201 
202         // create multiplexed estimation
203         final EstimatedMeasurementBase<MultiplexedMeasurement> multiplexed =
204                         new EstimatedMeasurementBase<>(this, iteration, evaluation,
205                                                        evaluationStates,
206                                                        new TimeStampedPVCoordinates[0]);
207 
208         // copy multiplexed value
209         multiplexed.setEstimatedValue(value);
210 
211         return multiplexed;
212 
213     }
214 
215     /** {@inheritDoc} */
216     @Override
217     protected EstimatedMeasurement<MultiplexedMeasurement> theoreticalEvaluation(final int iteration, final int evaluation,
218                                                                                  final SpacecraftState[] states) {
219 
220         final SpacecraftState[]              evaluationStates = new SpacecraftState[nbSat];
221         final double[]                       value            = new double[dimension];
222 
223         // loop over all multiplexed measurements
224         estimatedMeasurements.clear();
225         int index = 0;
226         for (int i = 0; i < observedMeasurements.size(); ++i) {
227 
228             // filter states involved in the current measurement
229             final SpacecraftState[] filteredStates = new SpacecraftState[multiplexedToUnderlying[i].length];
230             for (int j = 0; j < multiplexedToUnderlying[i].length; ++j) {
231                 filteredStates[j] = states[multiplexedToUnderlying[i][j]];
232             }
233 
234             // perform evaluation
235             final EstimatedMeasurement<?> eI = observedMeasurements.get(i).estimate(iteration, evaluation, filteredStates);
236             estimatedMeasurements.add(eI);
237 
238             // extract results
239             final double[] valueI = eI.getEstimatedValue();
240             System.arraycopy(valueI, 0, value, index, valueI.length);
241             index += valueI.length;
242 
243             // extract states
244             final SpacecraftState[] statesI = eI.getStates();
245             for (int j = 0; j < multiplexedToUnderlying[i].length; ++j) {
246                 evaluationStates[multiplexedToUnderlying[i][j]] = statesI[j];
247             }
248 
249         }
250 
251         // create multiplexed estimation
252         final EstimatedMeasurement<MultiplexedMeasurement> multiplexed =
253                         new EstimatedMeasurement<>(this, iteration, evaluation,
254                                                    evaluationStates,
255                                                    new TimeStampedPVCoordinates[0]);
256 
257         // copy multiplexed value
258         multiplexed.setEstimatedValue(value);
259 
260         // combine derivatives
261         final int                            stateSize             = estimatedMeasurements.getFirst().getStateSize();
262         final double[]                       zeroDerivative        = new double[stateSize];
263         final double[][][]                   stateDerivatives      = new double[nbSat][dimension][];
264         for (final double[][] m : stateDerivatives) {
265             Arrays.fill(m, zeroDerivative);
266         }
267 
268         final Map<ParameterDriver, double[]> parametersDerivatives = new IdentityHashMap<>();
269         index = 0;
270         for (int i = 0; i < observedMeasurements.size(); ++i) {
271 
272             final EstimatedMeasurement<?> eI   = estimatedMeasurements.get(i);
273             final int                     idx  = index;
274             final int                     dimI = eI.getObservedMeasurement().getDimension();
275 
276             // state derivatives
277             for (int j = 0; j < multiplexedToUnderlying[i].length; ++j) {
278                 System.arraycopy(eI.getStateDerivatives(j), 0,
279                                  stateDerivatives[multiplexedToUnderlying[i][j]], index,
280                                  dimI);
281             }
282 
283             // parameters derivatives
284             eI.getDerivativesDrivers().forEach(driver -> {
285                 final ParameterDriversList.DelegatingDriver delegating = parametersDrivers.findByName(driver.getName());
286                 final double[] derivatives = parametersDerivatives.computeIfAbsent(delegating, d -> new double[dimension]);
287                 System.arraycopy(eI.getParameterDerivatives(driver), 0, derivatives, idx, dimI);
288             });
289 
290             index += dimI;
291 
292         }
293 
294         // set states derivatives
295         for (int i = 0; i < nbSat; ++i) {
296             multiplexed.setStateDerivatives(i, stateDerivatives[i]);
297         }
298 
299         // set parameters derivatives
300         parametersDerivatives.forEach(multiplexed::setParameterDerivatives);
301 
302         return multiplexed;
303 
304     }
305 
306     /** Multiplex measurements data.
307      * @param measurements measurements to multiplex
308      * @return multiplexed data
309      */
310     private static MeasurementQuality multiplexMeasurementQuality(final List<ObservedMeasurement<?>> measurements) {
311 
312         final int totalSize = measurements.stream().mapToInt(ObservedMeasurement::getDimension).sum();
313         final double[] weights = new double[totalSize];
314         final RealMatrix covarianceMatrix = MatrixUtils.createRealMatrix(totalSize, totalSize);
315         int n = 0;
316         for (final ObservedMeasurement<?> measurement : measurements) {
317             System.arraycopy(measurement.getBaseWeight(), 0, weights, n, measurement.getDimension());
318             covarianceMatrix.setSubMatrix(measurement.getMeasurementQuality().getCovarianceMatrix().getData(), n, n);
319             n += measurement.getDimension();
320         }
321 
322         return new MeasurementQuality(covarianceMatrix.getData(), weights);
323 
324     }
325 
326     /** Multiplex measurements value.
327      * @param measurements measurements to multiplex
328      * @param extractor data extraction function
329      * @return multiplexed data
330      */
331     private static double[] multiplex(final List<ObservedMeasurement<?>> measurements,
332                                       final Function<ObservedMeasurement<?>, double[]> extractor) {
333 
334         // gather individual parts
335         final List<double[]> parts = new ArrayList<> (measurements.size());
336         int n = 0;
337         for (final ObservedMeasurement<?> measurement : measurements) {
338             final double[] p = extractor.apply(measurement);
339             parts.add(p);
340             n += p.length;
341         }
342 
343         // create multiplexed data
344         final double[] multiplexed = new double[n];
345         int index = 0;
346         for (final double[] p : parts) {
347             System.arraycopy(p, 0, multiplexed, index, p.length);
348             index += p.length;
349         }
350 
351         return multiplexed;
352 
353     }
354 
355     /** Multiplex satellites data.
356      * @param measurements measurements to multiplex
357      * @return multiplexed satellites data
358      */
359     private static List<ObservableSatellite> multiplex(final List<ObservedMeasurement<?>> measurements) {
360 
361         final List<ObservableSatellite> satellites = new ArrayList<>();
362 
363         // gather all satellites, removing duplicates
364         for (final ObservedMeasurement<?> measurement : measurements) {
365             for (final ObservableSatellite satellite : measurement.getSatellites()) {
366                 boolean searching = true;
367                 for (int i = 0; i < satellites.size() && searching; ++i) {
368                     // check if we already know this satellite
369                     searching = satellite.getPropagatorIndex() != satellites.get(i).getPropagatorIndex();
370                 }
371                 if (searching) {
372                     // this is a new satellite, add it to the global list
373                     satellites.add(satellite);
374                 }
375             }
376         }
377 
378         return satellites;
379 
380     }
381 
382 }