src/jalview/math/SparseMatrix.java

   1 /*
   2  * Jalview - A Sequence Alignment Editor and Viewer ($$Version-Rel$$)
   3  * Copyright (C) $$Year-Rel$$ The Jalview Authors
   4  *
   5  * This file is part of Jalview.
   6  *
   7  * Jalview is free software: you can redistribute it and/or
   8  * modify it under the terms of the GNU General Public License
   9  * as published by the Free Software Foundation, either version 3
  10  * of the License, or (at your option) any later version.
  11  *
  12  * Jalview is distributed in the hope that it will be useful, but
  13  * WITHOUT ANY WARRANTY; without even the implied warranty
  14  * of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  15  * PURPOSE.  See the GNU General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU General Public License
  18  * along with Jalview.  If not, see <http://www.gnu.org/licenses/>.
  19  * The Jalview Authors are detailed in the 'AUTHORS' file.
  20  */
  21 package jalview.math;
  22
  23 import jalview.ext.android.SparseDoubleArray;
  24
  25 /**
  26  * A variant of Matrix intended for use for sparse (mostly zero) matrices. This
  27  * class uses a SparseDoubleArray to hold each row of the matrix. The sparse
  28  * array only stores non-zero values. This gives a smaller memory footprint, and
  29  * fewer matrix calculation operations, for mostly zero matrices.
  30  *
  31  * @author gmcarstairs
  32  */
  33 public class SparseMatrix extends Matrix
  34 {
  35   /*
  36    * we choose columns for the sparse arrays as this allows
  37    * optimisation of the preMultiply() method used in PCA.run()
  38    */
  39   SparseDoubleArray[] sparseColumns;
  40
  41   /**
  42    * Constructor given data in [row][column] order
  43    *
  44    * @param v
  45    */
  46   public SparseMatrix(double[][] v)
  47   {
  48     super(v.length, v.length > 0 ? v[0].length : 0);
  49
  50     sparseColumns = new SparseDoubleArray[cols];
  51
  52     /*
  53      * transpose v[row][col] into [col][row] order
  54      */
  55     for (int col = 0; col < cols; col++)
  56     {
  57       SparseDoubleArray sparseColumn = new SparseDoubleArray();
  58       sparseColumns[col] = sparseColumn;
  59       for (int row = 0; row < rows; row++)
  60       {
  61         double value = v[row][col];
  62         if (value != 0d)
  63         {
  64           sparseColumn.put(row, value);
  65         }
  66       }
  67     }
  68   }
  69
  70   /**
  71    * Answers the value at row i, column j
  72    */
  73   @Override
  74   public double getValue(int i, int j)
  75   {
  76     return sparseColumns[j].get(i);
  77   }
  78
  79   /**
  80    * Sets the value at row i, column j to val
  81    */
  82   @Override
  83   public void setValue(int i, int j, double val)
  84   {
  85     if (val == 0d)
  86     {
  87       sparseColumns[j].delete(i);
  88     }
  89     else
  90     {
  91       sparseColumns[j].put(i, val);
  92     }
  93   }
  94
  95   @Override
  96   public double[] getColumn(int i)
  97   {
  98     double[] col = new double[height()];
  99
 100     SparseDoubleArray vals = sparseColumns[i];
 101     for (int nonZero = 0; nonZero < vals.size(); nonZero++)
 102     {
 103       col[vals.keyAt(nonZero)] = vals.valueAt(nonZero);
 104     }
 105     return col;
 106   }
 107
 108   @Override
 109   public MatrixI copy()
 110   {
 111     double[][] vals = new double[height()][width()];
 112     for (int i = 0; i < height(); i++)
 113     {
 114       vals[i] = getRow(i);
 115     }
 116     return new SparseMatrix(vals);
 117   }
 118
 119   @Override
 120   public MatrixI transpose()
 121   {
 122     double[][] out = new double[cols][rows];
 123
 124     /*
 125      * for each column...
 126      */
 127     for (int i = 0; i < cols; i++)
 128     {
 129       /*
 130        * put non-zero values into the corresponding row
 131        * of the transposed matrix
 132        */
 133       SparseDoubleArray vals = sparseColumns[i];
 134       for (int nonZero = 0; nonZero < vals.size(); nonZero++)
 135       {
 136         out[i][vals.keyAt(nonZero)] = vals.valueAt(nonZero);
 137       }
 138     }
 139
 140     return new SparseMatrix(out);
 141   }
 142
 143   /**
 144    * Answers a new matrix which is the product in.this. If the product contains
 145    * less than 20% non-zero values, it is returned as a SparseMatrix, else as a
 146    * Matrix.
 147    * <p>
 148    * This method is optimised for the sparse arrays which store column values
 149    * for a SparseMatrix. Note that postMultiply is not so optimised. That would
 150    * require redundantly also storing sparse arrays for the rows, which has not
 151    * been done. Currently only preMultiply is used in Jalview.
 152    */
 153   @Override
 154   public MatrixI preMultiply(MatrixI in)
 155   {
 156     if (in.width() != rows)
 157     {
 158       throw new IllegalArgumentException("Can't pre-multiply " + this.rows
 159               + " rows by " + in.width() + " columns");
 160     }
 161     double[][] tmp = new double[in.height()][this.cols];
 162
 163     long count = 0L;
 164     for (int i = 0; i < in.height(); i++)
 165     {
 166       for (int j = 0; j < this.cols; j++)
 167       {
 168         /*
 169          * result[i][j] is the vector product of
 170          * in.row[i] and this.column[j]
 171          * we only need to use non-zero values from the column
 172          */
 173         SparseDoubleArray vals = sparseColumns[j];
 174         boolean added = false;
 175         for (int nonZero = 0; nonZero < vals.size(); nonZero++)
 176         {
 177           int myRow = vals.keyAt(nonZero);
 178           double myValue = vals.valueAt(nonZero);
 179           tmp[i][j] += (in.getValue(i, myRow) * myValue);
 180           added = true;
 181         }
 182         if (added && tmp[i][j] != 0d)
 183         {
 184           count++; // non-zero entry in product
 185         }
 186       }
 187     }
 188
 189     /*
 190      * heuristic rule - if product is more than 80% zero
 191      * then construct a SparseMatrix, else a Matrix
 192      */
 193     if (count * 5 < in.height() * cols)
 194     {
 195       return new SparseMatrix(tmp);
 196     }
 197     else
 198     {
 199       return new Matrix(tmp);
 200     }
 201   }
 202
 203   @Override
 204   protected double divideValue(int i, int j, double divisor)
 205   {
 206     if (divisor == 0d)
 207     {
 208       return getValue(i, j);
 209     }
 210     double v = sparseColumns[j].divide(i, divisor);
 211     return v;
 212   }
 213
 214   @Override
 215   protected double addValue(int i, int j, double addend)
 216   {
 217     double v = sparseColumns[j].add(i, addend);
 218     return v;
 219   }
 220
 221   /**
 222    * Returns the fraction of the whole matrix size that is actually modelled in
 223    * sparse arrays (normally, the non-zero values)
 224    *
 225    * @return
 226    */
 227   public float getFillRatio()
 228   {
 229     long count = 0L;
 230     for (SparseDoubleArray col : sparseColumns)
 231     {
 232       count += col.size();
 233     }
 234     return count / (float) (height() * width());
 235   }
 236 }