Search in sources :

Example 1 with ArrayRealVector

use of org.hipparchus.linear.ArrayRealVector in project Orekit by CS-SI.

the class JacobianPropagatorConverter method getModel.

/**
 * {@inheritDoc}
 */
protected MultivariateJacobianFunction getModel() {
    return new MultivariateJacobianFunction() {

        /**
         * {@inheritDoc}
         */
        public Pair<RealVector, RealMatrix> value(final RealVector point) throws IllegalArgumentException, OrekitExceptionWrapper {
            try {
                final RealVector value = new ArrayRealVector(getTargetSize());
                final RealMatrix jacobian = MatrixUtils.createRealMatrix(getTargetSize(), point.getDimension());
                final NumericalPropagator prop = builder.buildPropagator(point.toArray());
                final int stateSize = isOnlyPosition() ? 3 : 6;
                final ParameterDriversList orbitalParameters = builder.getOrbitalParametersDrivers();
                final PartialDerivativesEquations pde = new PartialDerivativesEquations("pde", prop);
                final ParameterDriversList propagationParameters = pde.getSelectedParameters();
                prop.setInitialState(pde.setInitialJacobians(prop.getInitialState()));
                final JacobiansMapper mapper = pde.getMapper();
                final List<SpacecraftState> sample = getSample();
                for (int i = 0; i < sample.size(); ++i) {
                    final int row = i * stateSize;
                    if (prop.getInitialState().getDate().equals(sample.get(i).getDate())) {
                        // use initial state and Jacobians
                        fillRows(value, jacobian, row, prop.getInitialState(), stateSize, orbitalParameters, propagationParameters, mapper);
                    } else {
                        // use a date detector to pick up state and Jacobians
                        prop.addEventDetector(new DateDetector(sample.get(i).getDate()).withHandler(new EventHandler<DateDetector>() {

                            /**
                             * {@inheritDoc}
                             */
                            @Override
                            public Action eventOccurred(final SpacecraftState state, final DateDetector detector, final boolean increasing) throws OrekitException {
                                fillRows(value, jacobian, row, state, stateSize, orbitalParameters, propagationParameters, mapper);
                                return row + stateSize >= getTargetSize() ? Action.STOP : Action.CONTINUE;
                            }
                        }));
                    }
                }
                prop.propagate(sample.get(sample.size() - 1).getDate().shiftedBy(10.0));
                return new Pair<RealVector, RealMatrix>(value, jacobian);
            } catch (OrekitException ex) {
                throw new OrekitExceptionWrapper(ex);
            }
        }
    };
}
Also used : DateDetector(org.orekit.propagation.events.DateDetector) OrekitExceptionWrapper(org.orekit.errors.OrekitExceptionWrapper) ArrayRealVector(org.hipparchus.linear.ArrayRealVector) EventHandler(org.orekit.propagation.events.handlers.EventHandler) MultivariateJacobianFunction(org.hipparchus.optim.nonlinear.vector.leastsquares.MultivariateJacobianFunction) SpacecraftState(org.orekit.propagation.SpacecraftState) RealMatrix(org.hipparchus.linear.RealMatrix) NumericalPropagator(org.orekit.propagation.numerical.NumericalPropagator) ParameterDriversList(org.orekit.utils.ParameterDriversList) PartialDerivativesEquations(org.orekit.propagation.numerical.PartialDerivativesEquations) ArrayRealVector(org.hipparchus.linear.ArrayRealVector) RealVector(org.hipparchus.linear.RealVector) OrekitException(org.orekit.errors.OrekitException) JacobiansMapper(org.orekit.propagation.numerical.JacobiansMapper) Pair(org.hipparchus.util.Pair)

Example 2 with ArrayRealVector

use of org.hipparchus.linear.ArrayRealVector in project Orekit by CS-SI.

the class Model method unNormalizeStateVector.

/**
 * Un-normalize a state vector.
 * A state vector S is of size M = nbOrb + nbPropag + nbMeas
 * For each parameter i the normalized value of the state vector is:
 * Sn[i] = S[i] / scale[i]
 * @param normalizedStateVector The normalized state vector in input
 * @return the "physical" state vector
 */
private RealVector unNormalizeStateVector(final RealVector normalizedStateVector) {
    // Initialize output matrix
    final int nbParams = normalizedStateVector.getDimension();
    final RealVector physicalStateVector = new ArrayRealVector(nbParams);
    // Retrieve the scaling factors
    final double[] scale = getParametersScale();
    // Normalize the state matrix
    for (int i = 0; i < nbParams; ++i) {
        physicalStateVector.setEntry(i, normalizedStateVector.getEntry(i) * scale[i]);
    }
    return physicalStateVector;
}
Also used : ArrayRealVector(org.hipparchus.linear.ArrayRealVector) RealVector(org.hipparchus.linear.RealVector) ArrayRealVector(org.hipparchus.linear.ArrayRealVector)

