1 package org.lcsim.recon.tracking.trfutil;
2 import java.util.Random;
3 /** This class smears a vector of quantities in
4 * accordance with a multivariate normal distribution.
5 * The constructor takes a covariance matrix as input;
6 * only a minimal amount of checking is done.
7 * The method 'smear' takes an array and smears its
8 * members according to a gaussian distribution, whose
9 * width corresponds to the square root of the diagonal
10 * element of the covariance matrix. Here is an example
11 * of how the covariance terms are handled:
12 *<pre>
13 * v11
14 * cov = v12 v22
15 * v13 v23 v33
16 *
17 * x1 = x1 + a1*rn1
18 * x2 = x2 + a2*rn1 + b2*rn2
19 * x3 = x3 + a3*rn1 + b3*rn2 + c3*rn3
20 *</pre>
21 *<p>
22 * where rn1, rn2 and rn3 are independently generated normal random numbers
23 *
24 *<pre>
25 * a1 = sqrt(v11)
26 * if a1 = 0, then a2 = a3 = 0
27 * otherwise a2 = v12/a1, a3 = v13/a1
28 *
29 * b2 = sqrt(v22 - a2*a2)
30 * if b2 = 0, then b3 = 0
31 * otherwise b3 = (v23 - a2*a3)/b2
32 *
33 * c3 = sqrt(v33 - a3*a3 - b3*b3)
34 *</pre>
35 *
36 *@author Norman A. Graf
37 *@version 1.0
38 *
39 */
40
41
42 public class ArraySmearer
43 {
44 //attributes
45 // store the matrix used later for smearing
46 private double[][] _smear;
47 //store a copy of the covariance matrix for printing
48 // could be eliminated for efficiency
49 private double[][] _cov;
50 // the size of the matrix and arrays to be smeared
51 private int _ndim;
52 // the random number generator
53 private Random _r = new Random();
54
55 // methods
56
57 /**
58 *constructor from the covariance matrix
59 *
60 * @param cov The covariance matrix as
61 */
62 public ArraySmearer( double[][] cov)
63 {
64 // sanity checks...
65 _ndim = cov.length;
66 // dimension
67 for(int i = 0; i< _ndim; ++i)
68 {
69 if(cov[i].length != _ndim) throw new IllegalArgumentException("Array Smearer: Bad covariance matrix!");
70 }
71 // symmetry
72 for(int i = 0; i < _ndim; ++i)
73 {
74 for (int j = 0; j< _ndim ; ++j )
75 {
76 if( cov[i][j] != cov[j][i] ) throw new IllegalArgumentException("Array Smearer: Bad covariance matrix!");
77 }
78 }
79 // should also check the determinant, punt for now
80
81 _smear = new double[_ndim][_ndim];
82 // store local copy, could be elimated
83 _cov = new double[_ndim][_ndim];
84 System.arraycopy(cov,0,_cov,0,_ndim);
85 // calculate the elements of the smearing matrix
86 double ck = 0.;
87 double sum;
88 for (int i=0; i<_ndim; i++)
89 {
90 // Diagonal terms...
91 sum = 0;
92 for (int j=0; j<i; j++)
93 {
94 sum += _smear[i][j]*_smear[i][j];
95 }
96 _smear[i][i] = Math.sqrt(Math.abs(cov[i][i] - sum));
97 // Off-Diagonal terms...
98 for (int k=i+1; k<_ndim; k++)
99 {
100 sum = 0;
101 for (int j=0; j<i; j++)
102 {
103 sum += _smear[k][j]*_smear[i][j];
104 }
105 _smear[k][i] = (cov[k][i] - sum)/_smear[i][i];
106 }
107 }
108 }
109
110 /**
111 * Smear a vector according to the covariance matrix.
112 * Note overwrite of input argument!
113 * @param vec The array to be smeared.
114 */
115 public void smear(double[] vec)
116 {
117 if(vec.length!=_ndim) throw new IllegalArgumentException("Array Smearer: Bad input vector!");
118 double smrfac;
119
120 for (int i=0; i<_ndim; i++)
121 {
122 smrfac = 0.0;
123 for (int j=0; j<=i; j++)
124 {
125 smrfac += _smear[i][j]*_r.nextGaussian();
126 }
127 vec[i] += smrfac;
128 }
129 }
130
131 // generate random vector according to the covariance matrix
132
133 /**
134 * Generate a vector of variables distributed according to the covariance matrix.
135 *
136 * @param vec Array to be filled
137 */
138 public void generate(double[] vec)
139 {
140 if(vec.length!=_ndim) throw new IllegalArgumentException("Array Generator: Bad input vector!");
141 double[] z = new double[_ndim];
142 // generate a vector of uncorrelated gaussians
143 for (int i=0; i<_ndim; i++)
144 {
145 z[i] = _r.nextGaussian();
146 }
147 // add correlations
148 for (int i=0; i<_ndim; i++)
149 {
150 vec[i]=0.;
151 for (int j=0; j<=i; j++)
152 {
153 vec[i] += _smear[i][j]*z[j];
154 }
155
156 }
157 }
158
159
160 /**
161 *Output Stream
162 *
163 * @return String representation of object
164 */
165 public String toString()
166 {
167 String className = getClass().getName();
168 int lastDot = className.lastIndexOf('.');
169 if(lastDot!=-1)className = className.substring(lastDot+1);
170 StringBuffer sb = new StringBuffer(className+" of dimension "+_ndim+"\n");
171 for ( int i = 0; i<_ndim ; ++i )
172 {
173 for (int j = 0 ; j<_ndim ; ++j )
174 {
175 sb.append( _cov[i][j] +" ");
176 }
177 sb.append("\n");
178 }
179 sb.append("\n smear= \n");
180 for ( int i = 0; i<_ndim ; ++i )
181 {
182 for (int j = 0 ; j<_ndim ; ++j )
183 {
184 sb.append( _smear[i][j] +" ");
185 }
186 sb.append("\n");
187 }
188 return sb.toString();
189 }
190 }