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