Example 3 with ArrayRealVector

use of org.hipparchus.linear.ArrayRealVector in project Orekit by CS-SI.

the class FiniteDifferencePropagatorConverterTest method testGetObjectiveFunctionParametersOnlyScaledOnce.

/**
 * Test case for bug #362. Check that scaling is only applied once.
 *
 * @throws OrekitException on error.
 */
@Test
public void testGetObjectiveFunctionParametersOnlyScaledOnce() throws OrekitException {
    // setup
    // create some arbitrary sample data to run with
    Frame eci = FramesFactory.getGCRF();
    double gm = Constants.EIGEN5C_EARTH_MU;
    AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
    Propagator source = new KeplerianPropagator(new KeplerianOrbit(6878137, 0, 0, 0, 0, 0, PositionAngle.TRUE, eci, date, gm));
    // Create a mock builder that allows us to check the values passed to it
    PropagatorBuilder builder = Mockito.mock(PropagatorBuilder.class);
    ParameterDriversList list = new ParameterDriversList();
    list.add(new ParameterDriver("p1", 0, 1e-3, Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY));
    Mockito.when(builder.getOrbitalParametersDrivers()).thenReturn(list);
    Mockito.when(builder.getPropagationParametersDrivers()).thenReturn(new ParameterDriversList());
    Mockito.when(builder.getFrame()).thenReturn(eci);
    Mockito.when(builder.getSelectedNormalizedParameters()).thenReturn(new double[1]);
    Mockito.when(builder.buildPropagator(Mockito.any(double[].class))).thenReturn(source);
    // subject under test
    FiniteDifferencePropagatorConverter converter = new FiniteDifferencePropagatorConverter(builder, 1, 100);
    // set some internal variables in FDPC that are not settable using another
    // interface
    converter.convert(source, 1, 2);
    // Forget all the side effect of the call to convert()
    Mockito.clearInvocations(builder);
    // action
    converter.getModel().value(new ArrayRealVector(1));
    // verify
    Mockito.verify(builder).buildPropagator(new double[] { 0 });
    Mockito.verify(builder).buildPropagator(new double[] { 1 });
}
Also used : KeplerianPropagator(org.orekit.propagation.analytical.KeplerianPropagator) Frame(org.orekit.frames.Frame) ParameterDriversList(org.orekit.utils.ParameterDriversList) Propagator(org.orekit.propagation.Propagator) KeplerianPropagator(org.orekit.propagation.analytical.KeplerianPropagator) ArrayRealVector(org.hipparchus.linear.ArrayRealVector) KeplerianOrbit(org.orekit.orbits.KeplerianOrbit) ParameterDriver(org.orekit.utils.ParameterDriver) AbsoluteDate(org.orekit.time.AbsoluteDate) Test(org.junit.Test)

Example 4 with ArrayRealVector

use of org.hipparchus.linear.ArrayRealVector in project Orekit by CS-SI.

the class ModelTest method testPerfectValue.

