Merge branch 'features/r2_11_2/JAL-3821_reinstate_patch' into develop
[jalview.git] / unused / javajs / util / Matrix.java
1 package javajs.util;
2
3 /**
4  * 
5  * streamlined and refined for Jmol by Bob Hanson
6  * 
7  * from http://math.nist.gov/javanumerics/jama/
8  * 
9  * Jama = Java Matrix class.
10  * 
11  * @author The MathWorks, Inc. and the National Institute of Standards and
12  *         Technology.
13  * @version 5 August 1998
14  */
15
16 public class Matrix implements Cloneable {
17
18   public double[][] a;
19   protected int m, n;
20
21   /**
22    * Construct a matrix quickly without checking arguments.
23    * 
24    * @param a
25    *        Two-dimensional array of doubles or null
26    * @param m
27    *        Number of rows.
28    * @param n
29    *        Number of colums.
30    */
31
32   public Matrix(double[][] a, int m, int n) {
33     this.a = (a == null ? new double[m][n] : a);
34     this.m = m;
35     this.n = n;
36   }
37
38   /**
39    * Get row dimension.
40    * 
41    * @return m, the number of rows.
42    */
43
44   public int getRowDimension() {
45     return m;
46   }
47
48   /**
49    * Get column dimension.
50    * 
51    * @return n, the number of columns.
52    */
53
54   public int getColumnDimension() {
55     return n;
56   }
57
58   /**
59    * Access the internal two-dimensional array.
60    * 
61    * @return Pointer to the two-dimensional array of matrix elements.
62    */
63
64   public double[][] getArray() {
65     return a;
66   }
67
68   /**
69    * Copy the internal two-dimensional array.
70    * 
71    * @return Two-dimensional array copy of matrix elements.
72    */
73
74   public double[][] getArrayCopy() {
75     double[][] x = new double[m][n];
76     for (int i = m; --i >= 0;)
77       for (int j = n; --j >= 0;)
78         x[i][j] = a[i][j];
79     return x;
80   }
81
82   /**
83    * Make a deep copy of a matrix
84    * 
85    * @return copy
86    */
87
88   public Matrix copy() {
89     Matrix x = new Matrix(null, m, n);
90     double[][] c = x.a;
91     for (int i = m; --i >= 0;)
92       for (int j = n; --j >= 0;)
93         c[i][j] = a[i][j];
94     return x;
95   }
96
97   /**
98    * Clone the Matrix object.
99    */
100
101   @Override
102   public Object clone() {
103     return copy();
104   }
105
106   /**
107    * Get a submatrix.
108    * 
109    * @param i0
110    *        Initial row index
111    * @param j0
112    *        Initial column index
113    * @param nrows
114    *        Number of rows
115    * @param ncols
116    *        Number of columns
117    * @return submatrix
118    * 
119    */
120
121   public Matrix getSubmatrix(int i0, int j0, int nrows, int ncols) {
122     Matrix x = new Matrix(null, nrows, ncols);
123     double[][] xa = x.a;
124     for (int i = nrows; --i >= 0;)
125       for (int j = ncols; --j >= 0;)
126         xa[i][j] = a[i0 + i][j0 + j];
127     return x;
128   }
129
130   /**
131    * Get a submatrix for a give number of columns and selected row set.
132    * 
133    * @param r
134    *        Array of row indices.
135    * @param n
136    *        number of rows 
137    * @return submatrix
138    */
139
140   public Matrix getMatrixSelected(int[] r, int n) {
141     Matrix x = new Matrix(null, r.length, n);
142     double[][] xa = x.a;
143     for (int i = r.length; --i >= 0;) {
144       double[] b = a[r[i]];
145       for (int j = n; --j >= 0;)
146         xa[i][j] = b[j];
147     }
148     return x;
149   }
150
151   /**
152    * Matrix transpose.
153    * 
154    * @return A'
155    */
156
157   public Matrix transpose() {
158     Matrix x = new Matrix(null, n, m);
159     double[][] c = x.a;
160     for (int i = m; --i >= 0;)
161       for (int j = n; --j >= 0;)
162         c[j][i] = a[i][j];
163     return x;
164   }
165
166   /**
167    * add two matrices
168    * @param b
169    * @return new Matrix this + b
170    */
171   public Matrix add(Matrix b) {
172     return scaleAdd(b, 1);
173   }
174
175   /**
176    * subtract two matrices
177    * @param b
178    * @return new Matrix this - b
179    */
180   public Matrix sub(Matrix b) {
181     return scaleAdd(b, -1);
182   }
183   
184   /**
185    * X = A + B*scale
186    * @param b 
187    * @param scale 
188    * @return X
189    * 
190    */
191   public Matrix scaleAdd(Matrix b, double scale) {
192     Matrix x = new Matrix(null, m, n);
193     double[][] xa = x.a;
194     double[][] ba = b.a;
195     for (int i = m; --i >= 0;)
196       for (int j = n; --j >= 0;)
197         xa[i][j] = ba[i][j] * scale + a[i][j];
198     return x;
199   }
200
201   /**
202    * Linear algebraic matrix multiplication, A * B
203    * 
204    * @param b
205    *        another matrix
206    * @return Matrix product, A * B or null for wrong dimension
207    */
208
209   public Matrix mul(Matrix b) {
210     if (b.m != n)
211       return null;
212     Matrix x = new Matrix(null, m, b.n);
213     double[][] xa = x.a;
214     double[][] ba = b.a;
215     for (int j = b.n; --j >= 0;)
216       for (int i = m; --i >= 0;) {
217         double[] arowi = a[i];
218         double s = 0;
219         for (int k = n; --k >= 0;)
220           s += arowi[k] * ba[k][j];
221         xa[i][j] = s;
222       }
223     return x;
224   }
225
226   /**
227    * Matrix inverse or pseudoinverse
228    * 
229    * @return inverse (m == n) or pseudoinverse (m != n)
230    */
231
232   public Matrix inverse() {
233     return new LUDecomp(m, n).solve(identity(m, m), n);
234   }
235
236   /**
237    * Matrix trace.
238    * 
239    * @return sum of the diagonal elements.
240    */
241
242   public double trace() {
243     double t = 0;
244     for (int i = Math.min(m, n); --i >= 0;)
245       t += a[i][i];
246     return t;
247   }
248
249   /**
250    * Generate identity matrix
251    * 
252    * @param m
253    *        Number of rows.
254    * @param n
255    *        Number of columns.
256    * @return An m-by-n matrix with ones on the diagonal and zeros elsewhere.
257    */
258
259   public static Matrix identity(int m, int n) {
260     Matrix x = new Matrix(null, m, n);
261     double[][] xa = x.a;
262     for (int i = Math.min(m, n); --i >= 0;)
263       xa[i][i] = 1;
264     return x;
265   }
266
267   /**
268    * similarly to M3/M4 standard rotation/translation matrix
269    * we set a rotationTranslation matrix to be:
270    * 
271    * [   nxn rot    nx1 trans
272    * 
273    *     1xn  0     1x1 1      ]
274    * 
275    * 
276    * @return rotation matrix
277    */
278   public Matrix getRotation() {
279     return getSubmatrix(0, 0, m - 1, n - 1);
280   }
281
282   public Matrix getTranslation() {
283     return getSubmatrix(0, n - 1, m - 1, 1);
284   }
285
286   public static Matrix newT(T3 r, boolean asColumn) {
287     return (asColumn ? new Matrix(new double[][] { new double[] { r.x },
288         new double[] { r.y }, new double[] { r.z } }, 3, 1) : new Matrix(
289         new double[][] { new double[] { r.x, r.y, r.z } }, 1, 3));
290   }
291
292   @Override
293   public String toString() {
294     String s = "[\n";
295     for (int i = 0; i < m; i++) {
296       s += "  [";
297       for (int j = 0; j < n; j++)
298         s += " " + a[i][j];
299       s += "]\n";
300     }
301     s += "]";
302     return s;
303   }
304
305   /**
306    * 
307    * Edited down by Bob Hanson for minimum needed by Jmol -- just constructor
308    * and solve
309    * 
310    * LU Decomposition.
311    * <P>
312    * For an m-by-n matrix A with m >= n, the LU decomposition is an m-by-n unit
313    * lower triangular matrix L, an n-by-n upper triangular matrix U, and a
314    * permutation vector piv of length m so that A(piv,:) = L*U. If m < n, then L
315    * is m-by-m and U is m-by-n.
316    * <P>
317    * The LU decompostion with pivoting always exists, even if the matrix is
318    * singular, so the constructor will never fail. The primary use of the LU
319    * decomposition is in the solution of square systems of simultaneous linear
320    * equations. This will fail if isNonsingular() returns false.
321    */
322
323   private class LUDecomp {
324
325     /* ------------------------
326        Class variables
327      * ------------------------ */
328
329     /**
330      * Array for internal storage of decomposition.
331      * 
332      */
333     private double[][] LU;
334
335     /**
336      * Internal storage of pivot vector.
337      * 
338      */
339     private int[] piv;
340
341     private int pivsign;
342
343     /* ------------------------
344        Constructor
345      * ------------------------ */
346
347     /**
348      * LU Decomposition Structure to access L, U and piv.
349      * @param m 
350      * @param n 
351      * 
352      */
353
354     protected LUDecomp(int m, int n) {
355       
356       // Use a "left-looking", dot-product, Crout/Doolittle algorithm.
357
358       LU = getArrayCopy();
359       piv = new int[m];
360       for (int i = m; --i >= 0;)
361         piv[i] = i;
362       pivsign = 1;
363       double[] LUrowi;
364       double[] LUcolj = new double[m];
365
366       // Outer loop.
367
368       for (int j = 0; j < n; j++) {
369
370         // Make a copy of the j-th column to localize references.
371
372         for (int i = m; --i >= 0;)
373           LUcolj[i] = LU[i][j];
374
375         // Apply previous transformations.
376
377         for (int i = m; --i >= 0;) {
378           LUrowi = LU[i];
379
380           // Most of the time is spent in the following dot product.
381
382           int kmax = Math.min(i, j);
383           double s = 0.0;
384           for (int k = kmax; --k >= 0;)
385             s += LUrowi[k] * LUcolj[k];
386
387           LUrowi[j] = LUcolj[i] -= s;
388         }
389
390         // Find pivot and exchange if necessary.
391
392         int p = j;
393         for (int i = m; --i > j;)
394           if (Math.abs(LUcolj[i]) > Math.abs(LUcolj[p]))
395             p = i;
396         if (p != j) {
397           for (int k = n; --k >= 0;) {
398             double t = LU[p][k];
399             LU[p][k] = LU[j][k];
400             LU[j][k] = t;
401           }
402           int k = piv[p];
403           piv[p] = piv[j];
404           piv[j] = k;
405           pivsign = -pivsign;
406         }
407
408         // Compute multipliers.
409
410         if (j < m & LU[j][j] != 0.0)
411           for (int i = m; --i > j;)
412             LU[i][j] /= LU[j][j];
413       }
414     }
415
416     /* ------------------------
417        default Methods
418      * ------------------------ */
419
420     /**
421      * Solve A*X = B
422      * 
423      * @param b
424      *        A Matrix with as many rows as A and any number of columns.
425      * @param n 
426      * @return X so that L*U*X = B(piv,:) or null for wrong size or singular matrix
427      */
428
429     protected Matrix solve(Matrix b, int n) {
430       for (int j = 0; j < n; j++)
431         if (LU[j][j] == 0)
432           return null; // matrix is singular
433
434       // Copy right hand side with pivoting
435       int nx = b.n;
436       Matrix x = b.getMatrixSelected(piv, nx);
437       double[][] a = x.a;
438
439       // Solve L*Y = B(piv,:)
440       for (int k = 0; k < n; k++)
441         for (int i = k + 1; i < n; i++)
442           for (int j = 0; j < nx; j++)
443             a[i][j] -= a[k][j] * LU[i][k];
444
445       // Solve U*X = Y;
446       for (int k = n; --k >= 0;) {
447         for (int j = nx; --j >= 0;)
448           a[k][j] /= LU[k][k];
449         for (int i = k; --i >= 0;)
450           for (int j = nx; --j >= 0;)
451             a[i][j] -= a[k][j] * LU[i][k];
452       }
453       return x;
454     }
455   }
456
457 }