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