From: kiramt Date: Tue, 7 Mar 2017 10:01:14 +0000 (+0000) Subject: Merge remote-tracking branch 'origin/develop' into features/JAL-2388OverviewWindow X-Git-Tag: Release_2_10_2~3^2~92^2~49 X-Git-Url: http://source.jalview.org/gitweb/?a=commitdiff_plain;h=14a3cc9277d30d14bec559a12d9af72bd2955aff;hp=366b1c4a421d8c8455668b06e24620c2101d8a4b;p=jalview.git Merge remote-tracking branch 'origin/develop' into features/JAL-2388OverviewWindow --- diff --git a/src/jalview/analysis/PCA.java b/src/jalview/analysis/PCA.java index eaea7bf..9babaee 100755 --- a/src/jalview/analysis/PCA.java +++ b/src/jalview/analysis/PCA.java @@ -20,9 +20,7 @@ */ package jalview.analysis; -import jalview.datamodel.BinarySequence; -import jalview.datamodel.BinarySequence.InvalidSequenceTypeException; -import jalview.math.Matrix; +import jalview.math.MatrixI; import jalview.schemes.ResidueProperties; import jalview.schemes.ScoreMatrix; @@ -30,23 +28,22 @@ import java.io.PrintStream; /** * Performs Principal Component Analysis on given sequences - * - * @author $author$ - * @version $Revision$ */ public class PCA implements Runnable { - Matrix m; - - Matrix symm; + boolean jvCalcMode = true; - Matrix m2; + MatrixI symm; double[] eigenvalue; - Matrix eigenvector; + MatrixI eigenvector; - StringBuffer details = new StringBuffer(); + StringBuilder details = new StringBuilder(1024); + + private String[] seqs; + + private ScoreMatrix scoreMatrix; /** * Creates a new PCA object. By default, uses blosum62 matrix to generate @@ -77,87 +74,23 @@ public class PCA implements Runnable public PCA(String[] s, boolean nucleotides, String s_m) { + this.seqs = s; - BinarySequence[] bs = new BinarySequence[s.length]; - int ii = 0; - - while ((ii < s.length) && (s[ii] != null)) - { - bs[ii] = new BinarySequence(s[ii], nucleotides); - bs[ii].encode(); - ii++; - } - - BinarySequence[] bs2 = new BinarySequence[s.length]; - ii = 0; - ScoreMatrix smtrx = null; + scoreMatrix = null; String sm = s_m; if (sm != null) { - smtrx = ResidueProperties.getScoreMatrix(sm); + scoreMatrix = ResidueProperties.getScoreMatrix(sm); } - if (smtrx == null) + if (scoreMatrix == null) { // either we were given a non-existent score matrix or a scoremodel that // isn't based on a pairwise symbol score matrix - smtrx = ResidueProperties.getScoreMatrix(sm = (nucleotides ? "DNA" - : "BLOSUM62")); + scoreMatrix = ResidueProperties + .getScoreMatrix(sm = (nucleotides ? "DNA" : "BLOSUM62")); } details.append("PCA calculation using " + sm + " sequence similarity matrix\n========\n\n"); - while ((ii < s.length) && (s[ii] != null)) - { - bs2[ii] = new BinarySequence(s[ii], nucleotides); - if (smtrx != null) - { - try - { - bs2[ii].matrixEncode(smtrx); - } catch (InvalidSequenceTypeException x) - { - details.append("Unexpected mismatch of sequence type and score matrix. Calculation will not be valid!\n\n"); - } - } - ii++; - } - - // System.out.println("Created binary encoding"); - // printMemory(rt); - int count = 0; - - while ((count < bs.length) && (bs[count] != null)) - { - count++; - } - - double[][] seqmat = new double[count][bs[0].getDBinary().length]; - double[][] seqmat2 = new double[count][bs2[0].getDBinary().length]; - int i = 0; - - while (i < count) - { - seqmat[i] = bs[i].getDBinary(); - seqmat2[i] = bs2[i].getDBinary(); - i++; - } - - // System.out.println("Created array"); - // printMemory(rt); - // System.out.println(" --- Original matrix ---- "); - m = new Matrix(seqmat); - m2 = new Matrix(seqmat2); - - } - - /** - * Returns the matrix used in PCA calculation - * - * @return java.math.Matrix object - */ - - public Matrix getM() - { - return m; } /** @@ -170,7 +103,7 @@ public class PCA implements Runnable */ public double getEigenvalue(int i) { - return eigenvector.d[i]; + return eigenvector.getD()[i]; } /** @@ -189,9 +122,9 @@ public class PCA implements Runnable */ public float[][] getComponents(int l, int n, int mm, float factor) { - float[][] out = new float[m.rows][3]; + float[][] out = new float[getHeight()][3]; - for (int i = 0; i < m.rows; i++) + for (int i = 0; i < getHeight(); i++) { out[i][0] = (float) component(i, l) * factor; out[i][1] = (float) component(i, n) * factor; @@ -212,9 +145,9 @@ public class PCA implements Runnable public double[] component(int n) { // n = index of eigenvector - double[] out = new double[m.rows]; + double[] out = new double[getHeight()]; - for (int i = 0; i < m.rows; i++) + for (int i = 0; i < out.length; i++) { out[i] = component(i, n); } @@ -236,12 +169,12 @@ public class PCA implements Runnable { double out = 0.0; - for (int i = 0; i < symm.cols; i++) + for (int i = 0; i < symm.width(); i++) { - out += (symm.value[row][i] * eigenvector.value[i][n]); + out += (symm.getValue(row, i) * eigenvector.getValue(i, n)); } - return out / eigenvector.d[n]; + return out / eigenvector.getD()[n]; } public String getDetails() @@ -270,25 +203,17 @@ public class PCA implements Runnable } }; + // long now = System.currentTimeMillis(); try { details.append("PCA Calculation Mode is " + (jvCalcMode ? "Jalview variant" : "Original SeqSpace") + "\n"); - Matrix mt = m.transpose(); - details.append(" --- OrigT * Orig ---- \n"); - if (!jvCalcMode) - { - eigenvector = mt.preMultiply(m); // standard seqspace comparison matrix - } - else - { - eigenvector = mt.preMultiply(m2); // jalview variation on seqsmace - // method - } + eigenvector = scoreMatrix.computePairwiseScores(seqs); - eigenvector.print(ps); + details.append(" --- OrigT * Orig ---- \n"); + eigenvector.print(ps, "%8.2f"); symm = eigenvector.copy(); @@ -296,10 +221,10 @@ public class PCA implements Runnable details.append(" ---Tridiag transform matrix ---\n"); details.append(" --- D vector ---\n"); - eigenvector.printD(ps); + eigenvector.printD(ps, "%15.4e"); ps.println(); details.append("--- E vector ---\n"); - eigenvector.printE(ps); + eigenvector.printE(ps, "%15.4e"); ps.println(); // Now produce the diagonalization matrix @@ -313,9 +238,9 @@ public class PCA implements Runnable } details.append(" --- New diagonalization matrix ---\n"); - eigenvector.print(ps); + eigenvector.print(ps, "%8.2f"); details.append(" --- Eigenvalues ---\n"); - eigenvector.printD(ps); + eigenvector.printD(ps, "%15.4e"); ps.println(); /* * for (int seq=0;seq + * Note that this container keeps its mappings in an array data structure, using + * a binary search to find keys. The implementation is not intended to be + * appropriate for data structures that may contain large numbers of items. It + * is generally slower than a traditional HashMap, since lookups require a + * binary search and adds and removes require inserting and deleting entries in + * the array. For containers holding up to hundreds of items, the performance + * difference is not significant, less than 50%. + *

+ * + *

