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);
}
}
};
}
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;
}
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 });
}
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());
}
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);
}
Aggregations