2 * Jalview - A Sequence Alignment Editor and Viewer ($$Version-Rel$$)
3 * Copyright (C) $$Year-Rel$$ The Jalview Authors
5 * This file is part of Jalview.
7 * Jalview is free software: you can redistribute it and/or
8 * modify it under the terms of the GNU General Public License
9 * as published by the Free Software Foundation, either version 3
10 * of the License, or (at your option) any later version.
12 * Jalview is distributed in the hope that it will be useful, but
13 * WITHOUT ANY WARRANTY; without even the implied warranty
14 * of MERCHANTABILITY or FITNESS FOR A PARTICULAR
15 * PURPOSE. See the GNU General Public License for more details.
17 * You should have received a copy of the GNU General Public License
18 * along with Jalview. If not, see <http://www.gnu.org/licenses/>.
19 * The Jalview Authors are detailed in the 'AUTHORS' file.
23 import jalview.util.Format;
24 import jalview.util.MessageManager;
26 import java.io.PrintStream;
27 import java.util.Arrays;
30 * A class to model rectangular matrices of double values and operations on them
32 public class Matrix implements MatrixI
35 * maximum number of iterations for tqli
37 private static final int MAX_ITER = 45;
38 // fudge - add 15 iterations, just in case
43 final protected int rows;
46 * the number of columns
48 final protected int cols;
51 * the cell values in row-major order
53 private double[][] value;
55 protected double[] d; // Diagonal
57 protected double[] e; // off diagonal
60 * Constructor given number of rows and columns
65 protected Matrix(int rowCount, int colCount)
72 * Creates a new Matrix object containing a copy of the supplied array values.
76 * new Matrix(new double[][] {{2, 3, 4}, {5, 6, 7})
82 * Note that ragged arrays (with not all rows, or columns, of the same
83 * length), are not supported by this class. They can be constructed, but
84 * results of operations on them are undefined and may throw exceptions.
87 * the matrix values in row-major order
89 public Matrix(double[][] values)
91 this.rows = values.length;
92 this.cols = this.rows == 0 ? 0 : values[0].length;
95 * make a copy of the values array, for immutability
97 this.value = new double[rows][];
99 for (double[] row : values)
103 value[i] = new double[row.length];
104 System.arraycopy(row, 0, value[i], 0, row.length);
111 public MatrixI transpose()
113 double[][] out = new double[cols][rows];
115 for (int i = 0; i < cols; i++)
117 for (int j = 0; j < rows; j++)
119 out[i][j] = value[j][i];
123 return new Matrix(out);
134 public void print(PrintStream ps, String format)
136 for (int i = 0; i < rows; i++)
138 for (int j = 0; j < cols; j++)
140 Format.print(ps, format, getValue(i, j));
148 public MatrixI preMultiply(MatrixI in)
150 if (in.width() != rows)
152 throw new IllegalArgumentException("Can't pre-multiply " + this.rows
153 + " rows by " + in.width() + " columns");
155 double[][] tmp = new double[in.height()][this.cols];
157 for (int i = 0; i < in.height(); i++)
159 for (int j = 0; j < this.cols; j++)
162 * result[i][j] is the vector product of
163 * in.row[i] and this.column[j]
165 for (int k = 0; k < in.width(); k++)
167 tmp[i][j] += (in.getValue(i, k) * this.value[k][j]);
172 return new Matrix(tmp);
181 public double[] vectorPostMultiply(double[] in)
183 double[] out = new double[in.length];
185 for (int i = 0; i < in.length; i++)
189 for (int k = 0; k < in.length; k++)
191 out[i] += (value[i][k] * in[k]);
199 public MatrixI postMultiply(MatrixI in)
201 if (in.height() != this.cols)
203 throw new IllegalArgumentException("Can't post-multiply " + this.cols
204 + " columns by " + in.height() + " rows");
206 return in.preMultiply(this);
210 public MatrixI copy()
212 double[][] newmat = new double[rows][cols];
214 for (int i = 0; i < rows; i++)
216 System.arraycopy(value[i], 0, newmat[i], 0, value[i].length);
219 Matrix m = new Matrix(newmat);
222 m.d = Arrays.copyOf(this.d, this.d.length);
226 m.e = Arrays.copyOf(this.e, this.e.length);
249 this.d = new double[rows];
250 this.e = new double[rows];
252 for (i = n; i >= 2; i--)
260 for (k = 1; k <= l; k++)
262 double v = Math.abs(getValue(i - 1, k - 1));
268 e[i - 1] = getValue(i - 1, l - 1);
272 for (k = 1; k <= l; k++)
274 double v = divideValue(i - 1, k - 1, scale);
278 f = getValue(i - 1, l - 1);
282 g = -1.0 * Math.sqrt(h);
289 e[i - 1] = scale * g;
291 setValue(i - 1, l - 1, f - g);
294 for (j = 1; j <= l; j++)
296 double val = getValue(i - 1, j - 1) / h;
297 setValue(j - 1, i - 1, val);
300 for (k = 1; k <= j; k++)
302 g += (getValue(j - 1, k - 1) * getValue(i - 1, k - 1));
305 for (k = j + 1; k <= l; k++)
307 g += (getValue(k - 1, j - 1) * getValue(i - 1, k - 1));
311 f += (e[j - 1] * getValue(i - 1, j - 1));
316 for (j = 1; j <= l; j++)
318 f = getValue(i - 1, j - 1);
319 g = e[j - 1] - (hh * f);
322 for (k = 1; k <= j; k++)
324 double val = (f * e[k - 1]) + (g * getValue(i - 1, k - 1));
325 addValue(j - 1, k - 1, -val);
332 e[i - 1] = getValue(i - 1, l - 1);
341 for (i = 1; i <= n; i++)
347 for (j = 1; j <= l; j++)
351 for (k = 1; k <= l; k++)
353 g += (getValue(i - 1, k - 1) * getValue(k - 1, j - 1));
356 for (k = 1; k <= l; k++)
358 addValue(k - 1, j - 1, -(g * getValue(k - 1, i - 1)));
363 d[i - 1] = getValue(i - 1, i - 1);
364 setValue(i - 1, i - 1, 1.0);
366 for (j = 1; j <= l; j++)
368 setValue(j - 1, i - 1, 0.0);
369 setValue(i - 1, j - 1, 0.0);
375 * Adds f to the value at [i, j] and returns the new value
381 protected double addValue(int i, int j, double f)
383 double v = value[i][j] + f;
389 * Divides the value at [i, j] by divisor and returns the new value. If d is
390 * zero, returns the unchanged value.
397 protected double divideValue(int i, int j, double divisor)
401 return getValue(i, j);
403 double v = value[i][j];
413 public void tqli() throws Exception
432 for (i = 2; i <= n; i++)
439 for (l = 1; l <= n; l++)
445 for (m = l; m <= (n - 1); m++)
447 dd = Math.abs(d[m - 1]) + Math.abs(d[m]);
449 if ((Math.abs(e[m - 1]) + dd) == dd)
459 if (iter == MAX_ITER)
461 throw new Exception(MessageManager.formatMessage(
462 "exception.matrix_too_many_iteration", new String[]
463 { "tqli", Integer.valueOf(MAX_ITER).toString() }));
467 // jalview.bin.Console.outPrintln("Iteration " + iter);
470 g = (d[l] - d[l - 1]) / (2.0 * e[l - 1]);
471 r = Math.sqrt((g * g) + 1.0);
472 g = d[m - 1] - d[l - 1] + (e[l - 1] / (g + sign(r, g)));
477 for (i = m - 1; i >= l; i--)
482 if (Math.abs(f) >= Math.abs(g))
485 r = Math.sqrt((c * c) + 1.0);
493 r = Math.sqrt((s * s) + 1.0);
500 r = ((d[i - 1] - g) * s) + (2.0 * c * b);
505 for (k = 1; k <= n; k++)
507 f = getValue(k - 1, i);
508 setValue(k - 1, i, (s * getValue(k - 1, i - 1)) + (c * f));
509 setValue(k - 1, i - 1,
510 (c * getValue(k - 1, i - 1)) - (s * f));
514 d[l - 1] = d[l - 1] - p;
523 public double getValue(int i, int j)
529 public void setValue(int i, int j, double val)
551 this.d = new double[rows];
552 this.e = new double[rows];
554 for (i = n - 1; i >= 1; i--)
562 for (k = 0; k < l; k++)
564 scale += Math.abs(value[i][k]);
573 for (k = 0; k < l; k++)
575 value[i][k] /= scale;
576 h += (value[i][k] * value[i][k]);
583 g = -1.0 * Math.sqrt(h);
595 for (j = 0; j < l; j++)
597 value[j][i] = value[i][j] / h;
600 for (k = 0; k < j; k++)
602 g += (value[j][k] * value[i][k]);
605 for (k = j; k < l; k++)
607 g += (value[k][j] * value[i][k]);
611 f += (e[j] * value[i][j]);
616 for (j = 0; j < l; j++)
622 for (k = 0; k < j; k++)
624 value[j][k] -= ((f * e[k]) + (g * value[i][k]));
640 for (i = 0; i < n; i++)
646 for (j = 0; j < l; j++)
650 for (k = 0; k < l; k++)
652 g += (value[i][k] * value[k][j]);
655 for (k = 0; k < l; k++)
657 value[k][j] -= (g * value[k][i]);
665 for (j = 0; j < l; j++)
676 public void tqli2() throws Exception
696 for (i = 2; i <= n; i++)
703 for (l = 1; l <= n; l++)
709 for (m = l; m <= (n - 1); m++)
711 dd = Math.abs(d[m - 1]) + Math.abs(d[m]);
713 if ((Math.abs(e[m - 1]) + dd) == dd)
723 if (iter == MAX_ITER)
725 throw new Exception(MessageManager.formatMessage(
726 "exception.matrix_too_many_iteration", new String[]
727 { "tqli2", Integer.valueOf(MAX_ITER).toString() }));
731 // jalview.bin.Console.outPrintln("Iteration " + iter);
734 g = (d[l] - d[l - 1]) / (2.0 * e[l - 1]);
735 r = Math.sqrt((g * g) + 1.0);
736 g = d[m - 1] - d[l - 1] + (e[l - 1] / (g + sign(r, g)));
741 for (i = m - 1; i >= l; i--)
746 if (Math.abs(f) >= Math.abs(g))
749 r = Math.sqrt((c * c) + 1.0);
757 r = Math.sqrt((s * s) + 1.0);
764 r = ((d[i - 1] - g) * s) + (2.0 * c * b);
769 for (k = 1; k <= n; k++)
772 value[k - 1][i] = (s * value[k - 1][i - 1]) + (c * f);
773 value[k - 1][i - 1] = (c * value[k - 1][i - 1]) - (s * f);
777 d[l - 1] = d[l - 1] - p;
786 * Answers the first argument with the sign of the second argument
793 static double sign(double a, double b)
806 * Returns an array containing the values in the specified column
812 public double[] getColumn(int col)
814 double[] out = new double[rows];
816 for (int i = 0; i < rows; i++)
818 out[i] = value[i][col];
832 public void printD(PrintStream ps, String format)
834 for (int j = 0; j < rows; j++)
836 Format.print(ps, format, d[j]);
849 public void printE(PrintStream ps, String format)
851 for (int j = 0; j < rows; j++)
853 Format.print(ps, format, e[j]);
858 public double[] getD()
864 public double[] getE()
882 public double[] getRow(int i)
884 double[] row = new double[cols];
885 System.arraycopy(value[i], 0, row, 0, cols);
890 * Returns a length 2 array of {minValue, maxValue} of all values in the
891 * matrix. Returns null if the matrix is null or empty.
895 double[] findMinMax()
901 double min = Double.MAX_VALUE;
902 double max = -Double.MAX_VALUE;
903 boolean empty = true;
904 for (double[] row : value)
922 return empty ? null : new double[] { min, max };
929 public void reverseRange(boolean maxToZero)
935 double[] minMax = findMinMax();
938 return; // empty matrix
940 double subtractFrom = maxToZero ? minMax[1] : minMax[0] + minMax[1];
942 for (double[] row : value)
949 row[j] = subtractFrom - x;
957 * Multiplies every entry in the matrix by the given value.
962 public void multiply(double by)
964 for (double[] row : value)
968 for (int i = 0; i < row.length; i++)
977 public void setD(double[] v)
983 public void setE(double[] v)
988 public double getTotal()
991 for (int i = 0; i < this.height(); i++)
993 for (int j = 0; j < this.width(); j++)
1005 public boolean equals(MatrixI m2, double delta)
1007 if (m2 == null || this.height() != m2.height()
1008 || this.width() != m2.width())
1012 for (int i = 0; i < this.height(); i++)
1014 for (int j = 0; j < this.width(); j++)
1016 double diff = this.getValue(i, j) - m2.getValue(i, j);
1017 if (Math.abs(diff) > delta)