source: src/main/java/agents/anac/y2019/harddealer/math3/ode/JacobianMatrices.java

Last change on this file was 204, checked in by Katsuhide Fujita, 5 years ago

Fixed errors of ANAC2019 agents
  • Property svn:executable set to *
File size: 19.2 KB
Line 
1/*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF 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 */
17package agents.anac.y2019.harddealer.math3.ode;
18
19import java.lang.reflect.Array;
20import java.util.ArrayList;
21import java.util.Arrays;
22import java.util.List;
23
24import agents.anac.y2019.harddealer.math3.exception.DimensionMismatchException;
25import agents.anac.y2019.harddealer.math3.exception.MathIllegalArgumentException;
26import agents.anac.y2019.harddealer.math3.exception.MaxCountExceededException;
27import agents.anac.y2019.harddealer.math3.exception.util.LocalizedFormats;
28
29/**
30 * This class defines a set of {@link SecondaryEquations secondary equations} to
31 * compute the Jacobian matrices with respect to the initial state vector and, if
32 * any, to some parameters of the primary ODE set.
33 * <p>
34 * It is intended to be packed into an {@link ExpandableStatefulODE}
35 * in conjunction with a primary set of ODE, which may be:
36 * <ul>
37 * <li>a {@link FirstOrderDifferentialEquations}</li>
38 * <li>a {@link MainStateJacobianProvider}</li>
39 * </ul>
40 * In order to compute Jacobian matrices with respect to some parameters of the
41 * primary ODE set, the following parameter Jacobian providers may be set:
42 * <ul>
43 * <li>a {@link ParameterJacobianProvider}</li>
44 * <li>a {@link ParameterizedODE}</li>
45 * </ul>
46 * </p>
47 *
48 * @see ExpandableStatefulODE
49 * @see FirstOrderDifferentialEquations
50 * @see MainStateJacobianProvider
51 * @see ParameterJacobianProvider
52 * @see ParameterizedODE
53 *
54 * @since 3.0
55 */
56public class JacobianMatrices {
57
58 /** Expandable first order differential equation. */
59 private ExpandableStatefulODE efode;
60
61 /** Index of the instance in the expandable set. */
62 private int index;
63
64 /** FODE with exact primary Jacobian computation skill. */
65 private MainStateJacobianProvider jode;
66
67 /** FODE without exact parameter Jacobian computation skill. */
68 private ParameterizedODE pode;
69
70 /** Main state vector dimension. */
71 private int stateDim;
72
73 /** Selected parameters for parameter Jacobian computation. */
74 private ParameterConfiguration[] selectedParameters;
75
76 /** FODE with exact parameter Jacobian computation skill. */
77 private List<ParameterJacobianProvider> jacobianProviders;
78
79 /** Parameters dimension. */
80 private int paramDim;
81
82 /** Boolean for selected parameters consistency. */
83 private boolean dirtyParameter;
84
85 /** State and parameters Jacobian matrices in a row. */
86 private double[] matricesData;
87
88 /** Simple constructor for a secondary equations set computing Jacobian matrices.
89 * <p>
90 * Parameters must belong to the supported ones given by {@link
91 * Parameterizable#getParametersNames()}, so the primary set of differential
92 * equations must be {@link Parameterizable}.
93 * </p>
94 * <p>Note that each selection clears the previous selected parameters.</p>
95 *
96 * @param fode the primary first order differential equations set to extend
97 * @param hY step used for finite difference computation with respect to state vector
98 * @param parameters parameters to consider for Jacobian matrices processing
99 * (may be null if parameters Jacobians is not desired)
100 * @exception DimensionMismatchException if there is a dimension mismatch between
101 * the steps array {@code hY} and the equation dimension
102 */
103 public JacobianMatrices(final FirstOrderDifferentialEquations fode, final double[] hY,
104 final String... parameters)
105 throws DimensionMismatchException {
106 this(new MainStateJacobianWrapper(fode, hY), parameters);
107 }
108
109 /** Simple constructor for a secondary equations set computing Jacobian matrices.
110 * <p>
111 * Parameters must belong to the supported ones given by {@link
112 * Parameterizable#getParametersNames()}, so the primary set of differential
113 * equations must be {@link Parameterizable}.
114 * </p>
115 * <p>Note that each selection clears the previous selected parameters.</p>
116 *
117 * @param jode the primary first order differential equations set to extend
118 * @param parameters parameters to consider for Jacobian matrices processing
119 * (may be null if parameters Jacobians is not desired)
120 */
121 public JacobianMatrices(final MainStateJacobianProvider jode,
122 final String... parameters) {
123
124 this.efode = null;
125 this.index = -1;
126
127 this.jode = jode;
128 this.pode = null;
129
130 this.stateDim = jode.getDimension();
131
132 if (parameters == null) {
133 selectedParameters = null;
134 paramDim = 0;
135 } else {
136 this.selectedParameters = new ParameterConfiguration[parameters.length];
137 for (int i = 0; i < parameters.length; ++i) {
138 selectedParameters[i] = new ParameterConfiguration(parameters[i], Double.NaN);
139 }
140 paramDim = parameters.length;
141 }
142 this.dirtyParameter = false;
143
144 this.jacobianProviders = new ArrayList<ParameterJacobianProvider>();
145
146 // set the default initial state Jacobian to the identity
147 // and the default initial parameters Jacobian to the null matrix
148 matricesData = new double[(stateDim + paramDim) * stateDim];
149 for (int i = 0; i < stateDim; ++i) {
150 matricesData[i * (stateDim + 1)] = 1.0;
151 }
152
153 }
154
155 /** Register the variational equations for the Jacobians matrices to the expandable set.
156 * @param expandable expandable set into which variational equations should be registered
157 * @throws DimensionMismatchException if the dimension of the partial state does not
158 * match the selected equations set dimension
159 * @exception MismatchedEquations if the primary set of the expandable set does
160 * not match the one used to build the instance
161 * @see ExpandableStatefulODE#addSecondaryEquations(SecondaryEquations)
162 */
163 public void registerVariationalEquations(final ExpandableStatefulODE expandable)
164 throws DimensionMismatchException, MismatchedEquations {
165
166 // safety checks
167 final FirstOrderDifferentialEquations ode = (jode instanceof MainStateJacobianWrapper) ?
168 ((MainStateJacobianWrapper) jode).ode :
169 jode;
170 if (expandable.getPrimary() != ode) {
171 throw new MismatchedEquations();
172 }
173
174 efode = expandable;
175 index = efode.addSecondaryEquations(new JacobiansSecondaryEquations());
176 efode.setSecondaryState(index, matricesData);
177
178 }
179
180 /** Add a parameter Jacobian provider.
181 * @param provider the parameter Jacobian provider to compute exactly the parameter Jacobian matrix
182 */
183 public void addParameterJacobianProvider(final ParameterJacobianProvider provider) {
184 jacobianProviders.add(provider);
185 }
186
187 /** Set a parameter Jacobian provider.
188 * @param parameterizedOde the parameterized ODE to compute the parameter Jacobian matrix using finite differences
189 */
190 public void setParameterizedODE(final ParameterizedODE parameterizedOde) {
191 this.pode = parameterizedOde;
192 dirtyParameter = true;
193 }
194
195 /** Set the step associated to a parameter in order to compute by finite
196 * difference the Jacobian matrix.
197 * <p>
198 * Needed if and only if the primary ODE set is a {@link ParameterizedODE}.
199 * </p>
200 * <p>
201 * Given a non zero parameter value pval for the parameter, a reasonable value
202 * for such a step is {@code pval * FastMath.sqrt(Precision.EPSILON)}.
203 * </p>
204 * <p>
205 * A zero value for such a step doesn't enable to compute the parameter Jacobian matrix.
206 * </p>
207 * @param parameter parameter to consider for Jacobian processing
208 * @param hP step for Jacobian finite difference computation w.r.t. the specified parameter
209 * @see ParameterizedODE
210 * @exception UnknownParameterException if the parameter is not supported
211 */
212 public void setParameterStep(final String parameter, final double hP)
213 throws UnknownParameterException {
214
215 for (ParameterConfiguration param: selectedParameters) {
216 if (parameter.equals(param.getParameterName())) {
217 param.setHP(hP);
218 dirtyParameter = true;
219 return;
220 }
221 }
222
223 throw new UnknownParameterException(parameter);
224
225 }
226
227 /** Set the initial value of the Jacobian matrix with respect to state.
228 * <p>
229 * If this method is not called, the initial value of the Jacobian
230 * matrix with respect to state is set to identity.
231 * </p>
232 * @param dYdY0 initial Jacobian matrix w.r.t. state
233 * @exception DimensionMismatchException if matrix dimensions are incorrect
234 */
235 public void setInitialMainStateJacobian(final double[][] dYdY0)
236 throws DimensionMismatchException {
237
238 // Check dimensions
239 checkDimension(stateDim, dYdY0);
240 checkDimension(stateDim, dYdY0[0]);
241
242 // store the matrix in row major order as a single dimension array
243 int i = 0;
244 for (final double[] row : dYdY0) {
245 System.arraycopy(row, 0, matricesData, i, stateDim);
246 i += stateDim;
247 }
248
249 if (efode != null) {
250 efode.setSecondaryState(index, matricesData);
251 }
252
253 }
254
255 /** Set the initial value of a column of the Jacobian matrix with respect to one parameter.
256 * <p>
257 * If this method is not called for some parameter, the initial value of
258 * the column of the Jacobian matrix with respect to this parameter is set to zero.
259 * </p>
260 * @param pName parameter name
261 * @param dYdP initial Jacobian column vector with respect to the parameter
262 * @exception UnknownParameterException if a parameter is not supported
263 * @throws DimensionMismatchException if the column vector does not match state dimension
264 */
265 public void setInitialParameterJacobian(final String pName, final double[] dYdP)
266 throws UnknownParameterException, DimensionMismatchException {
267
268 // Check dimensions
269 checkDimension(stateDim, dYdP);
270
271 // store the column in a global single dimension array
272 int i = stateDim * stateDim;
273 for (ParameterConfiguration param: selectedParameters) {
274 if (pName.equals(param.getParameterName())) {
275 System.arraycopy(dYdP, 0, matricesData, i, stateDim);
276 if (efode != null) {
277 efode.setSecondaryState(index, matricesData);
278 }
279 return;
280 }
281 i += stateDim;
282 }
283
284 throw new UnknownParameterException(pName);
285
286 }
287
288 /** Get the current value of the Jacobian matrix with respect to state.
289 * @param dYdY0 current Jacobian matrix with respect to state.
290 */
291 public void getCurrentMainSetJacobian(final double[][] dYdY0) {
292
293 // get current state for this set of equations from the expandable fode
294 double[] p = efode.getSecondaryState(index);
295
296 int j = 0;
297 for (int i = 0; i < stateDim; i++) {
298 System.arraycopy(p, j, dYdY0[i], 0, stateDim);
299 j += stateDim;
300 }
301
302 }
303
304 /** Get the current value of the Jacobian matrix with respect to one parameter.
305 * @param pName name of the parameter for the computed Jacobian matrix
306 * @param dYdP current Jacobian matrix with respect to the named parameter
307 */
308 public void getCurrentParameterJacobian(String pName, final double[] dYdP) {
309
310 // get current state for this set of equations from the expandable fode
311 double[] p = efode.getSecondaryState(index);
312
313 int i = stateDim * stateDim;
314 for (ParameterConfiguration param: selectedParameters) {
315 if (param.getParameterName().equals(pName)) {
316 System.arraycopy(p, i, dYdP, 0, stateDim);
317 return;
318 }
319 i += stateDim;
320 }
321
322 }
323
324 /** Check array dimensions.
325 * @param expected expected dimension
326 * @param array (may be null if expected is 0)
327 * @throws DimensionMismatchException if the array dimension does not match the expected one
328 */
329 private void checkDimension(final int expected, final Object array)
330 throws DimensionMismatchException {
331 int arrayDimension = (array == null) ? 0 : Array.getLength(array);
332 if (arrayDimension != expected) {
333 throw new DimensionMismatchException(arrayDimension, expected);
334 }
335 }
336
337 /** Local implementation of secondary equations.
338 * <p>
339 * This class is an inner class to ensure proper scheduling of calls
340 * by forcing the use of {@link JacobianMatrices#registerVariationalEquations(ExpandableStatefulODE)}.
341 * </p>
342 */
343 private class JacobiansSecondaryEquations implements SecondaryEquations {
344
345 /** {@inheritDoc} */
346 public int getDimension() {
347 return stateDim * (stateDim + paramDim);
348 }
349
350 /** {@inheritDoc} */
351 public void computeDerivatives(final double t, final double[] y, final double[] yDot,
352 final double[] z, final double[] zDot)
353 throws MaxCountExceededException, DimensionMismatchException {
354
355 // Lazy initialization
356 if (dirtyParameter && (paramDim != 0)) {
357 jacobianProviders.add(new ParameterJacobianWrapper(jode, pode, selectedParameters));
358 dirtyParameter = false;
359 }
360
361 // variational equations:
362 // from d[dy/dt]/dy0 and d[dy/dt]/dp to d[dy/dy0]/dt and d[dy/dp]/dt
363
364 // compute Jacobian matrix with respect to primary state
365 double[][] dFdY = new double[stateDim][stateDim];
366 jode.computeMainStateJacobian(t, y, yDot, dFdY);
367
368 // Dispatch Jacobian matrix in the compound secondary state vector
369 for (int i = 0; i < stateDim; ++i) {
370 final double[] dFdYi = dFdY[i];
371 for (int j = 0; j < stateDim; ++j) {
372 double s = 0;
373 final int startIndex = j;
374 int zIndex = startIndex;
375 for (int l = 0; l < stateDim; ++l) {
376 s += dFdYi[l] * z[zIndex];
377 zIndex += stateDim;
378 }
379 zDot[startIndex + i * stateDim] = s;
380 }
381 }
382
383 if (paramDim != 0) {
384 // compute Jacobian matrices with respect to parameters
385 double[] dFdP = new double[stateDim];
386 int startIndex = stateDim * stateDim;
387 for (ParameterConfiguration param: selectedParameters) {
388 boolean found = false;
389 for (int k = 0 ; (!found) && (k < jacobianProviders.size()); ++k) {
390 final ParameterJacobianProvider provider = jacobianProviders.get(k);
391 if (provider.isSupported(param.getParameterName())) {
392 provider.computeParameterJacobian(t, y, yDot,
393 param.getParameterName(), dFdP);
394 for (int i = 0; i < stateDim; ++i) {
395 final double[] dFdYi = dFdY[i];
396 int zIndex = startIndex;
397 double s = dFdP[i];
398 for (int l = 0; l < stateDim; ++l) {
399 s += dFdYi[l] * z[zIndex];
400 zIndex++;
401 }
402 zDot[startIndex + i] = s;
403 }
404 found = true;
405 }
406 }
407 if (! found) {
408 Arrays.fill(zDot, startIndex, startIndex + stateDim, 0.0);
409 }
410 startIndex += stateDim;
411 }
412 }
413
414 }
415 }
416
417 /** Wrapper class to compute jacobian matrices by finite differences for ODE
418 * which do not compute them by themselves.
419 */
420 private static class MainStateJacobianWrapper implements MainStateJacobianProvider {
421
422 /** Raw ODE without jacobians computation skill to be wrapped into a MainStateJacobianProvider. */
423 private final FirstOrderDifferentialEquations ode;
424
425 /** Steps for finite difference computation of the jacobian df/dy w.r.t. state. */
426 private final double[] hY;
427
428 /** Wrap a {@link FirstOrderDifferentialEquations} into a {@link MainStateJacobianProvider}.
429 * @param ode original ODE problem, without jacobians computation skill
430 * @param hY step sizes to compute the jacobian df/dy
431 * @exception DimensionMismatchException if there is a dimension mismatch between
432 * the steps array {@code hY} and the equation dimension
433 */
434 MainStateJacobianWrapper(final FirstOrderDifferentialEquations ode,
435 final double[] hY)
436 throws DimensionMismatchException {
437 this.ode = ode;
438 this.hY = hY.clone();
439 if (hY.length != ode.getDimension()) {
440 throw new DimensionMismatchException(ode.getDimension(), hY.length);
441 }
442 }
443
444 /** {@inheritDoc} */
445 public int getDimension() {
446 return ode.getDimension();
447 }
448
449 /** {@inheritDoc} */
450 public void computeDerivatives(double t, double[] y, double[] yDot)
451 throws MaxCountExceededException, DimensionMismatchException {
452 ode.computeDerivatives(t, y, yDot);
453 }
454
455 /** {@inheritDoc} */
456 public void computeMainStateJacobian(double t, double[] y, double[] yDot, double[][] dFdY)
457 throws MaxCountExceededException, DimensionMismatchException {
458
459 final int n = ode.getDimension();
460 final double[] tmpDot = new double[n];
461
462 for (int j = 0; j < n; ++j) {
463 final double savedYj = y[j];
464 y[j] += hY[j];
465 ode.computeDerivatives(t, y, tmpDot);
466 for (int i = 0; i < n; ++i) {
467 dFdY[i][j] = (tmpDot[i] - yDot[i]) / hY[j];
468 }
469 y[j] = savedYj;
470 }
471 }
472
473 }
474
475 /**
476 * Special exception for equations mismatch.
477 * @since 3.1
478 */
479 public static class MismatchedEquations extends MathIllegalArgumentException {
480
481 /** Serializable UID. */
482 private static final long serialVersionUID = 20120902L;
483
484 /** Simple constructor. */
485 public MismatchedEquations() {
486 super(LocalizedFormats.UNMATCHED_ODE_IN_EXPANDED_SET);
487 }
488
489 }
490
491}
492
Note: See TracBrowser for help on using the repository browser.