1 package org.lcsim.recon.emid.hmatrix;
2 import java.io.*;
3
4
5
6
7
8
9 public class HMatrixBuilder
10 {
11 private double[][] _cov;
12 private double[] _vec;
13 private double[] _vals;
14 private double[][] _valsq;
15 private int _dim;
16 private int _key;
17 private int _ndata;
18 private boolean _isValid;
19 private HMatrix _hmx;
20
21
22
23
24
25
26
27
28
29 public HMatrixBuilder(int n, int key)
30 {
31 _dim = n;
32 _key = key;
33 _cov = new double[_dim][_dim];
34 _vec = new double[_dim];
35 _vals = new double[_dim];
36 _valsq = new double[_dim][_dim];
37 _ndata = 0;
38 _isValid = false;
39 }
40
41
42
43
44
45
46 public void accumulate(double[] dat)
47 {
48
49 for(int i=0; i<_dim; ++i)
50 {
51 for(int j=0; j<_dim; ++j)
52 {
53 _valsq[i][j]+=dat[i]*dat[j];
54 if(i==j) _vals[i]+=dat[i];
55
56 }
57 }
58 _ndata++;
59 }
60
61
62
63
64 public void validate()
65 {
66 double[] _covDiag = new double[_dim];
67 for(int i=0; i<_dim; ++i)
68 {
69 _vec[i] = _vals[i]/_ndata;
70 for(int j=0; j<_dim; ++j)
71 {
72 double data = (double)_ndata;
73 _cov[i][j] = _valsq[i][j]/(data) - (_vals[i]/(data)*(_vals[j]/(data)));
74 }
75 _covDiag[i] = _cov[i][i];
76 }
77
78 _hmx = new HMatrix(_dim, _key, _vec, _covDiag, invert(_cov,_dim));
79 _isValid = true;
80 }
81
82
83
84
85
86
87 public String toString()
88 {
89 StringBuffer sb = new StringBuffer("HMatrixBuilder: dimension "+_dim+" ");
90 if(!_isValid) return (sb.append("INVALID").toString());
91 sb.append(_hmx);
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107 return sb.toString();
108 }
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147 public HMatrix read(String filename)
148 {
149 return HMatrix.read(filename);
150 }
151
152
153
154
155
156
157
158 public void write(String filename, String comment)
159 {
160 _hmx.write(filename, comment);
161 }
162
163
164
165
166 private double[][] invert(double[][] mat, int dim)
167 {
168
169 int i = dim;
170 double[][] matrix = new double[i][2 * i];
171 double[][] matrix1 = new double[i][1];
172 double[][] matrix2 = new double[i][1];
173 for(int j = 0; j < i; j++)
174 {
175 for(int i2 = 0; i2 < i; i2++)
176 matrix[j][i2]= mat[j][i2];
177
178 }
179
180 for(int k = 0; k < i; k++)
181 {
182 for(int j2 = i; j2 < 2 * i; j2++)
183 matrix[k][j2]= 0.0D;
184
185 matrix[k][i + k]=1.0D;
186 }
187
188 for(int i4 = 0; i4 < i; i4++)
189 {
190 for(int l = i4; l < i; l++)
191 {
192 matrix2[l][0]=0.0D;
193 for(int k2 = i4; k2 < i; k2++)
194 matrix2[l][0]=matrix2[l][0] + matrix[l][k2];
195
196 matrix1[l][0]= Math.abs(matrix[l][i4]) / matrix2[l][0];
197 }
198
199 int k4 = i4;
200 for(int i1 = i4 + 1; i1 < i; i1++)
201 if(matrix1[i1][0] > matrix1[i1 - 1][0])
202 k4 = i1;
203
204 if(matrix1[k4][0] == 0.0D)
205 {
206 System.err.println("Matrix::inverse -> Matrix is singular.");
207 System.exit(0);
208 }
209 if(k4 != i4)
210 {
211 for(int l2 = 0; l2 < 2 * i; l2++)
212 {
213 double d = matrix[i4][l2];
214 matrix[i4][l2]=matrix[k4][l2];
215 matrix[k4][l2]=d;
216 }
217
218 }
219 for(int i3 = 2 * i - 1; i3 >= i4; i3--)
220 matrix[i4][i3]=matrix[i4][i3] / matrix[i4][i4];
221
222 for(int j1 = i4 + 1; j1 < i; j1++)
223 {
224 for(int j3 = 2 * i - 1; j3 >= i4 + 1; j3--)
225 matrix[j1][j3] =matrix[j1][j3] - matrix[i4][j3] * matrix[j1][i4];
226
227 }
228
229 }
230
231 for(int j4 = i - 1; j4 >= 0; j4--)
232 {
233 for(int k1 = j4 - 1; k1 >= 0; k1--)
234 {
235 for(int k3 = i; k3 < 2 * i; k3++)
236 matrix[k1][k3] =matrix[k1][k3] - matrix[j4][k3] * matrix[k1][j4];
237
238 }
239
240 }
241
242 double[][] matrix3 = new double[i][i];
243 for(int l1 = 0; l1 < i; l1++)
244 {
245 for(int l3 = 0; l3 < i; l3++)
246 matrix3[l1][l3] =matrix[l1][l3 + i];
247
248 }
249
250 return matrix3;
251 }
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308 }