Source code

001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017package org.apache.commons.math4.legacy.ode;
018
019import java.lang.reflect.Array;
020import java.util.ArrayList;
021import java.util.Arrays;
022import java.util.List;
023
024import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
025import org.apache.commons.math4.legacy.exception.MathIllegalArgumentException;
026import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
027import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
028
029/**
030 * This class defines a set of {@link SecondaryEquations secondary equations} to
031 * compute the Jacobian matrices with respect to the initial state vector and, if
032 * any, to some parameters of the primary ODE set.
033 * <p>
034 * It is intended to be packed into an {@link ExpandableStatefulODE}
035 * in conjunction with a primary set of ODE, which may be:
036 * <ul>
037 * <li>a {@link FirstOrderDifferentialEquations}</li>
038 * <li>a {@link MainStateJacobianProvider}</li>
039 * </ul>
040 * In order to compute Jacobian matrices with respect to some parameters of the
041 * primary ODE set, the following parameter Jacobian providers may be set:
042 * <ul>
043 * <li>a {@link ParameterJacobianProvider}</li>
044 * <li>a {@link ParameterizedODE}</li>
045 * </ul>
046 *
047 * @see ExpandableStatefulODE
048 * @see FirstOrderDifferentialEquations
049 * @see MainStateJacobianProvider
050 * @see ParameterJacobianProvider
051 * @see ParameterizedODE
052 *
053 * @since 3.0
054 */
055public class JacobianMatrices {
056
057    /** Expandable first order differential equation. */
058    private ExpandableStatefulODE efode;
059
060    /** Index of the instance in the expandable set. */
061    private int index;
062
063    /** FODE with exact primary Jacobian computation skill. */
064    private MainStateJacobianProvider jode;
065
066    /** FODE without exact parameter Jacobian computation skill. */
067    private ParameterizedODE pode;
068
069    /** Main state vector dimension. */
070    private int stateDim;
071
072    /** Selected parameters for parameter Jacobian computation. */
073    private ParameterConfiguration[] selectedParameters;
074
075    /** FODE with exact parameter Jacobian computation skill. */
076    private List<ParameterJacobianProvider> jacobianProviders;
077
078    /** Parameters dimension. */
079    private int paramDim;
080
081    /** Boolean for selected parameters consistency. */
082    private boolean dirtyParameter;
083
084    /** State and parameters Jacobian matrices in a row. */
085    private double[] matricesData;
086
087    /** Simple constructor for a secondary equations set computing Jacobian matrices.
088     * <p>
089     * Parameters must belong to the supported ones given by {@link
090     * Parameterizable#getParametersNames()}, so the primary set of differential
091     * equations must be {@link Parameterizable}.
092     * </p>
093     * <p>Note that each selection clears the previous selected parameters.</p>
094     *
095     * @param fode the primary first order differential equations set to extend
096     * @param hY step used for finite difference computation with respect to state vector
097     * @param parameters parameters to consider for Jacobian matrices processing
098     * (may be null if parameters Jacobians is not desired)
099     * @exception DimensionMismatchException if there is a dimension mismatch between
100     * the steps array {@code hY} and the equation dimension
101     */
102    public JacobianMatrices(final FirstOrderDifferentialEquations fode, final double[] hY,
103                            final String... parameters)
104        throws DimensionMismatchException {
105        this(new MainStateJacobianWrapper(fode, hY), parameters);
106    }
107
108    /** Simple constructor for a secondary equations set computing Jacobian matrices.
109     * <p>
110     * Parameters must belong to the supported ones given by {@link
111     * Parameterizable#getParametersNames()}, so the primary set of differential
112     * equations must be {@link Parameterizable}.
113     * </p>
114     * <p>Note that each selection clears the previous selected parameters.</p>
115     *
116     * @param jode the primary first order differential equations set to extend
117     * @param parameters parameters to consider for Jacobian matrices processing
118     * (may be null if parameters Jacobians is not desired)
119     */
120    public JacobianMatrices(final MainStateJacobianProvider jode,
121                            final String... parameters) {
122
123        this.efode = null;
124        this.index = -1;
125
126        this.jode = jode;
127        this.pode = null;
128
129        this.stateDim = jode.getDimension();
130
131        if (parameters == null) {
132            selectedParameters = null;
133            paramDim = 0;
134        } else {
135            this.selectedParameters = new ParameterConfiguration[parameters.length];
136            for (int i = 0; i < parameters.length; ++i) {
137                selectedParameters[i] = new ParameterConfiguration(parameters[i], Double.NaN);
138            }
139            paramDim = parameters.length;
140        }
141        this.dirtyParameter = false;
142
143        this.jacobianProviders = new ArrayList<>();
144
145        // set the default initial state Jacobian to the identity
146        // and the default initial parameters Jacobian to the null matrix
147        matricesData = new double[(stateDim + paramDim) * stateDim];
148        for (int i = 0; i < stateDim; ++i) {
149            matricesData[i * (stateDim + 1)] = 1.0;
150        }
151    }
152
153    /** Register the variational equations for the Jacobians matrices to the expandable set.
154     * @param expandable expandable set into which variational equations should be registered
155     * @throws DimensionMismatchException if the dimension of the partial state does not
156     * match the selected equations set dimension
157     * @exception MismatchedEquations if the primary set of the expandable set does
158     * not match the one used to build the instance
159     * @see ExpandableStatefulODE#addSecondaryEquations(SecondaryEquations)
160     */
161    public void registerVariationalEquations(final ExpandableStatefulODE expandable)
162        throws DimensionMismatchException, MismatchedEquations {
163
164        // safety checks
165        final FirstOrderDifferentialEquations ode = (jode instanceof MainStateJacobianWrapper) ?
166                                                    ((MainStateJacobianWrapper) jode).ode :
167                                                    jode;
168        if (expandable.getPrimary() != ode) {
169            throw new MismatchedEquations();
170        }
171
172        efode = expandable;
173        index = efode.addSecondaryEquations(new JacobiansSecondaryEquations());
174        efode.setSecondaryState(index, matricesData);
175    }
176
177    /** Add a parameter Jacobian provider.
178     * @param provider the parameter Jacobian provider to compute exactly the parameter Jacobian matrix
179     */
180    public void addParameterJacobianProvider(final ParameterJacobianProvider provider) {
181        jacobianProviders.add(provider);
182    }
183
184    /** Set a parameter Jacobian provider.
185     * @param parameterizedOde the parameterized ODE to compute the parameter Jacobian matrix using finite differences
186     */
187    public void setParameterizedODE(final ParameterizedODE parameterizedOde) {
188        this.pode = parameterizedOde;
189        dirtyParameter = true;
190    }
191
192    /** Set the step associated to a parameter in order to compute by finite
193     *  difference the Jacobian matrix.
194     * <p>
195     * Needed if and only if the primary ODE set is a {@link ParameterizedODE}.
196     * </p>
197     * <p>
198     * Given a non zero parameter value pval for the parameter, a reasonable value
199     * for such a step is {@code pval * JdkMath.sqrt(Precision.EPSILON)}.
200     * </p>
201     * <p>
202     * A zero value for such a step doesn't enable to compute the parameter Jacobian matrix.
203     * </p>
204     * @param parameter parameter to consider for Jacobian processing
205     * @param hP step for Jacobian finite difference computation w.r.t. the specified parameter
206     * @see ParameterizedODE
207     * @exception UnknownParameterException if the parameter is not supported
208     */
209    public void setParameterStep(final String parameter, final double hP)
210        throws UnknownParameterException {
211
212        for (ParameterConfiguration param: selectedParameters) {
213            if (parameter.equals(param.getParameterName())) {
214                param.setHP(hP);
215                dirtyParameter = true;
216                return;
217            }
218        }
219
220        throw new UnknownParameterException(parameter);
221    }
222
223    /** Set the initial value of the Jacobian matrix with respect to state.
224     * <p>
225     * If this method is not called, the initial value of the Jacobian
226     * matrix with respect to state is set to identity.
227     * </p>
228     * @param dYdY0 initial Jacobian matrix w.r.t. state
229     * @exception DimensionMismatchException if matrix dimensions are incorrect
230     */
231    public void setInitialMainStateJacobian(final double[][] dYdY0)
232        throws DimensionMismatchException {
233
234        // Check dimensions
235        checkDimension(stateDim, dYdY0);
236        checkDimension(stateDim, dYdY0[0]);
237
238        // store the matrix in row major order as a single dimension array
239        int i = 0;
240        for (final double[] row : dYdY0) {
241            System.arraycopy(row, 0, matricesData, i, stateDim);
242            i += stateDim;
243        }
244
245        if (efode != null) {
246            efode.setSecondaryState(index, matricesData);
247        }
248    }
249
250    /** Set the initial value of a column of the Jacobian matrix with respect to one parameter.
251     * <p>
252     * If this method is not called for some parameter, the initial value of
253     * the column of the Jacobian matrix with respect to this parameter is set to zero.
254     * </p>
255     * @param pName parameter name
256     * @param dYdP initial Jacobian column vector with respect to the parameter
257     * @exception UnknownParameterException if a parameter is not supported
258     * @throws DimensionMismatchException if the column vector does not match state dimension
259     */
260    public void setInitialParameterJacobian(final String pName, final double[] dYdP)
261        throws UnknownParameterException, DimensionMismatchException {
262
263        // Check dimensions
264        checkDimension(stateDim, dYdP);
265
266        // store the column in a global single dimension array
267        int i = stateDim * stateDim;
268        for (ParameterConfiguration param: selectedParameters) {
269            if (pName.equals(param.getParameterName())) {
270                System.arraycopy(dYdP, 0, matricesData, i, stateDim);
271                if (efode != null) {
272                    efode.setSecondaryState(index, matricesData);
273                }
274                return;
275            }
276            i += stateDim;
277        }
278
279        throw new UnknownParameterException(pName);
280    }
281
282    /** Get the current value of the Jacobian matrix with respect to state.
283     * @param dYdY0 current Jacobian matrix with respect to state.
284     */
285    public void getCurrentMainSetJacobian(final double[][] dYdY0) {
286
287        // get current state for this set of equations from the expandable fode
288        double[] p = efode.getSecondaryState(index);
289
290        int j = 0;
291        for (int i = 0; i < stateDim; i++) {
292            System.arraycopy(p, j, dYdY0[i], 0, stateDim);
293            j += stateDim;
294        }
295    }
296
297    /** Get the current value of the Jacobian matrix with respect to one parameter.
298     * @param pName name of the parameter for the computed Jacobian matrix
299     * @param dYdP current Jacobian matrix with respect to the named parameter
300     */
301    public void getCurrentParameterJacobian(String pName, final double[] dYdP) {
302
303        // get current state for this set of equations from the expandable fode
304        double[] p = efode.getSecondaryState(index);
305
306        int i = stateDim * stateDim;
307        for (ParameterConfiguration param: selectedParameters) {
308            if (param.getParameterName().equals(pName)) {
309                System.arraycopy(p, i, dYdP, 0, stateDim);
310                return;
311            }
312            i += stateDim;
313        }
314    }
315
316    /** Check array dimensions.
317     * @param expected expected dimension
318     * @param array (may be null if expected is 0)
319     * @throws DimensionMismatchException if the array dimension does not match the expected one
320     */
321    private void checkDimension(final int expected, final Object array)
322        throws DimensionMismatchException {
323        int arrayDimension = (array == null) ? 0 : Array.getLength(array);
324        if (arrayDimension != expected) {
325            throw new DimensionMismatchException(arrayDimension, expected);
326        }
327    }
328
329    /** Local implementation of secondary equations.
330     * <p>
331     * This class is an inner class to ensure proper scheduling of calls
332     * by forcing the use of {@link JacobianMatrices#registerVariationalEquations(ExpandableStatefulODE)}.
333     * </p>
334     */
335    private class JacobiansSecondaryEquations implements SecondaryEquations {
336
337        /** {@inheritDoc} */
338        @Override
339        public int getDimension() {
340            return stateDim * (stateDim + paramDim);
341        }
342
343        /** {@inheritDoc} */
344        @Override
345        public void computeDerivatives(final double t, final double[] y, final double[] yDot,
346                                       final double[] z, final double[] zDot)
347            throws MaxCountExceededException, DimensionMismatchException {
348
349            // Lazy initialization
350            if (dirtyParameter && paramDim != 0) {
351                jacobianProviders.add(new ParameterJacobianWrapper(jode, pode, selectedParameters));
352                dirtyParameter = false;
353            }
354
355            // variational equations:
356            // from d[dy/dt]/dy0 and d[dy/dt]/dp to d[dy/dy0]/dt and d[dy/dp]/dt
357
358            // compute Jacobian matrix with respect to primary state
359            double[][] dFdY = new double[stateDim][stateDim];
360            jode.computeMainStateJacobian(t, y, yDot, dFdY);
361
362            // Dispatch Jacobian matrix in the compound secondary state vector
363            for (int i = 0; i < stateDim; ++i) {
364                final double[] dFdYi = dFdY[i];
365                for (int j = 0; j < stateDim; ++j) {
366                    double s = 0;
367                    final int startIndex = j;
368                    int zIndex = startIndex;
369                    for (int l = 0; l < stateDim; ++l) {
370                        s += dFdYi[l] * z[zIndex];
371                        zIndex += stateDim;
372                    }
373                    zDot[startIndex + i * stateDim] = s;
374                }
375            }
376
377            if (paramDim != 0) {
378                // compute Jacobian matrices with respect to parameters
379                double[] dFdP = new double[stateDim];
380                int startIndex = stateDim * stateDim;
381                for (ParameterConfiguration param: selectedParameters) {
382                    boolean found = false;
383                    for (int k = 0 ; !found && k < jacobianProviders.size(); ++k) {
384                        final ParameterJacobianProvider provider = jacobianProviders.get(k);
385                        if (provider.isSupported(param.getParameterName())) {
386                            provider.computeParameterJacobian(t, y, yDot,
387                                                              param.getParameterName(), dFdP);
388                            for (int i = 0; i < stateDim; ++i) {
389                                final double[] dFdYi = dFdY[i];
390                                int zIndex = startIndex;
391                                double s = dFdP[i];
392                                for (int l = 0; l < stateDim; ++l) {
393                                    s += dFdYi[l] * z[zIndex];
394                                    zIndex++;
395                                }
396                                zDot[startIndex + i] = s;
397                            }
398                            found = true;
399                        }
400                    }
401                    if (! found) {
402                        Arrays.fill(zDot, startIndex, startIndex + stateDim, 0.0);
403                    }
404                    startIndex += stateDim;
405                }
406            }
407        }
408    }
409
410    /** Wrapper class to compute jacobian matrices by finite differences for ODE
411     *  which do not compute them by themselves.
412     */
413    private static class MainStateJacobianWrapper implements MainStateJacobianProvider {
414
415        /** Raw ODE without jacobians computation skill to be wrapped into a MainStateJacobianProvider. */
416        private final FirstOrderDifferentialEquations ode;
417
418        /** Steps for finite difference computation of the jacobian df/dy w.r.t. state. */
419        private final double[] hY;
420
421        /** Wrap a {@link FirstOrderDifferentialEquations} into a {@link MainStateJacobianProvider}.
422         * @param ode original ODE problem, without jacobians computation skill
423         * @param hY step sizes to compute the jacobian df/dy
424         * @exception DimensionMismatchException if there is a dimension mismatch between
425         * the steps array {@code hY} and the equation dimension
426         */
427        MainStateJacobianWrapper(final FirstOrderDifferentialEquations ode,
428                                 final double[] hY)
429            throws DimensionMismatchException {
430            this.ode = ode;
431            this.hY = hY.clone();
432            if (hY.length != ode.getDimension()) {
433                throw new DimensionMismatchException(ode.getDimension(), hY.length);
434            }
435        }
436
437        /** {@inheritDoc} */
438        @Override
439        public int getDimension() {
440            return ode.getDimension();
441        }
442
443        /** {@inheritDoc} */
444        @Override
445        public void computeDerivatives(double t, double[] y, double[] yDot)
446            throws MaxCountExceededException, DimensionMismatchException {
447            ode.computeDerivatives(t, y, yDot);
448        }
449
450        /** {@inheritDoc} */
451        @Override
452        public void computeMainStateJacobian(double t, double[] y, double[] yDot, double[][] dFdY)
453            throws MaxCountExceededException, DimensionMismatchException {
454
455            final int n = ode.getDimension();
456            final double[] tmpDot = new double[n];
457
458            for (int j = 0; j < n; ++j) {
459                final double savedYj = y[j];
460                y[j] += hY[j];
461                ode.computeDerivatives(t, y, tmpDot);
462                for (int i = 0; i < n; ++i) {
463                    dFdY[i][j] = (tmpDot[i] - yDot[i]) / hY[j];
464                }
465                y[j] = savedYj;
466            }
467        }
468    }
469
470    /**
471     * Special exception for equations mismatch.
472     * @since 3.1
473     */
474    public static class MismatchedEquations extends MathIllegalArgumentException {
475
476        /** Serializable UID. */
477        private static final long serialVersionUID = 20120902L;
478
479        /** Simple constructor. */
480        public MismatchedEquations() {
481            super(LocalizedFormats.UNMATCHED_ODE_IN_EXPANDED_SET);
482        }
483    }
484}
485