@Test
public void testPerfectValue() throws OrekitException {
    final Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
    final NumericalPropagatorBuilder propagatorBuilder = context.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true, 1.0e-6, 60.0, 0.001);
    final NumericalPropagatorBuilder[] builders = { propagatorBuilder };
    // create perfect PV measurements
    final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder);
    final List<ObservedMeasurement<?>> measurements = EstimationTestUtils.createMeasurements(propagator, new PVMeasurementCreator(), 0.0, 1.0, 300.0);
    final ParameterDriversList estimatedMeasurementsParameters = new ParameterDriversList();
    for (ObservedMeasurement<?> measurement : measurements) {
        for (final ParameterDriver driver : measurement.getParametersDrivers()) {
            if (driver.isSelected()) {
                estimatedMeasurementsParameters.add(driver);
            }
        }
    }
    // create model
    final ModelObserver modelObserver = new ModelObserver() {

        /**
         * {@inheritDoc}
         */
        @Override
        public void modelCalled(final Orbit[] newOrbits, final Map<ObservedMeasurement<?>, EstimatedMeasurement<?>> newEvaluations) {
            Assert.assertEquals(1, newOrbits.length);
            Assert.assertEquals(0, context.initialOrbit.getDate().durationFrom(newOrbits[0].getDate()), 1.0e-15);
            Assert.assertEquals(0, Vector3D.distance(context.initialOrbit.getPVCoordinates().getPosition(), newOrbits[0].getPVCoordinates().getPosition()), 1.0e-15);
            Assert.assertEquals(measurements.size(), newEvaluations.size());
        }
    };
    final Model model = new Model(builders, measurements, estimatedMeasurementsParameters, modelObserver);
    model.setIterationsCounter(new Incrementor(100));
    model.setEvaluationsCounter(new Incrementor(100));
    // Test forward propagation flag to true
    assertEquals(true, model.isForwardPropagation());
    // evaluate model on perfect start point
    final double[] normalizedProp = propagatorBuilder.getSelectedNormalizedParameters();
    final double[] normalized = new double[normalizedProp.length + estimatedMeasurementsParameters.getNbParams()];
    System.arraycopy(normalizedProp, 0, normalized, 0, normalizedProp.length);
    int i = normalizedProp.length;
    for (final ParameterDriver driver : estimatedMeasurementsParameters.getDrivers()) {
        normalized[i++] = driver.getNormalizedValue();
    }
    Pair<RealVector, RealMatrix> value = model.value(new ArrayRealVector(normalized));
    int index = 0;
    for (ObservedMeasurement<?> measurement : measurements) {
        for (int k = 0; k < measurement.getDimension(); ++k) {
            // the value is already a weighted residual
            Assert.assertEquals(0.0, value.getFirst().getEntry(index++), 1.6e-7);
        }
    }
    Assert.assertEquals(index, value.getFirst().getDimension());
}
Also used : Context(org.orekit.estimation.Context) ArrayRealVector(org.hipparchus.linear.ArrayRealVector) Incrementor(org.hipparchus.util.Incrementor) ParameterDriver(org.orekit.utils.ParameterDriver) RealMatrix(org.hipparchus.linear.RealMatrix) ParameterDriversList(org.orekit.utils.ParameterDriversList) NumericalPropagatorBuilder(org.orekit.propagation.conversion.NumericalPropagatorBuilder) ArrayRealVector(org.hipparchus.linear.ArrayRealVector) RealVector(org.hipparchus.linear.RealVector) Propagator(org.orekit.propagation.Propagator) Map(java.util.Map) ObservedMeasurement(org.orekit.estimation.measurements.ObservedMeasurement) PVMeasurementCreator(org.orekit.estimation.measurements.PVMeasurementCreator) Test(org.junit.Test)

Example 5 with ArrayRealVector

use of org.hipparchus.linear.ArrayRealVector in project Orekit by CS-SI.

the class JacobianPropagatorConverterTest method doTestDerivatives.

