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