*/
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;
/**
* 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
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;
}
/**
*/
public double getEigenvalue(int i)
{
- return eigenvector.d[i];
+ return eigenvector.getD()[i];
}
/**
*/
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;
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);
}
{
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()
}
};
+ // 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();
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
}
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<symm.rows;seq++) { ps.print("\"Seq"+seq+"\""); for
*
* ps.print(","+component(seq, ev)); } ps.println(); }
*/
+ // System.out.println(("PCA.run() took "
+ // + (System.currentTimeMillis() - now) + "ms"));
}
- boolean jvCalcMode = true;
-
public void setJvCalcMode(boolean calcMode)
{
this.jvCalcMode = calcMode;
}
+
+ /**
+ * Answers the N dimensions of the NxN PCA matrix. This is the number of
+ * sequences involved in the pairwise score calculation.
+ *
+ * @return
+ */
+ public int getHeight()
+ {
+ // TODO can any of seqs[] be null?
+ return seqs.length;
+ }
}
{
int nores = (isNa) ? ResidueProperties.maxNucleotideIndex
: ResidueProperties.maxProteinIndex;
- // Set all matrix to 0
+
dbinary = new double[getSequence().length * nores];
- for (int i = 0; i < dbinary.length; i++)
- {
- dbinary[i] = 0.0;
- }
return nores;
}
private void matrixEncode(final int[] aaIndex, final int[][] matrix)
{
- // Set all matrix to 0
- // dbinary = new double[getSequence().length * 21];
-
int nores = initMatrixGetNoRes();
- // for (int i = 0; i < dbinary.length; i++) {
- // dbinary[i] = 0.0;
- // }
for (int i = 0, iSize = getSequence().length; i < iSize; i++)
{
int aanum = nores - 1;
/*
* Copied to Jalview September 2016.
* Only the members of this class required for SparseIntArray were copied.
- * Method binarySearch(short[] array, int size, short value) added to support
+ * Change Log:
+ * Sep 2016: Method binarySearch(short[] array, int size, short value) added to support
* SparseShortArray.
+ * Jan 2017: EMPTY_DOUBLES added
*/
class ContainerHelpers
{
static final int[] EMPTY_INTS = new int[0];
+ static final double[] EMPTY_DOUBLES = new double[0];
+
static final long[] EMPTY_LONGS = new long[0];
static final Object[] EMPTY_OBJECTS = new Object[0];
--- /dev/null
+package jalview.ext.android;
+
+/*
+ * Copyright (C) 2006 The Android Open Source Project
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+/**
+ * SparseDoubleArray map integers to doubles. Unlike a normal array of integers,
+ * there can be gaps in the indices. It is intended to be more memory efficient
+ * than using a HashMap to map Integer to Double, both because it avoids
+ * auto-boxing keys and values and its data structure doesn't rely on an extra
+ * entry object for each mapping.
+ *
+ * <p>
+ * 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%.
+ * </p>
+ *
+ * <p>
+ * It is possible to iterate over the items in this container using
+ * {@link #keyAt(int)} and {@link #valueAt(int)}. Iterating over the keys using
+ * <code>keyAt(int)</code> 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 <code>valueAt(int)<code>.
+ * </p>
+ */
+
+/*
+ * 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 <code>0</code> 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 <code>0...size()-1</code>, returns the key from
+ * the <code>index</code>th key-value mapping that this SparseDoubleArray
+ * stores.
+ *
+ * <p>
+ * The keys corresponding to indices in ascending order are guaranteed to be
+ * in ascending order, e.g., <code>keyAt(0)</code> will return the smallest
+ * key and <code>keyAt(size()-1)</code> will return the largest key.
+ * </p>
+ */
+ public int keyAt(int index)
+ {
+ return mKeys[index];
+ }
+
+ /**
+ * Given an index in the range <code>0...size()-1</code>, returns the value
+ * from the <code>index</code>th key-value mapping that this SparseDoubleArray
+ * stores.
+ *
+ * <p>
+ * The values corresponding to indices in ascending order are guaranteed to be
+ * associated with keys in ascending order, e.g., <code>valueAt(0)</code> will
+ * return the value associated with the smallest key and
+ * <code>valueAt(size()-1)</code> will return the value associated with the
+ * largest key.
+ * </p>
+ */
+ 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}
+ *
+ * <p>
+ * 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;
+ }
+}
/**
* 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
*
* <pre>
}
/**
- * 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];
*
* @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();
* 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]);
}
}
}
* 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);
}
*
* @return
*/
- public Matrix copy()
+ @Override
+ public MatrixI copy()
{
double[][] newmat = new double[rows][cols];
/**
* DOCUMENT ME!
*/
+ @Override
public void tred()
{
int n = rows;
- int l;
int k;
int j;
int i;
for (i = n; i >= 2; i--)
{
- l = i - 1;
+ final int l = i - 1;
h = 0.0;
scale = 0.0;
{
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)
{
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;
for (i = 1; i <= n; i++)
{
- l = i - 1;
+ final int l = i - 1;
if (d[i - 1] != 0.0)
{
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;
double s;
double r;
double p;
- ;
double g;
double f;
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));
}
}
}
}
+ @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!
*/
}
/**
- * 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)
{
*
* @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]);
}
}
*
* @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;
+ }
}
--- /dev/null
+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();
+
+}
--- /dev/null
+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.
+ * <p>
+ * 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());
+ }
+}
-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 },
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;
}
/**
+ * 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
{
: 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);
}
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
+ * <ul>
+ * <li>Sequences:</li>
+ * <li>FKL</li>
+ * <li>R-D</li>
+ * <li>QIA</li>
+ * <li>GWC</li>
+ * <li>Score matrix is BLOSUM62</li>
+ * <li>Gaps treated same as X (unknown)</li>
+ * <li>product [0, 0] = F.F + K.K + L.L = 6 + 5 + 4 = 15</li>
+ * <li>product [1, 1] = R.R + -.- + D.D = 5 + -1 + 6 = 10</li>
+ * <li>product [2, 2] = Q.Q + I.I + A.A = 5 + 4 + 4 = 13</li>
+ * <li>product [3, 3] = G.G + W.W + C.C = 6 + 11 + 9 = 26</li>
+ * <li>product[0, 1] = F.R + K.- + L.D = -3 + -1 + -3 = -8
+ * <li>and so on</li>
+ * </ul>
+ */
+ 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);
+ }
}
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)
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();
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<SequencePoint>();
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()
// 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];
}
}
--- /dev/null
+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); // - / -
+ }
+}
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()
{
* 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(
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")
* (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
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
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")
}
}
Matrix m1 = new Matrix(in);
- Matrix m2 = m1.copy();
+ Matrix m2 = (Matrix) m1.copy();
assertTrue(matrixEquals(m1, m2));
}
// / 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();
// }
// 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);
+ }
}
--- /dev/null
+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
import java.util.Map;
import org.testng.annotations.BeforeClass;
-import org.testng.annotations.Test;
public class ScoreMatrixPrinter
{
JvOptionPane.setMockResponse(JvOptionPane.CANCEL_OPTION);
}
- @Test(groups = { "Functional" })
public void printAllMatrices()
{
for (Map.Entry<String, ScoreModelI> sm : ResidueProperties.scoreMatrices
}
}
- @Test(groups = { "Functional" })
public void printHTMLMatrices()
{
for (Map.Entry<String, ScoreModelI> _sm : ResidueProperties.scoreMatrices
--- /dev/null
+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);
+ }
+}
git://source.jalview.org / jalview.git / commitdiff