View Javadoc

1   package org.lcsim.recon.emid.hmatrix;
2   import java.io.*;
3   /** A class to construct an HMatrix.
4    *
5    *@author Norman A. Graf
6    *@version 1.0
7    *@see HMatrix
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      /** Construct an <b>n</b> x <b>n</b> HMatrixBuilder.
22       * Leaves HMatrix in an invalid state.
23       * One must either accumulate entries to build up the HMatrix
24       * or read in the HMatrix values from a file.
25       *
26       * @param   n  The dimension of the measurement space.
27       * @param   key The key by which to index this HMatrix.
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       * Add the measurement vector to the accumulated HMatrix
43       *
44       * @param   dat  The array of measurements.
45       */
46      public void accumulate(double[] dat)
47      {
48          //if(dat.length !=_dim) throw new InvalidArgumentException();
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                  //System.out.println("i= "+i+", j= "+j+"val= "+_vals[i]+", valsq= "+_valsq[i][j]);
56              }
57          }
58          _ndata++;
59      }
60      
61      /**
62       * Validates the HMatrix by performing the averages
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       * Output Stream
84       *
85       * @return  String representation of the object
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  /*		sb.append("\n\nvector: \n");
93                  for(int i = 0; i<_dim; ++i)
94                  {
95                          sb.append(_vec[i]+" ");
96                  }
97                  sb.append("\n\ncovariance matrix: \n");
98                  for(int i = 0; i<_dim; ++i)
99                  {
100                         for(int j = 0; j<_dim; ++j)
101                         {
102                                 sb.append(_cov[i][j]+" ");
103                         }
104                         sb.append("\n");
105                 }
106  */
107         return sb.toString();
108     }
109     
110     /**
111      * Writes the HMatrix to a Serialized file
112      *
113      *  @param   filename  The file to which to write.
114      */
115     /*
116     public void writeSerialized(String filename)
117     {
118         if (_isValid)
119         {
120             _hmx.writeSerialized(filename);
121         }
122         else
123         {
124             //		throw new Exception();
125             System.out.println("HMatrix not valid! Cannot write!");
126         }
127     }
128     */
129     /**
130      * Reads the HMatrix from a Serialized file
131      *
132      * @param   filename The Serialized file from which to read.
133      * @return  The HMatrix
134      */
135     /*
136     public HMatrix readSerialized(String filename)
137     {
138         return HMatrix.readSerialized(filename);
139     }
140     */
141     /**
142      * Reads the HMatrix from an ASCII file
143      *
144      * @param   filename The ASCII file from which to read.
145      * @return  The HMatrix
146      */
147     public HMatrix read(String filename)
148     {
149         return HMatrix.read(filename);
150     }
151     
152     /**
153      *  Writes out the HMatrix to an ASCII file
154      *
155      * @param   filename  The file to which to write.
156      * @param   comment   A comment for the header.
157      */
158     public void write(String filename, String comment)
159     {
160         _hmx.write(filename, comment);
161     }
162     
163     /**
164      * Invert the matrix
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 // Inversion in place requires larger matrix.
255         private double[][] invert( double[][] D, int n)
256         {
257                 double alpha;
258                 double beta;
259                 int i;
260                 int j;
261                 int k;
262                 int error;
263  
264                 error = 0;
265                 int n2 = 2*n;
266  
267                 // initialize the reduction matrix
268                 for( i = 1; i <= n; i++ )
269                 {
270                         for( j = 1; j <= n; j++ )
271                         {
272                                 D[i][j+n] = 0.;
273                         }
274                         D[i][i+n] = 1.0;
275                 }
276  
277                 // perform the reductions
278                 for( i = 1; i <= n; i++ )
279                 {
280                         alpha = D[i][i];
281                         if( alpha == 0.0 ) // error - singular matrix
282                         {
283                                 error = 1;
284                                 break;
285                         }
286                         else
287                         {
288                                 for( j = 1; j <= n2; j++ )
289                                 {
290                                         D[i][j] = D[i][j]/alpha;
291                                 }
292                                 for( k = 1; k <= n; k++ )
293                                 {
294                                         if( (k-i) != 0 )
295                                         {
296                                                 beta = D[k][i];
297                                                 for( j = 1; j <= n2; j++ )
298                                                 {
299                                                         D[k][j] = D[k][j] - beta*D[i][j];
300                                                 }
301                                         }
302                                 }
303                         }
304                 }
305                 return D;
306         }
307  */
308 }