private void doTestDerivatives(double tolP, double tolV, String... names) throws OrekitException {
    // we use a fixed step integrator on purpose
    // as the test is based on external differentiation using finite differences,
    // an adaptive step size integrator would introduce *lots* of numerical noise
    NumericalPropagatorBuilder builder = new NumericalPropagatorBuilder(OrbitType.CARTESIAN.convertType(orbit), new LutherIntegratorBuilder(10.0), PositionAngle.TRUE, dP);
    builder.setMass(200.0);
    builder.addForceModel(drag);
    builder.addForceModel(gravity);
    // retrieve a state slightly different from the initial state,
    // using normalized values different from 0.0 for the sake of generality
    RandomGenerator random = new Well19937a(0xe67f19c1a678d037l);
    List<ParameterDriver> all = new ArrayList<ParameterDriver>();
    for (final ParameterDriver driver : builder.getOrbitalParametersDrivers().getDrivers()) {
        all.add(driver);
    }
    for (final ParameterDriver driver : builder.getPropagationParametersDrivers().getDrivers()) {
        all.add(driver);
    }
    double[] normalized = new double[names.length];
    List<ParameterDriver> selected = new ArrayList<ParameterDriver>(names.length);
    int index = 0;
    for (final ParameterDriver driver : all) {
        boolean found = false;
        for (final String name : names) {
            if (name.equals(driver.getName())) {
                found = true;
                normalized[index++] = driver.getNormalizedValue() + (2 * random.nextDouble() - 1);
                selected.add(driver);
            }
        }
        driver.setSelected(found);
    }
    // create a one hour sample that starts 10 minutes after initial state
    // the 10 minutes offset implies even the first point is influenced by model parameters
    final List<SpacecraftState> sample = new ArrayList<SpacecraftState>();
    Propagator propagator = builder.buildPropagator(normalized);
    propagator.setMasterMode(60.0, new OrekitFixedStepHandler() {

        @Override
        public void handleStep(SpacecraftState currentState, boolean isLast) {
            sample.add(currentState);
        }
    });
    propagator.propagate(orbit.getDate().shiftedBy(600.0), orbit.getDate().shiftedBy(4200.0));
    JacobianPropagatorConverter fitter = new JacobianPropagatorConverter(builder, 1.0e-3, 5000);
    try {
        Method setSample = AbstractPropagatorConverter.class.getDeclaredMethod("setSample", List.class);
        setSample.setAccessible(true);
        setSample.invoke(fitter, sample);
    } catch (NoSuchMethodException | SecurityException | IllegalAccessException | IllegalArgumentException | InvocationTargetException e) {
        Assert.fail(e.getLocalizedMessage());
    }
    MultivariateVectorFunction f = fitter.getObjectiveFunction();
    Pair<RealVector, RealMatrix> p = fitter.getModel().value(new ArrayRealVector(normalized));
    // check derivatives
    // a h offset on normalized parameter represents a physical offset of h * scale
    RealMatrix m = p.getSecond();
    double h = 10.0;
    double[] shifted = normalized.clone();
    double maxErrorP = 0;
    double maxErrorV = 0;
    for (int j = 0; j < selected.size(); ++j) {
        shifted[j] = normalized[j] + 2.0 * h;
        double[] valueP2 = f.value(shifted);
        shifted[j] = normalized[j] + 1.0 * h;
        double[] valueP1 = f.value(shifted);
        shifted[j] = normalized[j] - 1.0 * h;
        double[] valueM1 = f.value(shifted);
        shifted[j] = normalized[j] - 2.0 * h;
        double[] valueM2 = f.value(shifted);
        shifted[j] = normalized[j];
        for (int i = 0; i < valueP2.length; ++i) {
            double d = (8 * (valueP1[i] - valueM1[i]) - (valueP2[i] - valueM2[i])) / (12 * h);
            if (i % 6 < 3) {
                // position
                maxErrorP = FastMath.max(maxErrorP, FastMath.abs(m.getEntry(i, j) - d));
            } else {
                // velocity
                maxErrorV = FastMath.max(maxErrorV, FastMath.abs(m.getEntry(i, j) - d));
            }
        }
    }
    Assert.assertEquals(0.0, maxErrorP, tolP);
    Assert.assertEquals(0.0, maxErrorV, tolV);
}
Also used : ArrayList(java.util.ArrayList) Well19937a(org.hipparchus.random.Well19937a) RandomGenerator(org.hipparchus.random.RandomGenerator) SpacecraftState(org.orekit.propagation.SpacecraftState) ArrayRealVector(org.hipparchus.linear.ArrayRealVector) RealVector(org.hipparchus.linear.RealVector) Propagator(org.orekit.propagation.Propagator) OrekitFixedStepHandler(org.orekit.propagation.sampling.OrekitFixedStepHandler) MathIllegalArgumentException(org.hipparchus.exception.MathIllegalArgumentException) ArrayRealVector(org.hipparchus.linear.ArrayRealVector) Method(java.lang.reflect.Method) ParameterDriver(org.orekit.utils.ParameterDriver) MultivariateVectorFunction(org.hipparchus.analysis.MultivariateVectorFunction) InvocationTargetException(java.lang.reflect.InvocationTargetException) RealMatrix(org.hipparchus.linear.RealMatrix)

Aggregations

ArrayRealVector (org.hipparchus.linear.ArrayRealVector)5 RealVector (org.hipparchus.linear.RealVector)4 RealMatrix (org.hipparchus.linear.RealMatrix)3 Propagator (org.orekit.propagation.Propagator)3 ParameterDriver (org.orekit.utils.ParameterDriver)3 ParameterDriversList (org.orekit.utils.ParameterDriversList)3 Test (org.junit.Test)2 SpacecraftState (org.orekit.propagation.SpacecraftState)2 InvocationTargetException (java.lang.reflect.InvocationTargetException)1 Method (java.lang.reflect.Method)1 ArrayList (java.util.ArrayList)1 Map (java.util.Map)1 MultivariateVectorFunction (org.hipparchus.analysis.MultivariateVectorFunction)1 MathIllegalArgumentException (org.hipparchus.exception.MathIllegalArgumentException)1 MultivariateJacobianFunction (org.hipparchus.optim.nonlinear.vector.leastsquares.MultivariateJacobianFunction)1 RandomGenerator (org.hipparchus.random.RandomGenerator)1 Well19937a (org.hipparchus.random.Well19937a)1 Incrementor (org.hipparchus.util.Incrementor)1 Pair (org.hipparchus.util.Pair)1 OrekitException (org.orekit.errors.OrekitException)1