1 /*
2 * Licensed to the Hipparchus project 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 Hipparchus project 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 * https://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 */
17 package org.hipparchus.analysis.differentiation;
18
19 import org.hipparchus.exception.LocalizedCoreFormats;
20 import org.hipparchus.exception.MathIllegalArgumentException;
21 import org.hipparchus.linear.MatrixDecomposer;
22 import org.hipparchus.linear.MatrixUtils;
23 import org.hipparchus.linear.RealMatrix;
24 import org.hipparchus.util.MathUtils;
25
26 /** Container for a Taylor map.
27 * <p>
28 * A Taylor map is a set of n {@link DerivativeStructure}
29 * \((f_1, f_2, \ldots, f_n)\) depending on m parameters \((p_1, p_2, \ldots, p_m)\),
30 * with positive n and m.
31 * </p>
32 * @since 2.2
33 */
34 public class TaylorMap implements DifferentialAlgebra {
35
36 /** Evaluation point. */
37 private final double[] point;
38
39 /** Mapping functions. */
40 private final DerivativeStructure[] functions;
41
42 /** Simple constructor.
43 * <p>
44 * The number of number of parameters and derivation orders of all
45 * functions must match.
46 * </p>
47 * @param point point at which map is evaluated
48 * @param functions functions composing the map (must contain at least one element)
49 */
50 public TaylorMap(final double[] point, final DerivativeStructure[] functions) {
51 if (point == null || point.length == 0) {
52 throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_OF_ELEMENTS_SHOULD_BE_POSITIVE,
53 point == null ? 0 : point.length);
54 }
55 if (functions == null || functions.length == 0) {
56 throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_OF_ELEMENTS_SHOULD_BE_POSITIVE,
57 functions == null ? 0 : functions.length);
58 }
59 this.point = point.clone();
60 this.functions = functions.clone();
61 final DSFactory factory0 = functions[0].getFactory();
62 MathUtils.checkDimension(point.length, factory0.getCompiler().getFreeParameters());
63 for (int i = 1; i < functions.length; ++i) {
64 factory0.checkCompatibility(functions[i].getFactory());
65 }
66 }
67
68 /** Constructor for identity map.
69 * <p>
70 * The identity is considered to be evaluated at origin.
71 * </p>
72 * @param parameters number of free parameters
73 * @param order derivation order
74 * @param nbFunctions number of functions
75 */
76 public TaylorMap(final int parameters, final int order, final int nbFunctions) {
77 this(parameters, nbFunctions);
78 final DSFactory factory = new DSFactory(parameters, order);
79 for (int i = 0; i < nbFunctions; ++i) {
80 functions[i] = factory.variable(i, 0.0);
81 }
82 }
83
84 /** Build an empty map evaluated at origin.
85 * @param parameters number of free parameters
86 * @param nbFunctions number of functions
87 */
88 private TaylorMap(final int parameters, final int nbFunctions) {
89 this.point = new double[parameters];
90 this.functions = new DerivativeStructure[nbFunctions];
91 }
92
93 /** {@inheritDoc} */
94 @Override
95 public int getFreeParameters() {
96 return point.length;
97 }
98
99 /** {@inheritDoc} */
100 @Override
101 public int getOrder() {
102 return functions[0].getOrder();
103 }
104
105 /** Get the number of parameters of the map.
106 * @return number of parameters of the map
107 */
108 @Deprecated
109 public int getNbParameters() {
110 return getFreeParameters();
111 }
112
113 /** Get the number of functions of the map.
114 * @return number of functions of the map
115 */
116 public int getNbFunctions() {
117 return functions.length;
118 }
119
120 /** Get the point at which map is evaluated.
121 * @return point at which map is evaluated
122 */
123 public double[] getPoint() {
124 return point.clone();
125 }
126
127 /** Get a function from the map.
128 * @param i index of the function (must be between 0 included and {@link #getNbFunctions()} excluded
129 * @return function at index i
130 */
131 public DerivativeStructure getFunction(final int i) {
132 return functions[i];
133 }
134
135 /** Subtract two maps.
136 * @param map map to subtract from instance
137 * @return this - map
138 */
139 private TaylorMap subtract(final TaylorMap map) {
140 final TaylorMap result = new TaylorMap(point.length, functions.length);
141 for (int i = 0; i < result.functions.length; ++i) {
142 result.functions[i] = functions[i].subtract(map.functions[i]);
143 }
144 return result;
145 }
146
147 /** Evaluate Taylor expansion of the map at some offset.
148 * @param deltaP parameters offsets \((\Delta p_1, \Delta p_2, \ldots, \Delta p_n)\)
149 * @return value of the Taylor expansion at \((p_1 + \Delta p_1, p_2 + \Delta p_2, \ldots, p_n + \Delta p_n)\)
150 */
151 public double[] value(final double... deltaP) {
152 final double[] value = new double[functions.length];
153 for (int i = 0; i < functions.length; ++i) {
154 value[i] = functions[i].taylor(deltaP);
155 }
156 return value;
157 }
158
159 /** Compose the instance with another Taylor map as \(\mathrm{this} \circ \mathrm{other}\).
160 * @param other map with which instance must be composed
161 * @return composed map \(\mathrm{this} \circ \mathrm{other}\)
162 */
163 public TaylorMap compose(final TaylorMap other) {
164
165 // safety check
166 MathUtils.checkDimension(getFreeParameters(), other.getNbFunctions());
167
168 final DerivativeStructure[] composed = new DerivativeStructure[functions.length];
169 for (int i = 0; i < functions.length; ++i) {
170 composed[i] = functions[i].rebase(other.functions);
171 }
172
173 return new TaylorMap(other.point, composed);
174
175 }
176
177 /** Invert the instance.
178 * <p>
179 * Consider {@link #value(double[]) Taylor expansion} of the map with
180 * small parameters offsets \((\Delta p_1, \Delta p_2, \ldots, \Delta p_n)\)
181 * which leads to evaluation offsets \((f_1 + df_1, f_2 + df_2, \ldots, f_n + df_n)\).
182 * The map inversion defines a Taylor map that computes \((\Delta p_1,
183 * \Delta p_2, \ldots, \Delta p_n)\) from \((df_1, df_2, \ldots, df_n)\).
184 * </p>
185 * <p>
186 * The map must be square to be invertible (i.e. the number of functions and the
187 * number of parameters in the functions must match)
188 * </p>
189 * @param decomposer matrix decomposer to user for inverting the linear part
190 * @return inverted map
191 * @see <a href="https://doi.org/10.1016/S1076-5670(08)70228-3">chapter
192 * 2 of Advances in Imaging and Electron Physics, vol 108
193 * by Martin Berz</a>
194 */
195 public TaylorMap invert(final MatrixDecomposer decomposer) {
196
197 final DSFactory factory = functions[0].getFactory();
198 final DSCompiler compiler = factory.getCompiler();
199 final int n = functions.length;
200
201 // safety check
202 MathUtils.checkDimension(n, functions[0].getFreeParameters());
203
204 // set up an indirection array between linear terms and complete derivatives arrays
205 final int[] indirection = new int[n];
206 int linearIndex = 0;
207 for (int k = 1; linearIndex < n; ++k) {
208 if (compiler.getPartialDerivativeOrdersSum(k) == 1) {
209 indirection[linearIndex++] = k;
210 }
211 }
212
213 // separate linear and non-linear terms
214 final RealMatrix linear = MatrixUtils.createRealMatrix(n, n);
215 final TaylorMap nonLinearTM = new TaylorMap(n, n);
216 for (int i = 0; i < n; ++i) {
217 nonLinearTM.functions[i] = factory.build(functions[i].getAllDerivatives());
218 nonLinearTM.functions[i].setDerivativeComponent(0, 0.0);
219 for (int j = 0; j < n; ++j) {
220 final int k = indirection[j];
221 linear.setEntry(i, j, functions[i].getDerivativeComponent(k));
222 nonLinearTM.functions[i].setDerivativeComponent(k, 0.0);
223 }
224 }
225
226 // invert the linear part
227 final RealMatrix linearInvert = decomposer.decompose(linear).getInverse();
228
229 // convert the invert of linear part back to a Taylor map
230 final TaylorMap linearInvertTM = new TaylorMap(n, n);
231 for (int i = 0; i < n; ++i) {
232 linearInvertTM.functions[i] = new DerivativeStructure(factory);
233 for (int j = 0; j < n; ++j) {
234 linearInvertTM.functions[i].setDerivativeComponent(indirection[j], linearInvert.getEntry(i, j));
235 }
236 }
237
238 // perform fixed-point evaluation of the inverse
239 // adding one derivation order at each iteration
240 final TaylorMap identity = new TaylorMap(n, compiler.getOrder(), n);
241 TaylorMap invertTM = linearInvertTM;
242 for (int k = 1; k < compiler.getOrder(); ++k) {
243 invertTM = linearInvertTM.compose(identity.subtract(nonLinearTM.compose(invertTM)));
244 }
245
246 // set the constants
247 for (int i = 0; i < n; ++i) {
248 invertTM.point[i] = functions[i].getValue();
249 invertTM.functions[i].setDerivativeComponent(0, point[i]);
250 }
251
252 return invertTM;
253
254 }
255
256 }