+ * It is possible to iterate over the items in this container using + * {@link #keyAt(int)} and {@link #valueAt(int)}. Iterating over the keys using + * keyAt(int) with ascending values of the index will return the + * keys in ascending order, or the values corresponding to the keys in ascending + * order in the case of valueAt(int). + *

+ */ + +/* + * Change log: + * Jan 2017 cloned from SparseIntArray for Jalview to support SparseMatrix + * - SparseDoubleArray(double[]) constructor added + * - multiply() added for more efficient multiply (or divide) of a value + */ +public class SparseDoubleArray implements Cloneable +{ + private int[] mKeys; + + private double[] mValues; + + private int mSize; + + /** + * Creates a new SparseDoubleArray containing no mappings. + */ + public SparseDoubleArray() + { + this(10); + } + + /** + * Creates a new SparseDoubleArray containing no mappings that will not + * require any additional memory allocation to store the specified number of + * mappings. If you supply an initial capacity of 0, the sparse array will be + * initialized with a light-weight representation not requiring any additional + * array allocations. + */ + public SparseDoubleArray(int initialCapacity) + { + if (initialCapacity == 0) + { + mKeys = ContainerHelpers.EMPTY_INTS; + mValues = ContainerHelpers.EMPTY_DOUBLES; + } + else + { + initialCapacity = idealDoubleArraySize(initialCapacity); + mKeys = new int[initialCapacity]; + mValues = new double[initialCapacity]; + } + mSize = 0; + } + + /** + * Constructor given an array of double values; stores the non-zero values + * + * @param row + */ + public SparseDoubleArray(double[] row) + { + this(); + for (int i = 0; i < row.length; i++) + { + if (row[i] != 0d) + { + put(i, row[i]); + } + } + } + + @Override + public SparseDoubleArray clone() + { + SparseDoubleArray clone = null; + try + { + clone = (SparseDoubleArray) super.clone(); + clone.mKeys = mKeys.clone(); + clone.mValues = mValues.clone(); + } catch (CloneNotSupportedException cnse) + { + /* ignore */ + } + return clone; + } + + /** + * Gets the value mapped from the specified key, or 0 if no such + * mapping has been made. + */ + public double get(int key) + { + return get(key, 0d); + } + + /** + * Gets the int mapped from the specified key, or the specified value if no + * such mapping has been made. + */ + public double get(int key, double valueIfKeyNotFound) + { + int i = ContainerHelpers.binarySearch(mKeys, mSize, key); + if (i < 0) + { + return valueIfKeyNotFound; + } + else + { + return mValues[i]; + } + } + + /** + * Removes the mapping from the specified key, if there was any. + */ + public void delete(int key) + { + int i = ContainerHelpers.binarySearch(mKeys, mSize, key); + if (i >= 0) + { + removeAt(i); + } + } + + /** + * Removes the mapping at the given index. + */ + public void removeAt(int index) + { + System.arraycopy(mKeys, index + 1, mKeys, index, mSize - (index + 1)); + System.arraycopy(mValues, index + 1, mValues, index, mSize + - (index + 1)); + mSize--; + } + + /** + * Adds a mapping from the specified key to the specified value, replacing the + * previous mapping from the specified key if there was one. + */ + public void put(int key, double value) + { + int i = ContainerHelpers.binarySearch(mKeys, mSize, key); + if (i >= 0) + { + mValues[i] = value; + } + else + { + i = ~i; + if (mSize >= mKeys.length) + { + int n = idealDoubleArraySize(mSize + 1); + int[] nkeys = new int[n]; + double[] nvalues = new double[n]; + // Log.e("SparseDoubleArray", "grow " + mKeys.length + " to " + n); + System.arraycopy(mKeys, 0, nkeys, 0, mKeys.length); + System.arraycopy(mValues, 0, nvalues, 0, mValues.length); + mKeys = nkeys; + mValues = nvalues; + } + if (mSize - i != 0) + { + // Log.e("SparseDoubleArray", "move " + (mSize - i)); + System.arraycopy(mKeys, i, mKeys, i + 1, mSize - i); + System.arraycopy(mValues, i, mValues, i + 1, mSize - i); + } + mKeys[i] = key; + mValues[i] = value; + mSize++; + } + } + + /** + * Returns the number of key-value mappings that this SparseDoubleArray + * currently stores. + */ + public int size() + { + return mSize; + } + + /** + * Given an index in the range 0...size()-1, returns the key from + * the indexth key-value mapping that this SparseDoubleArray + * stores. + * + *

+ * The keys corresponding to indices in ascending order are guaranteed to be + * in ascending order, e.g., keyAt(0) will return the smallest + * key and keyAt(size()-1) will return the largest key. + *

+ */ + public int keyAt(int index) + { + return mKeys[index]; + } + + /** + * Given an index in the range 0...size()-1, returns the value + * from the indexth key-value mapping that this SparseDoubleArray + * stores. + * + *

+ * The values corresponding to indices in ascending order are guaranteed to be + * associated with keys in ascending order, e.g., valueAt(0) will + * return the value associated with the smallest key and + * valueAt(size()-1) will return the value associated with the + * largest key. + *

+ */ + public double valueAt(int index) + { + return mValues[index]; + } + + /** + * Returns the index for which {@link #keyAt} would return the specified key, + * or a negative number if the specified key is not mapped. + */ + public int indexOfKey(int key) + { + return ContainerHelpers.binarySearch(mKeys, mSize, key); + } + + /** + * Returns an index for which {@link #valueAt} would return the specified key, + * or a negative number if no keys map to the specified value. Beware that + * this is a linear search, unlike lookups by key, and that multiple keys can + * map to the same value and this will find only one of them. + */ + public int indexOfValue(double value) + { + for (int i = 0; i < mSize; i++) + { + if (mValues[i] == value) + { + return i; + } + } + return -1; + } + + /** + * Removes all key-value mappings from this SparseDoubleArray. + */ + public void clear() + { + mSize = 0; + } + + /** + * Puts a key/value pair into the array, optimizing for the case where the key + * is greater than all existing keys in the array. + */ + public void append(int key, double value) + { + if (mSize != 0 && key <= mKeys[mSize - 1]) + { + put(key, value); + return; + } + int pos = mSize; + if (pos >= mKeys.length) + { + int n = idealDoubleArraySize(pos + 1); + int[] nkeys = new int[n]; + double[] nvalues = new double[n]; + // Log.e("SparseDoubleArray", "grow " + mKeys.length + " to " + n); + System.arraycopy(mKeys, 0, nkeys, 0, mKeys.length); + System.arraycopy(mValues, 0, nvalues, 0, mValues.length); + mKeys = nkeys; + mValues = nvalues; + } + mKeys[pos] = key; + mValues[pos] = value; + mSize = pos + 1; + } + + /** + * Created by analogy with + * com.android.internal.util.ArrayUtils#idealLongArraySize + * + * @param i + * @return + */ + public static int idealDoubleArraySize(int need) + { + return idealByteArraySize(need * 8) / 8; + } + + /** + * Inlined here by copying from com.android.internal.util.ArrayUtils + * + * @param i + * @return + */ + public static int idealByteArraySize(int need) + { + for (int i = 4; i < 32; i++) + { + if (need <= (1 << i) - 12) + { + return (1 << i) - 12; + } + } + + return need; + } + + /** + * {@inheritDoc} + * + *

+ * This implementation composes a string by iterating over its mappings. + */ + @Override + public String toString() + { + if (size() <= 0) + { + return "{}"; + } + StringBuilder buffer = new StringBuilder(mSize * 28); + buffer.append('{'); + for (int i = 0; i < mSize; i++) + { + if (i > 0) + { + buffer.append(", "); + } + int key = keyAt(i); + buffer.append(key); + buffer.append('='); + double value = valueAt(i); + buffer.append(value); + } + buffer.append('}'); + return buffer.toString(); + } + + /** + * Method (copied from put) added for Jalview to efficiently increment a key's + * value if present, else add it with the given value. This avoids a double + * binary search (once to get the value, again to put the updated value). + * + * @param key + * @oparam toAdd + * @return the new value for the key + */ + public double add(int key, double toAdd) + { + double newValue = toAdd; + int i = ContainerHelpers.binarySearch(mKeys, mSize, key); + if (i >= 0) + { + mValues[i] += toAdd; + newValue = mValues[i]; + } + else + { + i = ~i; + if (mSize >= mKeys.length) + { + int n = idealDoubleArraySize(mSize + 1); + int[] nkeys = new int[n]; + double[] nvalues = new double[n]; + System.arraycopy(mKeys, 0, nkeys, 0, mKeys.length); + System.arraycopy(mValues, 0, nvalues, 0, mValues.length); + mKeys = nkeys; + mValues = nvalues; + } + if (mSize - i != 0) + { + System.arraycopy(mKeys, i, mKeys, i + 1, mSize - i); + System.arraycopy(mValues, i, mValues, i + 1, mSize - i); + } + mKeys[i] = key; + mValues[i] = toAdd; + mSize++; + } + return newValue; + } + + /** + * Method added for Jalview to efficiently multiply a key's value if present, + * else do nothing. This avoids a double binary search (once to get the value, + * again to put the updated value). + * + * @param key + * @oparam toAdd + * @return the new value for the key + */ + public double divide(int key, double divisor) + { + double newValue = 0d; + if (divisor == 0d) + { + return newValue; + } + int i = ContainerHelpers.binarySearch(mKeys, mSize, key); + if (i >= 0) + { + mValues[i] /= divisor; + newValue = mValues[i]; + } + return newValue; + } +} diff --git a/src/jalview/math/Matrix.java b/src/jalview/math/Matrix.java index 647fc3a..de0bf77 100755 --- a/src/jalview/math/Matrix.java +++ b/src/jalview/math/Matrix.java @@ -28,36 +28,43 @@ import java.io.PrintStream; /** * A class to model rectangular matrices of double values and operations on them */ -public class Matrix +public class Matrix implements MatrixI { /* * the cell values in row-major order */ - public double[][] value; + private double[][] value; /* * the number of rows */ - public int rows; + protected int rows; /* * the number of columns */ - public int cols; + protected int cols; - /** DOCUMENT ME!! */ - public double[] d; // Diagonal + protected double[] d; // Diagonal - /** DOCUMENT ME!! */ - public double[] e; // off diagonal + protected double[] e; // off diagonal /** * maximum number of iterations for tqli * */ - int maxIter = 45; // fudge - add 15 iterations, just in case + private static final int maxIter = 45; // fudge - add 15 iterations, just in + // case /** + * Default constructor + */ + public Matrix() + { + + } + + /** * Creates a new Matrix object. For example * * 

@@ -82,11 +89,12 @@ public class Matrix
   }
 
   /**
-   * Returns a new matrix which is the transposes of this one
+   * Returns a new matrix which is the transpose of this one
    * 
    * @return DOCUMENT ME!
    */
-  public Matrix transpose()
+  @Override
+  public MatrixI transpose()
   {
     double[][] out = new double[cols][rows];
 
@@ -106,14 +114,16 @@ public class Matrix
    * 
    * @param ps
    *          DOCUMENT ME!
+   * @param format
    */
-  public void print(PrintStream ps)
+  @Override
+  public void print(PrintStream ps, String format)
   {
     for (int i = 0; i < rows; i++)
     {
       for (int j = 0; j < cols; j++)
       {
-        Format.print(ps, "%8.2f", value[i][j]);
+        Format.print(ps, format, getValue(i, j));
       }
 
       ps.println();
@@ -132,24 +142,27 @@ public class Matrix
    *           if the number of columns in the pre-multiplier is not equal to
    *           the number of rows in the multiplicand (this)
    */
-  public Matrix preMultiply(Matrix in)
+  @Override
+  public MatrixI preMultiply(MatrixI in)
   {
-    if (in.cols != this.rows)
+    if (in.width() != rows)
     {
       throw new IllegalArgumentException("Can't pre-multiply " + this.rows
-              + " rows by " + in.cols + " columns");
+              + " rows by " + in.width() + " columns");
     }
-    double[][] tmp = new double[in.rows][this.cols];
+    double[][] tmp = new double[in.height()][this.cols];
 
-    for (int i = 0; i < in.rows; i++)
+    for (int i = 0; i < in.height(); i++)
     {
       for (int j = 0; j < this.cols; j++)
       {
-        tmp[i][j] = 0.0;
-
-        for (int k = 0; k < in.cols; k++)
+        /*
+         * result[i][j] is the vector product of 
+         * in.row[i] and this.column[j]
+         */
+        for (int k = 0; k < in.width(); k++)
         {
-          tmp[i][j] += (in.value[i][k] * this.value[k][j]);
+          tmp[i][j] += (in.getValue(i, k) * this.value[k][j]);
         }
       }
     }
@@ -195,12 +208,13 @@ public class Matrix
    *           number of columns in the multiplicand (this)
    * @see #preMultiply(Matrix)
    */
-  public Matrix postMultiply(Matrix in)
+  @Override
+  public MatrixI postMultiply(MatrixI in)
   {
-    if (in.rows != this.cols)
+    if (in.height() != this.cols)
     {
       throw new IllegalArgumentException("Can't post-multiply " + this.cols
-              + " columns by " + in.rows + " rows");
+              + " columns by " + in.height() + " rows");
     }
     return in.preMultiply(this);
   }
@@ -210,7 +224,8 @@ public class Matrix
    * 
    * @return
    */
-  public Matrix copy()
+  @Override
+  public MatrixI copy()
   {
     double[][] newmat = new double[rows][cols];
 
@@ -225,10 +240,10 @@ public class Matrix
   /**
    * DOCUMENT ME!
    */
+  @Override
   public void tred()
   {
     int n = rows;
-    int l;
     int k;
     int j;
     int i;
@@ -244,7 +259,7 @@ public class Matrix
 
     for (i = n; i >= 2; i--)
     {
-      l = i - 1;
+      final int l = i - 1;
       h = 0.0;
       scale = 0.0;
 
@@ -252,22 +267,23 @@ public class Matrix
       {
         for (k = 1; k <= l; k++)
         {
-          scale += Math.abs(value[i - 1][k - 1]);
+          double v = Math.abs(getValue(i - 1, k - 1));
+          scale += v;
         }
 
         if (scale == 0.0)
         {
-          e[i - 1] = value[i - 1][l - 1];
+          e[i - 1] = getValue(i - 1, l - 1);
         }
         else
         {
           for (k = 1; k <= l; k++)
           {
-            value[i - 1][k - 1] /= scale;
-            h += (value[i - 1][k - 1] * value[i - 1][k - 1]);
+            double v = divideValue(i - 1, k - 1, scale);
+            h += v * v;
           }
 
-          f = value[i - 1][l - 1];
+          f = getValue(i - 1, l - 1);
 
           if (f > 0)
           {
@@ -280,46 +296,48 @@ public class Matrix
 
           e[i - 1] = scale * g;
           h -= (f * g);
-          value[i - 1][l - 1] = f - g;
+          setValue(i - 1, l - 1, f - g);
           f = 0.0;
 
           for (j = 1; j <= l; j++)
           {
-            value[j - 1][i - 1] = value[i - 1][j - 1] / h;
+            double val = getValue(i - 1, j - 1) / h;
+            setValue(j - 1, i - 1, val);
             g = 0.0;
 
             for (k = 1; k <= j; k++)
             {
-              g += (value[j - 1][k - 1] * value[i - 1][k - 1]);
+              g += (getValue(j - 1, k - 1) * getValue(i - 1, k - 1));
             }
 
             for (k = j + 1; k <= l; k++)
             {
-              g += (value[k - 1][j - 1] * value[i - 1][k - 1]);
+              g += (getValue(k - 1, j - 1) * getValue(i - 1, k - 1));
             }
 
             e[j - 1] = g / h;
-            f += (e[j - 1] * value[i - 1][j - 1]);
+            f += (e[j - 1] * getValue(i - 1, j - 1));
           }
 
           hh = f / (h + h);
 
           for (j = 1; j <= l; j++)
           {
-            f = value[i - 1][j - 1];
+            f = getValue(i - 1, j - 1);
             g = e[j - 1] - (hh * f);
             e[j - 1] = g;
 
             for (k = 1; k <= j; k++)
             {
-              value[j - 1][k - 1] -= ((f * e[k - 1]) + (g * value[i - 1][k - 1]));
+              double val = (f * e[k - 1]) + (g * getValue(i - 1, k - 1));
+              addValue(j - 1, k - 1, -val);
             }
           }
         }
       }
       else
       {
-        e[i - 1] = value[i - 1][l - 1];
+        e[i - 1] = getValue(i - 1, l - 1);
       }
 
       d[i - 1] = h;
@@ -330,7 +348,7 @@ public class Matrix
 
     for (i = 1; i <= n; i++)
     {
-      l = i - 1;
+      final int l = i - 1;
 
       if (d[i - 1] != 0.0)
       {
@@ -340,30 +358,66 @@ public class Matrix
 
           for (k = 1; k <= l; k++)
           {
-            g += (value[i - 1][k - 1] * value[k - 1][j - 1]);
+            g += (getValue(i - 1, k - 1) * getValue(k - 1, j - 1));
           }
 
           for (k = 1; k <= l; k++)
           {
-            value[k - 1][j - 1] -= (g * value[k - 1][i - 1]);
+            addValue(k - 1, j - 1, -(g * getValue(k - 1, i - 1)));
           }
         }
       }
 
-      d[i - 1] = value[i - 1][i - 1];
-      value[i - 1][i - 1] = 1.0;
+      d[i - 1] = getValue(i - 1, i - 1);
+      setValue(i - 1, i - 1, 1.0);
 
       for (j = 1; j <= l; j++)
       {
-        value[j - 1][i - 1] = 0.0;
-        value[i - 1][j - 1] = 0.0;
+        setValue(j - 1, i - 1, 0.0);
+        setValue(i - 1, j - 1, 0.0);
       }
     }
   }
 
   /**
+   * Adds f to the value at [i, j] and returns the new value
+   * 
+   * @param i
+   * @param j
+   * @param f
+   */
+  protected double addValue(int i, int j, double f)
+  {
+    double v = value[i][j] + f;
+    value[i][j] = v;
+    return v;
+  }
+
+  /**
+   * Divides the value at [i, j] by divisor and returns the new value. If d is
+   * zero, returns the unchanged value.
+   * 
+   * @param i
+   * @param j
+   * @param divisor
+   * @return
+   */
+  protected double divideValue(int i, int j, double divisor)
+  {
+    if (divisor == 0d)
+    {
+      return getValue(i, j);
+    }
+    double v = value[i][j];
+    v = v / divisor;
+    value[i][j] = v;
+    return v;
+  }
+
+  /**
    * DOCUMENT ME!
    */
+  @Override
   public void tqli() throws Exception
   {
     int n = rows;
@@ -376,7 +430,6 @@ public class Matrix
     double s;
     double r;
     double p;
-    ;
 
     double g;
     double f;
@@ -459,9 +512,9 @@ public class Matrix
 
             for (k = 1; k <= n; k++)
             {
-              f = value[k - 1][i];
-              value[k - 1][i] = (s * value[k - 1][i - 1]) + (c * f);
-              value[k - 1][i - 1] = (c * value[k - 1][i - 1]) - (s * f);
+              f = getValue(k - 1, i);
+              setValue(k - 1, i, (s * getValue(k - 1, i - 1)) + (c * f));
+              setValue(k - 1, i - 1, (c * getValue(k - 1, i - 1)) - (s * f));
             }
           }
 
@@ -473,6 +526,17 @@ public class Matrix
     }
   }
 
+  @Override
+  public double getValue(int i, int j)
+  {
+    return value[i][j];
+  }
+
+  public void setValue(int i, int j, double val)
+  {
+    value[i][j] = val;
+  }
+
   /**
    * DOCUMENT ME!
    */
@@ -725,16 +789,14 @@ public class Matrix
   }
 
   /**
-   * DOCUMENT ME!
+   * Answers the first argument with the sign of the second argument
    * 
    * @param a
-   *          DOCUMENT ME!
    * @param b
-   *          DOCUMENT ME!
    * 
-   * @return DOCUMENT ME!
+   * @return
    */
-  public double sign(double a, double b)
+  static double sign(double a, double b)
   {
     if (b < 0)
     {
@@ -770,12 +832,14 @@ public class Matrix
    * 
    * @param ps
    *          DOCUMENT ME!
+   * @param format
    */
-  public void printD(PrintStream ps)
+  @Override
+  public void printD(PrintStream ps, String format)
   {
     for (int j = 0; j < rows; j++)
     {
-      Format.print(ps, "%15.4e", d[j]);
+      Format.print(ps, format, d[j]);
     }
   }
 
@@ -784,12 +848,45 @@ public class Matrix
    * 
    * @param ps
    *          DOCUMENT ME!
+   * @param format TODO
    */
-  public void printE(PrintStream ps)
+  @Override
+  public void printE(PrintStream ps, String format)
   {
     for (int j = 0; j < rows; j++)
     {
-      Format.print(ps, "%15.4e", e[j]);
+      Format.print(ps, format, e[j]);
     }
   }
+
+  @Override
+  public double[] getD()
+  {
+    return d;
+  }
+
+  @Override
+  public double[] getE()
+  {
+    return e;
+  }
+  
+  @Override
+  public int height() {
+    return rows;
+  }
+
+  @Override
+  public int width()
+  {
+    return cols;
+  }
+
+  @Override
+  public double[] getRow(int i)
+  {
+    double[] row = new double[cols];
+    System.arraycopy(value[i], 0, row, 0, cols);
+    return row;
+  }
 }
diff --git a/src/jalview/math/MatrixI.java b/src/jalview/math/MatrixI.java
new file mode 100644
index 0000000..d74a98b
--- /dev/null
+++ b/src/jalview/math/MatrixI.java
@@ -0,0 +1,59 @@
+package jalview.math;
+
+import java.io.PrintStream;
+
+public interface MatrixI
+{
+  /**
+   * Answers the number of columns
+   * 
+   * @return
+   */
+  int width();
+
+  /**
+   * Answers the number of rows
+   * 
+   * @return
+   */
+  int height();
+
+  /**
+   * Answers the value at row i, column j
+   * 
+   * @param i
+   * @param j
+   * @return
+   */
+  double getValue(int i, int j);
+
+  /**
+   * Answers a copy of the values in the i'th row
+   * 
+   * @return
+   */
+  double[] getRow(int i);
+  
+  MatrixI copy();
+
+  MatrixI transpose();
+
+  MatrixI preMultiply(MatrixI m);
+
+  MatrixI postMultiply(MatrixI m);
+
+  double[] getD();
+
+  double[] getE();
+
+  void print(PrintStream ps, String format);
+
+  void printD(PrintStream ps, String format);
+
+  void printE(PrintStream ps, String format);
+
+  void tqli() throws Exception;
+
+  void tred();
+
+}
diff --git a/src/jalview/math/SparseMatrix.java b/src/jalview/math/SparseMatrix.java
new file mode 100644
index 0000000..72f0963
--- /dev/null
+++ b/src/jalview/math/SparseMatrix.java
@@ -0,0 +1,218 @@
+package jalview.math;
+
+import jalview.ext.android.SparseDoubleArray;
+
+/**
+ * A variant of Matrix intended for use for sparse (mostly zero) matrices. This
+ * class uses a SparseDoubleArray to hold each row of the matrix. The sparse
+ * array only stores non-zero values. This gives a smaller memory footprint, and
+ * fewer matrix calculation operations, for mostly zero matrices.
+ * 
+ * @author gmcarstairs
+ */
+public class SparseMatrix extends Matrix
+{
+  /*
+   * we choose columns for the sparse arrays as this allows
+   * optimisation of the preMultiply() method used in PCA.run()
+   */
+  SparseDoubleArray[] sparseColumns;
+
+  /**
+   * Constructor given data in [row][column] order
+   * 
+   * @param v
+   */
+  public SparseMatrix(double[][] v)
+  {
+    rows = v.length;
+    if (rows > 0) {
+      cols = v[0].length;
+    }
+    sparseColumns = new SparseDoubleArray[cols];
+
+    /*
+     * transpose v[row][col] into [col][row] order
+     */
+    for (int col = 0; col < cols; col++)
+    {
+      SparseDoubleArray sparseColumn = new SparseDoubleArray();
+      sparseColumns[col] = sparseColumn;
+      for (int row = 0; row < rows; row++)
+      {
+        double value = v[row][col];
+        if (value != 0d)
+        {
+          sparseColumn.put(row, value);
+        }
+      }
+    }
+  }
+
+  /**
+   * Answers the value at row i, column j
+   */
+  @Override
+  public double getValue(int i, int j)
+  {
+    return sparseColumns[j].get(i);
+  }
+
+  /**
+   * Sets the value at row i, column j to val
+   */
+  @Override
+  public void setValue(int i, int j, double val)
+  {
+    if (val == 0d)
+    {
+      sparseColumns[j].delete(i);
+    }
+    else
+    {
+      sparseColumns[j].put(i, val);
+    }
+  }
+
+  @Override
+  public double[] getColumn(int i)
+  {
+    double[] col = new double[height()];
+
+    SparseDoubleArray vals = sparseColumns[i];
+    for (int nonZero = 0; nonZero < vals.size(); nonZero++)
+    {
+      col[vals.keyAt(nonZero)] = vals.valueAt(nonZero);
+    }
+    return col;
+  }
+
+  @Override
+  public MatrixI copy()
+  {
+    double[][] vals = new double[height()][width()];
+    for (int i = 0; i < height(); i++)
+    {
+      vals[i] = getRow(i);
+    }
+    return new SparseMatrix(vals);
+  }
+
+  @Override
+  public MatrixI transpose()
+  {
+    double[][] out = new double[cols][rows];
+
+    /*
+     * for each column...
+     */
+    for (int i = 0; i < cols; i++)
+    {
+      /*
+       * put non-zero values into the corresponding row
+       * of the transposed matrix
+       */
+      SparseDoubleArray vals = sparseColumns[i];
+      for (int nonZero = 0; nonZero < vals.size(); nonZero++)
+      {
+        out[i][vals.keyAt(nonZero)] = vals.valueAt(nonZero);
+      }
+    }
+
+    return new SparseMatrix(out);
+  }
+
+  /**
+   * Answers a new matrix which is the product in.this. If the product contains
+   * less than 20% non-zero values, it is returned as a SparseMatrix, else as a
+   * Matrix.
+   * 

+   * This method is optimised for the sparse arrays which store column values
+   * for a SparseMatrix. Note that postMultiply is not so optimised. That would
+   * require redundantly also storing sparse arrays for the rows, which has not
+   * been done. Currently only preMultiply is used in Jalview.
+   */
+  @Override
+  public MatrixI preMultiply(MatrixI in)
+  {
+    if (in.width() != rows)
+    {
+      throw new IllegalArgumentException("Can't pre-multiply " + this.rows
+              + " rows by " + in.width() + " columns");
+    }
+    double[][] tmp = new double[in.height()][this.cols];
+
+    long count = 0L;
+    for (int i = 0; i < in.height(); i++)
+    {
+      for (int j = 0; j < this.cols; j++)
+      {
+        /*
+         * result[i][j] is the vector product of 
+         * in.row[i] and this.column[j]
+         * we only need to use non-zero values from the column
+         */
+        SparseDoubleArray vals = sparseColumns[j];
+        boolean added = false;
+        for (int nonZero = 0; nonZero < vals.size(); nonZero++)
+        {
+          int myRow = vals.keyAt(nonZero);
+          double myValue = vals.valueAt(nonZero);
+          tmp[i][j] += (in.getValue(i, myRow) * myValue);
+          added = true;
+        }
+        if (added && tmp[i][j] != 0d)
+        {
+          count++; // non-zero entry in product
+        }
+      }
+    }
+
+    /*
+     * heuristic rule - if product is more than 80% zero
+     * then construct a SparseMatrix, else a Matrix
+     */
+    if (count * 5 < in.height() * cols)
+    {
+      return new SparseMatrix(tmp);
+    }
+    else
+    {
+      return new Matrix(tmp);
+    }
+  }
+
+  @Override
+  protected double divideValue(int i, int j, double divisor)
+  {
+    if (divisor == 0d)
+    {
+      return getValue(i, j);
+    }
+    double v = sparseColumns[j].divide(i, divisor);
+    return v;
+  }
+
+  @Override
+  protected double addValue(int i, int j, double addend)
+  {
+    double v = sparseColumns[j].add(i, addend);
+    return v;
+  }
+
+  /**
+   * Returns the fraction of the whole matrix size that is actually modelled in
+   * sparse arrays (normally, the non-zero values)
+   * 
+   * @return
+   */
+  public float getFillRatio()
+  {
+    long count = 0L;
+    for (SparseDoubleArray col : sparseColumns)
+    {
+      count += col.size();
+    }
+    return count / (float) (height() * width());
+  }
+}
diff --git a/src/jalview/schemes/ResidueProperties.java b/src/jalview/schemes/ResidueProperties.java
index 406814a..1e6142d 100755
--- a/src/jalview/schemes/ResidueProperties.java
+++ b/src/jalview/schemes/ResidueProperties.java
@@ -486,7 +486,7 @@ public class ResidueProperties
           -3, 3, 0, -1, -4 },
       { -2, -2, 1, 6, -3, 0, 2, -1, -1, -3, -4, -1, -3, -3, -1, 0, -1, -4,
           -3, -3, 4, 1, -1, -4 },
-      { 0, 3, -3, -3, 9, -3, -4, -3, -3, -1, -1, -3, -1, -2, -3, -1, -1,
+      { 0, -3, -3, -3, 9, -3, -4, -3, -3, -1, -1, -3, -1, -2, -3, -1, -1,
           -2, -2, -1, -3, -3, -2, -4 },
       { -1, 1, 0, 0, -3, 5, 2, -2, 0, -3, -2, 1, 0, -3, -1, 0, -1, -2, -1,
           -2, 0, 3, -1, -4 },
diff --git a/src/jalview/schemes/ScoreMatrix.java b/src/jalview/schemes/ScoreMatrix.java
index 5e5fa8d..d82f54c 100644
--- a/src/jalview/schemes/ScoreMatrix.java
+++ b/src/jalview/schemes/ScoreMatrix.java
@@ -21,10 +21,10 @@
 package jalview.schemes;
 
 import jalview.analysis.scoremodels.PairwiseSeqScoreModel;
-import jalview.api.analysis.ScoreModelI;
+import jalview.math.Matrix;
+import jalview.math.MatrixI;
 
-public class ScoreMatrix extends PairwiseSeqScoreModel implements
-        ScoreModelI
+public class ScoreMatrix extends PairwiseSeqScoreModel
 {
   String name;
 
@@ -79,19 +79,21 @@ public class ScoreMatrix extends PairwiseSeqScoreModel implements
   }
 
   /**
+   * Answers the score for substituting first char in A1 with first char in A2
    * 
    * @param A1
    * @param A2
-   * @return score for substituting first char in A1 with first char in A2
+   * @return
    */
   public int getPairwiseScore(String A1, String A2)
   {
     return getPairwiseScore(A1.charAt(0), A2.charAt(0));
   }
 
+  @Override
   public int getPairwiseScore(char c, char d)
   {
-    int pog = 0;
+    int score = 0;
 
     try
     {
@@ -99,19 +101,19 @@ public class ScoreMatrix extends PairwiseSeqScoreModel implements
               : ResidueProperties.nucleotideIndex[c];
       int b = (type == 0) ? ResidueProperties.aaIndex[d]
               : ResidueProperties.nucleotideIndex[d];
-
-      pog = matrix[a][b];
+      score = matrix[a][b];
     } catch (Exception e)
     {
       // System.out.println("Unknown residue in " + A1 + " " + A2);
     }
 
-    return pog;
+    return score;
   }
 
   /**
    * pretty print the matrix
    */
+  @Override
   public String toString()
   {
     return outputMatrix(false);
@@ -170,4 +172,47 @@ public class ScoreMatrix extends PairwiseSeqScoreModel implements
     }
     return sb.toString();
   }
+
+  /**
+   * Computes an NxN matrix where N is the number of sequences, and entry [i, j]
+   * is sequence[i] pairwise multiplied with sequence[j], as a sum of scores
+   * computed using the current score matrix. For example
+   * 
    
+   * 
  • Sequences:
    • 
+   * 
  • FKL
    • 
+   * 
  • R-D
    • 
+   * 
  • QIA
    • 
+   * 
  • GWC
    • 
+   * 
  • Score matrix is BLOSUM62
    • 
+   * 
  • Gaps treated same as X (unknown)
    • 
+   * 
  • product [0, 0] = F.F + K.K + L.L = 6 + 5 + 4 = 15
    • 
+   * 
  • product [1, 1] = R.R + -.- + D.D = 5 + -1 + 6 = 10
    • 
+   * 
  • product [2, 2] = Q.Q + I.I + A.A = 5 + 4 + 4 = 13
    • 
+   * 
  • product [3, 3] = G.G + W.W + C.C = 6 + 11 + 9 = 26
    • 
+   * 
  • product[0, 1] = F.R + K.- + L.D = -3 + -1 + -3 = -8
+   * 
  • and so on
    • 
+   * 

+   */
+  public MatrixI computePairwiseScores(String[] seqs)
+  {
+    double[][] values = new double[seqs.length][];
+    for (int row = 0; row < seqs.length; row++)
+    {
+      values[row] = new double[seqs.length];
+      for (int col = 0; col < seqs.length; col++)
+      {
+        int total = 0;
+        int width = Math.min(seqs[row].length(), seqs[col].length());
+        for (int i = 0; i < width; i++)
+        {
+          char c1 = seqs[row].charAt(i);
+          char c2 = seqs[col].charAt(i);
+          int score = getPairwiseScore(c1, c2);
+          total += score;
+        }
+        values[row][col] = total;
+      }
+    }
+    return new Matrix(values);
+  }
 }
diff --git a/src/jalview/viewmodel/PCAModel.java b/src/jalview/viewmodel/PCAModel.java
index b0af302..0623dab 100644
--- a/src/jalview/viewmodel/PCAModel.java
+++ b/src/jalview/viewmodel/PCAModel.java
@@ -30,6 +30,13 @@ import java.util.Vector;
 
 public class PCAModel
 {
+  /*
+   * Jalview 2.10.1 treated gaps as X (peptide) or N (nucleotide)
+   * for pairwise scoring; 2.10.2 uses gap score (last column) in
+   * score matrix (JAL-2397)
+   * Set this flag to true (via Groovy) for 2.10.1 behaviour
+   */
+  private static boolean scoreGapAsAny = false;
 
   public PCAModel(AlignmentView seqstrings2, SequenceI[] seqs2,
           boolean nucleotide2)
@@ -69,8 +76,9 @@ public class PCAModel
 
   public void run()
   {
-
-    pca = new PCA(seqstrings.getSequenceStrings(' '), nucleotide,
+    char gapChar = scoreGapAsAny ? (nucleotide ? 'N' : 'X') : ' ';
+    String[] sequenceStrings = seqstrings.getSequenceStrings(gapChar);
+    pca = new PCA(sequenceStrings, nucleotide,
             score_matrix);
     pca.setJvCalcMode(jvCalcMode);
     pca.run();
@@ -83,32 +91,23 @@ public class PCAModel
       ii++;
     }
 
-    double[][] comps = new double[ii][ii];
-
-    for (int i = 0; i < ii; i++)
-    {
-      if (pca.getEigenvalue(i) > 1e-4)
-      {
-        comps[i] = pca.component(i);
-      }
-    }
-
-    top = pca.getM().rows - 1;
+    int height = pca.getHeight();
+    // top = pca.getM().height() - 1;
+    top = height - 1;
 
     points = new Vector();
     float[][] scores = pca.getComponents(top - 1, top - 2, top - 3, 100);
 
-    for (int i = 0; i < pca.getM().rows; i++)
+    for (int i = 0; i < height; i++)
     {
       SequencePoint sp = new SequencePoint(seqs[i], scores[i]);
       points.addElement(sp);
     }
-
   }
 
   public void updateRc(RotatableCanvasI rc)
   {
-    rc.setPoints(points, pca.getM().rows);
+    rc.setPoints(points, pca.getHeight());
   }
 
   public boolean isNucleotide()
@@ -146,9 +145,9 @@ public class PCAModel
     // note: actual indices for components are dim1-1, etc (patch for JAL-1123)
     float[][] scores = pca.getComponents(dim1 - 1, dim2 - 1, dim3 - 1, 100);
 
-    for (int i = 0; i < pca.getM().rows; i++)
+    for (int i = 0; i < pca.getHeight(); i++)
     {
-      ((SequencePoint) points.elementAt(i)).coord = scores[i];
+      points.elementAt(i).coord = scores[i];
     }
   }
 
diff --git a/test/jalview/ext/android/SparseDoubleArrayTest.java b/test/jalview/ext/android/SparseDoubleArrayTest.java
new file mode 100644
index 0000000..7d64a28
--- /dev/null
+++ b/test/jalview/ext/android/SparseDoubleArrayTest.java
@@ -0,0 +1,53 @@
+package jalview.ext.android;
+
+import static org.testng.Assert.assertEquals;
+
+import org.testng.annotations.Test;
+
+public class SparseDoubleArrayTest
+{
+
+  @Test
+  public void testConstructor()
+  {
+    double[] d = new double[] { 0d, 0d, 1.2d, 0d, 0d, 3.4d };
+    SparseDoubleArray s = new SparseDoubleArray(d);
+    for (int i = 0; i < d.length; i++)
+    {
+      assertEquals(s.get(i), d[i], "At [" + i + "]");
+    }
+  }
+
+  @Test
+  public void testAdd()
+  {
+    double[] d = new double[] { 0d, 0d, 1.2d, 0d, 0d, 3.4d };
+    SparseDoubleArray s = new SparseDoubleArray(d);
+    // add to zero (absent)
+    s.add(0, 3.2d);
+    assertEquals(s.get(0), 3.2d);
+    // add to non-zero
+    s.add(0, 2.5d);
+    assertEquals(s.get(0), 5.7d);
+    // add negative value
+    s.add(2, -5.3d);
+    assertEquals(s.get(2), -4.1d);
+    // add to unset value
+    s.add(12, 9.8d);
+    assertEquals(s.get(12), 9.8d);
+  }
+
+  @Test
+  public void testDivide()
+  {
+    double delta = 1.0e-10;
+    double[] d = new double[] { 0d, 2.4d, 1.2d, 0d, -4.8d, -3.6d };
+    SparseDoubleArray s = new SparseDoubleArray(d);
+    assertEquals(s.divide(0, 1d), 0d); // no such entry
+    assertEquals(s.divide(2, 0d), 0d); // zero divisor
+    assertEquals(s.divide(1, 2d), 1.2d, delta); // + / +
+    assertEquals(s.divide(2, -2d), -0.6d, delta); // + / -
+    assertEquals(s.divide(4, 3d), -1.6d, delta); // - / +
+    assertEquals(s.divide(5, -3d), 1.2d, delta); // - / -
+  }
+}
diff --git a/test/jalview/math/MatrixTest.java b/test/jalview/math/MatrixTest.java
index a4acbd0..961602d 100644
--- a/test/jalview/math/MatrixTest.java
+++ b/test/jalview/math/MatrixTest.java
@@ -8,9 +8,12 @@ import java.util.Arrays;
 import java.util.Random;
 
 import org.testng.annotations.Test;
+import org.testng.internal.junit.ArrayAsserts;
 
 public class MatrixTest
 {
+  final static double DELTA = 0.0001d;
+
   @Test(groups = "Timing")
   public void testPreMultiply_timing()
   {
@@ -37,20 +40,20 @@ public class MatrixTest
      * 1x3 times 3x1 is 1x1
      * 2x5 + 3x6 + 4*7 =  56
      */
-    Matrix m3 = m2.preMultiply(m1);
-    assertEquals(m3.rows, 1);
-    assertEquals(m3.cols, 1);
-    assertEquals(m3.value[0][0], 56d);
+    MatrixI m3 = m2.preMultiply(m1);
+    assertEquals(m3.height(), 1);
+    assertEquals(m3.width(), 1);
+    assertEquals(m3.getValue(0, 0), 56d);
 
     /*
      * 3x1 times 1x3 is 3x3
      */
     m3 = m1.preMultiply(m2);
-    assertEquals(m3.rows, 3);
-    assertEquals(m3.cols, 3);
-    assertEquals(Arrays.toString(m3.value[0]), "[10.0, 15.0, 20.0]");
-    assertEquals(Arrays.toString(m3.value[1]), "[12.0, 18.0, 24.0]");
-    assertEquals(Arrays.toString(m3.value[2]), "[14.0, 21.0, 28.0]");
+    assertEquals(m3.height(), 3);
+    assertEquals(m3.width(), 3);
+    assertEquals(Arrays.toString(m3.getRow(0)), "[10.0, 15.0, 20.0]");
+    assertEquals(Arrays.toString(m3.getRow(1)), "[12.0, 18.0, 24.0]");
+    assertEquals(Arrays.toString(m3.getRow(2)), "[14.0, 21.0, 28.0]");
   }
 
   @Test(
@@ -85,7 +88,18 @@ public class MatrixTest
   
   
   private boolean matrixEquals(Matrix m1, Matrix m2) {
-    return Arrays.deepEquals(m1.value, m2.value);
+    if (m1.width() != m2.width() || m1.height() != m2.height())
+    {
+      return false;
+    }
+    for (int i = 0; i < m1.height(); i++)
+    {
+      if (!Arrays.equals(m1.getRow(i), m2.getRow(i)))
+      {
+        return false;
+      }
+    }
+    return true;
   }
 
   @Test(groups = "Functional")
@@ -99,18 +113,18 @@ public class MatrixTest
      * (3020 30200)
      * (5040 50400)
      */
-    Matrix m1 = new Matrix(new double[][] { { 2, 3 }, { 4, 5 } });
-    Matrix m2 = new Matrix(new double[][] { { 10, 100 }, { 1000, 10000 } });
-    Matrix m3 = m1.postMultiply(m2);
-    assertEquals(Arrays.toString(m3.value[0]), "[3020.0, 30200.0]");
-    assertEquals(Arrays.toString(m3.value[1]), "[5040.0, 50400.0]");
+    MatrixI m1 = new Matrix(new double[][] { { 2, 3 }, { 4, 5 } });
+    MatrixI m2 = new Matrix(new double[][] { { 10, 100 }, { 1000, 10000 } });
+    MatrixI m3 = m1.postMultiply(m2);
+    assertEquals(Arrays.toString(m3.getRow(0)), "[3020.0, 30200.0]");
+    assertEquals(Arrays.toString(m3.getRow(1)), "[5040.0, 50400.0]");
 
     /*
      * also check m2.preMultiply(m1) - should be same as m1.postMultiply(m2) 
      */
     m3 = m2.preMultiply(m1);
-    assertEquals(Arrays.toString(m3.value[0]), "[3020.0, 30200.0]");
-    assertEquals(Arrays.toString(m3.value[1]), "[5040.0, 50400.0]");
+    assertEquals(Arrays.toString(m3.getRow(0)), "[3020.0, 30200.0]");
+    assertEquals(Arrays.toString(m3.getRow(1)), "[5040.0, 50400.0]");
 
     /*
      * m1 has more rows than columns
@@ -120,15 +134,15 @@ public class MatrixTest
     m1 = new Matrix(new double[][] { { 2 }, { 3 } });
     m2 = new Matrix(new double[][] { { 10, 100, 1000 } });
     m3 = m1.postMultiply(m2);
-    assertEquals(m3.rows, 2);
-    assertEquals(m3.cols, 3);
-    assertEquals(Arrays.toString(m3.value[0]), "[20.0, 200.0, 2000.0]");
-    assertEquals(Arrays.toString(m3.value[1]), "[30.0, 300.0, 3000.0]");
+    assertEquals(m3.height(), 2);
+    assertEquals(m3.width(), 3);
+    assertEquals(Arrays.toString(m3.getRow(0)), "[20.0, 200.0, 2000.0]");
+    assertEquals(Arrays.toString(m3.getRow(1)), "[30.0, 300.0, 3000.0]");
     m3 = m2.preMultiply(m1);
-    assertEquals(m3.rows, 2);
-    assertEquals(m3.cols, 3);
-    assertEquals(Arrays.toString(m3.value[0]), "[20.0, 200.0, 2000.0]");
-    assertEquals(Arrays.toString(m3.value[1]), "[30.0, 300.0, 3000.0]");
+    assertEquals(m3.height(), 2);
+    assertEquals(m3.width(), 3);
+    assertEquals(Arrays.toString(m3.getRow(0)), "[20.0, 200.0, 2000.0]");
+    assertEquals(Arrays.toString(m3.getRow(1)), "[30.0, 300.0, 3000.0]");
 
     /*
      * m1 has more columns than rows
@@ -141,19 +155,19 @@ public class MatrixTest
     m1 = new Matrix(new double[][] { { 2, 3, 4 } });
     m2 = new Matrix(new double[][] { { 5, 4 }, { 6, 3 }, { 7, 2 } });
     m3 = m1.postMultiply(m2);
-    assertEquals(m3.rows, 1);
-    assertEquals(m3.cols, 2);
-    assertEquals(m3.value[0][0], 56d);
-    assertEquals(m3.value[0][1], 25d);
+    assertEquals(m3.height(), 1);
+    assertEquals(m3.width(), 2);
+    assertEquals(m3.getRow(0)[0], 56d);
+    assertEquals(m3.getRow(0)[1], 25d);
 
     /*
      * and check premultiply equivalent
      */
     m3 = m2.preMultiply(m1);
-    assertEquals(m3.rows, 1);
-    assertEquals(m3.cols, 2);
-    assertEquals(m3.value[0][0], 56d);
-    assertEquals(m3.value[0][1], 25d);
+    assertEquals(m3.height(), 1);
+    assertEquals(m3.width(), 2);
+    assertEquals(m3.getRow(0)[0], 56d);
+    assertEquals(m3.getRow(0)[1], 25d);
   }
 
   @Test(groups = "Functional")
@@ -172,7 +186,7 @@ public class MatrixTest
       }
     }
     Matrix m1 = new Matrix(in);
-    Matrix m2 = m1.copy();
+    Matrix m2 = (Matrix) m1.copy();
     assertTrue(matrixEquals(m1, m2));
   }
 
@@ -200,12 +214,12 @@ public class MatrixTest
     // / origmat.print(System.out);
     // System.out.println();
     // System.out.println(" --- transpose matrix ---- ");
-    Matrix trans = origmat.transpose();
+    MatrixI trans = origmat.transpose();
   
     // trans.print(System.out);
     // System.out.println();
     // System.out.println(" --- OrigT * Orig ---- ");
-    Matrix symm = trans.postMultiply(origmat);
+    MatrixI symm = trans.postMultiply(origmat);
   
     // symm.print(System.out);
     // System.out.println();
@@ -255,4 +269,113 @@ public class MatrixTest
     // }
     // System.out.println();
   }
+
+  @Test(groups = "Timing")
+  public void testSign()
+  {
+    assertEquals(Matrix.sign(-1, -2), -1d);
+    assertEquals(Matrix.sign(-1, 2), 1d);
+    assertEquals(Matrix.sign(-1, 0), 1d);
+    assertEquals(Matrix.sign(1, -2), -1d);
+    assertEquals(Matrix.sign(1, 2), 1d);
+    assertEquals(Matrix.sign(1, 0), 1d);
+  }
+
+  /**
+   * Helper method to make values for a sparse, pseudo-random symmetric matrix
+   * 
+   * @param rows
+   * @param cols
+   * @param occupancy
+   *          one in 'occupancy' entries will be non-zero
+   * @return
+   */
+  public double[][] getSparseValues(int rows, int cols, int occupancy)
+  {
+    Random r = new Random(1729);
+
+    /*
+     * generate whole number values between -12 and +12
+     * (to mimic score matrices used in Jalview)
+     */
+    double[][] d = new double[rows][cols];
+    int m = 0;
+    for (int i = 0; i < rows; i++)
+    {
+      if (++m % occupancy == 0)
+      {
+        d[i][i] = r.nextInt() % 13; // diagonal
+      }
+      for (int j = 0; j < i; j++)
+      {
+        if (++m % occupancy == 0)
+        {
+          d[i][j] = r.nextInt() % 13;
+          d[j][i] = d[i][j];
+        }
+      }
+    }
+    return d;
+  
+  }
+
+  /**
+   * Verify that the results of method tred() are the same if the calculation is
+   * redone
+   */
+  @Test(groups = "Functional")
+  public void testTred_reproducible()
+  {
+    /*
+     * make a pseudo-random symmetric matrix as required for tred/tqli
+     */
+    int rows = 10;
+    int cols = rows;
+    double[][] d = getSparseValues(rows, cols, 3);
+  
+    /*
+     * make a copy of the values so m1, m2 are not
+     * sharing arrays!
+     */
+    double[][] d1 = new double[rows][cols];
+    for (int row = 0; row < rows; row++)
+    {
+      for (int col = 0; col < cols; col++)
+      {
+        d1[row][col] = d[row][col];
+      }
+    }
+    Matrix m1 = new Matrix(d);
+    Matrix m2 = new Matrix(d1);
+    assertMatricesMatch(m1, m2); // sanity check
+    m1.tred();
+    m2.tred();
+    assertMatricesMatch(m1, m2);
+  }
+
+  private void assertMatricesMatch(MatrixI m1, MatrixI m2)
+  {
+    if (m1.height() != m2.height())
+    {
+      fail("height mismatch");
+    }
+    if (m1.width() != m2.width())
+    {
+      fail("width mismatch");
+    }
+    for (int row = 0; row < m1.height(); row++)
+    {
+      for (int col = 0; col < m1.width(); col++)
+      {
+        double v2 = m2.getValue(row, col);
+        double v1 = m1.getValue(row, col);
+        if (Math.abs(v1 - v2) > DELTA)
+        {
+          fail(String.format("At [%d, %d] %f != %f", row, col, v1, v2));
+        }
+      }
+    }
+    ArrayAsserts.assertArrayEquals(m1.getD(), m2.getD(), 0.00001d);
+    ArrayAsserts.assertArrayEquals(m1.getE(), m2.getE(), 0.00001d);
+  }
 }
diff --git a/test/jalview/math/SparseMatrixTest.java b/test/jalview/math/SparseMatrixTest.java
new file mode 100644
index 0000000..3c2ccaa
--- /dev/null
+++ b/test/jalview/math/SparseMatrixTest.java
@@ -0,0 +1,416 @@
+package jalview.math;
+
+import static org.testng.Assert.assertEquals;
+import static org.testng.Assert.assertFalse;
+import static org.testng.Assert.assertTrue;
+import static org.testng.Assert.fail;
+
+import java.util.Random;
+
+import org.testng.annotations.Test;
+import org.testng.internal.junit.ArrayAsserts;
+
+public class SparseMatrixTest
+{
+  final static double DELTA = 0.0001d;
+
+  Random r = new Random(1729);
+
+  @Test(groups = "Functional")
+  public void testConstructor()
+  {
+    MatrixI m1 = new SparseMatrix(
+            new double[][] { { 2, 0, 4 }, { 0, 6, 0 } });
+    assertEquals(m1.getValue(0, 0), 2d);
+    assertEquals(m1.getValue(0, 1), 0d);
+    assertEquals(m1.getValue(0, 2), 4d);
+    assertEquals(m1.getValue(1, 0), 0d);
+    assertEquals(m1.getValue(1, 1), 6d);
+    assertEquals(m1.getValue(1, 2), 0d);
+  }
+
+  @Test(groups = "Functional")
+  public void testTranspose()
+  {
+    MatrixI m1 = new SparseMatrix(
+            new double[][] { { 2, 0, 4 }, { 5, 6, 0 } });
+    MatrixI m2 = m1.transpose();
+    assertTrue(m2 instanceof SparseMatrix);
+    assertEquals(m2.height(), 3);
+    assertEquals(m2.width(), 2);
+    assertEquals(m2.getValue(0, 0), 2d);
+    assertEquals(m2.getValue(0, 1), 5d);
+    assertEquals(m2.getValue(1, 0), 0d);
+    assertEquals(m2.getValue(1, 1), 6d);
+    assertEquals(m2.getValue(2, 0), 4d);
+    assertEquals(m2.getValue(2, 1), 0d);
+  }
+  @Test(groups = "Functional")
+  public void testPreMultiply()
+  {
+    MatrixI m1 = new SparseMatrix(new double[][] { { 2, 3, 4 } }); // 1x3
+    MatrixI m2 = new SparseMatrix(new double[][] { { 5 }, { 6 }, { 7 } }); // 3x1
+
+    /*
+     * 1x3 times 3x1 is 1x1
+     * 2x5 + 3x6 + 4*7 =  56
+     */
+    MatrixI m3 = m2.preMultiply(m1);
+    assertFalse(m3 instanceof SparseMatrix);
+    assertEquals(m3.height(), 1);
+    assertEquals(m3.width(), 1);
+    assertEquals(m3.getValue(0, 0), 56d);
+
+    /*
+     * 3x1 times 1x3 is 3x3
+     */
+    m3 = m1.preMultiply(m2);
+    assertEquals(m3.height(), 3);
+    assertEquals(m3.width(), 3);
+    assertEquals(m3.getValue(0, 0), 10d);
+    assertEquals(m3.getValue(0, 1), 15d);
+    assertEquals(m3.getValue(0, 2), 20d);
+    assertEquals(m3.getValue(1, 0), 12d);
+    assertEquals(m3.getValue(1, 1), 18d);
+    assertEquals(m3.getValue(1, 2), 24d);
+    assertEquals(m3.getValue(2, 0), 14d);
+    assertEquals(m3.getValue(2, 1), 21d);
+    assertEquals(m3.getValue(2, 2), 28d);
+  }
+
+  @Test(
+    groups = "Functional",
+    expectedExceptions = { IllegalArgumentException.class })
+  public void testPreMultiply_tooManyColumns()
+  {
+    Matrix m1 = new SparseMatrix(
+            new double[][] { { 2, 3, 4 }, { 3, 4, 5 } }); // 2x3
+
+    /*
+     * 2x3 times 2x3 invalid operation - 
+     * multiplier has more columns than multiplicand has rows
+     */
+    m1.preMultiply(m1);
+    fail("Expected exception");
+  }
+
+  @Test(
+    groups = "Functional",
+    expectedExceptions = { IllegalArgumentException.class })
+  public void testPreMultiply_tooFewColumns()
+  {
+    Matrix m1 = new SparseMatrix(
+            new double[][] { { 2, 3, 4 }, { 3, 4, 5 } }); // 2x3
+
+    /*
+     * 3x2 times 3x2 invalid operation - 
+     * multiplier has more columns than multiplicand has row
+     */
+    m1.preMultiply(m1);
+    fail("Expected exception");
+  }
+  
+  @Test(groups = "Functional")
+  public void testPostMultiply()
+  {
+    /*
+     * Square matrices
+     * (2 3) . (10   100)
+     * (4 5)   (1000 10000)
+     * =
+     * (3020 30200)
+     * (5040 50400)
+     */
+    MatrixI m1 = new SparseMatrix(new double[][] { { 2, 3 }, { 4, 5 } });
+    MatrixI m2 = new SparseMatrix(new double[][] { { 10, 100 },
+        { 1000, 10000 } });
+    MatrixI m3 = m1.postMultiply(m2);
+    assertEquals(m3.getValue(0, 0), 3020d);
+    assertEquals(m3.getValue(0, 1), 30200d);
+    assertEquals(m3.getValue(1, 0), 5040d);
+    assertEquals(m3.getValue(1, 1), 50400d);
+
+    /*
+     * also check m2.preMultiply(m1) - should be same as m1.postMultiply(m2) 
+     */
+    MatrixI m4 = m2.preMultiply(m1);
+    assertMatricesMatch(m3, m4, 0.00001d);
+
+    /*
+     * m1 has more rows than columns
+     * (2).(10 100 1000) = (20 200 2000)
+     * (3)                 (30 300 3000)
+     */
+    m1 = new SparseMatrix(new double[][] { { 2 }, { 3 } });
+    m2 = new SparseMatrix(new double[][] { { 10, 100, 1000 } });
+    m3 = m1.postMultiply(m2);
+    assertEquals(m3.height(), 2);
+    assertEquals(m3.width(), 3);
+    assertEquals(m3.getValue(0, 0), 20d);
+    assertEquals(m3.getValue(0, 1), 200d);
+    assertEquals(m3.getValue(0, 2), 2000d);
+    assertEquals(m3.getValue(1, 0), 30d);
+    assertEquals(m3.getValue(1, 1), 300d);
+    assertEquals(m3.getValue(1, 2), 3000d);
+
+    m4 = m2.preMultiply(m1);
+    assertMatricesMatch(m3, m4, 0.00001d);
+
+    /*
+     * m1 has more columns than rows
+     * (2 3 4) . (5 4) = (56 25)
+     *           (6 3) 
+     *           (7 2)
+     * [0, 0] = 2*5 + 3*6 + 4*7 = 56
+     * [0, 1] = 2*4 + 3*3 + 4*2 = 25  
+     */
+    m1 = new SparseMatrix(new double[][] { { 2, 3, 4 } });
+    m2 = new SparseMatrix(new double[][] { { 5, 4 }, { 6, 3 }, { 7, 2 } });
+    m3 = m1.postMultiply(m2);
+    assertEquals(m3.height(), 1);
+    assertEquals(m3.width(), 2);
+    assertEquals(m3.getValue(0, 0), 56d);
+    assertEquals(m3.getValue(0, 1), 25d);
+
+    /*
+     * and check premultiply equivalent
+     */
+    m4 = m2.preMultiply(m1);
+    assertMatricesMatch(m3, m4, 0.00001d);
+  }
+
+  @Test(groups = "Timing")
+  public void testSign()
+  {
+    assertEquals(Matrix.sign(-1, -2), -1d);
+    assertEquals(Matrix.sign(-1, 2), 1d);
+    assertEquals(Matrix.sign(-1, 0), 1d);
+    assertEquals(Matrix.sign(1, -2), -1d);
+    assertEquals(Matrix.sign(1, 2), 1d);
+    assertEquals(Matrix.sign(1, 0), 1d);
+  }
+
+  /**
+   * Verify that the results of method tred() are the same for SparseMatrix as
+   * they are for Matrix (i.e. a regression test rather than an absolute test of
+   * correctness of results)
+   */
+  @Test(groups = "Functional")
+  public void testTred_matchesMatrix()
+  {
+    /*
+     * make a pseudo-random symmetric matrix as required for tred/tqli
+     */
+    int rows = 10;
+    int cols = rows;
+    double[][] d = getSparseValues(rows, cols, 3);
+
+    /*
+     * make a copy of the values so m1, m2 are not
+     * sharing arrays!
+     */
+    double[][] d1 = new double[rows][cols];
+    for (int row = 0; row < rows; row++)
+    {
+      for (int col = 0; col < cols; col++)
+      {
+        d1[row][col] = d[row][col];
+      }
+    }
+    Matrix m1 = new Matrix(d);
+    Matrix m2 = new SparseMatrix(d1);
+    assertMatricesMatch(m1, m2, 0.00001d); // sanity check
+    m1.tred();
+    m2.tred();
+    assertMatricesMatch(m1, m2, 0.00001d);
+  }
+
+  private void assertMatricesMatch(MatrixI m1, MatrixI m2, double delta)
+  {
+    if (m1.height() != m2.height())
+    {
+      fail("height mismatch");
+    }
+    if (m1.width() != m2.width())
+    {
+      fail("width mismatch");
+    }
+    for (int row = 0; row < m1.height(); row++)
+    {
+      for (int col = 0; col < m1.width(); col++)
+      {
+        double v2 = m2.getValue(row, col);
+        double v1 = m1.getValue(row, col);
+        if (Math.abs(v1 - v2) > DELTA)
+        {
+          fail(String.format("At [%d, %d] %f != %f", row, col, v1, v2));
+        }
+      }
+    }
+    ArrayAsserts.assertArrayEquals(m1.getD(), m2.getD(), delta);
+    ArrayAsserts.assertArrayEquals(m1.getE(), m2.getE(), 0.00001d);
+  }
+
+  @Test
+  public void testGetValue()
+  {
+    double[][] d = new double[][] { { 0, 0, 1, 0, 0 }, { 2, 3, 0, 0, 0 },
+        { 4, 0, 0, 0, 5 } };
+    MatrixI m = new SparseMatrix(d);
+    for (int row = 0; row < 3; row++)
+    {
+      for (int col = 0; col < 5; col++)
+      {
+        assertEquals(m.getValue(row, col), d[row][col],
+                String.format("At [%d, %d]", row, col));
+      }
+    }
+  }
+
+  /**
+   * Verify that the results of method tqli() are the same for SparseMatrix as
+   * they are for Matrix (i.e. a regression test rather than an absolute test of
+   * correctness of results)
+   * 
+   * @throws Exception
+   */
+  @Test(groups = "Functional")
+  public void testTqli_matchesMatrix() throws Exception
+  {
+    /*
+     * make a pseudo-random symmetric matrix as required for tred
+     */
+    int rows = 6;
+    int cols = rows;
+    double[][] d = getSparseValues(rows, cols, 3);
+  
+    /*
+     * make a copy of the values so m1, m2 are not
+     * sharing arrays!
+     */
+    double[][] d1 = new double[rows][cols];
+    for (int row = 0; row < rows; row++)
+    {
+      for (int col = 0; col < cols; col++)
+      {
+        d1[row][col] = d[row][col];
+      }
+    }
+    Matrix m1 = new Matrix(d);
+    Matrix m2 = new SparseMatrix(d1);
+
+    // have to do tred() before doing tqli()
+    m1.tred();
+    m2.tred();
+    assertMatricesMatch(m1, m2, 0.00001d);
+
+    m1.tqli();
+    m2.tqli();
+    assertMatricesMatch(m1, m2, 0.00001d);
+  }
+
+  /**
+   * Helper method to make values for a sparse, pseudo-random symmetric matrix
+   * 
+   * @param rows
+   * @param cols
+   * @param occupancy
+   *          one in 'occupancy' entries will be non-zero
+   * @return
+   */
+  public double[][] getSparseValues(int rows, int cols, int occupancy)
+  {
+    /*
+     * generate whole number values between -12 and +12
+     * (to mimic score matrices used in Jalview)
+     */
+    double[][] d = new double[rows][cols];
+    int m = 0;
+    for (int i = 0; i < rows; i++)
+    {
+      if (++m % occupancy == 0)
+      {
+        d[i][i] = r.nextInt() % 13; // diagonal
+      }
+      for (int j = 0; j < i; j++)
+      {
+        if (++m % occupancy == 0)
+        {
+          d[i][j] = r.nextInt() % 13;
+          d[j][i] = d[i][j];
+        }
+      }
+    }
+    return d;
+
+  }
+
+  /**
+   * Test that verifies that the result of preMultiply is a SparseMatrix if more
+   * than 80% zeroes, else a Matrix
+   */
+  @Test(groups = "Functional")
+  public void testPreMultiply_sparseProduct()
+  {
+    MatrixI m1 = new SparseMatrix(new double[][] { { 1 }, { 0 }, { 0 },
+        { 0 }, { 0 } }); // 5x1
+    MatrixI m2 = new SparseMatrix(new double[][] { { 1, 1, 1, 1 } }); // 1x4
+  
+    /*
+     * m1.m2 makes a row of 4 1's, and 4 rows of zeros
+     * 20% non-zero so not 'sparse'
+     */
+    MatrixI m3 = m2.preMultiply(m1);
+    assertFalse(m3 instanceof SparseMatrix);
+
+    /*
+     * replace a 1 with a 0 in the product:
+     * it is now > 80% zero so 'sparse'
+     */
+    m2 = new SparseMatrix(new double[][] { { 1, 1, 1, 0 } });
+    m3 = m2.preMultiply(m1);
+    assertTrue(m3 instanceof SparseMatrix);
+  }
+
+  @Test(groups = "Functional")
+  public void testFillRatio()
+  {
+    SparseMatrix m1 = new SparseMatrix(new double[][] { { 2, 0, 4, 1, 0 },
+    { 0, 6, 0, 0, 0 } });
+    assertEquals(m1.getFillRatio(), 0.4f);
+  }
+
+  /**
+   * Verify that the results of method tred() are the same if the calculation is
+   * redone
+   */
+  @Test(groups = "Functional")
+  public void testTred_reproducible()
+  {
+    /*
+     * make a pseudo-random symmetric matrix as required for tred/tqli
+     */
+    int rows = 10;
+    int cols = rows;
+    double[][] d = getSparseValues(rows, cols, 3);
+  
+    /*
+     * make a copy of the values so m1, m2 are not
+     * sharing arrays!
+     */
+    double[][] d1 = new double[rows][cols];
+    for (int row = 0; row < rows; row++)
+    {
+      for (int col = 0; col < cols; col++)
+      {
+        d1[row][col] = d[row][col];
+      }
+    }
+    Matrix m1 = new SparseMatrix(d);
+    Matrix m2 = new SparseMatrix(d1);
+    assertMatricesMatch(m1, m2, 1.0e16); // sanity check
+    m1.tred();
+    m2.tred();
+    assertMatricesMatch(m1, m2, 0.00001d);
+  }
+}
\ No newline at end of file
diff --git a/test/jalview/schemes/ScoreMatrixPrinter.java b/test/jalview/schemes/ScoreMatrixPrinter.java
index a743163..80241fb 100644
--- a/test/jalview/schemes/ScoreMatrixPrinter.java
+++ b/test/jalview/schemes/ScoreMatrixPrinter.java
@@ -26,7 +26,6 @@ import jalview.gui.JvOptionPane;
 import java.util.Map;
 
 import org.testng.annotations.BeforeClass;
-import org.testng.annotations.Test;
 
 public class ScoreMatrixPrinter
 {
@@ -38,7 +37,6 @@ public class ScoreMatrixPrinter
     JvOptionPane.setMockResponse(JvOptionPane.CANCEL_OPTION);
   }
 
-  @Test(groups = { "Functional" })
   public void printAllMatrices()
   {
     for (Map.Entry sm : ResidueProperties.scoreMatrices
@@ -49,7 +47,6 @@ public class ScoreMatrixPrinter
     }
   }
 
-  @Test(groups = { "Functional" })
   public void printHTMLMatrices()
   {
     for (Map.Entry _sm : ResidueProperties.scoreMatrices
diff --git a/test/jalview/schemes/ScoreMatrixTest.java b/test/jalview/schemes/ScoreMatrixTest.java
new file mode 100644
index 0000000..e15dd41
--- /dev/null
+++ b/test/jalview/schemes/ScoreMatrixTest.java
@@ -0,0 +1,168 @@
+package jalview.schemes;
+
+import static org.testng.Assert.assertEquals;
+
+import jalview.math.MatrixI;
+
+import org.testng.annotations.Test;
+
+public class ScoreMatrixTest
+{
+  @Test(groups = "Functional")
+  public void testSymmetric()
+  {
+    verifySymmetric(ResidueProperties.getScoreMatrix("BLOSUM62"));
+    verifySymmetric(ResidueProperties.getScoreMatrix("PAM250"));
+    verifySymmetric(ResidueProperties.getScoreMatrix("DNA"));
+  }
+
+  private void verifySymmetric(ScoreMatrix sm)
+  {
+    int[][] m = sm.getMatrix();
+    int rows = m.length;
+    for (int row = 0; row < rows; row++)
+    {
+      assertEquals(m[row].length, rows);
+      for (int col = 0; col < rows; col++)
+      {
+        assertEquals(m[row][col], m[col][row], String.format("%s [%s, %s]",
+                sm.getName(), ResidueProperties.aa[row],
+                ResidueProperties.aa[col]));
+      }
+    }
+
+    /*
+     * also check the score matrix is sized for 
+     * the number of symbols scored, plus gap
+     */
+    assertEquals(rows, (sm.isDNA() ? ResidueProperties.maxNucleotideIndex
+            : ResidueProperties.maxProteinIndex) + 1);
+  }
+
+  /**
+   * A test that just asserts the expected values in the Blosum62 score matrix
+   */
+  @Test(groups = "Functional")
+  public void testBlosum62_values()
+  {
+    ScoreMatrix sm = ResidueProperties.getScoreMatrix("BLOSUM62");
+
+    /*
+     * verify expected scores against ARNDCQEGHILKMFPSTWYVBZX
+     * scraped from https://www.ncbi.nlm.nih.gov/Class/FieldGuide/BLOSUM62.txt
+     */
+    verifyValues(sm, 'A', new int[] { 4, -1, -2, -2, 0, -1, -1, 0, -2, -1,
+        -1, -1, -1, -2, -1, 1, 0, -3, -2, 0, -2, -1, 0 });
+    verifyValues(sm, 'R', new int[] { -1, 5, 0, -2, -3, 1, 0, -2, 0, -3,
+        -2, 2, -1, -3, -2, -1, -1, -3, -2, -3, -1, 0, -1 });
+    verifyValues(sm, 'N', new int[] { -2, 0, 6, 1, -3, 0, 0, 0, 1, -3, -3,
+        0, -2, -3, -2, 1, 0, -4, -2, -3, 3, 0, -1 });
+    verifyValues(sm, 'D', new int[] { -2, -2, 1, 6, -3, 0, 2, -1, -1, -3,
+        -4, -1, -3, -3, -1, 0, -1, -4, -3, -3, 4, 1, -1 });
+    verifyValues(sm, 'C', new int[] { 0, -3, -3, -3, 9, -3, -4, -3, -3, -1,
+        -1, -3, -1, -2, -3, -1, -1, -2, -2, -1, -3, -3, -2 });
+    verifyValues(sm, 'Q', new int[] { -1, 1, 0, 0, -3, 5, 2, -2, 0, -3, -2,
+        1, 0, -3, -1, 0, -1, -2, -1, -2, 0, 3, -1 });
+    verifyValues(sm, 'E', new int[] { -1, 0, 0, 2, -4, 2, 5, -2, 0, -3, -3,
+        1, -2, -3, -1, 0, -1, -3, -2, -2, 1, 4, -1 });
+    verifyValues(sm, 'G', new int[] { 0, -2, 0, -1, -3, -2, -2, 6, -2, -4,
+        -4, -2, -3, -3, -2, 0, -2, -2, -3, -3, -1, -2, -1 });
+    verifyValues(sm, 'H', new int[] { -2, 0, 1, -1, -3, 0, 0, -2, 8, -3,
+        -3, -1, -2, -1, -2, -1, -2, -2, 2, -3, 0, 0, -1 });
+    verifyValues(sm, 'I', new int[] { -1, -3, -3, -3, -1, -3, -3, -4, -3,
+        4, 2, -3, 1, 0, -3, -2, -1, -3, -1, 3, -3, -3, -1 });
+    verifyValues(sm, 'L', new int[] { -1, -2, -3, -4, -1, -2, -3, -4, -3,
+        2, 4, -2, 2, 0, -3, -2, -1, -2, -1, 1, -4, -3, -1 });
+    verifyValues(sm, 'K', new int[] { -1, 2, 0, -1, -3, 1, 1, -2, -1, -3,
+        -2, 5, -1, -3, -1, 0, -1, -3, -2, -2, 0, 1, -1 });
+    verifyValues(sm, 'M', new int[] { -1, -1, -2, -3, -1, 0, -2, -3, -2, 1,
+        2, -1, 5, 0, -2, -1, -1, -1, -1, 1, -3, -1, -1 });
+    verifyValues(sm, 'F', new int[] { -2, -3, -3, -3, -2, -3, -3, -3, -1,
+        0, 0, -3, 0, 6, -4, -2, -2, 1, 3, -1, -3, -3, -1 });
+    verifyValues(sm, 'P', new int[] { -1, -2, -2, -1, -3, -1, -1, -2, -2,
+        -3, -3, -1, -2, -4, 7, -1, -1, -4, -3, -2, -2, -1, -2 });
+    verifyValues(sm, 'S', new int[] { 1, -1, 1, 0, -1, 0, 0, 0, -1, -2, -2,
+        0, -1, -2, -1, 4, 1, -3, -2, -2, 0, 0, 0 });
+    verifyValues(sm, 'T', new int[] { 0, -1, 0, -1, -1, -1, -1, -2, -2, -1,
+        -1, -1, -1, -2, -1, 1, 5, -2, -2, 0, -1, -1, 0 });
+    verifyValues(sm, 'W', new int[] { -3, -3, -4, -4, -2, -2, -3, -2, -2,
+        -3, -2, -3, -1, 1, -4, -3, -2, 11, 2, -3, -4, -3, -2 });
+    verifyValues(sm, 'Y', new int[] { -2, -2, -2, -3, -2, -1, -2, -3, 2,
+        -1, -1, -2, -1, 3, -3, -2, -2, 2, 7, -1, -3, -2, -1 });
+    verifyValues(sm, 'V', new int[] { 0, -3, -3, -3, -1, -2, -2, -3, -3, 3,
+        1, -2, 1, -1, -2, -2, 0, -3, -1, 4, -3, -2, -1 });
+    verifyValues(sm, 'B', new int[] { -2, -1, 3, 4, -3, 0, 1, -1, 0, -3,
+        -4, 0, -3, -3, -2, 0, -1, -4, -3, -3, 4, 1, -1 });
+    verifyValues(sm, 'Z', new int[] { -1, 0, 0, 1, -3, 3, 4, -2, 0, -3, -3,
+        1, -1, -3, -1, 0, -1, -3, -2, -2, 1, 4, -1 });
+    verifyValues(sm, 'X', new int[] { 0, -1, -1, -1, -2, -1, -1, -1, -1,
+        -1, -1, -1, -1, -1, -2, 0, 0, -2, -1, -1, -1, -1, -1 });
+  }
+  /**
+   * Helper method to check pairwise scores for one residue
+   * 
+   * @param sm
+   * @param res
+   * @param expected
+   *          score values against 'res', in ResidueProperties.aaIndex order
+   */
+  private void verifyValues(ScoreMatrix sm, char res, int[] expected)
+  {
+    for (int j = 0; j < expected.length; j++)
+    {
+      char c2 = ResidueProperties.aa[j].charAt(0);
+      assertEquals(sm.getPairwiseScore(res, c2), expected[j],
+              String.format("%s->%s", res, c2));
+    }
+  }
+
+  @Test(groups = "Functional")
+  public void testComputePairwiseScores()
+  {
+    String[] seqs = new String[] { "FKL", "R-D", "QIA", "GWC" };
+    ScoreMatrix sm = ResidueProperties.getScoreMatrix("BLOSUM62");
+  
+    MatrixI pairwise = sm.computePairwiseScores(seqs);
+  
+    /*
+     * should be NxN where N = number of sequences
+     */
+    assertEquals(pairwise.height(), 4);
+    assertEquals(pairwise.width(), 4);
+  
+    /*
+     * should be symmetrical (because BLOSUM62 is)
+     */
+    for (int i = 0; i < pairwise.height(); i++)
+    {
+      for (int j = 0; j < pairwise.width(); j++)
+      {
+        assertEquals(pairwise.getValue(i, j), pairwise.getValue(j, i),
+                "Not symmetric");
+      }
+    }
+    /*
+     * verify expected BLOSUM dot product scores
+     */
+    // F.F + K.K + L.L = 6 + 5 + 4 = 15
+    assertEquals(pairwise.getValue(0, 0), 15d);
+    // R.R + -.- + D.D = 5 + 1 + 6 = 12
+    assertEquals(pairwise.getValue(1, 1), 12d);
+    // Q.Q + I.I + A.A = 5 + 4 + 4 = 13
+    assertEquals(pairwise.getValue(2, 2), 13d);
+    // G.G + W.W + C.C = 6 + 11 + 9 = 26
+    assertEquals(pairwise.getValue(3, 3), 26d);
+    // F.R + K.- + L.D = -3 + -4 + -4 = -11
+    assertEquals(pairwise.getValue(0, 1), -11d);
+    // F.Q + K.I + L.A = -3 + -3 + -1 = -7
+    assertEquals(pairwise.getValue(0, 2), -7d);
+    // F.G + K.W + L.C = -3 + -3 + -1 = -7
+    assertEquals(pairwise.getValue(0, 3), -7d);
+    // R.Q + -.I + D.A = 1 + -4 + -2 = -5
+    assertEquals(pairwise.getValue(1, 2), -5d);
+    // R.G + -.W + D.C = -2 + -4 + -3 = -9
+    assertEquals(pairwise.getValue(1, 3), -9d);
+    // Q.G + I.W + A.C = -2 + -3 + 0 = -5
+    assertEquals(pairwise.getValue(2, 3), -5d);
+  }
+}