src/jalview/math/Matrix.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.util.Format;
  24 import jalview.util.MessageManager;
  25
  26 import java.io.PrintStream;
  27 import java.lang.Math;
  28 import java.util.Arrays;
  29
  30 /**
  31  * A class to model rectangular matrices of double values and operations on them
  32  */
  33 public class Matrix implements MatrixI
  34 {
  35   /*
  36    * maximum number of iterations for tqli
  37    */
  38   private static final int MAX_ITER = 45;
  39   // fudge - add 15 iterations, just in case
  40
  41   /*
  42    * the number of rows
  43    */
  44   final protected int rows;
  45
  46   /*
  47    * the number of columns
  48    */
  49   final protected int cols;
  50
  51   /*
  52    * the cell values in row-major order
  53    */
  54   private double[][] value;
  55
  56   protected double[] d; // Diagonal
  57
  58   protected double[] e; // off diagonal
  59
  60   /**
  61    * Constructor given number of rows and columns
  62    *
  63    * @param colCount
  64    * @param rowCount
  65    */
  66   protected Matrix(int rowCount, int colCount)
  67   {
  68     rows = rowCount;
  69     cols = colCount;
  70   }
  71
  72   /**
  73    * Creates a new Matrix object containing a copy of the supplied array values.
  74    * For example
  75    *
  76    * <pre>
  77    *   new Matrix(new double[][] {{2, 3, 4}, {5, 6, 7})
  78    * constructs
  79    *   (2 3 4)
  80    *   (5 6 7)
  81    * </pre>
  82    *
  83    * Note that ragged arrays (with not all rows, or columns, of the same
  84    * length), are not supported by this class. They can be constructed, but
  85    * results of operations on them are undefined and may throw exceptions.
  86    *
  87    * @param values
  88    *          the matrix values in row-major order
  89    */
  90   public Matrix(double[][] values)
  91   {
  92     this.rows = values.length;
  93     this.cols = this.rows == 0 ? 0 : values[0].length;
  94
  95     /*
  96      * make a copy of the values array, for immutability
  97      */
  98     this.value = new double[rows][];
  99     int i = 0;
 100     for (double[] row : values)
 101     {
 102       if (row != null)
 103       {
 104         value[i] = new double[row.length];
 105         System.arraycopy(row, 0, value[i], 0, row.length);
 106       }
 107       i++;
 108     }
 109   }
 110
 111   public Matrix(float[][] values)
 112   {
 113     this.rows = values.length;
 114     this.cols = this.rows == 0 ? 0 : values[0].length;
 115
 116     /*
 117      * make a copy of the values array, for immutability
 118      */
 119     this.value = new double[rows][];
 120     int i = 0;
 121     for (float[] row : values)
 122     {
 123       if (row != null)
 124       {
 125         value[i] = new double[row.length];
 126         int j = 0;
 127         for (float oldValue : row)
 128         {
 129           value[i][j] = oldValue;
 130           j++;
 131         }
 132       }
 133       i++;
 134     }
 135   }
 136
 137   @Override
 138   public MatrixI transpose()
 139   {
 140     double[][] out = new double[cols][rows];
 141
 142     for (int i = 0; i < cols; i++)
 143     {
 144       for (int j = 0; j < rows; j++)
 145       {
 146         out[i][j] = value[j][i];
 147       }
 148     }
 149
 150     return new Matrix(out);
 151   }
 152
 153   /**
 154    * DOCUMENT ME!
 155    *
 156    * @param ps
 157    *          DOCUMENT ME!
 158    * @param format
 159    */
 160   @Override
 161   public void print(PrintStream ps, String format)
 162   {
 163     for (int i = 0; i < rows; i++)
 164     {
 165       for (int j = 0; j < cols; j++)
 166       {
 167         Format.print(ps, format, getValue(i, j));
 168       }
 169
 170       ps.println();
 171     }
 172   }
 173
 174   @Override
 175   public MatrixI preMultiply(MatrixI in)
 176   {
 177     if (in.width() != rows)
 178     {
 179       throw new IllegalArgumentException("Can't pre-multiply " + this.rows
 180               + " rows by " + in.width() + " columns");
 181     }
 182     double[][] tmp = new double[in.height()][this.cols];
 183
 184     for (int i = 0; i < in.height(); i++)
 185     {
 186       for (int j = 0; j < this.cols; j++)
 187       {
 188         /*
 189          * result[i][j] is the vector product of
 190          * in.row[i] and this.column[j]
 191          */
 192         for (int k = 0; k < in.width(); k++)
 193         {
 194           if (!Double.isNaN(in.getValue(i,k)) && !Double.isNaN(this.value[k][j]))
 195           {
 196             tmp[i][j] += (in.getValue(i, k) * this.value[k][j]);
 197           }
 198         }
 199       }
 200     }
 201
 202     return new Matrix(tmp);
 203   }
 204
 205   /**
 206    *
 207    * @param in
 208    *
 209    * @return
 210    */
 211   public double[] vectorPostMultiply(double[] in)
 212   {
 213     double[] out = new double[in.length];
 214
 215     for (int i = 0; i < in.length; i++)
 216     {
 217       out[i] = 0.0;
 218
 219       for (int k = 0; k < in.length; k++)
 220       {
 221         out[i] += (value[i][k] * in[k]);
 222       }
 223     }
 224
 225     return out;
 226   }
 227
 228   @Override
 229   public MatrixI postMultiply(MatrixI in)
 230   {
 231     if (in.height() != this.cols)
 232     {
 233       throw new IllegalArgumentException("Can't post-multiply " + this.cols
 234               + " columns by " + in.height() + " rows");
 235     }
 236     return in.preMultiply(this);
 237   }
 238
 239   @Override
 240   public MatrixI copy()
 241   {
 242     double[][] newmat = new double[rows][cols];
 243
 244     for (int i = 0; i < rows; i++)
 245     {
 246       System.arraycopy(value[i], 0, newmat[i], 0, value[i].length);
 247     }
 248
 249     Matrix m = new Matrix(newmat);
 250     if (this.d != null)
 251     {
 252       m.d = Arrays.copyOf(this.d, this.d.length);
 253     }
 254     if (this.e != null)
 255     {
 256       m.e = Arrays.copyOf(this.e, this.e.length);
 257     }
 258
 259     return m;
 260   }
 261
 262   /**
 263    * DOCUMENT ME!
 264    */
 265   @Override
 266   public void tred()
 267   {
 268     int n = rows;
 269     int k;
 270     int j;
 271     int i;
 272
 273     double scale;
 274     double hh;
 275     double h;
 276     double g;
 277     double f;
 278
 279     this.d = new double[rows];
 280     this.e = new double[rows];
 281
 282     for (i = n; i >= 2; i--)
 283     {
 284       final int l = i - 1;
 285       h = 0.0;
 286       scale = 0.0;
 287
 288       if (l > 1)
 289       {
 290         for (k = 1; k <= l; k++)
 291         {
 292           double v = Math.abs(getValue(i - 1, k - 1));
 293           scale += v;
 294         }
 295
 296         if (scale == 0.0)
 297         {
 298           e[i - 1] = getValue(i - 1, l - 1);
 299         }
 300         else
 301         {
 302           for (k = 1; k <= l; k++)
 303           {
 304             double v = divideValue(i - 1, k - 1, scale);
 305             h += v * v;
 306           }
 307
 308           f = getValue(i - 1, l - 1);
 309
 310           if (f > 0)
 311           {
 312             g = -1.0 * Math.sqrt(h);
 313           }
 314           else
 315           {
 316             g = Math.sqrt(h);
 317           }
 318
 319           e[i - 1] = scale * g;
 320           h -= (f * g);
 321           setValue(i - 1, l - 1, f - g);
 322           f = 0.0;
 323
 324           for (j = 1; j <= l; j++)
 325           {
 326             double val = getValue(i - 1, j - 1) / h;
 327             setValue(j - 1, i - 1, val);
 328             g = 0.0;
 329
 330             for (k = 1; k <= j; k++)
 331             {
 332               g += (getValue(j - 1, k - 1) * getValue(i - 1, k - 1));
 333             }
 334
 335             for (k = j + 1; k <= l; k++)
 336             {
 337               g += (getValue(k - 1, j - 1) * getValue(i - 1, k - 1));
 338             }
 339
 340             e[j - 1] = g / h;
 341             f += (e[j - 1] * getValue(i - 1, j - 1));
 342           }
 343
 344           hh = f / (h + h);
 345
 346           for (j = 1; j <= l; j++)
 347           {
 348             f = getValue(i - 1, j - 1);
 349             g = e[j - 1] - (hh * f);
 350             e[j - 1] = g;
 351
 352             for (k = 1; k <= j; k++)
 353             {
 354               double val = (f * e[k - 1]) + (g * getValue(i - 1, k - 1));
 355               addValue(j - 1, k - 1, -val);
 356             }
 357           }
 358         }
 359       }
 360       else
 361       {
 362         e[i - 1] = getValue(i - 1, l - 1);
 363       }
 364
 365       d[i - 1] = h;
 366     }
 367
 368     d[0] = 0.0;
 369     e[0] = 0.0;
 370
 371     for (i = 1; i <= n; i++)
 372     {
 373       final int l = i - 1;
 374
 375       if (d[i - 1] != 0.0)
 376       {
 377         for (j = 1; j <= l; j++)
 378         {
 379           g = 0.0;
 380
 381           for (k = 1; k <= l; k++)
 382           {
 383             g += (getValue(i - 1, k - 1) * getValue(k - 1, j - 1));
 384           }
 385
 386           for (k = 1; k <= l; k++)
 387           {
 388             addValue(k - 1, j - 1, -(g * getValue(k - 1, i - 1)));
 389           }
 390         }
 391       }
 392
 393       d[i - 1] = getValue(i - 1, i - 1);
 394       setValue(i - 1, i - 1, 1.0);
 395
 396       for (j = 1; j <= l; j++)
 397       {
 398         setValue(j - 1, i - 1, 0.0);
 399         setValue(i - 1, j - 1, 0.0);
 400       }
 401     }
 402   }
 403
 404   /**
 405    * Adds f to the value at [i, j] and returns the new value
 406    *
 407    * @param i
 408    * @param j
 409    * @param f
 410    */
 411   protected double addValue(int i, int j, double f)
 412   {
 413     double v = value[i][j] + f;
 414     value[i][j] = v;
 415     return v;
 416   }
 417
 418   /**
 419    * Divides the value at [i, j] by divisor and returns the new value. If d is
 420    * zero, returns the unchanged value.
 421    *
 422    * @param i
 423    * @param j
 424    * @param divisor
 425    * @return
 426    */
 427   protected double divideValue(int i, int j, double divisor)
 428   {
 429     if (divisor == 0d)
 430     {
 431       return getValue(i, j);
 432     }
 433     double v = value[i][j];
 434     v = v / divisor;
 435     value[i][j] = v;
 436     return v;
 437   }
 438
 439   /**
 440    * DOCUMENT ME!
 441    */
 442   @Override
 443   public void tqli() throws Exception
 444   {
 445     int n = rows;
 446
 447     int m;
 448     int l;
 449     int iter;
 450     int i;
 451     int k;
 452     double s;
 453     double r;
 454     double p;
 455
 456     double g;
 457     double f;
 458     double dd;
 459     double c;
 460     double b;
 461
 462     for (i = 2; i <= n; i++)
 463     {
 464       e[i - 2] = e[i - 1];
 465     }
 466
 467     e[n - 1] = 0.0;
 468
 469     for (l = 1; l <= n; l++)
 470     {
 471       iter = 0;
 472
 473       do
 474       {
 475         for (m = l; m <= (n - 1); m++)
 476         {
 477           dd = Math.abs(d[m - 1]) + Math.abs(d[m]);
 478
 479           if ((Math.abs(e[m - 1]) + dd) == dd)
 480           {
 481             break;
 482           }
 483         }
 484
 485         if (m != l)
 486         {
 487           iter++;
 488
 489           if (iter == MAX_ITER)
 490           {
 491             throw new Exception(MessageManager.formatMessage(
 492                     "exception.matrix_too_many_iteration", new String[]
 493                     { "tqli", Integer.valueOf(MAX_ITER).toString() }));
 494           }
 495           else
 496           {
 497             // jalview.bin.Console.outPrintln("Iteration " + iter);
 498           }
 499
 500           g = (d[l] - d[l - 1]) / (2.0 * e[l - 1]);
 501           r = Math.sqrt((g * g) + 1.0);
 502           g = d[m - 1] - d[l - 1] + (e[l - 1] / (g + sign(r, g)));
 503           c = 1.0;
 504           s = c;
 505           p = 0.0;
 506
 507           for (i = m - 1; i >= l; i--)
 508           {
 509             f = s * e[i - 1];
 510             b = c * e[i - 1];
 511
 512             if (Math.abs(f) >= Math.abs(g))
 513             {
 514               c = g / f;
 515               r = Math.sqrt((c * c) + 1.0);
 516               e[i] = f * r;
 517               s = 1.0 / r;
 518               c *= s;
 519             }
 520             else
 521             {
 522               s = f / g;
 523               r = Math.sqrt((s * s) + 1.0);
 524               e[i] = g * r;
 525               c = 1.0 / r;
 526               s *= c;
 527             }
 528
 529             g = d[i] - p;
 530             r = ((d[i - 1] - g) * s) + (2.0 * c * b);
 531             p = s * r;
 532             d[i] = g + p;
 533             g = (c * r) - b;
 534
 535             for (k = 1; k <= n; k++)
 536             {
 537               f = getValue(k - 1, i);
 538               setValue(k - 1, i, (s * getValue(k - 1, i - 1)) + (c * f));
 539               setValue(k - 1, i - 1,
 540                       (c * getValue(k - 1, i - 1)) - (s * f));
 541             }
 542           }
 543
 544           d[l - 1] = d[l - 1] - p;
 545           e[l - 1] = g;
 546           e[m - 1] = 0.0;
 547         }
 548       } while (m != l);
 549     }
 550   }
 551
 552   @Override
 553   public double getValue(int i, int j)
 554   {
 555     return value[i][j];
 556   }
 557
 558   @Override
 559   public void setValue(int i, int j, double val)
 560   {
 561     value[i][j] = val;
 562   }
 563
 564   /**
 565    * DOCUMENT ME!
 566    */
 567   public void tred2()
 568   {
 569     int n = rows;
 570     int l;
 571     int k;
 572     int j;
 573     int i;
 574
 575     double scale;
 576     double hh;
 577     double h;
 578     double g;
 579     double f;
 580
 581     this.d = new double[rows];
 582     this.e = new double[rows];
 583
 584     for (i = n - 1; i >= 1; i--)
 585     {
 586       l = i - 1;
 587       h = 0.0;
 588       scale = 0.0;
 589
 590       if (l > 0)
 591       {
 592         for (k = 0; k < l; k++)
 593         {
 594           scale += Math.abs(value[i][k]);
 595         }
 596
 597         if (scale == 0.0)
 598         {
 599           e[i] = value[i][l];
 600         }
 601         else
 602         {
 603           for (k = 0; k < l; k++)
 604           {
 605             value[i][k] /= scale;
 606             h += (value[i][k] * value[i][k]);
 607           }
 608
 609           f = value[i][l];
 610
 611           if (f > 0)
 612           {
 613             g = -1.0 * Math.sqrt(h);
 614           }
 615           else
 616           {
 617             g = Math.sqrt(h);
 618           }
 619
 620           e[i] = scale * g;
 621           h -= (f * g);
 622           value[i][l] = f - g;
 623           f = 0.0;
 624
 625           for (j = 0; j < l; j++)
 626           {
 627             value[j][i] = value[i][j] / h;
 628             g = 0.0;
 629
 630             for (k = 0; k < j; k++)
 631             {
 632               g += (value[j][k] * value[i][k]);
 633             }
 634
 635             for (k = j; k < l; k++)
 636             {
 637               g += (value[k][j] * value[i][k]);
 638             }
 639
 640             e[j] = g / h;
 641             f += (e[j] * value[i][j]);
 642           }
 643
 644           hh = f / (h + h);
 645
 646           for (j = 0; j < l; j++)
 647           {
 648             f = value[i][j];
 649             g = e[j] - (hh * f);
 650             e[j] = g;
 651
 652             for (k = 0; k < j; k++)
 653             {
 654               value[j][k] -= ((f * e[k]) + (g * value[i][k]));
 655             }
 656           }
 657         }
 658       }
 659       else
 660       {
 661         e[i] = value[i][l];
 662       }
 663
 664       d[i] = h;
 665     }
 666
 667     d[0] = 0.0;
 668     e[0] = 0.0;
 669
 670     for (i = 0; i < n; i++)
 671     {
 672       l = i - 1;
 673
 674       if (d[i] != 0.0)
 675       {
 676         for (j = 0; j < l; j++)
 677         {
 678           g = 0.0;
 679
 680           for (k = 0; k < l; k++)
 681           {
 682             g += (value[i][k] * value[k][j]);
 683           }
 684
 685           for (k = 0; k < l; k++)
 686           {
 687             value[k][j] -= (g * value[k][i]);
 688           }
 689         }
 690       }
 691
 692       d[i] = value[i][i];
 693       value[i][i] = 1.0;
 694
 695       for (j = 0; j < l; j++)
 696       {
 697         value[j][i] = 0.0;
 698         value[i][j] = 0.0;
 699       }
 700     }
 701   }
 702
 703   /**
 704    * DOCUMENT ME!
 705    */
 706   public void tqli2() throws Exception
 707   {
 708     int n = rows;
 709
 710     int m;
 711     int l;
 712     int iter;
 713     int i;
 714     int k;
 715     double s;
 716     double r;
 717     double p;
 718     ;
 719
 720     double g;
 721     double f;
 722     double dd;
 723     double c;
 724     double b;
 725
 726     for (i = 2; i <= n; i++)
 727     {
 728       e[i - 2] = e[i - 1];
 729     }
 730
 731     e[n - 1] = 0.0;
 732
 733     for (l = 1; l <= n; l++)
 734     {
 735       iter = 0;
 736
 737       do
 738       {
 739         for (m = l; m <= (n - 1); m++)
 740         {
 741           dd = Math.abs(d[m - 1]) + Math.abs(d[m]);
 742
 743           if ((Math.abs(e[m - 1]) + dd) == dd)
 744           {
 745             break;
 746           }
 747         }
 748
 749         if (m != l)
 750         {
 751           iter++;
 752
 753           if (iter == MAX_ITER)
 754           {
 755             throw new Exception(MessageManager.formatMessage(
 756                     "exception.matrix_too_many_iteration", new String[]
 757                     { "tqli2", Integer.valueOf(MAX_ITER).toString() }));
 758           }
 759           else
 760           {
 761             // jalview.bin.Console.outPrintln("Iteration " + iter);
 762           }
 763
 764           g = (d[l] - d[l - 1]) / (2.0 * e[l - 1]);
 765           r = Math.sqrt((g * g) + 1.0);
 766           g = d[m - 1] - d[l - 1] + (e[l - 1] / (g + sign(r, g)));
 767           c = 1.0;
 768           s = c;
 769           p = 0.0;
 770
 771           for (i = m - 1; i >= l; i--)
 772           {
 773             f = s * e[i - 1];
 774             b = c * e[i - 1];
 775
 776             if (Math.abs(f) >= Math.abs(g))
 777             {
 778               c = g / f;
 779               r = Math.sqrt((c * c) + 1.0);
 780               e[i] = f * r;
 781               s = 1.0 / r;
 782               c *= s;
 783             }
 784             else
 785             {
 786               s = f / g;
 787               r = Math.sqrt((s * s) + 1.0);
 788               e[i] = g * r;
 789               c = 1.0 / r;
 790               s *= c;
 791             }
 792
 793             g = d[i] - p;
 794             r = ((d[i - 1] - g) * s) + (2.0 * c * b);
 795             p = s * r;
 796             d[i] = g + p;
 797             g = (c * r) - b;
 798
 799             for (k = 1; k <= n; k++)
 800             {
 801               f = value[k - 1][i];
 802               value[k - 1][i] = (s * value[k - 1][i - 1]) + (c * f);
 803               value[k - 1][i - 1] = (c * value[k - 1][i - 1]) - (s * f);
 804             }
 805           }
 806
 807           d[l - 1] = d[l - 1] - p;
 808           e[l - 1] = g;
 809           e[m - 1] = 0.0;
 810         }
 811       } while (m != l);
 812     }
 813   }
 814
 815   /**
 816    * Answers the first argument with the sign of the second argument
 817    *
 818    * @param a
 819    * @param b
 820    *
 821    * @return
 822    */
 823   static double sign(double a, double b)
 824   {
 825     if (b < 0)
 826     {
 827       return -Math.abs(a);
 828     }
 829     else
 830     {
 831       return Math.abs(a);
 832     }
 833   }
 834
 835   /**
 836   * returns the matrix as a double[][] array
 837   *
 838   * @return
 839   */
 840   @Override
 841   public double[][] asArray()
 842   {
 843     return value;
 844   }
 845
 846   /**
 847    * Returns an array containing the values in the specified column
 848    *
 849    * @param col
 850    *
 851    * @return
 852    */
 853   @Override
 854   public double[] getColumn(int col)
 855   {
 856     double[] out = new double[rows];
 857
 858     for (int i = 0; i < rows; i++)
 859     {
 860       out[i] = value[i][col];
 861     }
 862
 863     return out;
 864   }
 865
 866   /**
 867    * DOCUMENT ME!
 868    *
 869    * @param ps
 870    *          DOCUMENT ME!
 871    * @param format
 872    */
 873   @Override
 874   public void printD(PrintStream ps, String format)
 875   {
 876     for (int j = 0; j < d.length; j++)
 877     {
 878       Format.print(ps, format, d[j]);
 879     }
 880   }
 881
 882   /**
 883    * DOCUMENT ME!
 884    *
 885    * @param ps
 886    *          DOCUMENT ME!
 887    * @param format
 888    *          TODO
 889    */
 890   @Override
 891   public void printE(PrintStream ps, String format)
 892   {
 893     for (int j = 0; j < rows; j++)
 894     {
 895       Format.print(ps, format, e[j]);
 896     }
 897   }
 898
 899   @Override
 900   public double[] getD()
 901   {
 902     return d;
 903   }
 904
 905   @Override
 906   public double[] getE()
 907   {
 908     return e;
 909   }
 910
 911   @Override
 912   public int height()
 913   {
 914     return rows;
 915   }
 916
 917   @Override
 918   public int width()
 919   {
 920     return cols;
 921   }
 922
 923   @Override
 924   public double[] getRow(int i)
 925   {
 926     double[] row = new double[cols];
 927     System.arraycopy(value[i], 0, row, 0, cols);
 928     return row;
 929   }
 930
 931   /**
 932    * Returns a length 2 array of {minValue, maxValue} of all values in the
 933    * matrix. Returns null if the matrix is null or empty.
 934    *
 935    * @return
 936    */
 937   double[] findMinMax()
 938   {
 939     if (value == null)
 940     {
 941       return null;
 942     }
 943     double min = Double.MAX_VALUE;
 944     double max = -Double.MAX_VALUE;
 945     boolean empty = true;
 946     for (double[] row : value)
 947     {
 948       if (row != null)
 949       {
 950         for (double x : row)
 951         {
 952           empty = false;
 953           if (x > max)
 954           {
 955             max = x;
 956           }
 957           if (x < min)
 958           {
 959             min = x;
 960           }
 961         }
 962       }
 963     }
 964     return empty ? null : new double[] { min, max };
 965   }
 966
 967   /**
 968    * {@inheritDoc}
 969    */
 970   @Override
 971   public void reverseRange(boolean maxToZero)
 972   {
 973     if (value == null)
 974     {
 975       return;
 976     }
 977     double[] minMax = findMinMax();
 978     if (minMax == null)
 979     {
 980       return; // empty matrix
 981     }
 982     double subtractFrom = maxToZero ? minMax[1] : minMax[0] + minMax[1];
 983
 984     for (double[] row : value)
 985     {
 986       if (row != null)
 987       {
 988         int j = 0;
 989         for (double x : row)
 990         {
 991           row[j] = subtractFrom - x;
 992           j++;
 993         }
 994       }
 995     }
 996   }
 997
 998   /**
 999    * Multiplies every entry in the matrix by the given value.
1000    *
1001    * @param
1002    */
1003   @Override
1004   public void multiply(double by)
1005   {
1006     for (double[] row : value)
1007     {
1008       if (row != null)
1009       {
1010         for (int i = 0; i < row.length; i++)
1011         {
1012           row[i] *= by;
1013         }
1014       }
1015     }
1016   }
1017
1018   /**
1019    * Add d to all entries of this matrix
1020    *
1021    * @param d ~ value to add
1022    */
1023   @Override
1024   public void add(double d)
1025   {
1026     for (double[] row : value)
1027     {
1028       if (row != null)
1029       {
1030         for (int i = 0; i < row.length; i++)
1031         {
1032           row[i] += d;
1033         }
1034       }
1035     }
1036   }
1037
1038   @Override
1039   public void setD(double[] v)
1040   {
1041     d = v;
1042   }
1043
1044   @Override
1045   public void setE(double[] v)
1046   {
1047     e = v;
1048   }
1049
1050   public double getTotal()
1051   {
1052     double d = 0d;
1053     for (int i = 0; i < this.height(); i++)
1054     {
1055       for (int j = 0; j < this.width(); j++)
1056       {
1057         d += value[i][j];
1058       }
1059     }
1060     return d;
1061   }
1062
1063   /**
1064    * {@inheritDoc}
1065    */
1066   @Override
1067   public boolean equals(MatrixI m2, double delta)
1068   {
1069     if (m2 == null || this.height() != m2.height()
1070             || this.width() != m2.width())
1071     {
1072       return false;
1073     }
1074     for (int i = 0; i < this.height(); i++)
1075     {
1076       for (int j = 0; j < this.width(); j++)
1077       {
1078         double diff = this.getValue(i, j) - m2.getValue(i, j);
1079         if (Math.abs(diff) > delta)
1080         {
1081           return false;
1082         }
1083       }
1084     }
1085     return true;
1086   }
1087
1088   /**
1089    * Returns a copy in which  every value in the matrix is its absolute
1090    *
1091    * @return
1092    */
1093   @Override
1094   public MatrixI absolute()
1095   {
1096     MatrixI copy = this.copy();
1097     for (int i = 0; i < copy.width(); i++)
1098     {
1099       double[] row = copy.getRow(i);
1100       if (row != null)
1101       {
1102         for (int j = 0; j < row.length; j++)
1103         {
1104           row[j] = Math.abs(row[j]);
1105         }
1106       }
1107     }
1108     return copy;
1109   }
1110
1111   /**
1112    * Returns the mean of each row
1113    *
1114    * @return
1115    */
1116   @Override
1117   public double[] meanRow()
1118   {
1119     double[] mean = new double[rows];
1120     int i = 0;
1121     for (double[] row : value)
1122     {
1123       if (row != null)
1124       {
1125         mean[i++] = MiscMath.mean(row);
1126       }
1127     }
1128     return mean;
1129   }
1130
1131   /**
1132    * Returns the mean of each column
1133    *
1134    * @return
1135    */
1136   @Override
1137   public double[] meanCol()
1138   {
1139     double[] mean = new double[cols];
1140     for (int j = 0; j < cols; j++)
1141     {
1142       double[] column = getColumn(j);
1143       if (column != null)
1144       {
1145         mean[j] = MiscMath.mean(column);
1146       }
1147     }
1148     return mean;
1149   }
1150
1151   /**
1152   * return a flattened matrix containing the sum of each column
1153   *
1154   * @return
1155   */
1156   @Override
1157   public double[] sumCol()
1158   {
1159     double[] sum = new double[cols];
1160     for (int j = 0; j < cols; j++)
1161     {
1162       double[] column = getColumn(j);
1163       if (column != null)
1164       {
1165         sum[j] = MiscMath.sum(column);
1166       }
1167     }
1168     return sum;
1169   }
1170
1171   /**
1172   * returns the mean value of the complete matrix
1173   *
1174   * @return
1175   */
1176   @Override
1177   public double mean()
1178   {
1179     double sum = 0;
1180     int nanCount = 0;
1181     for (double[] row : value)
1182     {
1183       for (double col : row)
1184       {
1185         if (!Double.isNaN(col))
1186         {
1187           sum += col;
1188         } else {
1189           nanCount++;
1190         }
1191       }
1192     }
1193     return sum / (double) (this.rows * this.cols - nanCount);
1194   }
1195
1196   /**
1197   * fills up a diagonal matrix with its transposed copy
1198   * !other side should be filled with 0
1199   * !keeps Double.NaN found in either side
1200   *
1201   * TODO check on which side it was diagonal and only do calculations for the other side
1202   */
1203   @Override
1204   public void fillDiagonal()
1205   {
1206     int n = this.rows;
1207     int m = this.cols;
1208     MatrixI copy = this.transpose();    // goes through each element in the matrix and
1209     for (int i = 0; i < n; i++)         // adds the value in the transposed copy to the original value
1210     {
1211       for (int j = 0; j < m; j++)
1212       {
1213         if (i != j)
1214         {
1215           this.addValue(i, j, copy.getValue(i,j));
1216         }
1217       }
1218     }
1219   }
1220
1221   /**
1222   * counts the number of Double.NaN in the matrix
1223   *
1224   * @return
1225   */
1226   @Override
1227   public int countNaN()
1228   {
1229     int NaN = 0;
1230     for (int i = 0; i < this.rows; i++)
1231     {
1232       for (int j = 0; j < this.cols; j++)
1233       {
1234         if (Double.isNaN(this.getValue(i,j)))
1235         {
1236           NaN++;
1237         }
1238       }
1239     }
1240     return NaN;
1241   }
1242
1243   /**
1244   * performs an element-wise addition of this matrix by another matrix ~ this - m
1245   * @param m ~ other matrix
1246   *
1247   * @return
1248   */
1249   @Override
1250   public MatrixI add(MatrixI m)
1251   {
1252     if (m.width() != cols || m.height() != rows)
1253     {
1254       throw new IllegalArgumentException("Can't add a " + m.height() + "x" + m.width() + " to a " + this.rows + "x" + this.cols + " matrix");
1255     }
1256     double[][] tmp = new double[this.rows][this.cols];
1257     for (int i = 0; i < this.rows; i++)
1258     {
1259       for (int j = 0; j < this.cols; j++)
1260       {
1261         tmp[i][j] = this.getValue(i,j) + m.getValue(i,j);
1262       }
1263     }
1264     return new Matrix(tmp);
1265   }
1266
1267   /**
1268   * performs an element-wise subtraction of this matrix by another matrix ~ this - m
1269   * @param m ~ other matrix
1270   *
1271   * @return
1272   */
1273   @Override
1274   public MatrixI subtract(MatrixI m)
1275   {
1276     if (m.width() != cols || m.height() != rows)
1277     {
1278       throw new IllegalArgumentException("Can't subtract a " + m.height() + "x" + m.width() + " from a " + this.rows + "x" + this.cols + " matrix");
1279     }
1280     double[][] tmp = new double[this.rows][this.cols];
1281     for (int i = 0; i < this.rows; i++)
1282     {
1283       for (int j = 0; j < this.cols; j++)
1284       {
1285         tmp[i][j] = this.getValue(i,j) -  m.getValue(i,j);
1286       }
1287     }
1288     return new Matrix(tmp);
1289   }
1290
1291   /**
1292   * performs an element-wise multiplication of this matrix by another matrix ~ this * m
1293   * @param m ~ other matrix
1294   *
1295   * @return
1296   */
1297   @Override
1298   public MatrixI elementwiseMultiply(MatrixI m)
1299   {
1300     if (m.width() != cols || m.height() != rows)
1301     {
1302       throw new IllegalArgumentException("Can't multiply a " + this.rows + "x" + this.cols + " by a " + m.height() + "x" + m.width() + " matrix");
1303     }
1304     double[][] tmp = new double[this.rows][this.cols];
1305     for (int i = 0; i < this.rows; i++)
1306     {
1307       for (int j = 0; j < this.cols; j++)
1308       {
1309         tmp[i][j] = this.getValue(i, j) * m.getValue(i,j);
1310       }
1311     }
1312     return new Matrix(tmp);
1313   }
1314
1315   /**
1316   * performs an element-wise division of this matrix by another matrix ~ this / m
1317   * @param m ~ other matrix
1318   *
1319   * @return
1320   */
1321   @Override
1322   public MatrixI elementwiseDivide(MatrixI m)
1323   {
1324     if (m.width() != cols || m.height() != rows)
1325     {
1326       throw new IllegalArgumentException("Can't divide a " + this.rows + "x" + this.cols + " by a " + m.height() + "x" + m.width() + " matrix");
1327     }
1328     double[][] tmp = new double[this.rows][this.cols];
1329     for (int i = 0; i < this.rows; i++)
1330     {
1331       for (int j = 0; j < this.cols; j++)
1332       {
1333         tmp[i][j] = this.getValue(i, j) / m.getValue(i,j);
1334       }
1335     }
1336     return new Matrix(tmp);
1337   }
1338
1339   /**
1340   * calculate the root-mean-square for tow matrices
1341   * @param m ~ other matrix
1342   *
1343   * @return
1344   */
1345   @Override
1346   public double rmsd(MatrixI m)
1347   {
1348     MatrixI squaredDeviates = this.subtract(m);
1349     squaredDeviates = squaredDeviates.preMultiply(squaredDeviates);
1350     return Math.sqrt(squaredDeviates.mean());
1351   }
1352
1353   /**
1354   * calculates the Frobenius norm of this matrix
1355   *
1356   * @return
1357   */
1358   @Override
1359   public double norm()
1360   {
1361     double result = 0;
1362     for (double[] row : value)
1363     {
1364       for (double val : row)
1365       {
1366         result += Math.pow(val, 2);
1367       }
1368     }
1369     return Math.sqrt(result);
1370   }
1371
1372   /**
1373   * returns the sum of all values in this matrix
1374   *
1375   * @return
1376   */
1377   @Override
1378   public double sum()
1379   {
1380     double sum = 0;
1381     for (double[] row : value)
1382     {
1383       for (double val : row)
1384       {
1385         sum += (Double.isNaN(val)) ? 0.0 : val;
1386       }
1387     }
1388     return sum;
1389   }
1390
1391   /**
1392   * returns the sum-product of this matrix with vector v
1393   * @param v ~ vector
1394   *
1395   * @return
1396   */
1397   @Override
1398   public double[] sumProduct(double[] v)
1399   {
1400     if (v.length != cols)
1401     {
1402       throw new IllegalArgumentException("Vector and matrix do not have the same dimension! (" + v.length + " != " + cols + ")");
1403     }
1404     double[] result = new double[rows];
1405     for (int i = 0; i < rows; i++)
1406     {
1407       double[] row = value[i];
1408       double sum = 0;
1409       for (int j = 0; j < row.length; j++)
1410       {
1411         sum += row[j] * v[j];
1412       }
1413       result[i] = sum;
1414     }
1415     return result;
1416   }
1417
1418   /**
1419   * mirrors columns of the matrix
1420   *
1421   * @return
1422   */
1423   @Override
1424   public MatrixI mirrorCol()
1425   {
1426     double[][] result = new double[rows][cols];
1427     for (int i = 0; i < rows; i++)
1428     {
1429       int k = cols - 1; // reverse col
1430       for (int j = 0; j < cols; j++)
1431       {
1432         result[i][k--] = this.getValue(i,j);
1433       }
1434     }
1435     MatrixI resultMatrix = new Matrix(result);
1436     if (d != null)
1437       resultMatrix.setD(d);
1438     if (e != null)
1439       resultMatrix.setE(e);
1440
1441     return resultMatrix;
1442   }
1443 }