Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
TMultiDimFit.cxx
Go to the documentation of this file.
1// @(#)root/hist:$Id$
2// Author: Christian Holm Christensen 07/11/2000
3
4/** \class TMultiDimFit
5 \ingroup Hist
6
7 Multidimensional Fits in ROOT.
8 ## Overview
9 A common problem encountered in different fields of applied science is
10 to find an expression for one physical quantity in terms of several
11 others, which are directly measurable.
12
13 An example in high energy physics is the evaluation of the momentum of
14 a charged particle from the observation of its trajectory in a magnetic
15 field. The problem is to relate the momentum of the particle to the
16 observations, which may consists of positional measurements at
17 intervals along the particle trajectory.
18
19 The exact functional relationship between the measured quantities
20 (e.g., the space-points) and the dependent quantity (e.g., the
21 momentum) is in general not known, but one possible way of solving the
22 problem, is to find an expression which reliably approximates the
23 dependence of the momentum on the observations.
24
25 This explicit function of the observations can be obtained by a
26 <I>least squares</I> fitting procedure applied to a representative
27 sample of the data, for which the dependent quantity (e.g., momentum)
28 and the independent observations are known. The function can then be
29 used to compute the quantity of interest for new observations of the
30 independent variables.
31
32 This class <TT>TMultiDimFit</TT> implements such a procedure in
33 ROOT. It is largely based on the CERNLIB MUDIFI package [2].
34 Though the basic concepts are still sound, and
35 therefore kept, a few implementation details have changed, and this
36 class can take advantage of MINUIT [4] to improve the errors
37 of the fitting, thanks to the class TMinuit.
38
39 In [5] and [6] H. Wind demonstrates the utility
40 of this procedure in the context of tracking, magnetic field
41 parameterisation, and so on. The outline of the method used in this
42 class is based on Winds discussion, and I refer these two excellents
43 text for more information.
44
45 And example of usage is given in multidimfit.C.
46
47 ## The Method
48 Let \f$ D \f$ by the dependent quantity of interest, which depends smoothly
49 on the observable quantities \f$ x_1, \ldots, x_N \f$ which we'll denote by
50 \f$\mathbf{x}\f$. Given a training sample of \f$ M\f$ tuples of the form, (TMultiDimFit::AddRow)
51
52 \f[
53 \left(\mathbf{x}_j, D_j, E_j\right)\quad,
54 \f]
55 where \f$\mathbf{x}_j = (x_{1,j},\ldots,x_{N,j})\f$ are \f$ N\f$ independent
56 variables, \f$ D_j\f$ is the known, quantity dependent at \f$\mathbf{x}_j\f$ and \f$ E_j\f$ is
57 the square error in \f$ D_j\f$, the class will try to find the parameterization
58 \f[
59 D_p(\mathbf{x}) = \sum_{l=1}^{L} c_l \prod_{i=1}^{N} p_{li}\left(x_i\right)
60 = \sum_{l=1}^{L} c_l F_l(\mathbf{x})
61 \f]
62 such that
63
64 \f[
65 S \equiv \sum_{j=1}^{M} \left(D_j - D_p\left(\mathbf{x}_j\right)\right)^2
66 \f]
67 is minimal. Here \f$p_{li}(x_i)\f$ are monomials, or Chebyshev or Legendre
68 polynomials, labelled \f$l = 1, \ldots, L\f$, in each variable \f$ x_i\f$,\f$ i=1, \ldots, N\f$.
69
70 So what TMultiDimFit does, is to determine the number of terms \f$ L\f$, and then \f$ L\f$ terms
71 (or functions) \f$ F_l\f$, and the \f$ L\f$ coefficients \f$ c_l\f$, so that \f$ S\f$ is minimal
72 (TMultiDimFit::FindParameterization).
73
74 Of course it's more than a little unlikely that \f$ S\f$ will ever become
75 exact zero as a result of the procedure outlined below. Therefore, the
76 user is asked to provide a minimum relative error \f$ \epsilon\f$ (TMultiDimFit::SetMinRelativeError),
77 and \f$ S\f$ will be considered minimized when
78
79 \f[
80 R = \frac{S}{\sum_{j=1}^M D_j^2} < \epsilon
81 \f]
82 Optionally, the user may impose a functional expression by specifying
83 the powers of each variable in \f$ L\f$ specified functions \f$ F_1, \ldots,F_L\f$ (TMultiDimFit::SetPowers).
84 In that case, only the coefficients \f$ c_l\f$ is calculated by the class.
85
86 ## Limiting the Number of Terms
87 As always when dealing with fits, there's a real chance of *over fitting*. As is well-known, it's
88 always possible to fit an \f$ N-1\f$ polynomial in \f$ x\f$ to \f$ N\f$ points \f$ (x,y)\f$ with
89 \f$\chi^2 = 0\f$, but the polynomial is not likely to fit new data at all [1].
90 Therefore, the user is asked to provide an upper limit, \f$ L_{max}\f$ to the number of terms in
91 \f$ D_p\f$ (TMultiDimFit::SetMaxTerms).
92
93 However, since there's an infinite number of \f$ F_l\f$ to choose from, the
94 user is asked to give the maximum power. \f$ P_{max,i}\f$, of each variable
95 \f$ x_i\f$ to be considered in the minimization of \f$ S\f$ (TMultiDimFit::SetMaxPowers).
96
97 One way of obtaining values for the maximum power in variable \f$ i\f$, is
98 to perform a regular fit to the dependent quantity \f$ D\f$, using a
99 polynomial only in \f$ x_i\f$. The maximum power is \f$ P_{max,i}\f$ is then the
100 power that does not significantly improve the one-dimensional
101 least-square fit over \f$ x_i\f$ to \f$ D\f$ [5].
102
103 There are still a huge amount of possible choices for \f$ F_l\f$; in fact
104 there are \f$\prod_{i=1}^{N} (P_{max,i} + 1)\f$ possible
105 choices. Obviously we need to limit this. To this end, the user is
106 asked to set a *power control limit*, \f$ Q\f$ (TMultiDimFit::SetPowerLimit), and a function
107 \f$ F_l\f$ is only accepted if
108 \f[
109 Q_l = \sum_{i=1}^{N} \frac{P_{li}}{P_{max,i}} < Q
110 \f]
111 where \f$ P_{li}\f$ is the leading power of variable \f$ x_i\f$ in function \f$ F_l\f$ (TMultiDimFit::MakeCandidates).
112 So the number of functions increase with \f$ Q\f$ (1, 2 is fine, 5 is way out).
113
114 ## Gram-Schmidt Orthogonalisation
115 To further reduce the number of functions in the final expression,
116 only those functions that significantly reduce \f$ S\f$ is chosen. What
117 `significant' means, is chosen by the user, and will be
118 discussed below (see [2.3](TMultiFimFit.html#sec:selectiondetail)).
119
120 The functions \f$ F_l\f$ are generally not orthogonal, which means one will
121 have to evaluate all possible \f$ F_l\f$'s over all data-points before
122 finding the most significant [1]. We can, however, do
123 better then that. By applying the *modified Gram-Schmidt
124 orthogonalisation* algorithm [5] [3] to the
125 functions \f$ F_l\f$, we can evaluate the contribution to the reduction of
126 \f$ S\f$ from each function in turn, and we may delay the actual inversion
127 of the curvature-matrix (TMultiDimFit::MakeGramSchmidt).
128
129 So we are let to consider an \f$ M\times L\f$ matrix \f$\mathsf{F}\f$, an
130 element of which is given by
131 \f[
132 f_{jl} = F_j\left(x_{1j} , x_{2j}, \ldots, x_{Nj}\right)
133 = F_l(\mathbf{x}_j)\, \quad\mbox{with}~j=1,2,\ldots,M,
134 \f]
135 where \f$ j\f$ labels the \f$ M\f$ rows in the training sample and \f$ l\f$ labels
136 \f$ L\f$ functions of \f$ N\f$ variables, and \f$ L \leq M\f$. That is, \f$ f_{jl}\f$ is
137 the term (or function) numbered \f$ l\f$ evaluated at the data point
138 \f$ j\f$. We have to normalise \f$\mathbf{x}_j\f$ to \f$ [-1,1]\f$ for this to
139 succeed [5] (TMultiDimFit::MakeNormalized). We then define a
140 matrix \f$\mathsf{W}\f$ of which the columns \f$\mathbf{w}_j\f$ are given by
141 \f{eqnarray*}{
142 \mathbf{w}_1 &=& \mathbf{f}_1 = F_1\left(\mathbf x_1\right)\\
143 \mathbf{w}_l &=& \mathbf{f}_l - \sum^{l-1}_{k=1} \frac{\mathbf{f}_l \bullet
144 \mathbf{w}_k}{\mathbf{w}_k^2}\mathbf{w}_k\,.
145 \f}
146 and \f$\mathbf{w}_{l}\f$ is the component of \f$\mathbf{f}_{l} \f$ orthogonal
147 to \f$\mathbf{w}_{1}, \ldots, \mathbf{w}_{l-1}\f$. Hence we obtain [3],
148 \f[
149 \mathbf{w}_k\bullet\mathbf{w}_l = 0\quad\mbox{if}~k \neq l\quad.
150 \f]
151 We now take as a new model \f$\mathsf{W}\mathbf{a}\f$. We thus want to
152 minimize
153 \f[
154 S\equiv \left(\mathbf{D} - \mathsf{W}\mathbf{a}\right)^2\quad,
155 \f]
156 where \f$\mathbf{D} = \left(D_1,\ldots,D_M\right)\f$ is a vector of the
157 dependent quantity in the sample. Differentiation with respect to
158 \f$ a_j\f$ gives, using [6], <a name="eq:dS2"></a>
159 \f[
160 \mathbf{D}\bullet\mathbf{w}_l - a_l\mathbf{w}_l^2 = 0
161 \f]
162 or
163 \f[
164 a_l = \frac{\mathbf{D}_l\bullet\mathbf{w}_l}{\mathbf{w}_l^2}
165 \f]
166 Let \f$ S_j\f$ be the sum of squares of residuals when taking \f$ j\f$ functions
167 into account. Then
168 \f[
169 S_l = \left[\mathbf{D} - \sum^l_{k=1} a_k\mathbf{w}_k\right]^2
170 = \mathbf{D}^2 - 2\mathbf{D} \sum^l_{k=1} a_k\mathbf{w}_k
171 + \sum^l_{k=1} a_k^2\mathbf{w}_k^2
172 \f]
173 Using [9], we see that
174 \f{eqnarray*}{
175 S_l &=& \mathbf{D}^2 - 2 \sum^l_{k=1} a_k^2\mathbf{w}_k^2 +
176 \sum^j_{k=1} a_k^2\mathbf{w}_k^2\nonumber\\
177 &=& \mathbf{D}^2 - \sum^l_{k=1} a_k^2\mathbf{w}_k^2\nonumber\\
178 &=& \mathbf{D}^2 - \sum^l_{k=1} \frac{\left(\mathbf D\bullet \mathbf
179 w_k\right)}{\mathbf w_k^2}
180 \f}
181 So for each new function \f$ F_l\f$ included in the model, we get a
182 reduction of the sum of squares of residuals of \f$a_l^2\mathbf{w}_l^2\f$,
183 where \f$\mathbf{w}_l\f$ is given by [4] and \f$ a_l\f$ by [9]. Thus, using
184 the Gram-Schmidt orthogonalisation, we
185 can decide if we want to include this function in the final model,
186 *before* the matrix inversion.
187
188 ## Function Selection Based on Residual
189 Supposing that \f$ L-1\f$ steps of the procedure have been performed, the
190 problem now is to consider the \f$L^{\mbox{th}}\f$ function.
191
192 The sum of squares of residuals can be written as
193 \f[
194 S_L = \textbf{D}^T\bullet\textbf{D} -
195 \sum^L_{l=1}a^2_l\left(\textbf{w}_l^T\bullet\textbf{w}_l\right)
196 \f]
197 where the relation [9] have been taken into account. The
198 contribution of the \f$L^{\mbox{th}}\f$ function to the reduction of S, is
199 given by
200 \f[
201 \Delta S_L = a^2_L\left(\textbf{w}_L^T\bullet\textbf{w}_L\right)
202 \f]
203 Two test are now applied to decide whether this \f$L^{\mbox{th}}\f$
204 function is to be included in the final expression, or not.
205
206 ## Test 1
207 Denoting by \f$ H_{L-1}\f$ the subspace spanned by \f$\textbf{w}_1,\ldots,\textbf{w}_{L-1}\f$
208 the function \f$\textbf{w}_L\f$ is by construction (see 4) the projection of the function
209 \f$ F_L\f$ onto the direction perpendicular to \f$ H_{L-1}\f$. Now, if the
210 length of \f$\textbf{w}_L\f$ (given by \f$\textbf{w}_L\bullet\textbf{w}_L\f$)
211 is very small compared to the length of \f$\textbf{f}_L\f$ this new
212 function can not contribute much to the reduction of the sum of
213 squares of residuals. The test consists then in calculating the angle
214 \f$ \theta \f$ between the two vectors \f$\textbf{w}_L\f$ \f$ \textbf {f}_L\f$
215 (see also figure 1) and requiring that it's
216 *greater* then a threshold value which the user must set (TMultiDimFit::SetMinAngle).
217
218 \image html multidimfit_img86.gif "Figure 1: (a) angle \\f$\\theta\\f$ between \\f$\\textbf{w}_l\\f$ and \\f$\\textbf{f}_L\\f$, (b) angle \\f$ \\phi \\f$ between \\f$\\textbf{w}_L\\f$ and \\f$\\textbf{D}\\f$"
219
220 ## Test 2
221 Let \f$\textbf{D}\f$ be the data vector to be fitted. As illustrated in
222 figure 1, the \f$L^{\mbox{th}}\f$ function \f$\textbf{w}_L\f$
223 will contribute significantly to the reduction of \f$ S\f$, if the angle
224 \f$\phi^\prime\f$ between \f$\textbf{w}_L\f$ and \f$\textbf{D}\f$ is smaller than
225 an upper limit \f$ \phi \f$, defined by the user (MultiDimFit::SetMaxAngle)
226
227 However, the method automatically readjusts the value of this angle
228 while fitting is in progress, in order to make the selection criteria
229 less and less difficult to be fulfilled. The result is that the
230 functions contributing most to the reduction of \f$ S\f$ are chosen first
231 (TMultiDimFit::TestFunction).
232
233 In case \f$ \phi \f$ isn't defined, an alternative method of
234 performing this second test is used: The \f$L^{\mbox{th}}\f$
235 function \f$\textbf{f}_L\f$ is accepted if (refer also to equation (13))
236 \f[
237 \Delta S_L > \frac{S_{L-1}}{L_{max}-L}
238 \f]
239 where \f$ S_{L-1}\f$ is the sum of the \f$ L-1\f$ first residuals from the
240 \f$ L-1\f$ functions previously accepted; and \f$ L_{max}\f$ is the total number
241 of functions allowed in the final expression of the fit (defined by
242 user).
243
244 From this we see, that by restricting \f$ L_{max}\f$ -- the number of
245 terms in the final model -- the fit is more difficult to perform,
246 since the above selection criteria is more limiting.
247
248 The more coefficients we evaluate, the more the sum of squares of
249 residuals \f$ S\f$ will be reduced. We can evaluate \f$ S\f$ before inverting
250 \f$\mathsf{B}\f$ as shown below.
251
252 ## Coefficients and Coefficient Errors
253 Having found a parameterization, that is the \f$ F_l\f$'s and \f$ L\f$, that
254 minimizes \f$ S\f$, we still need to determine the coefficients
255 \f$ c_l\f$. However, it's a feature of how we choose the significant
256 functions, that the evaluation of the \f$ c_l\f$'s becomes trivial [5]. To derive
257 \f$\mathbf{c}\f$, we first note that
258 equation (4) can be written as
259 \f[
260 \mathsf{F} = \mathsf{W}\mathsf{B}
261 \f]
262 where
263 \f{eqnarray*}{
264 b_{ij} = \frac{\mathbf{f}_j \bullet \mathbf{w}_i}{\mathbf{w}_i^2}
265 & \mbox{if} & i < j\\
266 1 & \mbox{if} & i = j\\
267 0 & \mbox{if} & i > j
268 \f}
269 Consequently, \f$\mathsf{B}\f$ is an upper triangle matrix, which can be
270 readily inverted. So we now evaluate
271 \f[
272 \mathsf{F}\mathsf{B}^{-1} = \mathsf{W}
273 \f]
274 The model \f$\mathsf{W}\mathbf{a}\f$ can therefore be written as
275 \f$(\mathsf{F}\mathsf{B}^{-1})\mathbf{a} = \mathsf{F}(\mathsf{B}^{-1}\mathbf{a})\,.\f$
276
277 The original model \f$\mathsf{F}\mathbf{c}\f$ is therefore identical with
278 this if
279 \f[
280 \mathbf{c} = \left(\mathsf{B}^{-1}\mathbf{a}\right) =
281 \left[\mathbf{a}^T\left(\mathsf{B}^{-1}\right)^T\right]^T\,.
282 \f]
283 The reason we use \f$\left(\mathsf{B}^{-1}\right)^T\f$ rather then
284 \f$\mathsf{B}^{-1}\f$ is to save storage, since \f$\left(\mathsf{B}^{-1}\right)^T\f$
285 can be stored in the same matrix as \f$\mathsf{B}\f$ (TMultiDimFit::MakeCoefficients).
286 The errors in the coefficients is calculated by inverting the curvature matrix
287 of the non-orthogonal functions \f$ f_{lj}\f$ [1] (TMultiDimFit::MakeCoefficientErrors).
288
289 ## Considerations
290 It's important to realize that the training sample should be
291 representative of the problem at hand, in particular along the borders
292 of the region of interest. This is because the algorithm presented
293 here, is a *interpolation*, rather then a *extrapolation* [5].
294
295 Also, the independent variables \f$ x_{i}\f$ need to be linear
296 independent, since the procedure will perform poorly if they are not
297 [5]. One can find an linear transformation from ones
298 original variables \f$ \xi_{i}\f$ to a set of linear independent variables
299 \f$ x_{i}\f$, using a *Principal Components Analysis* (see TPrincipal), and
300 then use the transformed variable as input to this class [5] [6].
301
302 H. Wind also outlines a method for parameterising a multidimensional
303 dependence over a multidimensional set of variables. An example
304 of the method from [5], is a follows (please refer to
305 [5] for a full discussion):
306
307 1. Define \f$\mathbf{P} = (P_1, \ldots, P_5)\f$ are the 5 dependent
308 quantities that define a track.
309 2. Compute, for \f$ M\f$ different values of \f$\mathbf{P}\f$, the tracks
310 through the magnetic field, and determine the corresponding
311 \f$\mathbf{x} = (x_1, \ldots, x_N)\f$.
312 3. Use the simulated observations to determine, with a simple
313 approximation, the values of \f$\mathbf{P}_j\f$. We call these values
314 \f$\mathbf{P}^\prime_j, j = 1, \ldots, M\f$.
315 4. Determine from \f$\mathbf{x}\f$ a set of at least five relevant
316 coordinates \f$\mathbf{x}^\prime\f$, using contrains, *or
317 alternative:*
318 5. Perform a Principal Component Analysis (using TPrincipal), and use
319 to get a linear transformation \f$\mathbf{x} \rightarrow \mathbf{x}^\prime\f$, so that
320 \f$\mathbf{x}^\prime\f$ are constrained and linear independent.
321 6. Perform a Principal Component Analysis on
322 \f$Q_i = P_i / P^\prime_i\, i = 1, \ldots, 5\f$, to get linear
323 indenpendent (among themselves, but not independent of \f$\mathbf{x}\f$) quantities
324 \f$\mathbf{Q}^\prime\f$
325 7. For each component \f$Q^\prime_i\f$ make a multidimensional fit,
326 using \f$\mathbf{x}^\prime\f$ as the variables, thus determining a set of
327 coefficients \f$\mathbf{c}_i\f$.
328
329 To process data, using this parameterisation, do
330 1. Test wether the observation \f$\mathbf{x}\f$ within the domain of
331 the parameterization, using the result from the Principal Component
332 Analysis.
333 2. Determine \f$\mathbf{P}^\prime\f$ as before.
334 3. Determine \f$\mathbf{x}^\prime\f$ as before.
335 4. Use the result of the fit to determine \f$\mathbf{Q}^\prime\f$.
336 5. Transform back to \f$\mathbf{P}\f$ from \f$\mathbf{Q}^\prime\f$, using
337 the result from the Principal Component Analysis.
338
339 ## Testing the parameterization
340 The class also provides functionality for testing the, over the
341 training sample, found parameterization (TMultiDimFit::Fit). This is done by passing
342 the class a test sample of \f$ M_t\f$ tuples of the form
343 \f$(\mathbf{x}_{t,j},D_{t,j}, E_{t,j})\f$, where \f$\mathbf{x}_{t,j}\f$ are the independent
344 variables, \f$ D_{t,j}\f$ the known, dependent quantity, and \f$ E_{t,j}\f$ is
345 the square error in \f$ D_{t,j}\f$ (TMultiDimFit::AddTestRow).
346
347 The parameterization is then evaluated at every \f$\mathbf{x}_t\f$ in the
348 test sample, and
349 \f[
350 S_t \equiv \sum_{j=1}^{M_t} \left(D_{t,j} -
351 D_p\left(\mathbf{x}_{t,j}\right)\right)^2
352 \f]
353 is evaluated. The relative error over the test sample
354 \f[
355 R_t = \frac{S_t}{\sum_{j=1}^{M_t} D_{t,j}^2}
356 \f]
357 should not be to low or high compared to \f$ R\f$ from the training
358 sample. Also, multiple correlation coefficient from both samples should
359 be fairly close, otherwise one of the samples is not representative of
360 the problem. A large difference in the reduced \f$ \chi^2\f$ over the two
361 samples indicate an over fit, and the maximum number of terms in the
362 parameterisation should be reduced.
363
364 It's possible to use [4] to further improve the fit, using the test sample.
365
366 Christian Holm
367
368 ## Bibliography
369 - <a name="bevington"></a> Philip R. Bevington and D. Keith Robinson. *Data Reduction and Error Analysis for
370 the Physical Sciences*. McGraw-Hill, 2 edition, 1992.
371 - <a name="mudifi"></a> R. Brun et al. *Long writeup DD/75-23*, CERN, 1980.
372 - Gene H. Golub and Charles F. van Loan. *Matrix Computations*.
373 John Hopkins University Press, Baltimore, 3 edition, 1996.
374 - <a name="minuit"></a>F. James. *Minuit*. Long writeup D506, CERN, 1998.
375 - <a name="wind72"></a>H. Wind. *Function parameterization*. Proceedings of the 1972 CERN Computing and Data Processing
376 School, volume 72-21 of Yellow report. CERN, 1972.
377 - <a name="wind81"></a>H. Wind. 1. principal component analysis, 2. pattern recognition for track
378 finding, 3. interpolation and functional representation. Yellow report EP/81-12, CERN, 1981.
379
380[1]: classTMultiDimFit.html#bevington
381[2]: classTMultiDimFit.html#mudifi
382[4]: classTMultiDimFit.html#minuit
383[5]: classTMultiDimFit.html#wind72
384[6]: classTMultiDimFit.html#wind81
385[9]: classTMultiDimFit.html#eq:dS2
386*/
387
388
389#include "Riostream.h"
390#include "TMultiDimFit.h"
391#include "TMath.h"
392#include "TH1.h"
393#include "TH2.h"
394#include "TROOT.h"
395#include "TList.h"
396#include "TBrowser.h"
397#include "TDecompChol.h"
398#include "TDatime.h"
399
400
401#define RADDEG (180. / TMath::Pi())
402#define DEGRAD (TMath::Pi() / 180.)
403#define HIST_XORIG 0
404#define HIST_DORIG 1
405#define HIST_XNORM 2
406#define HIST_DSHIF 3
407#define HIST_RX 4
408#define HIST_RD 5
409#define HIST_RTRAI 6
410#define HIST_RTEST 7
411#define PARAM_MAXSTUDY 1
412#define PARAM_SEVERAL 2
413#define PARAM_RELERR 3
414#define PARAM_MAXTERMS 4
415
416
417////////////////////////////////////////////////////////////////////////////////
418
419static void mdfHelper(int&, double*, double&, double*, int);
420
421////////////////////////////////////////////////////////////////////////////////
422
424
425//____________________________________________________________________
426// Static instance. Used with mdfHelper and TMinuit
428
429
430////////////////////////////////////////////////////////////////////////////////
431/// Empty CTOR. Do not use
432
434{
435 fMeanQuantity = 0;
436 fMaxQuantity = 0;
437 fMinQuantity = 0;
438 fSumSqQuantity = 0;
440
441 fNVariables = 0;
442 fSampleSize = 0;
443 fTestSampleSize = 0;
444
445 fMinAngle = 1;
446 fMaxAngle = 0;
447 fMaxTerms = 0;
449 fMaxPowers = 0;
450 fPowerLimit = 0;
451
452 fMaxFunctions = 0;
453 fFunctionCodes = 0;
454 fMaxStudy = 0;
455 fMaxFuncNV = 0;
456
457 fMaxPowersFinal = 0;
458 fPowers = 0;
459 fPowerIndex = 0;
460
461 fMaxResidual = 0;
462 fMinResidual = 0;
463 fMaxResidualRow = 0;
464 fMinResidualRow = 0;
465 fSumSqResidual = 0;
466
467 fNCoefficients = 0;
468 fRMS = 0;
469 fChi2 = 0;
471
472 fError = 0;
473 fTestError = 0;
474 fPrecision = 0;
475 fTestPrecision = 0;
478
479 fHistograms = 0;
480 fHistogramMask = 0;
481 fBinVarX = 100;
482 fBinVarY = 100;
483
484 fFitter = 0;
489
490}
491
492
493////////////////////////////////////////////////////////////////////////////////
494/// Constructor
495/// Second argument is the type of polynomials to use in
496/// parameterisation, one of:
497/// TMultiDimFit::kMonomials
498/// TMultiDimFit::kChebyshev
499/// TMultiDimFit::kLegendre
500///
501/// Options:
502/// K Compute (k)correlation matrix
503/// V Be verbose
504///
505/// Default is no options.
506///
507
511: TNamed("multidimfit","Multi-dimensional fit object"),
512fQuantity(dimension),
513fSqError(dimension),
514fVariables(dimension*100),
515fMeanVariables(dimension),
516fMaxVariables(dimension),
517fMinVariables(dimension)
518{
519 fgInstance = this;
520
521 fMeanQuantity = 0;
522 fMaxQuantity = 0;
523 fMinQuantity = 0;
524 fSumSqQuantity = 0;
526
527 fNVariables = dimension;
528 fSampleSize = 0;
529 fTestSampleSize = 0;
530
531 fMinAngle = 1;
532 fMaxAngle = 0;
533 fMaxTerms = 0;
534 fMinRelativeError = 0.01;
535 fMaxPowers = new Int_t[dimension];
536 fPowerLimit = 1;
537
538 fMaxFunctions = 0;
539 fFunctionCodes = 0;
540 fMaxStudy = 0;
541 fMaxFuncNV = 0;
542
543 fMaxPowersFinal = new Int_t[dimension];
544 fPowers = 0;
545 fPowerIndex = 0;
546
547 fMaxResidual = 0;
548 fMinResidual = 0;
549 fMaxResidualRow = 0;
550 fMinResidualRow = 0;
551 fSumSqResidual = 0;
552
553 fNCoefficients = 0;
554 fRMS = 0;
555 fChi2 = 0;
557
558 fError = 0;
559 fTestError = 0;
560 fPrecision = 0;
561 fTestPrecision = 0;
564
565 fHistograms = 0;
566 fHistogramMask = 0;
567 fBinVarX = 100;
568 fBinVarY = 100;
569
570 fFitter = 0;
571 fPolyType = type;
575 TString opt = option;
576 opt.ToLower();
577
578 if (opt.Contains("k")) fShowCorrelation = kTRUE;
579 if (opt.Contains("v")) fIsVerbose = kTRUE;
580}
581
582
583////////////////////////////////////////////////////////////////////////////////
584/// Destructor
585
587{
588 delete [] fPowers;
589 delete [] fMaxPowers;
590 delete [] fMaxPowersFinal;
591 delete [] fPowerIndex;
592 delete [] fFunctionCodes;
593 if (fHistograms) fHistograms->Clear("nodelete");
594 delete fHistograms;
595}
596
597
598////////////////////////////////////////////////////////////////////////////////
599/// Add a row consisting of fNVariables independent variables, the
600/// known, dependent quantity, and optionally, the square error in
601/// the dependent quantity, to the training sample to be used for the
602/// parameterization.
603/// The mean of the variables and quantity is calculated on the fly,
604/// as outlined in TPrincipal::AddRow.
605/// This sample should be representative of the problem at hand.
606/// Please note, that if no error is given Poisson statistics is
607/// assumed and the square error is set to the value of dependent
608/// quantity. See also the
609/// <a href="#TMultiDimFit:description">class description</a>
610
612{
613 if (!x)
614 return;
615
616 if (++fSampleSize == 1) {
617 fMeanQuantity = D;
618 fMaxQuantity = D;
619 fMinQuantity = D;
620 fSumSqQuantity = D * D;// G.Q. erratum on August 15th, 2008
621 }
622 else {
625 fSumSqQuantity += D * D;
626
627 if (D >= fMaxQuantity) fMaxQuantity = D;
628 if (D <= fMinQuantity) fMinQuantity = D;
629 }
630
631
632 // If the vector isn't big enough to hold the new data, then
633 // expand the vector by half it's size.
635 if (fSampleSize > size) {
638 }
639
640 // Store the value
641 fQuantity(fSampleSize-1) = D;
642 fSqError(fSampleSize-1) = (E == 0 ? D : E);
643
644 // Store data point in internal vector
645 // If the vector isn't big enough to hold the new data, then
646 // expand the vector by half it's size
650
651
652 // Increment the data point counter
653 Int_t i,j;
654 for (i = 0; i < fNVariables; i++) {
655 if (fSampleSize == 1) {
656 fMeanVariables(i) = x[i];
657 fMaxVariables(i) = x[i];
658 fMinVariables(i) = x[i];
659 }
660 else {
661 fMeanVariables(i) *= 1 - 1./Double_t(fSampleSize);
663
664 // Update the maximum value for this component
665 if (x[i] >= fMaxVariables(i)) fMaxVariables(i) = x[i];
666
667 // Update the minimum value for this component
668 if (x[i] <= fMinVariables(i)) fMinVariables(i) = x[i];
669
670 }
671
672 // Store the data.
673 j = (fSampleSize-1) * fNVariables + i;
674 fVariables(j) = x[i];
675 }
676}
677
678
679////////////////////////////////////////////////////////////////////////////////
680/// Add a row consisting of fNVariables independent variables, the
681/// known, dependent quantity, and optionally, the square error in
682/// the dependent quantity, to the test sample to be used for the
683/// test of the parameterization.
684/// This sample needn't be representative of the problem at hand.
685/// Please note, that if no error is given Poisson statistics is
686/// assumed and the square error is set to the value of dependent
687/// quantity. See also the
688/// <a href="#TMultiDimFit:description">class description</a>
689
691{
692 if (fTestSampleSize++ == 0) {
696 }
697
698 // If the vector isn't big enough to hold the new data, then
699 // expand the vector by half it's size.
701 if (fTestSampleSize > size) {
704 }
705
706 // Store the value
708 fTestSqError(fTestSampleSize-1) = (E == 0 ? D : E);
709
710 // Store data point in internal vector
711 // If the vector isn't big enough to hold the new data, then
712 // expand the vector by half it's size
716
717
718 // Increment the data point counter
719 Int_t i,j;
720 for (i = 0; i < fNVariables; i++) {
721 j = fNVariables * (fTestSampleSize - 1) + i;
722 fTestVariables(j) = x[i];
723
724 if (x[i] > fMaxVariables(i))
725 Warning("AddTestRow", "variable %d (row: %d) too large: %f > %f",
726 i, fTestSampleSize, x[i], fMaxVariables(i));
727 if (x[i] < fMinVariables(i))
728 Warning("AddTestRow", "variable %d (row: %d) too small: %f < %f",
729 i, fTestSampleSize, x[i], fMinVariables(i));
730 }
731}
732
733
734////////////////////////////////////////////////////////////////////////////////
735/// Browse the TMultiDimFit object in the TBrowser.
736
738{
739 if (fHistograms) {
740 TIter next(fHistograms);
741 TH1* h = 0;
742 while ((h = (TH1*)next()))
743 b->Add(h,h->GetName());
744 }
745 if (fVariables.IsValid())
746 b->Add(&fVariables, "Variables (Training)");
747 if (fQuantity.IsValid())
748 b->Add(&fQuantity, "Quantity (Training)");
749 if (fSqError.IsValid())
750 b->Add(&fSqError, "Error (Training)");
752 b->Add(&fMeanVariables, "Mean of Variables (Training)");
754 b->Add(&fMaxVariables, "Mean of Variables (Training)");
756 b->Add(&fMinVariables, "Min of Variables (Training)");
758 b->Add(&fTestVariables, "Variables (Test)");
760 b->Add(&fTestQuantity, "Quantity (Test)");
761 if (fTestSqError.IsValid())
762 b->Add(&fTestSqError, "Error (Test)");
763 if (fFunctions.IsValid())
764 b->Add(&fFunctions, "Functions");
766 b->Add(&fCoefficients,"Coefficients");
768 b->Add(&fCoefficientsRMS,"Coefficients Errors");
770 b->Add(&fOrthFunctions, "Orthogonal Functions");
772 b->Add(&fOrthFunctionNorms, "Orthogonal Functions Norms");
773 if (fResiduals.IsValid())
774 b->Add(&fResiduals, "Residuals");
776 b->Add(&fOrthCoefficients,"Orthogonal Coefficients");
778 b->Add(&fOrthCurvatureMatrix,"Orthogonal curvature matrix");
780 b->Add(&fCorrelationMatrix,"Correlation Matrix");
781 if (fFitter)
782 b->Add(fFitter, fFitter->GetName());
783}
784
785
786////////////////////////////////////////////////////////////////////////////////
787/// Clear internal structures and variables
788
790{
791 Int_t i, j, n = fNVariables, m = fMaxFunctions;
792
793 // Training sample, dependent quantity
794 fQuantity.Zero();
795 fSqError.Zero();
796 fMeanQuantity = 0;
797 fMaxQuantity = 0;
798 fMinQuantity = 0;
799 fSumSqQuantity = 0;
801
802 // Training sample, independent variables
804 fNVariables = 0;
805 fSampleSize = 0;
809
810 // Test sample
814 fTestSampleSize = 0;
815
816 // Functions
818 //for (i = 0; i < fMaxTerms; i++) fPowerIndex[i] = 0;
819 //for (i = 0; i < fMaxTerms; i++) fFunctionCodes[i] = 0;
820 fMaxFunctions = 0;
821 fMaxStudy = 0;
824
825 // Control parameters
827 fMinAngle = 0;
828 fMaxAngle = 0;
829 fMaxTerms = 0;
830
831 // Powers
832 for (i = 0; i < n; i++) {
833 fMaxPowers[i] = 0;
834 fMaxPowersFinal[i] = 0;
835 for (j = 0; j < m; j++)
836 fPowers[i * n + j] = 0;
837 }
838 fPowerLimit = 0;
839
840 // Residuals
841 fMaxResidual = 0;
842 fMinResidual = 0;
843 fMaxResidualRow = 0;
844 fMinResidualRow = 0;
845 fSumSqResidual = 0;
846
847 // Fit
848 fNCoefficients = 0;
851 fRMS = 0;
853 fError = 0;
854 fTestError = 0;
855 fPrecision = 0;
856 fTestPrecision = 0;
857
858 // Coefficients
863
864 // Options
868}
869
870
871////////////////////////////////////////////////////////////////////////////////
872/// Evaluate parameterization at point x. Optional argument coeff is
873/// a vector of coefficients for the parameterisation, fNCoefficients
874/// elements long.
875
876Double_t TMultiDimFit::Eval(const Double_t *x, const Double_t* coeff) const
877{
878 Double_t returnValue = fMeanQuantity;
879 Double_t term = 0;
880 Int_t i, j;
881
882 for (i = 0; i < fNCoefficients; i++) {
883 // Evaluate the ith term in the expansion
884 term = (coeff ? coeff[i] : fCoefficients(i));
885 for (j = 0; j < fNVariables; j++) {
886 // Evaluate the factor (polynomial) in the j-th variable.
888 Double_t y = 1 + 2. / (fMaxVariables(j) - fMinVariables(j))
889 * (x[j] - fMaxVariables(j));
890 term *= EvalFactor(p,y);
891 }
892 // Add this term to the final result
893 returnValue += term;
894 }
895 return returnValue;
896}
897
898
899////////////////////////////////////////////////////////////////////////////////
900/// Evaluate parameterization error at point x. Optional argument coeff is
901/// a vector of coefficients for the parameterisation, fNCoefficients
902/// elements long.
903
905{
906 Double_t returnValue = 0;
907 Double_t term = 0;
908 Int_t i, j;
909
910 for (i = 0; i < fNCoefficients; i++) {
911 // std::cout << "Error coef " << i << " -> " << fCoefficientsRMS(i) << std::endl;
912 }
913 for (i = 0; i < fNCoefficients; i++) {
914 // Evaluate the ith term in the expansion
915 term = (coeff ? coeff[i] : fCoefficientsRMS(i));
916 for (j = 0; j < fNVariables; j++) {
917 // Evaluate the factor (polynomial) in the j-th variable.
919 Double_t y = 1 + 2. / (fMaxVariables(j) - fMinVariables(j))
920 * (x[j] - fMaxVariables(j));
921 term *= EvalFactor(p,y);
922 // std::cout << "i,j " << i << ", " << j << " " << p << " " << y << " " << EvalFactor(p,y) << " " << term << std::endl;
923 }
924 // Add this term to the final result
925 returnValue += term*term;
926 // std::cout << " i = " << i << " value = " << returnValue << std::endl;
927 }
928 returnValue = sqrt(returnValue);
929 return returnValue;
930}
931
932
933////////////////////////////////////////////////////////////////////////////////
934/// PRIVATE METHOD:
935/// Calculate the control parameter from the passed powers
936
938{
939 Double_t s = 0;
940 Double_t epsilon = 1e-6; // a small number
941 for (Int_t i = 0; i < fNVariables; i++) {
942 if (fMaxPowers[i] != 1)
943 s += (epsilon + iv[i] - 1) / (epsilon + fMaxPowers[i] - 1);
944 }
945 return s;
946}
947
948////////////////////////////////////////////////////////////////////////////////
949/// PRIVATE METHOD:
950/// Evaluate function with power p at variable value x
951
953{
954 Int_t i = 0;
955 Double_t p1 = 1;
956 Double_t p2 = 0;
957 Double_t p3 = 0;
958 Double_t r = 0;
959
960 switch(p) {
961 case 1:
962 r = 1;
963 break;
964 case 2:
965 r = x;
966 break;
967 default:
968 p2 = x;
969 for (i = 3; i <= p; i++) {
970 p3 = p2 * x;
971 if (fPolyType == kLegendre)
972 p3 = ((2 * i - 3) * p2 * x - (i - 2) * p1) / (i - 1);
973 else if (fPolyType == kChebyshev)
974 p3 = 2 * x * p2 - p1;
975 p1 = p2;
976 p2 = p3;
977 }
978 r = p3;
979 }
980
981 return r;
982}
983
984
985////////////////////////////////////////////////////////////////////////////////
986/// Find the parameterization
987///
988/// Options:
989/// None so far
990///
991/// For detailed description of what this entails, please refer to the
992/// <a href="#TMultiDimFit:description">class description</a>
993
995{
1002}
1003
1004////////////////////////////////////////////////////////////////////////////////
1005/// Try to fit the found parameterisation to the test sample.
1006///
1007/// Options
1008/// M use Minuit to improve coefficients
1009///
1010/// Also, refer to
1011/// <a href="#TMultiDimFit:description">class description</a>
1012
1014{
1015 Int_t i, j;
1017 Double_t sumSqD = 0;
1018 Double_t sumD = 0;
1019 Double_t sumSqR = 0;
1020 Double_t sumR = 0;
1021
1022 // Calculate the residuals over the test sample
1023 for (i = 0; i < fTestSampleSize; i++) {
1024 for (j = 0; j < fNVariables; j++)
1025 x[j] = fTestVariables(i * fNVariables + j);
1026 Double_t res = fTestQuantity(i) - Eval(x);
1027 sumD += fTestQuantity(i);
1028 sumSqD += fTestQuantity(i) * fTestQuantity(i);
1029 sumR += res;
1030 sumSqR += res * res;
1032 ((TH1D*)fHistograms->FindObject("res_test"))->Fill(res);
1033 }
1034 Double_t dAvg = sumSqD - (sumD * sumD) / fTestSampleSize;
1035 Double_t rAvg = sumSqR - (sumR * sumR) / fTestSampleSize;
1036 fTestCorrelationCoeff = (dAvg - rAvg) / dAvg;
1037 fTestError = sumSqR;
1038 fTestPrecision = sumSqR / sumSqD;
1039
1040 TString opt(option);
1041 opt.ToLower();
1042
1043 if (!opt.Contains("m"))
1044 MakeChi2();
1045
1047 Warning("Fit", "test sample is very small");
1048
1049 if (!opt.Contains("m")) {
1050 Error("Fit", "invalid option");
1051 delete [] x;
1052 return;
1053 }
1054
1056 if (!fFitter) {
1057 Error("Fit", "Cannot create Fitter");
1058 delete [] x;
1059 return;
1060 }
1062
1063 const Int_t maxArgs = 16;
1064 Int_t args = 1;
1065 Double_t* arglist = new Double_t[maxArgs];
1066 arglist[0] = -1;
1067 fFitter->ExecuteCommand("SET PRINT",arglist,args);
1068
1069 for (i = 0; i < fNCoefficients; i++) {
1070 Double_t startVal = fCoefficients(i);
1071 Double_t startErr = fCoefficientsRMS(i);
1072 fFitter->SetParameter(i, Form("coeff%02d",i),
1073 startVal, startErr, 0, 0);
1074 }
1075
1076 // arglist[0] = 0;
1077 args = 1;
1078 // fFitter->ExecuteCommand("SET PRINT",arglist,args);
1079 fFitter->ExecuteCommand("MIGRAD",arglist,args);
1080
1081 for (i = 0; i < fNCoefficients; i++) {
1082 Double_t val = 0, err = 0, low = 0, high = 0;
1083 fFitter->GetParameter(i, Form("coeff%02d",i),
1084 val, err, low, high);
1085 fCoefficients(i) = val;
1086 fCoefficientsRMS(i) = err;
1087 }
1088 delete [] x;
1089}
1090
1091////////////////////////////////////////////////////////////////////////////////
1092/// Return the static instance.
1093
1095{
1096 return fgInstance;
1097}
1098
1099////////////////////////////////////////////////////////////////////////////////
1100/// PRIVATE METHOD:
1101/// Create list of candidate functions for the parameterisation. See
1102/// also
1103/// <a href="#TMultiDimFit:description">class description</a>
1104
1106{
1107 Int_t i = 0;
1108 Int_t j = 0;
1109 Int_t k = 0;
1110
1111 // The temporary array to store the powers in. We don't need to
1112 // initialize this array however.
1114 Int_t *powers = new Int_t[fMaxFuncNV];
1115
1116 // store of `control variables'
1117 Double_t* control = new Double_t[fMaxFunctions];
1118
1119 // We've better initialize the variables
1120 Int_t *iv = new Int_t[fNVariables];
1121 for (i = 0; i < fNVariables; i++)
1122 iv[i] = 1;
1123
1124 if (!fIsUserFunction) {
1125
1126 // Number of funcs selected
1127 Int_t numberFunctions = 0;
1128
1129 while (kTRUE) {
1130 // Get the control value for this function
1131 Double_t s = EvalControl(iv);
1132
1133 if (s <= fPowerLimit) {
1134
1135 // Call over-loadable method Select, as to allow the user to
1136 // interfere with the selection of functions.
1137 if (Select(iv)) {
1138 numberFunctions++;
1139
1140 // If we've reached the user defined limit of how many
1141 // functions we can consider, break out of the loop
1142 if (numberFunctions > fMaxFunctions)
1143 break;
1144
1145 // Store the control value, so we can sort array of powers
1146 // later on
1147 control[numberFunctions-1] = Int_t(1.0e+6*s);
1148
1149 // Store the powers in powers array.
1150 for (i = 0; i < fNVariables; i++) {
1151 j = (numberFunctions - 1) * fNVariables + i;
1152 powers[j] = iv[i];
1153 }
1154 } // if (Select())
1155 } // if (s <= fPowerLimit)
1156
1157 for (i = 0; i < fNVariables; i++)
1158 if (iv[i] < fMaxPowers[i])
1159 break;
1160
1161 // If all variables have reached their maximum power, then we
1162 // break out of the loop
1163 if (i == fNVariables) {
1164 fMaxFunctions = numberFunctions;
1165 break;
1166 }
1167
1168 // Next power in variable i
1169 if (i < fNVariables) iv[i]++;
1170
1171 for (j = 0; j < i; j++)
1172 iv[j] = 1;
1173 } // while (kTRUE)
1174 }
1175 else {
1176 // In case the user gave an explicit function
1177 for (i = 0; i < fMaxFunctions; i++) {
1178 // Copy the powers to working arrays
1179 for (j = 0; j < fNVariables; j++) {
1180 powers[i * fNVariables + j] = fPowers[i * fNVariables + j];
1181 iv[j] = fPowers[i * fNVariables + j];
1182 }
1183
1184 control[i] = Int_t(1.0e+6*EvalControl(iv));
1185 }
1186 }
1187
1188 // Now we need to sort the powers according to least `control
1189 // variable'
1190 Int_t *order = new Int_t[fMaxFunctions];
1191 for (i = 0; i < fMaxFunctions; i++)
1192 order[i] = i;
1194 fPowers = new Int_t[fMaxFuncNV];
1195
1196 for (i = 0; i < fMaxFunctions; i++) {
1197 Double_t x = control[i];
1198 Int_t l = order[i];
1199 k = i;
1200
1201 for (j = i; j < fMaxFunctions; j++) {
1202 if (control[j] <= x) {
1203 x = control[j];
1204 l = order[j];
1205 k = j;
1206 }
1207 }
1208
1209 if (k != i) {
1210 control[k] = control[i];
1211 control[i] = x;
1212 order[k] = order[i];
1213 order[i] = l;
1214 }
1215 }
1216
1217 for (i = 0; i < fMaxFunctions; i++)
1218 for (j = 0; j < fNVariables; j++)
1219 fPowers[i * fNVariables + j] = powers[order[i] * fNVariables + j];
1220
1221 delete [] control;
1222 delete [] powers;
1223 delete [] order;
1224 delete [] iv;
1225}
1226
1227
1228////////////////////////////////////////////////////////////////////////////////
1229/// Calculate Chi square over either the test sample. The optional
1230/// argument coeff is a vector of coefficients to use in the
1231/// evaluation of the parameterisation. If coeff == 0, then the found
1232/// coefficients is used.
1233/// Used my MINUIT for fit (see TMultDimFit::Fit)
1234
1236{
1237 fChi2 = 0;
1238 Int_t i, j;
1240 for (i = 0; i < fTestSampleSize; i++) {
1241 // Get the stored point
1242 for (j = 0; j < fNVariables; j++)
1243 x[j] = fTestVariables(i * fNVariables + j);
1244
1245 // Evaluate function. Scale to shifted values
1246 Double_t f = Eval(x,coeff);
1247
1248 // Calculate contribution to Chic square
1249 fChi2 += 1. / TMath::Max(fTestSqError(i),1e-20)
1250 * (fTestQuantity(i) - f) * (fTestQuantity(i) - f);
1251 }
1252
1253 // Clean up
1254 delete [] x;
1255
1256 return fChi2;
1257}
1258
1259
1260////////////////////////////////////////////////////////////////////////////////
1261/// Generate the file `<filename>` with .C appended if argument doesn't
1262/// end in .cxx or .C. The contains the implementation of the
1263/// function:
1264///
1265/// `Double_t <funcname>(Double_t *x)`
1266///
1267/// which does the same as TMultiDimFit::Eval. Please refer to this
1268/// method.
1269///
1270/// Further, the static variables:
1271///
1272/// Int_t gNVariables
1273/// Int_t gNCoefficients
1274/// Double_t gDMean
1275/// Double_t gXMean[]
1276/// Double_t gXMin[]
1277/// Double_t gXMax[]
1278/// Double_t gCoefficient[]
1279/// Int_t gPower[]
1280///
1281/// are initialized. The only ROOT header file needed is Rtypes.h
1282///
1283/// See TMultiDimFit::MakeRealCode for a list of options
1284
1286{
1287
1288 TString outName(filename);
1289 if (!outName.EndsWith(".C") && !outName.EndsWith(".cxx"))
1290 outName += ".C";
1291
1292 MakeRealCode(outName.Data(),"",option);
1293}
1294
1295
1296
1297////////////////////////////////////////////////////////////////////////////////
1298/// PRIVATE METHOD:
1299/// Compute the errors on the coefficients. For this to be done, the
1300/// curvature matrix of the non-orthogonal functions, is computed.
1301
1303{
1304 Int_t i = 0;
1305 Int_t j = 0;
1306 Int_t k = 0;
1310
1311 TMatrixDSym curvatureMatrix(fNCoefficients);
1312
1313 // Build the curvature matrix
1314 for (i = 0; i < fNCoefficients; i++) {
1315 iF = TMatrixDRow(fFunctions,i);
1316 for (j = 0; j <= i; j++) {
1317 jF = TMatrixDRow(fFunctions,j);
1318 for (k = 0; k < fSampleSize; k++)
1319 curvatureMatrix(i,j) +=
1320 1 / TMath::Max(fSqError(k), 1e-20) * iF(k) * jF(k);
1321 curvatureMatrix(j,i) = curvatureMatrix(i,j);
1322 }
1323 }
1324
1325 // Calculate Chi Square
1326 fChi2 = 0;
1327 for (i = 0; i < fSampleSize; i++) {
1328 Double_t f = 0;
1329 for (j = 0; j < fNCoefficients; j++)
1330 f += fCoefficients(j) * fFunctions(j,i);
1331 fChi2 += 1. / TMath::Max(fSqError(i),1e-20) * (fQuantity(i) - f)
1332 * (fQuantity(i) - f);
1333 }
1334
1335 // Invert the curvature matrix
1336 const TVectorD diag = TMatrixDDiag_const(curvatureMatrix);
1337 curvatureMatrix.NormByDiag(diag);
1338
1339 TDecompChol chol(curvatureMatrix);
1340 if (!chol.Decompose())
1341 Error("MakeCoefficientErrors", "curvature matrix is singular");
1342 chol.Invert(curvatureMatrix);
1343
1344 curvatureMatrix.NormByDiag(diag);
1345
1346 for (i = 0; i < fNCoefficients; i++)
1347 fCoefficientsRMS(i) = TMath::Sqrt(curvatureMatrix(i,i));
1348}
1349
1350
1351////////////////////////////////////////////////////////////////////////////////
1352/// PRIVATE METHOD:
1353/// Invert the model matrix B, and compute final coefficients. For a
1354/// more thorough discussion of what this means, please refer to the
1355/// <a href="#TMultiDimFit:description">class description</a>
1356///
1357/// First we invert the lower triangle matrix fOrthCurvatureMatrix
1358/// and store the inverted matrix in the upper triangle.
1359
1361{
1362 Int_t i = 0, j = 0;
1363 Int_t col = 0, row = 0;
1364
1365 // Invert the B matrix
1366 for (col = 1; col < fNCoefficients; col++) {
1367 for (row = col - 1; row > -1; row--) {
1368 fOrthCurvatureMatrix(row,col) = 0;
1369 for (i = row; i <= col ; i++)
1370 fOrthCurvatureMatrix(row,col) -=
1372 * fOrthCurvatureMatrix(i,col);
1373 }
1374 }
1375
1376 // Compute the final coefficients
1378
1379 for (i = 0; i < fNCoefficients; i++) {
1380 Double_t sum = 0;
1381 for (j = i; j < fNCoefficients; j++)
1383 fCoefficients(i) = sum;
1384 }
1385
1386 // Compute the final residuals
1388 for (i = 0; i < fSampleSize; i++)
1389 fResiduals(i) = fQuantity(i);
1390
1391 for (i = 0; i < fNCoefficients; i++)
1392 for (j = 0; j < fSampleSize; j++)
1393 fResiduals(j) -= fCoefficients(i) * fFunctions(i,j);
1394
1395 // Compute the max and minimum, and squared sum of the evaluated
1396 // residuals
1397 fMinResidual = 10e10;
1398 fMaxResidual = -10e10;
1399 Double_t sqRes = 0;
1400 for (i = 0; i < fSampleSize; i++){
1401 sqRes += fResiduals(i) * fResiduals(i);
1402 if (fResiduals(i) <= fMinResidual) {
1404 fMinResidualRow = i;
1405 }
1406 if (fResiduals(i) >= fMaxResidual) {
1408 fMaxResidualRow = i;
1409 }
1410 }
1411
1414
1415 // If we use histograms, fill some more
1419 for (i = 0; i < fSampleSize; i++) {
1421 ((TH2D*)fHistograms->FindObject("res_d"))->Fill(fQuantity(i),
1422 fResiduals(i));
1424 ((TH1D*)fHistograms->FindObject("res_train"))->Fill(fResiduals(i));
1425
1427 for (j = 0; j < fNVariables; j++)
1428 ((TH2D*)fHistograms->FindObject(Form("res_x_%d",j)))
1429 ->Fill(fVariables(i * fNVariables + j),fResiduals(i));
1430 }
1431 } // If histograms
1432
1433}
1434
1435
1436////////////////////////////////////////////////////////////////////////////////
1437/// PRIVATE METHOD:
1438/// Compute the correlation matrix
1439
1441{
1442 if (!fShowCorrelation)
1443 return;
1444
1446
1447 Double_t d2 = 0;
1448 Double_t ddotXi = 0; // G.Q. needs to be reinitialized in the loop over i fNVariables
1449 Double_t xiNorm = 0; // G.Q. needs to be reinitialized in the loop over i fNVariables
1450 Double_t xidotXj = 0; // G.Q. needs to be reinitialized in the loop over j fNVariables
1451 Double_t xjNorm = 0; // G.Q. needs to be reinitialized in the loop over j fNVariables
1452
1453 Int_t i, j, k, l, m; // G.Q. added m variable
1454 for (i = 0; i < fSampleSize; i++)
1455 d2 += fQuantity(i) * fQuantity(i);
1456
1457 for (i = 0; i < fNVariables; i++) {
1458 ddotXi = 0.; // G.Q. reinitialisation
1459 xiNorm = 0.; // G.Q. reinitialisation
1460 for (j = 0; j< fSampleSize; j++) {
1461 // Index of sample j of variable i
1462 k = j * fNVariables + i;
1463 ddotXi += fQuantity(j) * (fVariables(k) - fMeanVariables(i));
1464 xiNorm += (fVariables(k) - fMeanVariables(i))
1465 * (fVariables(k) - fMeanVariables(i));
1466 }
1467 fCorrelationMatrix(i,0) = ddotXi / TMath::Sqrt(d2 * xiNorm);
1468
1469 for (j = 0; j < i; j++) {
1470 xidotXj = 0.; // G.Q. reinitialisation
1471 xjNorm = 0.; // G.Q. reinitialisation
1472 for (k = 0; k < fSampleSize; k++) {
1473 // Index of sample j of variable i
1474 // l = j * fNVariables + k; // G.Q.
1475 l = k * fNVariables + j; // G.Q.
1476 m = k * fNVariables + i; // G.Q.
1477 // G.Q. xidotXj += (fVariables(i) - fMeanVariables(i))
1478 // G.Q. * (fVariables(l) - fMeanVariables(j));
1479 xidotXj += (fVariables(m) - fMeanVariables(i))
1480 * (fVariables(l) - fMeanVariables(j)); // G.Q. modified index for Xi
1481 xjNorm += (fVariables(l) - fMeanVariables(j))
1482 * (fVariables(l) - fMeanVariables(j));
1483 }
1484 //fCorrelationMatrix(i+1,j) = xidotXj / TMath::Sqrt(xiNorm * xjNorm);
1485 fCorrelationMatrix(i,j+1) = xidotXj / TMath::Sqrt(xiNorm * xjNorm);
1486 }
1487 }
1488}
1489
1490
1491
1492////////////////////////////////////////////////////////////////////////////////
1493/// PRIVATE METHOD:
1494/// Make Gram-Schmidt orthogonalisation. The class description gives
1495/// a thorough account of this algorithm, as well as
1496/// references. Please refer to the
1497/// <a href="#TMultiDimFit:description">class description</a>
1498
1500{
1501
1502 // calculate w_i, that is, evaluate the current function at data
1503 // point i
1504 Double_t f2 = 0;
1507 Int_t j = 0;
1508 Int_t k = 0;
1509
1510 for (j = 0; j < fSampleSize; j++) {
1511 fFunctions(fNCoefficients, j) = 1;
1513 // First, however, we need to calculate f_fNCoefficients
1514 for (k = 0; k < fNVariables; k++) {
1515 Int_t p = fPowers[function * fNVariables + k];
1516 Double_t x = fVariables(j * fNVariables + k);
1518 }
1519
1520 // Calculate f dot f in f2
1522 // Assign to w_fNCoefficients f_fNCoefficients
1524 }
1525
1526 // the first column of w is equal to f
1527 for (j = 0; j < fNCoefficients; j++) {
1528 Double_t fdw = 0;
1529 // Calculate (f_fNCoefficients dot w_j) / w_j^2
1530 for (k = 0; k < fSampleSize; k++) {
1531 fdw += fFunctions(fNCoefficients, k) * fOrthFunctions(j,k)
1532 / fOrthFunctionNorms(j);
1533 }
1534
1536 // and subtract it from the current value of w_ij
1537 for (k = 0; k < fSampleSize; k++)
1539 }
1540
1541 for (j = 0; j < fSampleSize; j++) {
1542 // calculate squared length of w_fNCoefficients
1546
1547 // calculate D dot w_fNCoefficients in A
1550 }
1551
1552 // First test, but only if didn't user specify
1553 if (!fIsUserFunction)
1556 return 0;
1557
1558 // The result found by this code for the first residual is always
1559 // much less then the one found be MUDIFI. That's because it's
1560 // supposed to be. The cause is the improved precision of Double_t
1561 // over DOUBLE PRECISION!
1565
1566 // Calculate the residual from including this fNCoefficients.
1568
1569 return dResidur;
1570}
1571
1572
1573////////////////////////////////////////////////////////////////////////////////
1574/// Make histograms of the result of the analysis. This message
1575/// should be sent after having read all data points, but before
1576/// finding the parameterization
1577///
1578/// Options:
1579/// A All the below
1580/// X Original independent variables
1581/// D Original dependent variables
1582/// N Normalised independent variables
1583/// S Shifted dependent variables
1584/// R1 Residuals versus normalised independent variables
1585/// R2 Residuals versus dependent variable
1586/// R3 Residuals computed on training sample
1587/// R4 Residuals computed on test sample
1588///
1589/// For a description of these quantities, refer to
1590/// <a href="#TMultiDimFit:description">class description</a>
1591
1593{
1594 TString opt(option);
1595 opt.ToLower();
1596
1597 if (opt.Length() < 1)
1598 return;
1599
1600 if (!fHistograms)
1601 fHistograms = new TList;
1602
1603 // Counter variable
1604 Int_t i = 0;
1605
1606 // Histogram of original variables
1607 if (opt.Contains("x") || opt.Contains("a")) {
1609 for (i = 0; i < fNVariables; i++)
1610 if (!fHistograms->FindObject(Form("x_%d_orig",i)))
1611 fHistograms->Add(new TH1D(Form("x_%d_orig",i),
1612 Form("Original variable # %d",i),
1614 fMaxVariables(i)));
1615 }
1616
1617 // Histogram of original dependent variable
1618 if (opt.Contains("d") || opt.Contains("a")) {
1620 if (!fHistograms->FindObject("d_orig"))
1621 fHistograms->Add(new TH1D("d_orig", "Original Quantity",
1623 }
1624
1625 // Histograms of normalized variables
1626 if (opt.Contains("n") || opt.Contains("a")) {
1628 for (i = 0; i < fNVariables; i++)
1629 if (!fHistograms->FindObject(Form("x_%d_norm",i)))
1630 fHistograms->Add(new TH1D(Form("x_%d_norm",i),
1631 Form("Normalized variable # %d",i),
1632 fBinVarX, -1,1));
1633 }
1634
1635 // Histogram of shifted dependent variable
1636 if (opt.Contains("s") || opt.Contains("a")) {
1638 if (!fHistograms->FindObject("d_shifted"))
1639 fHistograms->Add(new TH1D("d_shifted", "Shifted Quantity",
1642 }
1643
1644 // Residual from training sample versus independent variables
1645 if (opt.Contains("r1") || opt.Contains("a")) {
1647 for (i = 0; i < fNVariables; i++)
1648 if (!fHistograms->FindObject(Form("res_x_%d",i)))
1649 fHistograms->Add(new TH2D(Form("res_x_%d",i),
1650 Form("Computed residual versus x_%d", i),
1651 fBinVarX, -1, 1,
1652 fBinVarY,
1655 }
1656
1657 // Residual from training sample versus. dependent variable
1658 if (opt.Contains("r2") || opt.Contains("a")) {
1660 if (!fHistograms->FindObject("res_d"))
1661 fHistograms->Add(new TH2D("res_d",
1662 "Computed residuals vs Quantity",
1663 fBinVarX,
1666 fBinVarY,
1669 }
1670
1671 // Residual from training sample
1672 if (opt.Contains("r3") || opt.Contains("a")) {
1674 if (!fHistograms->FindObject("res_train"))
1675 fHistograms->Add(new TH1D("res_train",
1676 "Computed residuals over training sample",
1679
1680 }
1681 if (opt.Contains("r4") || opt.Contains("a")) {
1683 if (!fHistograms->FindObject("res_test"))
1684 fHistograms->Add(new TH1D("res_test",
1685 "Distribution of residuals from test",
1688 }
1689}
1690
1691
1692////////////////////////////////////////////////////////////////////////////////
1693/// Generate the file `<classname>MDF.cxx` which contains the
1694/// implementation of the method:
1695///
1696/// `Double_t <classname>::%MDF(Double_t *x)`
1697///
1698/// which does the same as TMultiDimFit::Eval. Please refer to this
1699/// method.
1700///
1701/// Further, the public static members:
1702/// \code{.cpp}
1703/// Int_t <classname>::fgNVariables
1704/// Int_t <classname>::fgNCoefficients
1705/// Double_t <classname>::fgDMean
1706/// Double_t <classname>::fgXMean[] //[fgNVariables]
1707/// Double_t <classname>::fgXMin[] //[fgNVariables]
1708/// Double_t <classname>::fgXMax[] //[fgNVariables]
1709/// Double_t <classname>::fgCoefficient[] //[fgNCoeffficents]
1710/// Int_t <classname>::fgPower[] //[fgNCoeffficents*fgNVariables]
1711/// \endcode
1712///
1713/// are initialized, and assumed to exist. The class declaration is
1714/// assumed to be in `<classname>.h` and assumed to be provided by the
1715/// user.
1716///
1717/// \see TMultiDimFit::MakeRealCode for a list of options
1718///
1719/// The minimal class definition is:
1720/// \code{.cpp}
1721/// class <classname> {
1722/// public:
1723/// Int_t <classname>::fgNVariables; // Number of variables
1724/// Int_t <classname>::fgNCoefficients; // Number of terms
1725/// Double_t <classname>::fgDMean; // Mean from training sample
1726/// Double_t <classname>::fgXMean[]; // Mean from training sample
1727/// Double_t <classname>::fgXMin[]; // Min from training sample
1728/// Double_t <classname>::fgXMax[]; // Max from training sample
1729/// Double_t <classname>::fgCoefficient[]; // Coefficients
1730/// Int_t <classname>::fgPower[]; // Function powers
1731///
1732/// Double_t Eval(Double_t *x);
1733/// };
1734/// \endcode
1735///
1736/// Whether the method `<classname>::%Eval` should be static or not, is
1737/// up to the user.
1738
1740{
1741 MakeRealCode(Form("%sMDF.cxx", classname), classname, option);
1742}
1743
1744
1745
1746////////////////////////////////////////////////////////////////////////////////
1747/// PRIVATE METHOD:
1748/// Normalize data to the interval [-1;1]. This is needed for the
1749/// classes method to work.
1750
1752{
1753 Int_t i = 0;
1754 Int_t j = 0;
1755 Int_t k = 0;
1756
1757 for (i = 0; i < fSampleSize; i++) {
1759 ((TH1D*)fHistograms->FindObject("d_orig"))->Fill(fQuantity(i));
1760
1763
1765 ((TH1D*)fHistograms->FindObject("d_shifted"))->Fill(fQuantity(i));
1766
1767 for (j = 0; j < fNVariables; j++) {
1768 Double_t range = 1. / (fMaxVariables(j) - fMinVariables(j));
1769 k = i * fNVariables + j;
1770
1771 // Fill histograms of original independent variables
1773 ((TH1D*)fHistograms->FindObject(Form("x_%d_orig",j)))
1774 ->Fill(fVariables(k));
1775
1776 // Normalise independent variables
1777 fVariables(k) = 1 + 2 * range * (fVariables(k) - fMaxVariables(j));
1778
1779 // Fill histograms of normalised independent variables
1781 ((TH1D*)fHistograms->FindObject(Form("x_%d_norm",j)))
1782 ->Fill(fVariables(k));
1783
1784 }
1785 }
1786 // Shift min and max of dependent variable
1789
1790 // Shift mean of independent variables
1791 for (i = 0; i < fNVariables; i++) {
1792 Double_t range = 1. / (fMaxVariables(i) - fMinVariables(i));
1793 fMeanVariables(i) = 1 + 2 * range * (fMeanVariables(i)
1794 - fMaxVariables(i));
1795 }
1796}
1797
1798
1799////////////////////////////////////////////////////////////////////////////////
1800/// PRIVATE METHOD:
1801/// Find the parameterization over the training sample. A full account
1802/// of the algorithm is given in the
1803/// <a href="#TMultiDimFit:description">class description</a>
1804
1806{
1807 Int_t i = -1;
1808 Int_t j = 0;
1809 Int_t k = 0;
1810 Int_t maxPass = 3;
1811 Int_t studied = 0;
1812 Double_t squareResidual = fSumSqAvgQuantity;
1813 fNCoefficients = 0;
1820 fFunctions = 1;
1821
1824 Int_t l;
1825 for (l=0;l<fMaxFunctions;l++) fFunctionCodes[l] = 0;
1826 for (l=0;l<fMaxTerms;l++) fPowerIndex[l] = 0;
1827
1828 if (fMaxAngle != 0) maxPass = 100;
1829 if (fIsUserFunction) maxPass = 1;
1830
1831 // Loop over the number of functions we want to study.
1832 // increment inspection counter
1833 while(kTRUE) {
1834
1835 // Reach user defined limit of studies
1836 if (studied++ >= fMaxStudy) {
1838 break;
1839 }
1840
1841 // Considered all functions several times
1842 if (k >= maxPass) {
1844 break;
1845 }
1846
1847 // increment function counter
1848 i++;
1849
1850 // If we've reached the end of the functions, restart pass
1851 if (i == fMaxFunctions) {
1852 if (fMaxAngle != 0)
1853 fMaxAngle += (90 - fMaxAngle) / 2;
1854 i = 0;
1855 studied--;
1856 k++;
1857 continue;
1858 }
1859 if (studied == 1)
1860 fFunctionCodes[i] = 0;
1861 else if (fFunctionCodes[i] >= 2)
1862 continue;
1863
1864 // Print a happy message
1865 if (fIsVerbose && studied == 1)
1866 std::cout << "Coeff SumSqRes Contrib Angle QM Func"
1867 << " Value W^2 Powers" << std::endl;
1868
1869 // Make the Gram-Schmidt
1870 Double_t dResidur = MakeGramSchmidt(i);
1871
1872 if (dResidur == 0) {
1873 // This function is no good!
1874 // First test is in MakeGramSchmidt
1875 fFunctionCodes[i] = 1;
1876 continue;
1877 }
1878
1879 // If user specified function, assume they know what they are doing
1880 if (!fIsUserFunction) {
1881 // Flag this function as considered
1882 fFunctionCodes[i] = 2;
1883
1884 // Test if this function contributes to the fit
1885 if (!TestFunction(squareResidual, dResidur)) {
1886 fFunctionCodes[i] = 1;
1887 continue;
1888 }
1889 }
1890
1891 // If we get to here, the function currently considered is
1892 // fNCoefficients, so we increment the counter
1893 // Flag this function as OK, and store and the number in the
1894 // index.
1895 fFunctionCodes[i] = 3;
1898
1899 // We add the current contribution to the sum of square of
1900 // residuals;
1901 squareResidual -= dResidur;
1902
1903
1904 // Calculate control parameter from this function
1905 for (j = 0; j < fNVariables; j++) {
1906 if (fNCoefficients == 1
1907 || fMaxPowersFinal[j] <= fPowers[i * fNVariables + j] - 1)
1908 fMaxPowersFinal[j] = fPowers[i * fNVariables + j] - 1;
1909 }
1911
1912 // Print the statistics about this function
1913 if (fIsVerbose) {
1914 std::cout << std::setw(5) << fNCoefficients << " "
1915 << std::setw(10) << std::setprecision(4) << squareResidual << " "
1916 << std::setw(10) << std::setprecision(4) << dResidur << " "
1917 << std::setw(7) << std::setprecision(3) << fMaxAngle << " "
1918 << std::setw(7) << std::setprecision(3) << s << " "
1919 << std::setw(5) << i << " "
1920 << std::setw(10) << std::setprecision(4)
1922 << std::setw(10) << std::setprecision(4)
1924 << std::flush;
1925 for (j = 0; j < fNVariables; j++)
1926 std::cout << " " << fPowers[i * fNVariables + j] - 1 << std::flush;
1927 std::cout << std::endl;
1928 }
1929
1930 if (fNCoefficients >= fMaxTerms /* && fIsVerbose */) {
1932 break;
1933 }
1934
1935 Double_t err = TMath::Sqrt(TMath::Max(1e-20,squareResidual) /
1937 if (err < fMinRelativeError) {
1939 break;
1940 }
1941
1942 }
1943
1944 fError = TMath::Max(1e-20,squareResidual);
1947}
1948
1949
1950////////////////////////////////////////////////////////////////////////////////
1951/// PRIVATE METHOD:
1952/// This is the method that actually generates the code for the
1953/// evaluation the parameterization on some point.
1954/// It's called by TMultiDimFit::MakeCode and TMultiDimFit::MakeMethod.
1955///
1956/// The options are: NONE so far
1957
1959 const char *classname,
1960 Option_t *)
1961{
1962 Int_t i, j;
1963
1964 Bool_t isMethod = (classname[0] == '\0' ? kFALSE : kTRUE);
1965 const char *prefix = (isMethod ? Form("%s::", classname) : "");
1966 const char *cv_qual = (isMethod ? "" : "static ");
1967
1968 std::ofstream outFile(filename,std::ios::out|std::ios::trunc);
1969 if (!outFile) {
1970 Error("MakeRealCode","couldn't open output file '%s'",filename);
1971 return;
1972 }
1973
1974 if (fIsVerbose)
1975 std::cout << "Writing on file \"" << filename << "\" ... " << std::flush;
1976 //
1977 // Write header of file
1978 //
1979 // Emacs mode line ;-)
1980 outFile << "// -*- mode: c++ -*-" << std::endl;
1981 // Info about creator
1982 outFile << "// " << std::endl
1983 << "// File " << filename
1984 << " generated by TMultiDimFit::MakeRealCode" << std::endl;
1985 // Time stamp
1986 TDatime date;
1987 outFile << "// on " << date.AsString() << std::endl;
1988 // ROOT version info
1989 outFile << "// ROOT version " << gROOT->GetVersion()
1990 << std::endl << "//" << std::endl;
1991 // General information on the code
1992 outFile << "// This file contains the function " << std::endl
1993 << "//" << std::endl
1994 << "// double " << prefix << "MDF(double *x); " << std::endl
1995 << "//" << std::endl
1996 << "// For evaluating the parameterization obtained" << std::endl
1997 << "// from TMultiDimFit and the point x" << std::endl
1998 << "// " << std::endl
1999 << "// See TMultiDimFit class documentation for more "
2000 << "information " << std::endl << "// " << std::endl;
2001 // Header files
2002 if (isMethod)
2003 // If these are methods, we need the class header
2004 outFile << "#include \"" << classname << ".h\"" << std::endl;
2005
2006 //
2007 // Now for the data
2008 //
2009 outFile << "//" << std::endl
2010 << "// Static data variables" << std::endl
2011 << "//" << std::endl;
2012 outFile << cv_qual << "int " << prefix << "gNVariables = "
2013 << fNVariables << ";" << std::endl;
2014 outFile << cv_qual << "int " << prefix << "gNCoefficients = "
2015 << fNCoefficients << ";" << std::endl;
2016 outFile << cv_qual << "double " << prefix << "gDMean = "
2017 << fMeanQuantity << ";" << std::endl;
2018
2019 // Assignment to mean vector.
2020 outFile << "// Assignment to mean vector." << std::endl;
2021 outFile << cv_qual << "double " << prefix
2022 << "gXMean[] = {" << std::endl;
2023 for (i = 0; i < fNVariables; i++)
2024 outFile << (i != 0 ? ", " : " ") << fMeanVariables(i) << std::flush;
2025 outFile << " };" << std::endl << std::endl;
2026
2027 // Assignment to minimum vector.
2028 outFile << "// Assignment to minimum vector." << std::endl;
2029 outFile << cv_qual << "double " << prefix
2030 << "gXMin[] = {" << std::endl;
2031 for (i = 0; i < fNVariables; i++)
2032 outFile << (i != 0 ? ", " : " ") << fMinVariables(i) << std::flush;
2033 outFile << " };" << std::endl << std::endl;
2034
2035 // Assignment to maximum vector.
2036 outFile << "// Assignment to maximum vector." << std::endl;
2037 outFile << cv_qual << "double " << prefix
2038 << "gXMax[] = {" << std::endl;
2039 for (i = 0; i < fNVariables; i++)
2040 outFile << (i != 0 ? ", " : " ") << fMaxVariables(i) << std::flush;
2041 outFile << " };" << std::endl << std::endl;
2042
2043 // Assignment to coefficients vector.
2044 outFile << "// Assignment to coefficients vector." << std::endl;
2045 outFile << cv_qual << "double " << prefix
2046 << "gCoefficient[] = {" << std::flush;
2047 for (i = 0; i < fNCoefficients; i++)
2048 outFile << (i != 0 ? "," : "") << std::endl
2049 << " " << fCoefficients(i) << std::flush;
2050 outFile << std::endl << " };" << std::endl << std::endl;
2051
2052 // Assignment to error coefficients vector.
2053 outFile << "// Assignment to error coefficients vector." << std::endl;
2054 outFile << cv_qual << "double " << prefix
2055 << "gCoefficientRMS[] = {" << std::flush;
2056 for (i = 0; i < fNCoefficients; i++)
2057 outFile << (i != 0 ? "," : "") << std::endl
2058 << " " << fCoefficientsRMS(i) << std::flush;
2059 outFile << std::endl << " };" << std::endl << std::endl;
2060
2061 // Assignment to powers vector.
2062 outFile << "// Assignment to powers vector." << std::endl
2063 << "// The powers are stored row-wise, that is" << std::endl
2064 << "// p_ij = " << prefix
2065 << "gPower[i * NVariables + j];" << std::endl;
2066 outFile << cv_qual << "int " << prefix
2067 << "gPower[] = {" << std::flush;
2068 for (i = 0; i < fNCoefficients; i++) {
2069 for (j = 0; j < fNVariables; j++) {
2070 if (j != 0) outFile << std::flush << " ";
2071 else outFile << std::endl << " ";
2072 outFile << fPowers[fPowerIndex[i] * fNVariables + j]
2073 << (i == fNCoefficients - 1 && j == fNVariables - 1 ? "" : ",")
2074 << std::flush;
2075 }
2076 }
2077 outFile << std::endl << "};" << std::endl << std::endl;
2078
2079
2080 //
2081 // Finally we reach the function itself
2082 //
2083 outFile << "// " << std::endl
2084 << "// The "
2085 << (isMethod ? "method " : "function ")
2086 << " double " << prefix
2087 << "MDF(double *x)"
2088 << std::endl << "// " << std::endl;
2089 outFile << "double " << prefix
2090 << "MDF(double *x) {" << std::endl
2091 << " double returnValue = " << prefix << "gDMean;" << std::endl
2092 << " int i = 0, j = 0, k = 0;" << std::endl
2093 << " for (i = 0; i < " << prefix << "gNCoefficients ; i++) {"
2094 << std::endl
2095 << " // Evaluate the ith term in the expansion" << std::endl
2096 << " double term = " << prefix << "gCoefficient[i];"
2097 << std::endl
2098 << " for (j = 0; j < " << prefix << "gNVariables; j++) {"
2099 << std::endl
2100 << " // Evaluate the polynomial in the jth variable." << std::endl
2101 << " int power = "<< prefix << "gPower["
2102 << prefix << "gNVariables * i + j]; " << std::endl
2103 << " double p1 = 1, p2 = 0, p3 = 0, r = 0;" << std::endl
2104 << " double v = 1 + 2. / ("
2105 << prefix << "gXMax[j] - " << prefix
2106 << "gXMin[j]) * (x[j] - " << prefix << "gXMax[j]);" << std::endl
2107 << " // what is the power to use!" << std::endl
2108 << " switch(power) {" << std::endl
2109 << " case 1: r = 1; break; " << std::endl
2110 << " case 2: r = v; break; " << std::endl
2111 << " default: " << std::endl
2112 << " p2 = v; " << std::endl
2113 << " for (k = 3; k <= power; k++) { " << std::endl
2114 << " p3 = p2 * v;" << std::endl;
2115 if (fPolyType == kLegendre)
2116 outFile << " p3 = ((2 * i - 3) * p2 * v - (i - 2) * p1)"
2117 << " / (i - 1);" << std::endl;
2118 if (fPolyType == kChebyshev)
2119 outFile << " p3 = 2 * v * p2 - p1; " << std::endl;
2120 outFile << " p1 = p2; p2 = p3; " << std::endl << " }" << std::endl
2121 << " r = p3;" << std::endl << " }" << std::endl
2122 << " // multiply this term by the poly in the jth var" << std::endl
2123 << " term *= r; " << std::endl << " }" << std::endl
2124 << " // Add this term to the final result" << std::endl
2125 << " returnValue += term;" << std::endl << " }" << std::endl
2126 << " return returnValue;" << std::endl << "}" << std::endl << std::endl;
2127
2128 // EOF
2129 outFile << "// EOF for " << filename << std::endl;
2130
2131 // Close the file
2132 outFile.close();
2133
2134 if (fIsVerbose)
2135 std::cout << "done" << std::endl;
2136}
2137
2138
2139////////////////////////////////////////////////////////////////////////////////
2140/// Print statistics etc.
2141/// Options are
2142/// P Parameters
2143/// S Statistics
2144/// C Coefficients
2145/// R Result of parameterisation
2146/// F Result of fit
2147/// K Correlation Matrix
2148/// M Pretty print formula
2149///
2150
2152{
2153 Int_t i = 0;
2154 Int_t j = 0;
2155
2156 TString opt(option);
2157 opt.ToLower();
2158
2159 if (opt.Contains("p")) {
2160 // Print basic parameters for this object
2161 std::cout << "User parameters:" << std::endl
2162 << "----------------" << std::endl
2163 << " Variables: " << fNVariables << std::endl
2164 << " Data points: " << fSampleSize << std::endl
2165 << " Max Terms: " << fMaxTerms << std::endl
2166 << " Power Limit Parameter: " << fPowerLimit << std::endl
2167 << " Max functions: " << fMaxFunctions << std::endl
2168 << " Max functions to study: " << fMaxStudy << std::endl
2169 << " Max angle (optional): " << fMaxAngle << std::endl
2170 << " Min angle: " << fMinAngle << std::endl
2171 << " Relative Error accepted: " << fMinRelativeError << std::endl
2172 << " Maximum Powers: " << std::flush;
2173 for (i = 0; i < fNVariables; i++)
2174 std::cout << " " << fMaxPowers[i] - 1 << std::flush;
2175 std::cout << std::endl << std::endl
2176 << " Parameterisation will be done using " << std::flush;
2177 if (fPolyType == kChebyshev)
2178 std::cout << "Chebyshev polynomials" << std::endl;
2179 else if (fPolyType == kLegendre)
2180 std::cout << "Legendre polynomials" << std::endl;
2181 else
2182 std::cout << "Monomials" << std::endl;
2183 std::cout << std::endl;
2184 }
2185
2186 if (opt.Contains("s")) {
2187 // Print statistics for read data
2188 std::cout << "Sample statistics:" << std::endl
2189 << "------------------" << std::endl
2190 << " D" << std::flush;
2191 for (i = 0; i < fNVariables; i++)
2192 std::cout << " " << std::setw(10) << i+1 << std::flush;
2193 std::cout << std::endl << " Max: " << std::setw(10) << std::setprecision(7)
2194 << fMaxQuantity << std::flush;
2195 for (i = 0; i < fNVariables; i++)
2196 std::cout << " " << std::setw(10) << std::setprecision(4)
2197 << fMaxVariables(i) << std::flush;
2198 std::cout << std::endl << " Min: " << std::setw(10) << std::setprecision(7)
2199 << fMinQuantity << std::flush;
2200 for (i = 0; i < fNVariables; i++)
2201 std::cout << " " << std::setw(10) << std::setprecision(4)
2202 << fMinVariables(i) << std::flush;
2203 std::cout << std::endl << " Mean: " << std::setw(10) << std::setprecision(7)
2204 << fMeanQuantity << std::flush;
2205 for (i = 0; i < fNVariables; i++)
2206 std::cout << " " << std::setw(10) << std::setprecision(4)
2207 << fMeanVariables(i) << std::flush;
2208 std::cout << std::endl << " Function Sum Squares: " << fSumSqQuantity
2209 << std::endl << std::endl;
2210 }
2211
2212 if (opt.Contains("r")) {
2213 std::cout << "Results of Parameterisation:" << std::endl
2214 << "----------------------------" << std::endl
2215 << " Total reduction of square residuals "
2216 << fSumSqResidual << std::endl
2217 << " Relative precision obtained: "
2218 << fPrecision << std::endl
2219 << " Error obtained: "
2220 << fError << std::endl
2221 << " Multiple correlation coefficient: "
2222 << fCorrelationCoeff << std::endl
2223 << " Reduced Chi square over sample: "
2224 << fChi2 / (fSampleSize - fNCoefficients) << std::endl
2225 << " Maximum residual value: "
2226 << fMaxResidual << std::endl
2227 << " Minimum residual value: "
2228 << fMinResidual << std::endl
2229 << " Estimated root mean square: "
2230 << fRMS << std::endl
2231 << " Maximum powers used: " << std::flush;
2232 for (j = 0; j < fNVariables; j++)
2233 std::cout << fMaxPowersFinal[j] << " " << std::flush;
2234 std::cout << std::endl
2235 << " Function codes of candidate functions." << std::endl
2236 << " 1: considered,"
2237 << " 2: too little contribution,"
2238 << " 3: accepted." << std::flush;
2239 for (i = 0; i < fMaxFunctions; i++) {
2240 if (i % 60 == 0)
2241 std::cout << std::endl << " " << std::flush;
2242 else if (i % 10 == 0)
2243 std::cout << " " << std::flush;
2244 std::cout << fFunctionCodes[i];
2245 }
2246 std::cout << std::endl << " Loop over candidates stopped because " << std::flush;
2247 switch(fParameterisationCode){
2248 case PARAM_MAXSTUDY:
2249 std::cout << "max allowed studies reached" << std::endl; break;
2250 case PARAM_SEVERAL:
2251 std::cout << "all candidates considered several times" << std::endl; break;
2252 case PARAM_RELERR:
2253 std::cout << "wanted relative error obtained" << std::endl; break;
2254 case PARAM_MAXTERMS:
2255 std::cout << "max number of terms reached" << std::endl; break;
2256 default:
2257 std::cout << "some unknown reason" << std::endl;
2258 break;
2259 }
2260 std::cout << std::endl;
2261 }
2262
2263 if (opt.Contains("f")) {
2264 std::cout << "Results of Fit:" << std::endl
2265 << "---------------" << std::endl
2266 << " Test sample size: "
2267 << fTestSampleSize << std::endl
2268 << " Multiple correlation coefficient: "
2269 << fTestCorrelationCoeff << std::endl
2270 << " Relative precision obtained: "
2271 << fTestPrecision << std::endl
2272 << " Error obtained: "
2273 << fTestError << std::endl
2274 << " Reduced Chi square over sample: "
2275 << fChi2 / (fSampleSize - fNCoefficients) << std::endl
2276 << std::endl;
2277 if (fFitter) {
2278 fFitter->PrintResults(1,1);
2279 std::cout << std::endl;
2280 }
2281 }
2282
2283 if (opt.Contains("c")){
2284 std::cout << "Coefficients:" << std::endl
2285 << "-------------" << std::endl
2286 << " # Value Error Powers" << std::endl
2287 << " ---------------------------------------" << std::endl;
2288 for (i = 0; i < fNCoefficients; i++) {
2289 std::cout << " " << std::setw(3) << i << " "
2290 << std::setw(12) << fCoefficients(i) << " "
2291 << std::setw(12) << fCoefficientsRMS(i) << " " << std::flush;
2292 for (j = 0; j < fNVariables; j++)
2293 std::cout << " " << std::setw(3)
2294 << fPowers[fPowerIndex[i] * fNVariables + j] - 1 << std::flush;
2295 std::cout << std::endl;
2296 }
2297 std::cout << std::endl;
2298 }
2299 if (opt.Contains("k") && fCorrelationMatrix.IsValid()) {
2300 std::cout << "Correlation Matrix:" << std::endl
2301 << "-------------------";
2303 }
2304
2305 if (opt.Contains("m")) {
2306 std::cout << "Parameterization:" << std::endl
2307 << "-----------------" << std::endl
2308 << " Normalised variables: " << std::endl;
2309 for (i = 0; i < fNVariables; i++)
2310 std::cout << "\ty_" << i << "\t= 1 + 2 * (x_" << i << " - "
2311 << fMaxVariables(i) << ") / ("
2312 << fMaxVariables(i) << " - " << fMinVariables(i) << ")"
2313 << std::endl;
2314 std::cout << std::endl
2315 << " f(";
2316 for (i = 0; i < fNVariables; i++) {
2317 std::cout << "y_" << i;
2318 if (i != fNVariables-1) std::cout << ", ";
2319 }
2320 std::cout << ") = ";
2321 for (i = 0; i < fNCoefficients; i++) {
2322 if (i != 0)
2323 std::cout << std::endl << "\t" << (fCoefficients(i) < 0 ? "- " : "+ ")
2325 else
2326 std::cout << fCoefficients(i);
2327 for (j = 0; j < fNVariables; j++) {
2329 switch (p) {
2330 case 1: break;
2331 case 2: std::cout << " * y_" << j; break;
2332 default:
2333 switch(fPolyType) {
2334 case kLegendre: std::cout << " * L_" << p-1 << "(y_" << j << ")"; break;
2335 case kChebyshev: std::cout << " * C_" << p-1 << "(y_" << j << ")"; break;
2336 default: std::cout << " * y_" << j << "^" << p-1; break;
2337 }
2338 }
2339
2340 }
2341 }
2342 std::cout << std::endl;
2343 }
2344}
2345
2346
2347////////////////////////////////////////////////////////////////////////////////
2348/// Selection method. User can override this method for specialized
2349/// selection of acceptable functions in fit. Default is to select
2350/// all. This message is sent during the build-up of the function
2351/// candidates table once for each set of powers in
2352/// variables. Notice, that the argument array contains the powers
2353/// PLUS ONE. For example, to De select the function
2354/// f = x1^2 * x2^4 * x3^5,
2355/// this method should return kFALSE if given the argument
2356/// { 3, 4, 6 }
2357
2359{
2360 return kTRUE;
2361}
2362
2363////////////////////////////////////////////////////////////////////////////////
2364/// Set the max angle (in degrees) between the initial data vector to
2365/// be fitted, and the new candidate function to be included in the
2366/// fit. By default it is 0, which automatically chooses another
2367/// selection criteria. See also
2368/// <a href="#TMultiDimFit:description">class description</a>
2369
2371{
2372 if (ang >= 90 || ang < 0) {
2373 Warning("SetMaxAngle", "angle must be in [0,90)");
2374 return;
2375 }
2376
2377 fMaxAngle = ang;
2378}
2379
2380////////////////////////////////////////////////////////////////////////////////
2381/// Set the min angle (in degrees) between a new candidate function
2382/// and the subspace spanned by the previously accepted
2383/// functions. See also
2384/// <a href="#TMultiDimFit:description">class description</a>
2385
2387{
2388 if (ang > 90 || ang <= 0) {
2389 Warning("SetMinAngle", "angle must be in [0,90)");
2390 return;
2391 }
2392
2393 fMinAngle = ang;
2394
2395}
2396
2397
2398////////////////////////////////////////////////////////////////////////////////
2399/// Define a user function. The input array must be of the form
2400/// (p11, ..., p1N, ... ,pL1, ..., pLN)
2401/// Where N is the dimension of the data sample, L is the number of
2402/// terms (given in terms) and the first number, labels the term, the
2403/// second the variable. More information is given in the
2404/// <a href="#TMultiDimFit:description">class description</a>
2405
2406void TMultiDimFit::SetPowers(const Int_t* powers, Int_t terms)
2407{
2409 fMaxFunctions = terms;
2410 fMaxTerms = terms;
2411 fMaxStudy = terms;
2413 fPowers = new Int_t[fMaxFuncNV];
2414 Int_t i, j;
2415 for (i = 0; i < fMaxFunctions; i++)
2416 for(j = 0; j < fNVariables; j++)
2417 fPowers[i * fNVariables + j] = powers[i * fNVariables + j] + 1;
2418}
2419
2420////////////////////////////////////////////////////////////////////////////////
2421/// Set the user parameter for the function selection. The bigger the
2422/// limit, the more functions are used. The meaning of this variable
2423/// is defined in the
2424/// <a href="#TMultiDimFit:description">class description</a>
2425
2427{
2428 fPowerLimit = limit;
2429}
2430
2431////////////////////////////////////////////////////////////////////////////////
2432/// Set the maximum power to be considered in the fit for each
2433/// variable. See also
2434/// <a href="#TMultiDimFit:description">class description</a>
2435
2437{
2438 if (!powers)
2439 return;
2440
2441 for (Int_t i = 0; i < fNVariables; i++)
2442 fMaxPowers[i] = powers[i]+1;
2443}
2444
2445////////////////////////////////////////////////////////////////////////////////
2446/// Set the acceptable relative error for when sum of square
2447/// residuals is considered minimized. For a full account, refer to
2448/// the
2449/// <a href="#TMultiDimFit:description">class description</a>
2450
2452{
2453 fMinRelativeError = error;
2454}
2455
2456
2457////////////////////////////////////////////////////////////////////////////////
2458/// PRIVATE METHOD:
2459/// Test whether the currently considered function contributes to the
2460/// fit. See also
2461/// <a href="#TMultiDimFit:description">class description</a>
2462
2464 Double_t dResidur)
2465{
2466 if (fNCoefficients != 0) {
2467 // Now for the second test:
2468 if (fMaxAngle == 0) {
2469 // If the user hasn't supplied a max angle do the test as,
2470 if (dResidur <
2471 squareResidual / (fMaxTerms - fNCoefficients + 1 + 1E-10)) {
2472 return kFALSE;
2473 }
2474 }
2475 else {
2476 // If the user has provided a max angle, test if the calculated
2477 // angle is less then the max angle.
2478 if (TMath::Sqrt(dResidur/fSumSqAvgQuantity) <
2480 return kFALSE;
2481 }
2482 }
2483 }
2484 // If we get here, the function is OK
2485 return kTRUE;
2486}
2487
2488
2489////////////////////////////////////////////////////////////////////////////////
2490/// Helper function for doing the minimisation of Chi2 using Minuit
2491
2492void mdfHelper(int& /*npar*/, double* /*divs*/, double& chi2,
2493 double* coeffs, int /*flag*/)
2494{
2495 // Get pointer to current TMultiDimFit object.
2497 chi2 = mdf->MakeChi2(coeffs);
2498}
#define b(i)
Definition RSha256.hxx:100
#define f(i)
Definition RSha256.hxx:104
#define h(i)
Definition RSha256.hxx:106
#define e(i)
Definition RSha256.hxx:103
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
int Int_t
Definition RtypesCore.h:45
char Char_t
Definition RtypesCore.h:37
constexpr Bool_t kFALSE
Definition RtypesCore.h:101
double Double_t
Definition RtypesCore.h:59
constexpr Bool_t kTRUE
Definition RtypesCore.h:100
const char Option_t
Definition RtypesCore.h:66
#define ClassImp(name)
Definition Rtypes.h:377
#define SETBIT(n, i)
Definition Rtypes.h:86
#define TESTBIT(n, i)
Definition Rtypes.h:88
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t option
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char filename
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t Atom_t Atom_t Time_t type
TMatrixTRow< Double_t > TMatrixDRow
TMatrixTDiag_const< Double_t > TMatrixDDiag_const
#define HIST_DSHIF
#define PARAM_RELERR
#define HIST_DORIG
#define HIST_XNORM
#define HIST_RX
#define HIST_XORIG
#define PARAM_SEVERAL
#define PARAM_MAXSTUDY
static void mdfHelper(int &, double *, double &, double *, int)
Helper function for doing the minimisation of Chi2 using Minuit.
#define HIST_RTEST
#define DEGRAD
#define HIST_RTRAI
#define HIST_RD
#define PARAM_MAXTERMS
#define gROOT
Definition TROOT.h:405
char * Form(const char *fmt,...)
Formats a string in a circular formatting buffer.
Definition TString.cxx:2467
Using a TBrowser one can browse all ROOT objects.
Definition TBrowser.h:37
This class stores the date and time with a precision of one second in an unsigned 32 bit word (950130...
Definition TDatime.h:37
const char * AsString() const
Return the date & time as a string (ctime() format).
Definition TDatime.cxx:102
Cholesky Decomposition class.
Definition TDecompChol.h:25
Bool_t Decompose() override
Matrix A is decomposed in component U so that A = U^T * U If the decomposition succeeds,...
Bool_t Invert(TMatrixDSym &inv)
For a symmetric matrix A(m,m), its inverse A_inv(m,m) is returned .
1-D histogram with a double per channel (see TH1 documentation)}
Definition TH1.h:620
TH1 is the base class of all histogram classes in ROOT.
Definition TH1.h:58
virtual Bool_t Add(TF1 *h1, Double_t c1=1, Option_t *option="")
Performs the operation: this = this + c1*f1 if errors are defined (see TH1::Sumw2),...
Definition TH1.cxx:807
2-D histogram with a double per channel (see TH1 documentation)}
Definition TH2.h:300
A doubly linked list.
Definition TList.h:38
void Clear(Option_t *option="") override
Remove all objects from the list.
Definition TList.cxx:402
TObject * FindObject(const char *name) const override
Find an object in this list using its name.
Definition TList.cxx:578
void Add(TObject *obj) override
Definition TList.h:81
virtual TMatrixTBase< Element > & Zero()
Set matrix elements to zero.
void Print(Option_t *name="") const override
Print the matrix as a table of elements.
Bool_t IsValid() const
virtual TMatrixTBase< Element > & NormByDiag(const TVectorT< Element > &v, Option_t *option="D")
option:
TMatrixTBase< Element > & ResizeTo(Int_t nrows, Int_t ncols, Int_t=-1) override
Set size of the matrix to nrows x ncols New dynamic elements are created, the overlapping part of the...
Multidimensional Fits in ROOT.
void Print(Option_t *option="ps") const override
Print statistics etc.
virtual void MakeCorrelation()
PRIVATE METHOD: Compute the correlation matrix.
TVectorD fCoefficients
Vector of the final coefficients.
Double_t fPrecision
Relative precision of param.
Byte_t fHistogramMask
Bit pattern of histograms used.
virtual void MakeCoefficientErrors()
PRIVATE METHOD: Compute the errors on the coefficients.
TMatrixD fOrthCurvatureMatrix
Model matrix.
Bool_t fShowCorrelation
print correlation matrix
Bool_t fIsUserFunction
Flag for user defined function.
Int_t fNVariables
Number of independent variables.
Int_t fMaxResidualRow
Row giving max residual.
Double_t fMinQuantity
Min value of dependent quantity.
virtual Double_t MakeChi2(const Double_t *coeff=nullptr)
Calculate Chi square over either the test sample.
Double_t fSumSqQuantity
SumSquare of dependent quantity.
TVectorD fMaxVariables
max value of independent variables
Int_t fSampleSize
Size of training sample.
Double_t fSumSqAvgQuantity
Sum of squares away from mean.
Int_t fMaxTerms
Max terms expected in final expr.
Int_t * fPowers
[fMaxFuncNV] where fMaxFuncNV = fMaxFunctions*fNVariables
virtual void MakeNormalized()
PRIVATE METHOD: Normalize data to the interval [-1;1].
Int_t * fMaxPowers
[fNVariables] maximum powers
Int_t fMaxStudy
max functions to study
TVectorD fTestQuantity
Test sample, dependent quantity.
TVectorD fTestSqError
Test sample, Error in quantity.
virtual Double_t Eval(const Double_t *x, const Double_t *coeff=nullptr) const
Evaluate parameterization at point x.
Int_t fMaxFunctions
max number of functions
Int_t fBinVarY
Number of bin in dependent variables.
static TMultiDimFit * Instance()
Return the static instance.
virtual void MakeCode(const char *functionName="MDF", Option_t *option="")
Generate the file <filename> with .C appended if argument doesn't end in .cxx or ....
virtual void MakeCandidates()
PRIVATE METHOD: Create list of candidate functions for the parameterisation.
Double_t fTestError
Error from test.
Double_t fMinResidual
Min residual value.
void SetPowerLimit(Double_t limit=1e-3)
Set the user parameter for the function selection.
virtual void FindParameterization(Option_t *option="")
Find the parameterization.
virtual Double_t MakeGramSchmidt(Int_t function)
PRIVATE METHOD: Make Gram-Schmidt orthogonalisation.
TVectorD fSqError
Training sample, error in quantity.
virtual void AddRow(const Double_t *x, Double_t D, Double_t E=0)
Add a row consisting of fNVariables independent variables, the known, dependent quantity,...
TVectorD fQuantity
Training sample, dependent quantity.
TVectorD fTestVariables
Test sample, independent variables.
TVectorD fOrthCoefficients
The model coefficients.
Int_t fTestSampleSize
Size of test sample.
TVectorD fResiduals
Vector of the final residuals.
virtual Bool_t TestFunction(Double_t squareResidual, Double_t dResidur)
PRIVATE METHOD: Test whether the currently considered function contributes to the fit.
TVectorD fVariables
Training sample, independent variables.
Double_t fCorrelationCoeff
Multi Correlation coefficient.
void SetMaxPowers(const Int_t *powers)
Set the maximum power to be considered in the fit for each variable.
Int_t * fMaxPowersFinal
[fNVariables] maximum powers from fit;
Int_t fMinResidualRow
Row giving min residual.
virtual void SetPowers(const Int_t *powers, Int_t terms)
Define a user function.
Double_t fMaxQuantity
Max value of dependent quantity.
virtual void MakeRealCode(const char *filename, const char *classname, Option_t *option="")
PRIVATE METHOD: This is the method that actually generates the code for the evaluation the parameteri...
Double_t fRMS
Root mean square of fit.
TVirtualFitter * fFitter
TMatrixD fFunctions
Functions evaluated over sample.
virtual void AddTestRow(const Double_t *x, Double_t D, Double_t E=0)
Add a row consisting of fNVariables independent variables, the known, dependent quantity,...
Int_t fParameterisationCode
Exit code of parameterisation.
Double_t fMaxAngle
Max angle for accepting new function.
Int_t fNCoefficients
Dimension of model coefficients.
virtual void Fit(Option_t *option="")
Try to fit the found parameterisation to the test sample.
TMultiDimFit()
Empty CTOR. Do not use.
Double_t fTestCorrelationCoeff
Multi Correlation coefficient.
Double_t fMinAngle
Min angle for accepting new function.
void SetMinRelativeError(Double_t error)
Set the acceptable relative error for when sum of square residuals is considered minimized.
Double_t fSumSqResidual
Sum of Square residuals.
TVectorD fMeanVariables
mean value of independent variables
virtual void MakeMethod(const Char_t *className="MDF", Option_t *option="")
Generate the file <classname>MDF.cxx which contains the implementation of the method:
void SetMaxAngle(Double_t angle=0)
Set the max angle (in degrees) between the initial data vector to be fitted, and the new candidate fu...
Double_t fTestPrecision
Relative precision of test.
virtual void MakeCoefficients()
PRIVATE METHOD: Invert the model matrix B, and compute final coefficients.
TMatrixD fCorrelationMatrix
Correlation matrix.
Int_t * fPowerIndex
[fMaxTerms] Index of accepted powers
TMatrixD fOrthFunctions
As above, but orthogonalised.
virtual Double_t EvalControl(const Int_t *powers) const
PRIVATE METHOD: Calculate the control parameter from the passed powers.
void Clear(Option_t *option="") override
Clear internal structures and variables.
virtual void MakeParameterization()
PRIVATE METHOD: Find the parameterization over the training sample.
EMDFPolyType fPolyType
Fit object (MINUIT)
void Browse(TBrowser *b) override
Browse the TMultiDimFit object in the TBrowser.
virtual Bool_t Select(const Int_t *iv)
Selection method.
TList * fHistograms
List of histograms.
Double_t fChi2
Chi square of fit.
Double_t fMaxResidual
Max residual value.
Double_t fError
Error from parametrization.
Int_t fBinVarX
Number of bin in independent variables.
virtual Double_t EvalFactor(Int_t p, Double_t x) const
PRIVATE METHOD: Evaluate function with power p at variable value x.
void SetMinAngle(Double_t angle=1)
Set the min angle (in degrees) between a new candidate function and the subspace spanned by the previ...
virtual Double_t EvalError(const Double_t *x, const Double_t *coeff=nullptr) const
Evaluate parameterization error at point x.
TVectorD fMinVariables
min value of independent variables
TVectorD fOrthFunctionNorms
Norm of the evaluated functions.
Double_t fMeanQuantity
Mean of dependent quantity.
TVectorD fCoefficientsRMS
Vector of RMS of coefficients.
Double_t fPowerLimit
Control parameter.
static TMultiDimFit * fgInstance
~TMultiDimFit() override
Destructor.
virtual void MakeHistograms(Option_t *option="A")
Make histograms of the result of the analysis.
Bool_t fIsVerbose
Int_t fMaxFuncNV
fMaxFunctions*fNVariables
Double_t fMinRelativeError
Min relative error accepted.
Int_t * fFunctionCodes
[fMaxFunctions] acceptance code
The TNamed class is the base class for all named ROOT classes.
Definition TNamed.h:29
const char * GetName() const override
Returns name of object.
Definition TNamed.h:47
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition TObject.cxx:956
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition TObject.cxx:970
Basic string class.
Definition TString.h:139
Ssiz_t Length() const
Definition TString.h:421
void ToLower()
Change string to lower-case.
Definition TString.cxx:1170
Bool_t EndsWith(const char *pat, ECaseCompare cmp=kExact) const
Return true if string ends with the specified string.
Definition TString.cxx:2222
const char * Data() const
Definition TString.h:380
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition TString.h:636
TVectorT< Element > & Zero()
Set vector elements to zero.
Definition TVectorT.cxx:453
TVectorT< Element > & ResizeTo(Int_t lwb, Int_t upb)
Resize the vector to [lwb:upb] .
Definition TVectorT.cxx:294
Int_t GetNrows() const
Definition TVectorT.h:75
Bool_t IsValid() const
Definition TVectorT.h:83
virtual void PrintResults(Int_t level, Double_t amin) const =0
virtual void SetFCN(void(*fcn)(Int_t &, Double_t *, Double_t &f, Double_t *, Int_t))
To set the address of the minimization objective function called by the native compiler (see function...
virtual Int_t SetParameter(Int_t ipar, const char *parname, Double_t value, Double_t verr, Double_t vlow, Double_t vhigh)=0
virtual Int_t ExecuteCommand(const char *command, Double_t *args, Int_t nargs)=0
virtual Double_t GetParameter(Int_t ipar) const =0
static TVirtualFitter * Fitter(TObject *obj, Int_t maxpar=25)
Static function returning a pointer to the current fitter.
Double_t y[n]
Definition legend1.C:17
Double_t x[n]
Definition legend1.C:17
const Int_t n
Definition legend1.C:16
Short_t Max(Short_t a, Short_t b)
Returns the largest of a and b.
Definition TMathBase.h:250
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition TMath.h:660
Double_t Cos(Double_t)
Returns the cosine of an angle of x radians.
Definition TMath.h:592
Double_t Sin(Double_t)
Returns the sine of an angle of x radians.
Definition TMath.h:586
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:123
TMarker m
Definition textangle.C:8
TLine l
Definition textangle.C:4
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2345
double epsilon
Definition triangle.c:618