001/** 002 * Copyright (c) 2011, The University of Southampton and the individual contributors. 003 * All rights reserved. 004 * 005 * Redistribution and use in source and binary forms, with or without modification, 006 * are permitted provided that the following conditions are met: 007 * 008 * * Redistributions of source code must retain the above copyright notice, 009 * this list of conditions and the following disclaimer. 010 * 011 * * Redistributions in binary form must reproduce the above copyright notice, 012 * this list of conditions and the following disclaimer in the documentation 013 * and/or other materials provided with the distribution. 014 * 015 * * Neither the name of the University of Southampton nor the names of its 016 * contributors may be used to endorse or promote products derived from this 017 * software without specific prior written permission. 018 * 019 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND 020 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 021 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 022 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR 023 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES 024 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; 025 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON 026 * ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 027 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 028 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 029 */ 030package org.openimaj.math.matrix; 031 032import gnu.trove.TIntCollection; 033import gnu.trove.list.array.TIntArrayList; 034import gnu.trove.map.hash.TIntIntHashMap; 035import gnu.trove.procedure.TIntProcedure; 036import gov.sandia.cognition.math.matrix.MatrixEntry; 037import no.uib.cipr.matrix.sparse.FlexCompRowMatrix; 038 039import org.openimaj.util.array.ArrayUtils; 040import org.openimaj.util.array.SparseBinSearchDoubleArray; 041import org.openimaj.util.array.SparseDoubleArray; 042 043import ch.akuhn.matrix.DenseMatrix; 044import ch.akuhn.matrix.DenseVector; 045import ch.akuhn.matrix.Matrix; 046import ch.akuhn.matrix.SparseMatrix; 047import ch.akuhn.matrix.SparseVector; 048import ch.akuhn.matrix.Vector; 049import ch.akuhn.matrix.Vector.Entry; 050 051import com.jmatio.types.MLArray; 052import com.jmatio.types.MLDouble; 053 054/** 055 * Some helpful operations on {@link Matrix} instances from Adrian Kuhn's 056 * library. 057 * 058 * @author Sina Samangooei (ss@ecs.soton.ac.uk) 059 */ 060public class MatlibMatrixUtils { 061 private static final double EPS = 1e-8; 062 063 /** 064 * Compute the matrix sparsity (i.e. proportion of zero elements) 065 * 066 * @param mat 067 * the matrix 068 * @return the sparsity 069 * @see SparseMatrix#density() 070 */ 071 public static double sparsity(Matrix mat) { 072 return 1 - mat.density(); 073 } 074 075 /** 076 * Raise each element to the power d, operating on the matrix itself 077 * 078 * @param matrix 079 * the matrix 080 * @param d 081 * the power 082 * @return the input matrix, with elements raised to the given power 083 */ 084 public static <T extends Matrix> T powInplace(T matrix, double d) { 085 int rowN = 0; 086 for (final Vector ent : matrix.rows()) { 087 for (final Entry row : ent.entries()) { 088 matrix.put(rowN, row.index, Math.pow(row.value, d)); 089 } 090 rowN++; 091 } 092 return matrix; 093 } 094 095 /** 096 * Left multiply two matrices: <code>R = D . A</code> 097 * 098 * @param D 099 * first matrix 100 * @param A 101 * second matrix 102 * @return result of multiplication 103 */ 104 public static SparseMatrix times(DiagonalMatrix D, SparseMatrix A) { 105 final SparseMatrix mat = new SparseMatrix(A.rowCount(), A.columnCount()); 106 final double[] Dvals = D.getVals(); 107 int rowIndex = 0; 108 109 for (final Vector row : A.rows()) { 110 for (final Entry ent : row.entries()) { 111 mat.put(rowIndex, ent.index, ent.value * Dvals[rowIndex]); 112 } 113 rowIndex++; 114 } 115 116 return mat; 117 } 118 119 /** 120 * Right multiply two matrices: <code>R = A . D</code> 121 * 122 * @param D 123 * first matrix 124 * @param A 125 * second matrix 126 * @return result of multiplication 127 */ 128 public static SparseMatrix times(SparseMatrix A, DiagonalMatrix D) { 129 final SparseMatrix mat = new SparseMatrix(A.rowCount(), A.columnCount()); 130 int rowIndex = 0; 131 final double[] Dvals = D.getVals(); 132 for (final Vector row : A.rows()) { 133 for (final Entry ent : row.entries()) { 134 mat.put(rowIndex, ent.index, ent.value * Dvals[ent.index]); 135 } 136 rowIndex++; 137 } 138 return mat; 139 } 140 141 /** 142 * Add two matrices, storing the results in the first: 143 * <code>A = A + B</code> 144 * 145 * @param A 146 * first matrix 147 * @param B 148 * matrix to add 149 * @return A first matrix 150 */ 151 public static SparseMatrix plusInplace(SparseMatrix A, Matrix B) { 152 for (int i = 0; i < A.rowCount(); i++) { 153 A.addToRow(i, B.row(i)); 154 } 155 return A; 156 } 157 158 /** 159 * Add two matrices, storing the results in the first: 160 * <code>A = A + B</code> 161 * 162 * @param A 163 * first matrix 164 * @param B 165 * matrix to add 166 * @return A first matrix 167 */ 168 @SuppressWarnings("unchecked") 169 public static <T extends Matrix> T plusInplace(T A, Matrix B) { 170 if (A instanceof SparseMatrix) 171 return (T) plusInplace((SparseMatrix) A, B); 172 for (int i = 0; i < A.rowCount(); i++) { 173 final Vector brow = B.row(i); 174 for (int j = 0; j < A.columnCount(); j++) { 175 A.row(i).add(j, brow.get(j)); 176 } 177 } 178 return A; 179 } 180 181 /** 182 * Add a constant inplace <code>A = A + d</code> 183 * 184 * @param A 185 * first matrix 186 * @param d 187 * the constant to add 188 * @return A first matrix 189 */ 190 public static <T extends Matrix> T plusInplace(T A, double d) { 191 for (int i = 0; i < A.rowCount(); i++) { 192 for (int j = 0; j < A.columnCount(); j++) { 193 A.row(i).add(j, d); 194 } 195 } 196 return A; 197 } 198 199 /** 200 * Subtract two matrices, storing the result in the second: 201 * <code>A = D - A</code> 202 * 203 * @param D 204 * first matrix 205 * @param A 206 * second matrix 207 * @return second matrix 208 * 209 */ 210 public static <T extends Matrix> T minusInplace(DiagonalMatrix D, T A) { 211 final double[] Dval = D.getVals(); 212 for (int i = 0; i < Dval.length; i++) { 213 final Iterable<Entry> rowents = A.row(i).entries(); 214 for (final Entry entry : rowents) { 215 A.put(i, entry.index, -entry.value); 216 } 217 A.put(i, i, Dval[i] - A.get(i, i)); 218 } 219 return A; 220 } 221 222 /** 223 * Subtract two matrices, storing the result in the first: 224 * <code>A = A - B</code> 225 * 226 * @param A 227 * first matrix 228 * @param B 229 * second matrix 230 * @return first matrix 231 * 232 */ 233 public static Matrix minusInplace(Matrix A, Matrix B) { 234 for (int i = 0; i < A.rowCount(); i++) { 235 final Iterable<Entry> rowents = A.row(i).entries(); 236 for (final Entry entry : rowents) { 237 A.put(i, entry.index, entry.value - B.get(i, entry.index)); 238 } 239 } 240 return A; 241 } 242 243 /** 244 * Subtract a vector from another vector <code>A = A - D</code> 245 * 246 * @param A 247 * first matrix 248 * @param D 249 * second matrix 250 * @return second matrix 251 * 252 */ 253 public static <T extends Vector> T minusInplace(T A, Vector D) { 254 for (int i = 0; i < A.size(); i++) { 255 A.put(i, A.get(i) - D.get(i)); 256 } 257 return A; 258 } 259 260 /** 261 * Subtract a vector from another vector <code>A = A + D</code> 262 * 263 * @param A 264 * first matrix 265 * @param D 266 * second matrix 267 * @return second matrix 268 * 269 */ 270 public static <T extends Vector> T plusInplace(T A, Vector D) { 271 for (int i = 0; i < A.size(); i++) { 272 A.put(i, A.get(i) + D.get(i)); 273 } 274 return A; 275 } 276 277 /** 278 * Add two matrices, storing the results in the second: 279 * <code>A = D + A</code> 280 * 281 * @param D 282 * first matrix 283 * @param A 284 * second matrix 285 * @return second matrix 286 * 287 */ 288 public static <T extends Matrix> T plusInplace(DiagonalMatrix D, T A) { 289 final double[] Dval = D.getVals(); 290 for (int i = 0; i < Dval.length; i++) { 291 A.put(i, i, A.get(i, i) + Dval[i]); 292 } 293 return A; 294 } 295 296 /** 297 * Compute A^T . B 298 * 299 * @param A 300 * @param B 301 * @return A^T . B 302 */ 303 public static Matrix transposeDotProduct(Matrix A, Matrix B) { 304 final int mA = A.columnCount(); 305 final int nB = B.columnCount(); 306 final Matrix ret = A.newInstance(mA, nB); 307 for (int i = 0; i < mA; i++) { 308 final Vector column = A.column(i); 309 for (int j = 0; j < nB; j++) { 310 final double dot = column.dot(B.column(j)); 311 if (Math.abs(dot) > EPS) 312 ret.put(i, j, dot); 313 } 314 } 315 return ret; 316 } 317 318 /** 319 * Compute Y = A . B^T 320 * 321 * @param A 322 * @param B 323 * @return Y 324 */ 325 public static Matrix dotProductTranspose(Matrix A, Matrix B) { 326 final Matrix ret = A.newInstance(A.rowCount(), B.rowCount()); 327 return dotProductTranspose(A, B, ret); 328 } 329 330 /** 331 * Perform: A.T.dot(B.T) without performing the transpose. This is fine for 332 * dense matricies but is very inefficient for sparse matrices, consider 333 * performing the transpose manually. 334 * 335 * @param A 336 * @param B 337 * @return A.T.dot(B.T) 338 */ 339 public static Matrix dotProductTransposeTranspose(Matrix A, Matrix B) { 340 final int mA = A.columnCount(); 341 final int nB = B.rowCount(); 342 final Matrix ret = A.newInstance(mA, nB); 343 for (int i = 0; i < mA; i++) { 344 final Vector column = A.column(i); 345 for (int j = 0; j < nB; j++) { 346 final double dot = column.dot(B.row(j)); 347 if (Math.abs(dot) > EPS) 348 ret.put(i, j, dot); 349 } 350 } 351 return ret; 352 } 353 354 /** 355 * Y = A . Bt 356 * 357 * @param A 358 * @param B 359 * @param Y 360 * @return Y 361 */ 362 public static <T extends Matrix> T dotProductTranspose(Matrix A, Matrix B, 363 T Y) 364 { 365 if (A.columnCount() != B.columnCount()) 366 throw new RuntimeException( 367 String.format("Matrix size mismatch, A.cols == %d and B.T.cols == %d", A.columnCount(), 368 B.columnCount())); 369 final int mA = A.rowCount(); 370 final int nB = B.rowCount(); 371 for (int i = 0; i < mA; i++) { 372 final Vector arow = A.row(i); 373 for (int j = 0; j < nB; j++) { 374 final Vector brow = B.row(j); 375 final double dot = arow.dot(brow); 376 if (Math.abs(dot) > EPS) 377 Y.put(i, j, dot); 378 } 379 } 380 return Y; 381 } 382 383 /** 384 * A = A . s 385 * 386 * @param A 387 * @param s 388 * @return A 389 */ 390 public static <T extends Matrix> T scaleInplace(T A, double s) { 391 for (final Vector row : A.rows()) { 392 row.timesEquals(s); 393 } 394 return A; 395 } 396 397 /** 398 * @param laplacian 399 * @return returns a dense jama matrix 400 */ 401 public static Jama.Matrix toJama(Matrix laplacian) { 402 double[][] asArray = null; 403 if (laplacian instanceof DenseMatrix) { 404 asArray = ((DenseMatrix) laplacian).unwrap(); 405 } else { 406 asArray = laplacian.asArray(); 407 } 408 final Jama.Matrix ret = new Jama.Matrix(asArray, laplacian.rowCount(), 409 laplacian.columnCount()); 410 return ret; 411 } 412 413 /** 414 * @param vector 415 * @return the vector as a column in a matrix 416 */ 417 public static Jama.Matrix toColJama(Vector vector) { 418 final double[] vec = new double[vector.size()]; 419 vector.storeOn(vec, 0); 420 final Jama.Matrix ret = new Jama.Matrix(vec.length, 1); 421 for (int i = 0; i < vec.length; i++) { 422 ret.set(i, 0, vec[i]); 423 } 424 425 return ret; 426 } 427 428 /** 429 * @param vector 430 * @return the vector as a row in a matrix 431 */ 432 public static Jama.Matrix toRowJama(Vector vector) { 433 final double[] vec = new double[vector.size()]; 434 vector.storeOn(vec, 0); 435 final Jama.Matrix ret = new Jama.Matrix(1, vec.length); 436 for (int i = 0; i < vec.length; i++) { 437 ret.set(0, i, vec[i]); 438 } 439 440 return ret; 441 } 442 443 /** 444 * @param sol 445 * @return Dense matrix from a {@link Jama.Matrix} 446 */ 447 public static Matrix fromJama(Jama.Matrix sol) { 448 // final DenseMatrix mat = new DenseMatrix(sol.getRowDimension(), 449 // sol.getColumnDimension()); 450 // for (int i = 0; i < mat.rowCount(); i++) { 451 // for (int j = 0; j < mat.columnCount(); j++) { 452 // mat.put(i, j, sol.get(i, j)); 453 // } 454 // } 455 final Matrix mat = new DenseMatrix(sol.getArray()); 456 return mat; 457 } 458 459 /** 460 * @param sol 461 * @return Dense matrix from a {@link Jama.Matrix} 462 */ 463 public static no.uib.cipr.matrix.Matrix toMTJ(Matrix sol) { 464 no.uib.cipr.matrix.Matrix mat; 465 if (sol instanceof SparseMatrix) { 466 final FlexCompRowMatrix fmat = new FlexCompRowMatrix( 467 sol.rowCount(), sol.columnCount()); 468 int i = 0; 469 for (final Vector vec : sol.rows()) { 470 471 final no.uib.cipr.matrix.sparse.SparseVector x = new no.uib.cipr.matrix.sparse.SparseVector( 472 vec.size(), vec.used()); 473 for (final Entry ve : vec.entries()) { 474 x.set(ve.index, ve.value); 475 } 476 fmat.setRow(i, x); 477 i++; 478 } 479 mat = fmat; 480 } else { 481 mat = new no.uib.cipr.matrix.DenseMatrix(sol.rowCount(), 482 sol.columnCount()); 483 for (int i = 0; i < sol.rowCount(); i++) { 484 for (int j = 0; j < sol.columnCount(); j++) { 485 mat.set(i, j, sol.get(i, j)); 486 } 487 } 488 } 489 return mat; 490 } 491 492 /** 493 * Extract the submatrix of the same type of mat 494 * 495 * @param mat 496 * @param rowstart 497 * @param rowend 498 * @param colstart 499 * @param colend 500 * @return new instance 501 */ 502 public static <T extends Matrix> T subMatrix(T mat, int rowstart, 503 int rowend, int colstart, int colend) 504 { 505 @SuppressWarnings("unchecked") 506 final T ret = (T) mat.newInstance(rowend - rowstart, colend - colstart); 507 508 for (int i = 0; i < ret.rowCount(); i++) { 509 final Vector row = mat.row(i + rowstart); 510 for (final Entry ent : row.entries()) { 511 if (ent.index >= colstart && ent.index < colend) { 512 ret.put(i, ent.index - colstart, ent.value); 513 } 514 } 515 } 516 517 return ret; 518 } 519 520 /** 521 * Extract a submatrix from the given rows and cols 522 * 523 * @param mat 524 * the matrix to extract from 525 * @param rows 526 * the rows to extract 527 * @param cols 528 * the columns to extract 529 * @return the extracted matrix 530 */ 531 public static <T extends Matrix> T subMatrix(final T mat, 532 TIntCollection rows, TIntCollection cols) 533 { 534 final TIntIntHashMap actualCols; 535 if (!(cols instanceof TIntIntHashMap)) { 536 actualCols = new TIntIntHashMap(); 537 cols.forEach(new TIntProcedure() { 538 int seen = 0; 539 540 @Override 541 public boolean execute(int value) { 542 actualCols.put(value, seen++); 543 return true; 544 } 545 }); 546 } else { 547 actualCols = (TIntIntHashMap) cols; 548 } 549 @SuppressWarnings("unchecked") 550 final T ret = (T) mat.newInstance(rows.size(), cols.size()); 551 rows.forEach(new TIntProcedure() { 552 int seenrows = 0; 553 554 @Override 555 public boolean execute(final int rowIndex) { 556 final Vector row = mat.row(rowIndex); 557 for (final Entry ent : row.entries()) { 558 if (actualCols.contains(ent.index)) { 559 ret.put(seenrows, actualCols.get(ent.index), ent.value); 560 } 561 } 562 seenrows++; 563 return true; 564 } 565 }); 566 567 return ret; 568 } 569 570 /** 571 * @param mat 572 * @param rows 573 * @param colstart 574 * @param colend 575 * @return the submatrix 576 */ 577 public static <T extends Matrix> T subMatrix(final T mat, 578 TIntArrayList rows, final int colstart, final int colend) 579 { 580 @SuppressWarnings("unchecked") 581 final T ret = (T) mat.newInstance(rows.size(), colend - colstart); 582 rows.forEach(new TIntProcedure() { 583 int seen = 0; 584 585 @Override 586 public boolean execute(int rowIndex) { 587 final Vector row = mat.row(rowIndex); 588 for (final Entry ent : row.entries()) { 589 if (ent.index >= colstart && ent.index < colend) { 590 ret.put(seen, ent.index - colstart, ent.value); 591 } 592 } 593 seen++; 594 return true; 595 } 596 }); 597 598 return ret; 599 600 } 601 602 /** 603 * @param mat 604 * @param rowstart 605 * @param rowend 606 * @param cols 607 * @return the submatrix 608 */ 609 public static <T extends Matrix> T subMatrix(final T mat, 610 final int rowstart, final int rowend, TIntArrayList cols) 611 { 612 @SuppressWarnings("unchecked") 613 final T ret = (T) mat.newInstance(rowend - rowstart, cols.size()); 614 cols.forEach(new TIntProcedure() { 615 int seen = 0; 616 617 @Override 618 public boolean execute(int colIndex) { 619 final Vector col = mat.column(colIndex); 620 for (final Entry ent : col.entries()) { 621 if (ent.index >= rowstart && ent.index < rowend) { 622 ret.put(ent.index - rowstart, seen, ent.value); 623 } 624 } 625 seen++; 626 return true; 627 } 628 }); 629 630 return ret; 631 632 } 633 634 /** 635 * @param m 636 * @return a {@link MLDouble} for matlab 637 */ 638 public static MLDouble asMatlab(Matrix m) { 639 final double[][] retArr = new double[m.rowCount()][m.columnCount()]; 640 for (int i = 0; i < retArr.length; i++) { 641 for (int j = 0; j < retArr[i].length; j++) { 642 retArr[i][j] = m.get(i, j); 643 } 644 } 645 final MLDouble ret = new MLDouble("out", retArr); 646 return ret; 647 } 648 649 /** 650 * Calculate all 3, used by {@link #min(Matrix)}, {@link #max(Matrix)} and 651 * {@link #mean(Matrix)} 652 * 653 * @param mat 654 * @return the min, max and mean of the provided matrix 655 */ 656 public static double[] minmaxmean(Matrix mat) { 657 double min = Double.MAX_VALUE, max = -Double.MAX_VALUE, mean = 0; 658 final double size = mat.rowCount() * mat.columnCount(); 659 for (final Vector v : mat.rows()) { 660 for (final Entry ent : v.entries()) { 661 min = Math.min(min, ent.value); 662 max = Math.max(max, ent.value); 663 mean += ent.value / size; 664 } 665 } 666 return new double[] { min, max, mean }; 667 } 668 669 /** 670 * uses the first value returned by {@link #minmaxmean(Matrix)} 671 * 672 * @param mat 673 * @return the min 674 */ 675 public static double min(Matrix mat) { 676 return minmaxmean(mat)[0]; 677 } 678 679 /** 680 * uses the second value returned by {@link #minmaxmean(Matrix)} 681 * 682 * @param mat 683 * @return the min 684 */ 685 public static double max(Matrix mat) { 686 return minmaxmean(mat)[1]; 687 } 688 689 /** 690 * uses the third value returned by {@link #minmaxmean(Matrix)} 691 * 692 * @param mat 693 * @return the min 694 */ 695 public static double mean(Matrix mat) { 696 return minmaxmean(mat)[2]; 697 } 698 699 /** 700 * @param l 701 * @param v 702 * @return performs l - v returning a matrix of type T 703 */ 704 public static <T extends Matrix> T minus(T l, double v) { 705 @SuppressWarnings("unchecked") 706 final T ret = (T) l.newInstance(l.rowCount(), l.columnCount()); 707 int r = 0; 708 for (final Vector vec : l.rows()) { 709 for (final Entry ent : vec.entries()) { 710 ret.put(r, ent.index, ent.value - v); 711 } 712 r++; 713 } 714 return ret; 715 } 716 717 /** 718 * @param l 719 * @param v 720 * @return performs l - v returning a matrix of type T 721 */ 722 public static Vector minus(Vector l, Vector v) { 723 final Vector ret = DenseVector.dense(l.size()); 724 for (int i = 0; i < l.size(); i++) { 725 ret.put(i, l.get(i) - v.get(i)); 726 } 727 return ret; 728 } 729 730 /** 731 * @param v 732 * @param l 733 * @return performs v - l returning a matrix of type T 734 */ 735 public static <T extends Matrix> T minus(double v, T l) { 736 @SuppressWarnings("unchecked") 737 final T ret = (T) l.newInstance(l.rowCount(), l.columnCount()); 738 for (int i = 0; i < l.rowCount(); i++) { 739 for (int j = 0; j < l.columnCount(); j++) { 740 ret.put(i, j, v - l.get(i, j)); 741 } 742 } 743 744 return ret; 745 } 746 747 /** 748 * Create a {@link Matrix} from a matlab {@link MLArray} 749 * 750 * @param mlArray 751 * the matlab array 752 * @return the matrix 753 */ 754 public static Matrix fromMatlab(MLArray mlArray) { 755 final Matrix mat = new DenseMatrix(mlArray.getM(), mlArray.getN()); 756 final MLDouble mlDouble = (MLDouble) mlArray; 757 for (int i = 0; i < mat.rowCount(); i++) { 758 for (int j = 0; j < mat.columnCount(); j++) { 759 mat.put(i, j, mlDouble.get(i, j)); 760 } 761 } 762 return mat; 763 } 764 765 /** 766 * Sum the diagonal of the given matrix 767 * 768 * @param d 769 * the matrix 770 * @return the sum along the diagonal 771 */ 772 public static double sum(DiagonalMatrix d) { 773 return ArrayUtils.sumValues(d.getVals()); 774 } 775 776 /** 777 * Set a submatrix of a larger matrix 778 * 779 * @param to 780 * the matrix to write into 781 * @param row 782 * the row to start inserting from 783 * @param col 784 * the column to start inserting from 785 * @param from 786 * the matrix to insert 787 */ 788 public static void setSubMatrix(Matrix to, int row, int col, Matrix from) { 789 for (int i = row; i < row + from.rowCount(); i++) { 790 for (int j = col; j < col + from.columnCount(); j++) { 791 to.put(i, j, from.get(i - row, j - col)); 792 } 793 } 794 } 795 796 /** 797 * Transpose a matrix, returning a new matrix. 798 * 799 * @param mat 800 * the matrix to transpose 801 * @return the transposed matrix 802 */ 803 public static <T extends Matrix> T transpose(T mat) { 804 @SuppressWarnings("unchecked") 805 final T ret = (T) mat.newInstance(mat.columnCount(), mat.rowCount()); 806 // for (int i = 0; i < mat.columnCount(); i++) { 807 // for (int j = 0; j < mat.rowCount(); j++) { 808 // ret.put(i, j, mat.get(j, i)); 809 // } 810 // } 811 for (int i = 0; i < mat.rowCount(); i++) { 812 final Vector v = mat.row(i); 813 for (final Entry ent : v.entries()) { 814 ret.put(ent.index, i, ent.value); 815 } 816 } 817 818 return ret; 819 } 820 821 /** 822 * @param A 823 * @param B 824 * @return A = MAX(A,B) 825 */ 826 public static SparseMatrix maxInplace(SparseMatrix A, SparseMatrix B) { 827 for (int i = 0; i < A.rowCount(); i++) { 828 final Vector arow = A.row(i); 829 final Vector brow = B.row(i); 830 for (final Entry br : brow.entries()) { 831 if (arow.get(br.index) < br.value) { 832 A.put(i, br.index, br.value); 833 } 834 } 835 } 836 return A; 837 } 838 839 /** 840 * @param A 841 * @param B 842 * @return A = MIN(A,B) 843 */ 844 public static SparseMatrix minInplace(SparseMatrix A, SparseMatrix B) { 845 for (int i = 0; i < A.rowCount(); i++) { 846 final Vector arow = A.row(i); 847 final Vector brow = B.row(i); 848 for (final Entry br : brow.entries()) { 849 if (arow.get(br.index) > br.value) { 850 A.put(i, br.index, br.value); 851 } 852 } 853 } 854 return A; 855 } 856 857 /** 858 * @param A 859 * @param B 860 * @return A = A.*B 861 */ 862 public static SparseMatrix timesInplace(SparseMatrix A, SparseMatrix B) { 863 for (int i = 0; i < A.rowCount(); i++) { 864 final Vector arow = A.row(i); 865 final Vector brow = B.row(i); 866 for (final Entry br : brow.entries()) { 867 A.put(i, br.index, br.value * arow.get(br.index)); 868 } 869 870 for (final Entry ar : arow.entries()) { 871 if (brow.get(ar.index) == 0) // All items in A not in B must be 872 // set to 0 873 A.put(i, ar.index, 0); 874 } 875 } 876 return A; 877 } 878 879 /** 880 * @param A 881 * @param B 882 * @return A = A.*B 883 */ 884 public static Matrix timesInplace(Matrix A, Matrix B) { 885 for (int i = 0; i < A.rowCount(); i++) { 886 final Vector arow = A.row(i); 887 final Vector brow = B.row(i); 888 for (final Entry br : brow.entries()) { 889 A.put(i, br.index, br.value * arow.get(br.index)); 890 } 891 892 for (final Entry ar : arow.entries()) { 893 if (brow.get(ar.index) == 0) // All items in A not in B must be 894 // set to 0 895 A.put(i, ar.index, 0); 896 } 897 } 898 return A; 899 } 900 901 /** 902 * Copy a matrix 903 * 904 * @param sparseMatrix 905 * the matrix to copy 906 * @return the copy 907 */ 908 public static <T extends Matrix> T copy(T sparseMatrix) { 909 @SuppressWarnings("unchecked") 910 final T t = (T) sparseMatrix.newInstance(); 911 912 for (int r = 0; r < sparseMatrix.rowCount(); r++) { 913 final Vector row = sparseMatrix.row(r); 914 for (final Entry ent : row.entries()) { 915 t.put(r, ent.index, ent.value); 916 } 917 } 918 return t; 919 } 920 921 /** 922 * Set values below the given threshold to zero in the output matrix. 923 * 924 * @param data 925 * the input matrix 926 * @param thresh 927 * the threshold 928 * @return a new matrix with values in the input matrix set to zero. 929 */ 930 public static SparseMatrix threshold(SparseMatrix data, double thresh) { 931 final SparseMatrix newdata = (SparseMatrix) data.newInstance(); 932 933 for (int r = 0; r < data.rowCount(); r++) { 934 final Vector vec = data.row(r); 935 for (final Entry ent : vec.entries()) { 936 if (ent.value > thresh) { 937 newdata.put(r, ent.index, 1); 938 } 939 } 940 } 941 return newdata; 942 } 943 944 /** 945 * Convert a {@link ch.akuhn.matrix.SparseVector} to a 946 * {@link SparseDoubleArray}. 947 * 948 * @param row 949 * the vector to convert 950 * @return the converted vector 951 */ 952 public static SparseDoubleArray sparseVectorToSparseArray( 953 ch.akuhn.matrix.SparseVector row) 954 { 955 final SparseDoubleArray sda = new SparseBinSearchDoubleArray( 956 row.size(), row.used(), row.keys(), row.values()); 957 return sda; 958 } 959 960 /** 961 * 962 * @param to 963 * add items to this 964 * @param startindex 965 * starting index in to 966 * @param from 967 * add items from this 968 */ 969 public static void setSubVector(Vector to, int startindex, Vector from) { 970 if (to instanceof DenseVector && from instanceof DenseVector) { 971 final double[] tod = ((DenseVector) to).unwrap(); 972 final double[] fromd = ((DenseVector) from).unwrap(); 973 System.arraycopy(fromd, 0, tod, startindex, fromd.length); 974 return; 975 } 976 for (int i = 0; i < from.size(); i++) { 977 to.put(i + startindex, from.get(i)); 978 } 979 } 980 981 /** 982 * Starting from a given column of a row, set the values of a matrix to the 983 * values of v 984 * 985 * @param to 986 * @param row 987 * @param col 988 * @param v 989 */ 990 public static void setSubMatrixRow(Matrix to, int row, int col, Vector v) { 991 for (int i = col, j = 0; i < col + v.size(); i++, j++) { 992 to.put(row, i, v.get(j)); 993 } 994 } 995 996 /** 997 * Starting from a given row of a column, set the values of a matrix to the 998 * values of v 999 * 1000 * @param to 1001 * @param row 1002 * @param col 1003 * @param v 1004 */ 1005 public static void setSubMatrixCol(Matrix to, int row, int col, Vector v) { 1006 for (int i = row, j = 0; i < row + v.size(); i++, j++) { 1007 to.put(i, col, v.get(j)); 1008 } 1009 } 1010 1011 /** 1012 * @param v 1013 * the value vector 1014 * @param d 1015 * the check value 1016 * @return for each item in the vector, returns 1 if the value is less than 1017 * the check value 1018 */ 1019 public static Vector lessThan(Vector v, double d) { 1020 final Vector out = new SparseVector(v.size(), 1); 1021 for (int i = 0; i < v.size(); i++) { 1022 if (v.get(i) < d) 1023 out.put(i, 1); 1024 } 1025 return out; 1026 } 1027 1028 /** 1029 * @param v 1030 * @return any values of the vector are not-zero 1031 */ 1032 public static boolean any(Vector v) { 1033 for (int i = 0; i < v.size(); i++) { 1034 if (v.get(i) != 0) 1035 return true; 1036 } 1037 return false; 1038 } 1039 1040 /** 1041 * @param m 1042 * @param col 1043 * @return m with the added column 1044 */ 1045 public static <T extends Matrix> T appendColumn(T m, Vector col) { 1046 @SuppressWarnings("unchecked") 1047 final T ret = (T) m.newInstance(m.rowCount(), m.columnCount() + 1); 1048 setSubMatrixCol(ret, 0, m.columnCount(), col); 1049 return ret; 1050 } 1051 1052 /** 1053 * @param m 1054 * @param row 1055 * @return m with the added column 1056 */ 1057 public static <T extends Matrix> T appendRow(T m, Vector row) { 1058 @SuppressWarnings("unchecked") 1059 final T ret = (T) m.newInstance(m.rowCount() + 1, m.columnCount()); 1060 setSubMatrixRow(ret, m.rowCount(), 0, row); 1061 return ret; 1062 } 1063 1064 /** 1065 * Create a {@link Matrix} from the Cognitive Foundry equivalent 1066 * 1067 * @param init 1068 * the matrix 1069 * @return the converted matrix 1070 */ 1071 public static Matrix fromCF(gov.sandia.cognition.math.matrix.Matrix init) { 1072 Matrix ret; 1073 final int m = init.getNumRows(); 1074 final int n = init.getNumColumns(); 1075 if (init instanceof gov.sandia.cognition.math.matrix.mtj.SparseMatrix) { 1076 ret = SparseMatrix.sparse(init.getNumRows(), init.getNumColumns()); 1077 } else { 1078 ret = DenseMatrix.dense(m, n); 1079 } 1080 for (final MatrixEntry ent : init) { 1081 ret.put(ent.getRowIndex(), ent.getColumnIndex(), ent.getValue()); 1082 } 1083 return ret; 1084 } 1085 1086 /** 1087 * Compute the dot product X.W 1088 * 1089 * @param X 1090 * @param W 1091 * @return dot product 1092 */ 1093 public static Matrix dotProduct(Matrix X, Matrix W) { 1094 final Matrix ret = X.newInstance(X.rowCount(), W.columnCount()); 1095 for (int j = 0; j < ret.columnCount(); j++) { 1096 final Vector column = W.column(j); 1097 for (int i = 0; i < ret.rowCount(); i++) { 1098 ret.put(i, j, X.row(i).dot(column)); 1099 } 1100 } 1101 return ret; 1102 } 1103 1104 /** 1105 * Compute the 2-norm (Euclidean norm) of the vector 1106 * 1107 * @param row 1108 * the vector 1109 * @return the Euclidean norm 1110 */ 1111 public static double norm2(Vector row) { 1112 double norm = 0; 1113 for (final Entry e : row.entries()) 1114 norm += e.value * e.value; 1115 return Math.sqrt(norm); 1116 } 1117 1118 /** 1119 * Subtract matrices A-B 1120 * 1121 * @param A 1122 * @param B 1123 * @return A-B 1124 */ 1125 public static Matrix minus(Matrix A, Matrix B) { 1126 final Matrix ret = copy(A); 1127 minusInplace(ret, B); 1128 return ret; 1129 } 1130 1131 /** 1132 * Compute the Frobenius norm 1133 * 1134 * @param A 1135 * the matrix 1136 * @return the F norm 1137 */ 1138 public static double normF(Matrix A) { 1139 double scale = 0, ssq = 1; 1140 for (final Vector v : A.rows()) { 1141 for (final Entry e : v.entries()) { 1142 final double Aval = e.value; 1143 if (Aval != 0) { 1144 final double absxi = Math.abs(Aval); 1145 if (scale < absxi) { 1146 ssq = 1 + ssq * Math.pow(scale / absxi, 2); 1147 scale = absxi; 1148 } else { 1149 ssq = ssq + Math.pow(absxi / scale, 2); 1150 } 1151 } 1152 } 1153 } 1154 return scale * Math.sqrt(ssq); 1155 } 1156 1157 /** 1158 * Stack matrices vertically 1159 * 1160 * @param matricies 1161 * matrices to stack 1162 * @return matrix created from the stacking 1163 */ 1164 public static Matrix vstack(Matrix... matricies) { 1165 int nrows = 0; 1166 int ncols = 0; 1167 for (final Matrix matrix : matricies) { 1168 nrows += matrix.rowCount(); 1169 ncols = matrix.columnCount(); 1170 } 1171 final Matrix ret = matricies[0].newInstance(nrows, ncols); 1172 int currentRow = 0; 1173 for (final Matrix matrix : matricies) { 1174 setSubMatrix(ret, currentRow, 0, matrix); 1175 currentRow += matrix.rowCount(); 1176 } 1177 return ret; 1178 } 1179 1180}