001// --- BEGIN LICENSE BLOCK ---
002/* 
003 * Copyright (c) 2009, Mikio L. Braun
004 * Copyright (c) 2008, Johannes Schaback
005 * Copyright (c) 2009, Jan Saputra M?ller
006 * All rights reserved.
007 * 
008 * Redistribution and use in source and binary forms, with or without
009 * modification, are permitted provided that the following conditions are
010 * met:
011 * 
012 *     * Redistributions of source code must retain the above copyright
013 *       notice, this list of conditions and the following disclaimer.
014 * 
015 *     * Redistributions in binary form must reproduce the above
016 *       copyright notice, this list of conditions and the following
017 *       disclaimer in the documentation and/or other materials provided
018 *       with the distribution.
019 * 
020 *     * Neither the name of the Technische Universit?t Berlin nor the
021 *       names of its contributors may be used to endorse or promote
022 *       products derived from this software without specific prior
023 *       written permission.
024 * 
025 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
026 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
027 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
028 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
029 * HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
030 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
031 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
032 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
033 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
034 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
035 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
036 */
037// --- END LICENSE BLOCK ---
038package org.jblas;
039
040import org.jblas.exceptions.SizeException;
041import org.jblas.ranges.Range;
042import java.io.BufferedReader;
043import java.io.DataInputStream;
044import java.io.DataOutputStream;
045import java.io.FileInputStream;
046import java.io.FileOutputStream;
047import java.io.IOException;
048import java.io.InputStreamReader;
049
050import java.io.ObjectInputStream;
051import java.io.ObjectOutputStream;
052import java.io.PrintWriter;
053import java.io.Serializable;
054import java.io.StringWriter;
055import java.util.AbstractList;
056import java.util.Arrays;
057import java.util.Comparator;
058import java.util.Iterator;
059import java.util.LinkedList;
060import java.util.List;
061
062/**
063 * A general matrix class for <tt>double</tt> typed values.
064 * 
065 * Don't be intimidated by the large number of methods this function defines. Most
066 * are overloads provided for ease of use. For example, for each arithmetic operation,
067 * up to six overloaded versions exist to handle in-place computations, and
068 * scalar arguments.
069 * 
070 * <h3>Construction</h3>
071 * 
072 * <p>To construct a two-dimensional matrices, you can use the following constructors
073 * and static methods.</p>
074 * 
075 * <table class="my">
076 * <tr><th>Method<th>Description
077 * <tr><td>DoubleMatrix(m,n, [value1, value2, value3...])<td>Values are filled in row by row.
078 * <tr><td>DoubleMatrix(new double[][] {{value1, value2, ...}, ...}<td>Inner arrays are columns.
079 * <tr><td>DoubleMatrix.zeros(m,n) <td>Initial values set to 0.0.
080 * <tr><td>DoubleMatrix.ones(m,n) <td>Initial values set to 1.0.
081 * <tr><td>DoubleMatrix.rand(m,n) <td>Values drawn at random between 0.0 and 1.0.
082 * <tr><td>DoubleMatrix.randn(m,n) <td>Values drawn from normal distribution.
083 * <tr><td>DoubleMatrix.eye(n) <td>Unit matrix (values 0.0 except for 1.0 on the diagonal).
084 * <tr><td>DoubleMatrix.diag(array) <td>Diagonal matrix with given diagonal elements.
085 * </table>
086 * 
087 * <p>Alternatively, you can construct (column) vectors, if you just supply the length
088 * using the following constructors and static methods.</p>
089 * 
090 * <table class="my">
091 * <tr><th>Method<th>Description
092 * <tr><td>DoubleMatrix(m)<td>Constructs a column vector.
093 * <tr><td>DoubleMatrix(new double[] {value1, value2, ...})<td>Constructs a column vector.
094 * <tr><td>DoubleMatrix.zeros(m) <td>Initial values set to 0.0.
095 * <tr><td>DoubleMatrix.ones(m) <td>Initial values set to 1.0.
096 * <tr><td>DoubleMatrix.rand(m) <td>Values drawn at random between 0.0 and 1.0.
097 * <tr><td>DoubleMatrix.randn(m) <td>Values drawn from normal distribution.
098 * </table>
099 * 
100 * <p>You can also construct new matrices by concatenating matrices either horziontally
101 * or vertically:</p>
102 * 
103 * <table class="my">
104 * <tr><th>Method<th>Description
105 * <tr><td>x.concatHorizontally(y)<td>New matrix will be x next to y.
106 * <tr><td>x.concatVertically(y)<td>New matrix will be x atop y.
107 * </table>
108 * 
109 * <h3>Element Access, Copying and Duplication</h3>
110 * 
111 * <p>To access individual elements, or whole rows and columns, use the following
112 * methods:<p>
113 * 
114 * <table class="my">
115 * <tr><th>x.Method<th>Description
116 * <tr><td>x.get(i,j)<td>Get element in row i and column j.
117 * <tr><td>x.put(i, j, v)<td>Set element in row i and column j to value v
118 * <tr><td>x.get(i)<td>Get the ith element of the matrix (traversing rows first).
119 * <tr><td>x.put(i, v)<td>Set the ith element of the matrix (traversing rows first).
120 * <tr><td>x.getColumn(i)<td>Get a copy of column i.
121 * <tr><td>x.putColumn(i, c)<td>Put matrix c into column i.
122 * <tr><td>x.getRow(i)<td>Get a copy of row i.
123 * <tr><td>x.putRow(i, c)<td>Put matrix c into row i.
124 * <tr><td>x.swapColumns(i, j)<td>Swap the contents of columns i and j.
125 * <tr><td>x.swapRows(i, j)<td>Swap the contents of columns i and j.
126 * </table>
127 * 
128 * <p>For <tt>get</tt> and <tt>put</tt>, you can also pass integer arrays,
129 * DoubleMatrix objects, or Range objects, which then specify the indices used 
130 * as follows:
131 * 
132 * <ul>
133 * <li><em>integer array:</em> the elements will be used as indices.
134 * <li><em>DoubleMatrix object:</em> non-zero entries specify the indices.
135 * <li><em>Range object:</em> see below.
136 * </ul>
137 * 
138 * <p>When using <tt>put</tt> with multiple indices, the assigned object must
139 * have the correct size or be a scalar.</p>
140 *
141 * <p>There exist the following Range objects. The Class <tt>RangeUtils</tt> also
142 * contains the a number of handy helper methods for constructing these ranges.</p>
143 * <table class="my">
144 * <tr><th>Class <th>RangeUtils method <th>Indices
145 * <tr><td>AllRange <td>all() <td>All legal indices.
146 * <tr><td>PointRange <td>point(i) <td> A single point.
147 * <tr><td>IntervalRange <td>interval(a, b)<td> All indices from a to b (inclusive)
148 * <tr><td rowspan=3>IndicesRange <td>indices(int[])<td> The specified indices.
149 * <tr><td>indices(DoubleMatrix)<td>The specified indices.
150 * <tr><td>find(DoubleMatrix)<td>The non-zero entries of the matrix.
151 * </table>
152 * 
153 * <p>The following methods can be used for duplicating and copying matrices.</p>
154 * 
155 * <table class="my">
156 * <tr><th>Method<th>Description
157 * <tr><td>x.dup()<td>Get a copy of x.
158 * <tr><td>x.copy(y)<td>Copy the contents of y to x (possible resizing x).
159 * </table>
160 *    
161 * <h3>Size and Shape</h3>
162 * 
163 * <p>The following methods permit to acces the size of a matrix and change its size or shape.</p>
164 * 
165 * <table class="my">
166 * <tr><th>x.Method<th>Description
167 * <tr><td>x.rows<td>Number of rows.
168 * <tr><td>x.columns<td>Number of columns.
169 * <tr><td>x.length<td>Total number of elements.
170 * <tr><td>x.isEmpty()<td>Checks whether rows == 0 and columns == 0.
171 * <tr><td>x.isRowVector()<td>Checks whether rows == 1.
172 * <tr><td>x.isColumnVector()<td>Checks whether columns == 1.
173 * <tr><td>x.isVector()<td>Checks whether rows == 1 or columns == 1.
174 * <tr><td>x.isSquare()<td>Checks whether rows == columns.
175 * <tr><td>x.isScalar()<td>Checks whether length == 1.
176 * <tr><td>x.resize(r, c)<td>Resize the matrix to r rows and c columns, discarding the content.
177 * <tr><td>x.reshape(r, c)<td>Resize the matrix to r rows and c columns.<br> Number of elements must not change.
178 * </table>
179 * 
180 * <p>The size is stored in the <tt>rows</tt> and <tt>columns</tt> member variables.
181 * The total number of elements is stored in <tt>length</tt>. Do not change these
182 * values unless you know what you're doing!</p>
183 * 
184 * <h3>Arithmetics</h3>
185 * 
186 * <p>The usual arithmetic operations are implemented. Each operation exists in a
187 * in-place version, recognizable by the suffix <tt>"i"</tt>, to which you can supply
188 * the result matrix (or <tt>this</tt> is used, if missing). Using in-place operations
189 * can also lead to a smaller memory footprint, as the number of temporary objects
190 * which are directly garbage collected again is reduced.</p>
191 * 
192 * <p>Whenever you specify a result vector, the result vector must already have the
193 * correct dimensions.</p>
194 * 
195 * <p>For example, you can add two matrices using the <tt>add</tt> method. If you want
196 * to store the result in of <tt>x + y</tt> in <tt>z</tt>, type
197 * <span class=code>
198 * x.addi(y, z)   // computes x = y + z.
199 * </span>
200 * Even in-place methods return the result, such that you can easily chain in-place methods,
201 * for example:
202 * <span class=code>
203 * x.addi(y).addi(z) // computes x += y; x += z
204 * </span></p> 
205 *
206 * <p>Methods which operate element-wise only make sure that the length of the matrices
207 * is correct. Therefore, you can add a 3 * 3 matrix to a 1 * 9 matrix, for example.</p>
208 * 
209 * <p>Finally, there exist versions which take doubles instead of DoubleMatrix Objects
210 * as arguments. These then compute the operation with the same value as the
211 * right-hand-side. The same effect can be achieved by passing a DoubleMatrix with
212 * exactly one element.</p>
213 * 
214 * <table class="my">
215 * <tr><th>Operation <th>Method <th>Comment
216 * <tr><td>x + y <td>x.add(y)                 <td>
217 * <tr><td>x - y <td>x.sub(y), y.rsub(x) <td>rsub subtracts left from right hand side
218 * <tr><td rowspan=3>x * y  <td>x.mul(y) <td>element-wise multiplication 
219 * <tr>                                           <td>x.mmul(y)<td>matrix-matrix multiplication
220 * <tr>                                           <td>x.dot(y) <td>scalar-product
221 * <tr><td>x / y <td>x.div(y), y.rdiv(x) <td>rdiv divides right hand side by left hand side.
222 * <tr><td>- x      <td>x.neg()                               <td>
223 * </table>
224 * 
225 * <p>There also exist operations which work on whole columns or rows.</p>
226 * 
227 * <table class="my">
228 * <tr><th>Method <th>Description
229 * <tr><td>x.addRowVector<td>adds a vector to each row (addiRowVector works in-place)
230 * <tr><td>x.addColumnVector<td>adds a vector to each column
231 * <tr><td>x.subRowVector<td>subtracts a vector from each row
232 * <tr><td>x.subColumnVector<td>subtracts a vector from each column
233 * <tr><td>x.mulRow<td>Multiplies a row by a scalar
234 * <tr><td>x.mulColumn<td>multiplies a row by a column
235 * </table>
236 * 
237 * <p>In principle, you could achieve the same result by first calling getColumn(), 
238 * adding, and then calling putColumn, but these methods are much faster.</p>
239 * 
240 * <p>The following comparison operations are available</p>
241 *  
242 * <table class="my">
243 * <tr><th>Operation <th>Method
244 * <tr><td>x &lt; y             <td>x.lt(y)
245 * <tr><td>x &lt;= y    <td>x.le(y)
246 * <tr><td>x &gt; y             <td>x.gt(y)
247 * <tr><td>x &gt;= y    <td>x.ge(y)
248 * <tr><td>x == y           <td>x.eq(y)
249 * <tr><td>x != y           <td>x.ne(y)
250 * </table>
251 *
252 * <p> Logical operations are also supported. For these operations, a value different from
253 * zero is treated as "true" and zero is treated as "false". All operations are carried
254 * out elementwise.</p>
255 * 
256 * <table class="my">
257 * <tr><th>Operation <th>Method
258 * <tr><td>x & y        <td>x.and(y)
259 * <tr><td>x | y    <td>x.or(y)
260 * <tr><td>x ^ y    <td>x.xor(y)
261 * <tr><td>! x              <td>x.not()
262 * </table>
263 * 
264 * <p>Finally, there are a few more methods to compute various things:</p>
265 * 
266 * <table class="my">
267 * <tr><th>Method <th>Description
268 * <tr><td>x.max() <td>Return maximal element
269 * <tr><td>x.argmax() <td>Return index of largest element
270 * <tr><td>x.min() <td>Return minimal element
271 * <tr><td>x.argmin() <td>Return index of largest element
272 * <tr><td>x.columnMins() <td>Return column-wise minima
273 * <tr><td>x.columnArgmins() <td>Return column-wise index of minima
274 * <tr><td>x.columnMaxs() <td>Return column-wise maxima
275 * <tr><td>x.columnArgmaxs() <td>Return column-wise index of maxima
276 * </table>
277 * 
278 * @author Mikio Braun, Johannes Schaback
279 */
280public class DoubleMatrix implements Serializable {
281
282    /** Number of rows. */
283    public int rows;
284    /** Number of columns. */
285    public int columns;
286    /** Total number of elements (for convenience). */
287    public int length;
288    /** The actual data stored by rows (that is, row 0, row 1...). */
289    public double[] data = null; // rows are contiguous
290    public static final DoubleMatrix EMPTY = new DoubleMatrix();
291    static final long serialVersionUID = -1249281332731183060L;
292
293    /**************************************************************************
294     *
295     * Constructors and factory functions
296     *
297     **************************************************************************/
298    /** Create a new matrix with <i>newRows</i> rows, <i>newColumns</i> columns
299     * using <i>newData></i> as the data. The length of the data is not checked!
300     */
301    public DoubleMatrix(int newRows, int newColumns, double... newData) {
302        rows = newRows;
303        columns = newColumns;
304        length = rows * columns;
305
306        if (newData != null && newData.length != newRows * newColumns) {
307            throw new IllegalArgumentException(
308                    "Passed data must match matrix dimensions.");
309        }
310
311        data = newData;
312        //System.err.printf("%d * %d matrix created\n", rows, columns);
313    }
314
315    /**
316     * Creates a new <i>n</i> times <i>m</i> <tt>DoubleMatrix</tt>.
317     * @param newRows the number of rows (<i>n</i>) of the new matrix.
318     * @param newColumns the number of columns (<i>m</i>) of the new matrix.
319     */
320    public DoubleMatrix(int newRows, int newColumns) {
321        this(newRows, newColumns, new double[newRows * newColumns]);
322    }
323
324    /**
325     * Creates a new <tt>DoubleMatrix</tt> of size 0 times 0.
326     */
327    public DoubleMatrix() {
328        this(0, 0, (double[]) null);
329    }
330
331    /**
332     * Create a Matrix of length <tt>len</tt>. By default, this creates a row vector.
333     * @param len
334     */
335    public DoubleMatrix(int len) {
336        this(len, 1, new double[len]);
337    }
338
339    public DoubleMatrix(double[] newData) {
340        this(newData.length);
341        data = newData;
342    }
343
344    /**
345     * Creates a new matrix by reading it from a file.
346     * @param filename the path and name of the file to read the matrix from
347     * @throws IOException
348     */
349    public DoubleMatrix(String filename) throws IOException {
350        load(filename);
351    }
352
353    /**
354     * Creates a new <i>n</i> times <i>m</i> <tt>DoubleMatrix</tt> from
355     * the given <i>n</i> times <i>m</i> 2D data array. The first dimension of the array makes the
356     * rows (<i>n</i>) and the second dimension the columns (<i>m</i>). For example, the
357     * given code <br/><br/>
358     * <code>new DoubleMatrix(new double[][]{{1d, 2d, 3d}, {4d, 5d, 6d}, {7d, 8d, 9d}}).print();</code><br/><br/>
359     * will constructs the following matrix:
360     * <pre>
361     * 1.0      2.0     3.0
362     * 4.0      5.0     6.0
363     * 7.0      8.0     9.0
364     * </pre>.
365     * @param data <i>n</i> times <i>m</i> data array
366     */
367    public DoubleMatrix(double[][] data) {
368        this(data.length, data[0].length);
369
370        for (int r = 0; r < rows; r++) {
371            assert (data[r].length == columns);
372        }
373
374        for (int r = 0; r < rows; r++) {
375            for (int c = 0; c < columns; c++) {
376                put(r, c, data[r][c]);
377            }
378        }
379    }
380
381    public DoubleMatrix(List<Double> data) {
382        this(data.size());
383
384        int c = 0;
385        for (java.lang.Double d : data) {
386            put(c++, d);
387        }
388    }
389
390    /**
391     * Construct DoubleMatrix from ASCII representation.
392     *
393     * This is not very fast, but can be quiet useful when
394     * you want to "just" construct a matrix, for example
395     * when testing.
396     *
397     * The format is semicolon separated rows of space separated values,
398     * for example "1 2 3; 4 5 6; 7 8 9".
399     */
400    public static DoubleMatrix valueOf(String text) {
401        String[] rowValues = text.split(";");
402
403        // process first line
404        String[] columnValues = rowValues[0].trim().split("\\s+");
405
406        DoubleMatrix result = null;
407
408        // process rest
409        for (int r = 0; r < rowValues.length; r++) {
410            columnValues = rowValues[r].trim().split("\\s+");
411
412            if (r == 0) {
413                result = new DoubleMatrix(rowValues.length, columnValues.length);
414            }
415
416            for (int c = 0; c < columnValues.length; c++) {
417                result.put(r, c, Double.valueOf(columnValues[c]));
418            }
419        }
420
421        return result;
422    }
423
424    /**
425     * Serialization
426     */
427    private void writeObject(ObjectOutputStream s) throws IOException {
428        s.defaultWriteObject();
429    }
430
431    private void readObject(ObjectInputStream s) throws IOException, ClassNotFoundException {
432        s.defaultReadObject();
433    }
434
435    /** Create matrix with random values uniformly in 0..1. */
436    public static DoubleMatrix rand(int rows, int columns) {
437        DoubleMatrix m = new DoubleMatrix(rows, columns);
438
439        java.util.Random r = new java.util.Random();
440        for (int i = 0; i < rows * columns; i++) {
441            m.data[i] = r.nextDouble();
442        }
443
444        return m;
445    }
446
447    /** Creates a row vector with random values uniformly in 0..1. */
448    public static DoubleMatrix rand(int len) {
449        return rand(len, 1);
450    }
451
452    /** Create matrix with normally distributed random values. */
453    public static DoubleMatrix randn(int rows, int columns) {
454        DoubleMatrix m = new DoubleMatrix(rows, columns);
455
456        java.util.Random r = new java.util.Random();
457        for (int i = 0; i < rows * columns; i++) {
458            m.data[i] = (double) r.nextGaussian();
459        }
460
461        return m;
462    }
463
464    /** Create row vector with normally distributed random values. */
465    public static DoubleMatrix randn(int len) {
466        return randn(len, 1);
467    }
468
469    /** Creates a new matrix in which all values are equal 0. */
470    public static DoubleMatrix zeros(int rows, int columns) {
471        return new DoubleMatrix(rows, columns);
472    }
473
474    /** Creates a row vector of given length. */
475    public static DoubleMatrix zeros(int length) {
476        return zeros(length, 1);
477    }
478
479    /** Creates a new matrix in which all values are equal 1. */
480    public static DoubleMatrix ones(int rows, int columns) {
481        DoubleMatrix m = new DoubleMatrix(rows, columns);
482
483        for (int i = 0; i < rows * columns; i++) {
484            m.put(i, 1.0);
485        }
486
487        return m;
488    }
489
490    /** Creates a row vector with all elements equal to 1. */
491    public static DoubleMatrix ones(int length) {
492        return ones(length, 1);
493    }
494
495    /** Construct a new n-by-n identity matrix. */
496    public static DoubleMatrix eye(int n) {
497        DoubleMatrix m = new DoubleMatrix(n, n);
498
499        for (int i = 0; i < n; i++) {
500            m.put(i, i, 1.0);
501        }
502
503        return m;
504    }
505
506    /**
507     * Creates a new matrix where the values of the given vector are the diagonal values of
508     * the matrix.
509     */
510    public static DoubleMatrix diag(DoubleMatrix x) {
511        DoubleMatrix m = new DoubleMatrix(x.length, x.length);
512
513        for (int i = 0; i < x.length; i++) {
514            m.put(i, i, x.get(i));
515        }
516
517        return m;
518    }
519
520    /**
521     * Create a 1-by-1 matrix. For many operations, this matrix functions like a
522     * normal double.
523     */
524    public static DoubleMatrix scalar(double s) {
525        DoubleMatrix m = new DoubleMatrix(1, 1);
526        m.put(0, 0, s);
527        return m;
528    }
529
530    /** Test whether a matrix is scalar. */
531    public boolean isScalar() {
532        return length == 1;
533    }
534
535    /** Return the first element of the matrix. */
536    public double scalar() {
537        return get(0);
538    }
539
540    public static DoubleMatrix logspace(double lower, double upper, int size) {
541        DoubleMatrix result = new DoubleMatrix(size);
542        for (int i = 0; i < size; i++) {
543            double t = (double) i / (size - 1);
544            double e = lower * (1 - t) + t * upper;
545            result.put(i, (double) Math.pow(10.0, e));
546        }
547        return result;
548    }
549
550    public static DoubleMatrix linspace(int lower, int upper, int size) {
551        DoubleMatrix result = new DoubleMatrix(size);
552        for (int i = 0; i < size; i++) {
553            double t = (double) i / (size - 1);
554            result.put(i, lower * (1 - t) + t * upper);
555        }
556        return result;
557    }
558
559    /**
560     * Concatenates two matrices horizontally. Matrices must have identical
561     * numbers of rows.
562     */
563    public static DoubleMatrix concatHorizontally(DoubleMatrix A, DoubleMatrix B) {
564        if (A.rows != B.rows) {
565            throw new SizeException("Matrices don't have same number of rows.");
566        }
567
568        DoubleMatrix result = new DoubleMatrix(A.rows, A.columns + B.columns);
569        SimpleBlas.copy(A, result);
570        JavaBlas.rcopy(B.length, B.data, 0, 1, result.data, A.length, 1);
571        return result;
572    }
573
574    /**
575     * Concatenates two matrices vertically. Matrices must have identical
576     * numbers of columns.
577     */
578    public static DoubleMatrix concatVertically(DoubleMatrix A, DoubleMatrix B) {
579        if (A.columns != B.columns) {
580            throw new SizeException("Matrices don't have same number of columns (" + A.columns + " != " + B.columns + ".");
581        }
582
583        DoubleMatrix result = new DoubleMatrix(A.rows + B.rows, A.columns);
584
585        for (int i = 0; i < A.columns; i++) {
586            JavaBlas.rcopy(A.rows, A.data, A.index(0, i), 1, result.data, result.index(0, i), 1);
587            JavaBlas.rcopy(B.rows, B.data, B.index(0, i), 1, result.data, result.index(A.rows, i), 1);
588        }
589
590        return result;
591    }
592
593    /**************************************************************************
594     * Working with slices (Man! 30+ methods just to make this a bit flexible...)
595     */
596    /** Get all elements specified by the linear indices. */
597    public DoubleMatrix get(int[] indices) {
598        DoubleMatrix result = new DoubleMatrix(indices.length);
599
600        for (int i = 0; i < indices.length; i++) {
601            result.put(i, get(indices[i]));
602        }
603
604        return result;
605    }
606
607    /** Get all elements for a given row and the specified columns. */
608    public DoubleMatrix get(int r, int[] indices) {
609        DoubleMatrix result = new DoubleMatrix(1, indices.length);
610
611        for (int i = 0; i < indices.length; i++) {
612            result.put(i, get(r, indices[i]));
613        }
614
615        return result;
616    }
617
618    /** Get all elements for a given column and the specified rows. */
619    public DoubleMatrix get(int[] indices, int c) {
620        DoubleMatrix result = new DoubleMatrix(indices.length, c);
621
622        for (int i = 0; i < indices.length; i++) {
623            result.put(i, get(indices[i], c));
624        }
625
626        return result;
627    }
628
629    /** Get all elements from the specified rows and columns. */
630    public DoubleMatrix get(int[] rindices, int[] cindices) {
631        DoubleMatrix result = new DoubleMatrix(rindices.length, cindices.length);
632
633        for (int i = 0; i < rindices.length; i++) {
634            for (int j = 0; j < cindices.length; j++) {
635                result.put(i, j, get(rindices[i], cindices[j]));
636            }
637        }
638
639        return result;
640    }
641
642    /** Get elements from specified rows and columns. */
643    public DoubleMatrix get(Range rs, Range cs) {
644        rs.init(0, rows);
645        cs.init(0, columns);
646        DoubleMatrix result = new DoubleMatrix(rs.length(), cs.length());
647
648        for (; rs.hasMore(); rs.next()) {
649            for (; cs.hasMore(); cs.next()) {
650                result.put(rs.index(), cs.index(), get(rs.value(), cs.value()));
651            }
652        }
653
654        return result;
655    }
656
657    public DoubleMatrix get(Range rs, int c) {
658        rs.init(0, rows);
659        DoubleMatrix result = new DoubleMatrix(rs.length(), 1);
660
661        for (; rs.hasMore(); rs.next()) {
662            result.put(rs.index(), 0, get(rs.value(), c));
663        }
664
665        return result;
666    }
667
668    public DoubleMatrix get(int r, Range cs) {
669        cs.init(0, columns);
670        DoubleMatrix result = new DoubleMatrix(1, cs.length());
671
672        for (; cs.hasMore(); cs.next()) {
673            result.put(0, cs.index(), get(r, cs.value()));
674        }
675
676        return result;
677
678    }
679
680    /** Get elements specified by the non-zero entries of the passed matrix. */
681    public DoubleMatrix get(DoubleMatrix indices) {
682        return get(indices.findIndices());
683    }
684
685    /**
686     * Get elements from a row and columns as specified by the non-zero entries of
687     * a matrix.
688     */
689    public DoubleMatrix get(int r, DoubleMatrix indices) {
690        return get(r, indices.findIndices());
691    }
692
693    /**
694     * Get elements from a column and rows as specified by the non-zero entries of
695     * a matrix.
696     */
697    public DoubleMatrix get(DoubleMatrix indices, int c) {
698        return get(indices.findIndices(), c);
699    }
700
701    /**
702     * Get elements from columns and rows as specified by the non-zero entries of
703     * the passed matrices.
704     */
705    public DoubleMatrix get(DoubleMatrix rindices, DoubleMatrix cindices) {
706        return get(rindices.findIndices(), cindices.findIndices());
707    }
708
709    /** Return all elements with linear index a, a + 1, ..., b - 1.*/
710    public DoubleMatrix getRange(int a, int b) {
711        DoubleMatrix result = new DoubleMatrix(b - a);
712
713        for (int k = 0; k < b - a; k++) {
714            result.put(k, get(a + k));
715        }
716
717        return result;
718    }
719
720    /** Get elements from a row and columns <tt>a</tt> to <tt>b</tt>. */
721    public DoubleMatrix getColumnRange(int r, int a, int b) {
722        DoubleMatrix result = new DoubleMatrix(1, b - a);
723
724        for (int k = 0; k < b - a; k++) {
725            result.put(k, get(r, a + k));
726        }
727
728        return result;
729    }
730
731    /** Get elements from a column and rows <tt>a/tt> to <tt>b</tt>. */
732    public DoubleMatrix getRowRange(int a, int b, int c) {
733        DoubleMatrix result = new DoubleMatrix(b - a);
734
735        for (int k = 0; k < b - a; k++) {
736            result.put(k, get(a + k, c));
737        }
738
739        return result;
740    }
741
742    /**
743     * Get elements from rows <tt>ra</tt> to <tt>rb</tt> and
744     * columns <tt>ca</tt> to <tt>cb</tt>.
745     */
746    public DoubleMatrix getRange(int ra, int rb, int ca, int cb) {
747        DoubleMatrix result = new DoubleMatrix(rb - ra, cb - ca);
748
749        for (int i = 0; i < rb - ra; i++) {
750            for (int j = 0; j < cb - ca; j++) {
751                result.put(i, j, get(ra + i, ca + j));
752            }
753        }
754
755        return result;
756    }
757
758    /** Get whole rows from the passed indices. */
759    public DoubleMatrix getRows(int[] rindices) {
760        DoubleMatrix result = new DoubleMatrix(rindices.length, columns);
761        for (int i = 0; i < rindices.length; i++) {
762            JavaBlas.rcopy(columns, data, index(rindices[i], 0), rows, result.data, result.index(i, 0), result.rows);
763        }
764        return result;
765    }
766
767    /** Get whole rows as specified by the non-zero entries of a matrix. */
768    public DoubleMatrix getRows(DoubleMatrix rindices) {
769        return getRows(rindices.findIndices());
770    }
771
772    public DoubleMatrix getRows(Range indices, DoubleMatrix result) {
773        indices.init(0, rows);
774        if (result.rows < indices.length()) {
775            throw new SizeException("Result matrix does not have enough rows (" + result.rows + " < " + indices.length() + ")");
776        }
777        result.checkColumns(columns);
778
779        for (int c = 0; c < columns; c++) {
780            indices.init(0, rows);
781            for (int r = 0; indices.hasMore(); indices.next(), r++) {
782                result.put(r, c, get(indices.index(), c));
783            }
784        }
785        return result;
786    }
787
788    public DoubleMatrix getRows(Range indices) {
789        indices.init(0, rows);
790        DoubleMatrix result = new DoubleMatrix(indices.length(), columns);
791        return getRows(indices, result);
792    }
793
794    /** Get whole columns from the passed indices. */
795    public DoubleMatrix getColumns(int[] cindices) {
796        DoubleMatrix result = new DoubleMatrix(rows, cindices.length);
797        for (int i = 0; i < cindices.length; i++) {
798            JavaBlas.rcopy(rows, data, index(0, cindices[i]), 1, result.data, result.index(0, i), 1);
799        }
800        return result;
801    }
802
803    /** Get whole columns as specified by the non-zero entries of a matrix. */
804    public DoubleMatrix getColumns(DoubleMatrix cindices) {
805        return getColumns(cindices.findIndices());
806    }
807
808    /**
809     * Assert that the matrix has a certain length.
810     * @throws SizeException
811     */
812    public void checkLength(int l) {
813        if (length != l) {
814            throw new SizeException("Matrix does not have the necessary length (" + length + " != " + l + ").");
815        }
816    }
817
818    /**
819     * Asserts that the matrix has a certain number of rows.
820     * @throws SizeException
821     */
822    public void checkRows(int r) {
823        if (rows != r) {
824            throw new SizeException("Matrix does not have the necessary number of rows (" + rows + " != " + r + ").");
825        }
826    }
827
828    /**
829     * Asserts that the amtrix has a certain number of columns.
830     * @throws SizeException
831     */
832    public void checkColumns(int c) {
833        if (columns != c) {
834            throw new SizeException("Matrix does not have the necessary number of columns (" + columns + " != " + c + ").");
835        }
836    }
837
838
839    /** Set elements in linear ordering in the specified indices. */
840    public DoubleMatrix put(int[] indices, DoubleMatrix x) {
841        if (x.isScalar()) {
842            return put(indices, x.scalar());
843        }
844        x.checkLength(indices.length);
845
846        for (int i = 0; i < indices.length; i++) {
847            put(indices[i], x.get(i));
848        }
849
850        return this;
851    }
852
853    /** Set multiple elements in a row. */
854    public DoubleMatrix put(int r, int[] indices, DoubleMatrix x) {
855        if (x.isScalar()) {
856            return put(r, indices, x.scalar());
857        }
858        x.checkColumns(indices.length);
859
860        for (int i = 0; i < indices.length; i++) {
861            put(r, indices[i], x.get(i));
862        }
863
864        return this;
865    }
866
867    /** Set multiple elements in a row. */
868    public DoubleMatrix put(int[] indices, int c, DoubleMatrix x) {
869        if (x.isScalar()) {
870            return put(indices, c, x.scalar());
871        }
872        x.checkRows(indices.length);
873
874        for (int i = 0; i < indices.length; i++) {
875            put(indices[i], c, x.get(i));
876        }
877
878        return this;
879    }
880
881    /** Put a sub-matrix as specified by the indices. */
882    public DoubleMatrix put(int[] rindices, int[] cindices, DoubleMatrix x) {
883        if (x.isScalar()) {
884            return put(rindices, cindices, x.scalar());
885        }
886        x.checkRows(rindices.length);
887        x.checkColumns(cindices.length);
888
889        for (int i = 0; i < rindices.length; i++) {
890            for (int j = 0; j < cindices.length; j++) {
891                put(rindices[i], cindices[j], x.get(i, j));
892            }
893        }
894
895        return this;
896    }
897
898    /** Put a matrix into specified indices. */
899    public DoubleMatrix put(Range rs, Range cs, DoubleMatrix x) {
900        rs.init(0, rows);
901        cs.init(0, columns);
902
903        x.checkRows(rs.length());
904        x.checkColumns(cs.length());
905
906        for (; rs.hasMore(); rs.next()) {
907            for (; cs.hasMore(); cs.next()) {
908                put(rs.value(), cs.value(), x.get(rs.index(), cs.index()));
909            }
910        }
911
912        return this;
913    }
914
915    /** Put a single value into the specified indices (linear adressing). */
916    public DoubleMatrix put(int[] indices, double v) {
917        for (int i = 0; i < indices.length; i++) {
918            put(indices[i], v);
919        }
920
921        return this;
922    }
923
924    /** Put a single value into a row and the specified columns. */
925    public DoubleMatrix put(int r, int[] indices, double v) {
926        for (int i = 0; i < indices.length; i++) {
927            put(r, indices[i], v);
928        }
929
930        return this;
931    }
932
933    /** Put a single value into the specified rows of a column. */
934    public DoubleMatrix put(int[] indices, int c, double v) {
935        for (int i = 0; i < indices.length; i++) {
936            put(indices[i], c, v);
937        }
938
939        return this;
940    }
941
942    /** Put a single value into the specified rows and columns. */
943    public DoubleMatrix put(int[] rindices, int[] cindices, double v) {
944        for (int i = 0; i < rindices.length; i++) {
945            for (int j = 0; j < cindices.length; j++) {
946                put(rindices[i], cindices[j], v);
947            }
948        }
949
950        return this;
951    }
952
953    /**
954     * Put a sub-matrix into the indices specified by the non-zero entries
955     * of <tt>indices</tt> (linear adressing).
956     */
957    public DoubleMatrix put(DoubleMatrix indices, DoubleMatrix v) {
958        return put(indices.findIndices(), v);
959    }
960
961    /** Put a sub-vector into the specified columns (non-zero entries of <tt>indices</tt>) of a row. */
962    public DoubleMatrix put(int r, DoubleMatrix indices, DoubleMatrix v) {
963        return put(r, indices.findIndices(), v);
964    }
965
966    /** Put a sub-vector into the specified rows (non-zero entries of <tt>indices</tt>) of a column. */
967    public DoubleMatrix put(DoubleMatrix indices, int c, DoubleMatrix v) {
968        return put(indices.findIndices(), c, v);
969    }
970
971    /**
972     * Put a sub-matrix into the specified rows and columns (non-zero entries of
973     * <tt>rindices</tt> and <tt>cindices</tt>.
974     */
975    public DoubleMatrix put(DoubleMatrix rindices, DoubleMatrix cindices, DoubleMatrix v) {
976        return put(rindices.findIndices(), cindices.findIndices(), v);
977    }
978
979    /**
980     * Put a single value into the elements specified by the non-zero
981     * entries of <tt>indices</tt> (linear adressing).
982     */
983    public DoubleMatrix put(DoubleMatrix indices, double v) {
984        return put(indices.findIndices(), v);
985    }
986
987    /**
988     * Put a single value into the specified columns (non-zero entries of
989     * <tt>indices</tt>) of a row.
990     */
991    public DoubleMatrix put(int r, DoubleMatrix indices, double v) {
992        return put(r, indices.findIndices(), v);
993    }
994
995    /**
996     * Put a single value into the specified rows (non-zero entries of
997     * <tt>indices</tt>) of a column.
998     */
999    public DoubleMatrix put(DoubleMatrix indices, int c, double v) {
1000        return put(indices.findIndices(), c, v);
1001    }
1002
1003    /**
1004     * Put a single value in the specified rows and columns (non-zero entries
1005     * of <tt>rindices</tt> and <tt>cindices</tt>.
1006     */
1007    public DoubleMatrix put(DoubleMatrix rindices, DoubleMatrix cindices, double v) {
1008        return put(rindices.findIndices(), cindices.findIndices(), v);
1009    }
1010
1011    /** Find the linear indices of all non-zero elements. */
1012    public int[] findIndices() {
1013        int len = 0;
1014        for (int i = 0; i < length; i++) {
1015            if (get(i) != 0.0) {
1016                len++;
1017            }
1018        }
1019
1020        int[] indices = new int[len];
1021        int c = 0;
1022
1023        for (int i = 0; i < length; i++) {
1024            if (get(i) != 0.0) {
1025                indices[c++] = i;
1026            }
1027        }
1028
1029        return indices;
1030    }
1031
1032    /**************************************************************************
1033     * Basic operations (copying, resizing, element access)
1034     */
1035    /** Return transposed copy of this matrix. */
1036    public DoubleMatrix transpose() {
1037        DoubleMatrix result = new DoubleMatrix(columns, rows);
1038
1039        for (int i = 0; i < rows; i++) {
1040            for (int j = 0; j < columns; j++) {
1041                result.put(j, i, get(i, j));
1042            }
1043        }
1044
1045        return result;
1046    }
1047
1048    /**
1049     * Compare two matrices. Returns true if and only if other is also a
1050     * DoubleMatrix which has the same size and the maximal absolute
1051     * difference in matrix elements is smaller thatn 1e-6.
1052     */
1053    @Override
1054    public boolean equals(Object o) {
1055        if (!(o instanceof DoubleMatrix)) {
1056            return false;
1057        }
1058
1059        DoubleMatrix other = (DoubleMatrix) o;
1060
1061        if (!sameSize(other)) {
1062            return false;
1063        }
1064
1065        DoubleMatrix diff = MatrixFunctions.absi(sub(other));
1066
1067        return diff.max() / (rows * columns) < 1e-6;
1068    }
1069
1070    @Override
1071    public int hashCode() {
1072        int hash = 7;
1073        hash = 83 * hash + this.rows;
1074        hash = 83 * hash + this.columns;
1075        hash = 83 * hash + Arrays.hashCode(this.data);
1076        return hash;
1077    }
1078    
1079    /** Resize the matrix. All elements will be set to zero. */
1080    public void resize(int newRows, int newColumns) {
1081        rows = newRows;
1082        columns = newColumns;
1083        length = newRows * newColumns;
1084        data = new double[rows * columns];
1085    }
1086
1087    /** Reshape the matrix. Number of elements must not change. */
1088    public DoubleMatrix reshape(int newRows, int newColumns) {
1089        if (length != newRows * newColumns) {
1090            throw new IllegalArgumentException(
1091                    "Number of elements must not change.");
1092        }
1093
1094        rows = newRows;
1095        columns = newColumns;
1096
1097        return this;
1098    }
1099
1100    /** Generate a new matrix which has the given number of replications of this. */
1101    public DoubleMatrix repmat(int rowMult, int columnMult) {
1102        DoubleMatrix result = new DoubleMatrix(rows * rowMult, columns * columnMult);
1103
1104        for (int c = 0; c < columnMult; c++) {
1105            for (int r = 0; r < rowMult; r++) {
1106                for (int i = 0; i < rows; i++) {
1107                    for (int j = 0; j < columns; j++) {
1108                        result.put(r * rows + i, c * columns + j, get(i, j));
1109                    }
1110                }
1111            }
1112        }
1113        return result;
1114    }
1115
1116    /** Checks whether two matrices have the same size. */
1117    public boolean sameSize(DoubleMatrix a) {
1118        return rows == a.rows && columns == a.columns;
1119    }
1120
1121    /** Throws SizeException unless two matrices have the same size. */
1122    public void assertSameSize(DoubleMatrix a) {
1123        if (!sameSize(a)) {
1124            throw new SizeException("Matrices must have the same size.");
1125        }
1126    }
1127
1128    /** Checks whether two matrices can be multiplied (that is, number of columns of
1129     * this must equal number of rows of a. */
1130    public boolean multipliesWith(DoubleMatrix a) {
1131        return columns == a.rows;
1132    }
1133
1134    /** Throws SizeException unless matrices can be multiplied with one another. */
1135    public void assertMultipliesWith(DoubleMatrix a) {
1136        if (!multipliesWith(a)) {
1137            throw new SizeException("Number of columns of left matrix must be equal to number of rows of right matrix.");
1138        }
1139    }
1140
1141    /** Checks whether two matrices have the same length. */
1142    public boolean sameLength(DoubleMatrix a) {
1143        return length == a.length;
1144    }
1145
1146    /** Throws SizeException unless matrices have the same length. */
1147    public void assertSameLength(DoubleMatrix a) {
1148        if (!sameLength(a)) {
1149            throw new SizeException("Matrices must have same length (is: " + length + " and " + a.length + ")");
1150        }
1151    }
1152
1153    /** Copy DoubleMatrix a to this. this a is resized if necessary. */
1154    public DoubleMatrix copy(DoubleMatrix a) {
1155        if (!sameSize(a)) {
1156            resize(a.rows, a.columns);
1157        }
1158
1159        System.arraycopy(a.data, 0, data, 0, length);
1160        return a;
1161    }
1162
1163    /**
1164     * Returns a duplicate of this matrix. Geometry is the same (including offsets, transpose, etc.),
1165     * but the buffer is not shared.
1166     */
1167    public DoubleMatrix dup() {
1168        DoubleMatrix out = new DoubleMatrix(rows, columns);
1169
1170        JavaBlas.rcopy(length, data, 0, 1, out.data, 0, 1);
1171
1172        return out;
1173    }
1174
1175    /** Swap two columns of a matrix. */
1176    public DoubleMatrix swapColumns(int i, int j) {
1177        NativeBlas.dswap(rows, data, index(0, i), 1, data, index(0, j), 1);
1178        return this;
1179    }
1180
1181    /** Swap two rows of a matrix. */
1182    public DoubleMatrix swapRows(int i, int j) {
1183        NativeBlas.dswap(columns, data, index(i, 0), rows, data, index(j, 0), rows);
1184        return this;
1185    }
1186
1187    /** Set matrix element */
1188    public DoubleMatrix put(int rowIndex, int columnIndex, double value) {
1189        data[index(rowIndex, columnIndex)] = value;
1190        return this;
1191    }
1192
1193    /** Retrieve matrix element */
1194    public double get(int rowIndex, int columnIndex) {
1195        return data[index(rowIndex, columnIndex)];
1196    }
1197
1198    /** Get index of an element */
1199    public int index(int rowIndex, int columnIndex) {
1200        return rowIndex + rows * columnIndex;
1201    }
1202
1203    /** Compute the row index of a linear index. */
1204    public int indexRows(int i) {
1205        return i / rows;
1206    }
1207
1208    /** Compute the column index of a linear index. */
1209    public int indexColumns(int i) {
1210        return i - indexRows(i) * rows;
1211    }
1212
1213    /** Get a matrix element (linear indexing). */
1214    public double get(int i) {
1215        return data[i];
1216    }
1217
1218    /** Set a matrix element (linear indexing). */
1219    public DoubleMatrix put(int i, double v) {
1220        data[i] = v;
1221        return this;
1222    }
1223
1224    /** Set all elements to a value. */
1225    public DoubleMatrix fill(double value) {
1226        for (int i = 0; i < length; i++) {
1227            put(i, value);
1228        }
1229        return this;
1230    }
1231
1232    /** Get number of rows. */
1233    public int getRows() {
1234        return rows;
1235    }
1236
1237    /** Get number of columns. */
1238    public int getColumns() {
1239        return columns;
1240    }
1241
1242    /** Get total number of elements. */
1243    public int getLength() {
1244        return length;
1245    }
1246
1247    /** Checks whether the matrix is empty. */
1248    public boolean isEmpty() {
1249        return columns == 0 || rows == 0;
1250    }
1251
1252    /** Checks whether the matrix is square. */
1253    public boolean isSquare() {
1254        return columns == rows;
1255    }
1256
1257    /** Throw SizeException unless matrix is square. */
1258    public void assertSquare() {
1259        if (!isSquare()) {
1260            throw new SizeException("Matrix must be square!");
1261        }
1262    }
1263
1264    /** Checks whether the matrix is a vector. */
1265    public boolean isVector() {
1266        return columns == 1 || rows == 1;
1267    }
1268
1269    /** Checks whether the matrix is a row vector. */
1270    public boolean isRowVector() {
1271        return rows == 1;
1272    }
1273
1274    /** Checks whether the matrix is a column vector. */
1275    public boolean isColumnVector() {
1276        return columns == 1;
1277    }
1278
1279    /** Returns the diagonal of the matrix. */
1280    public DoubleMatrix diag() {
1281        assertSquare();
1282        DoubleMatrix d = new DoubleMatrix(rows);
1283        JavaBlas.rcopy(rows, data, 0, rows + 1, d.data, 0, 1);
1284        return d;
1285    }
1286
1287    /** Pretty-print this matrix to <tt>System.out</tt>. */
1288    public void print() {
1289        System.out.println(toString());
1290    }
1291
1292    /** Generate string representation of the matrix. */
1293    @Override
1294    public String toString() {
1295        StringBuilder s = new StringBuilder();
1296
1297        s.append("[");
1298
1299        for (int i = 0; i < rows; i++) {
1300            for (int j = 0; j < columns; j++) {
1301                s.append(get(i, j));
1302                if (j < columns - 1) {
1303                    s.append(", ");
1304                }
1305            }
1306            if (i < rows - 1) {
1307                s.append("; ");
1308            }
1309        }
1310
1311        s.append("]");
1312
1313        return s.toString();
1314    }
1315
1316    /**
1317     * Generate string representation of the matrix, with specified
1318     * format for the entries. For example, <code>x.toString("%.1f")</code>
1319     * generates a string representations having only one position after the
1320     * decimal point.
1321     */
1322    public String toString(String fmt) {
1323        StringWriter s = new StringWriter();
1324        PrintWriter p = new PrintWriter(s);
1325
1326        p.print("[");
1327
1328        for (int r = 0; r < rows; r++) {
1329            for (int c = 0; c < columns; c++) {
1330                p.printf(fmt, get(r, c));
1331                if (c < columns - 1) {
1332                    p.print(", ");
1333                }
1334            }
1335            if (r < rows - 1) {
1336                p.print("; ");
1337            }
1338        }
1339
1340        p.print("]");
1341
1342        return s.toString();
1343    }
1344
1345    /** Converts the matrix to a one-dimensional array of doubles. */
1346    public double[] toArray() {
1347        double[] array = new double[length];
1348
1349        System.arraycopy(data, 0, array, 0, length);
1350
1351        return array;
1352    }
1353
1354    /** Converts the matrix to a two-dimensional array of doubles. */
1355    public double[][] toArray2() {
1356        double[][] array = new double[rows][columns];
1357
1358        for (int r = 0; r < rows; r++) {
1359            for (int c = 0; c < columns; c++) {
1360                array[r][c] = get(r, c);
1361            }
1362        }
1363
1364        return array;
1365    }
1366
1367    /** Converts the matrix to a one-dimensional array of integers. */
1368    public int[] toIntArray() {
1369        int[] array = new int[length];
1370
1371        for (int i = 0; i < length; i++) {
1372            array[i] = (int) Math.rint(get(i));
1373        }
1374
1375        return array;
1376    }
1377
1378    /** Convert the matrix to a two-dimensional array of integers. */
1379    public int[][] toIntArray2() {
1380        int[][] array = new int[rows][columns];
1381
1382        for (int r = 0; r < rows; r++) {
1383            for (int c = 0; c < columns; c++) {
1384                array[r][c] = (int) Math.rint(get(r, c));
1385            }
1386        }
1387
1388        return array;
1389    }
1390
1391    /** Convert the matrix to a one-dimensional array of boolean values. */
1392    public boolean[] toBooleanArray() {
1393        boolean[] array = new boolean[length];
1394
1395        for (int i = 0; i < length; i++) {
1396            array[i] = get(i) != 0.0 ? true : false;
1397        }
1398
1399        return array;
1400    }
1401
1402    /** Convert the matrix to a two-dimensional array of boolean values. */
1403    public boolean[][] toBooleanArray2() {
1404        boolean[][] array = new boolean[rows][columns];
1405
1406        for (int r = 0; r < rows; r++) {
1407            for (int c = 0; c < columns; c++) {
1408                array[r][c] = get(r, c) != 0.0 ? true : false;
1409            }
1410        }
1411
1412        return array;
1413    }
1414
1415    public FloatMatrix toFloat() {
1416         FloatMatrix result = new FloatMatrix(rows, columns);
1417         for (int i = 0; i < length; i++) {
1418            result.put(i, (float) get(i));
1419         }
1420         return result;
1421    }
1422
1423    /**
1424     * A wrapper which allows to view a matrix as a List of Doubles (read-only!).
1425     * Also implements the {@link ConvertsToDoubleMatrix} interface.
1426     */
1427    public class ElementsAsListView extends AbstractList<Double> implements ConvertsToDoubleMatrix {
1428
1429        private DoubleMatrix me;
1430
1431        public ElementsAsListView(DoubleMatrix me) {
1432            this.me = me;
1433        }
1434
1435        @Override
1436        public Double get(int index) {
1437            return me.get(index);
1438        }
1439
1440        @Override
1441        public int size() {
1442            return me.length;
1443        }
1444
1445        public DoubleMatrix convertToDoubleMatrix() {
1446            return me;
1447        }
1448    }
1449
1450    public class RowsAsListView extends AbstractList<DoubleMatrix> implements ConvertsToDoubleMatrix {
1451
1452        private DoubleMatrix me;
1453
1454        public RowsAsListView(DoubleMatrix me) {
1455            this.me = me;
1456        }
1457
1458        @Override
1459        public DoubleMatrix get(int index) {
1460            return getRow(index);
1461        }
1462
1463        @Override
1464        public int size() {
1465            return rows;
1466        }
1467
1468        public DoubleMatrix convertToDoubleMatrix() {
1469            return me;
1470        }
1471    }
1472
1473    public class ColumnsAsListView extends AbstractList<DoubleMatrix> implements ConvertsToDoubleMatrix {
1474
1475        private DoubleMatrix me;
1476
1477        public ColumnsAsListView(DoubleMatrix me) {
1478            this.me = me;
1479        }
1480
1481        @Override
1482        public DoubleMatrix get(int index) {
1483            return getColumn(index);
1484        }
1485
1486        @Override
1487        public int size() {
1488            return columns;
1489        }
1490
1491        public DoubleMatrix convertToDoubleMatrix() {
1492            return me;
1493        }
1494    }
1495
1496    public List<Double> elementsAsList() {
1497        return new ElementsAsListView(this);
1498    }
1499
1500    public List<DoubleMatrix> rowsAsList() {
1501        return new RowsAsListView(this);
1502    }
1503
1504    public List<DoubleMatrix> columnsAsList() {
1505        return new ColumnsAsListView(this);
1506    }
1507
1508    /**************************************************************************
1509     * Arithmetic Operations
1510     */
1511    /**
1512     * Ensures that the result vector has the same length as this. If not,
1513     * resizing result is tried, which fails if result == this or result == other.
1514     */
1515    private void ensureResultLength(DoubleMatrix other, DoubleMatrix result) {
1516        if (!sameLength(result)) {
1517            if (result == this || result == other) {
1518                throw new SizeException("Cannot resize result matrix because it is used in-place.");
1519            }
1520            result.resize(rows, columns);
1521        }
1522    }
1523
1524    /** Add two matrices (in-place). */
1525    public DoubleMatrix addi(DoubleMatrix other, DoubleMatrix result) {
1526        if (other.isScalar()) {
1527            return addi(other.scalar(), result);
1528        }
1529        if (isScalar()) {
1530            return other.addi(scalar(), result);
1531        }
1532
1533        assertSameLength(other);
1534        ensureResultLength(other, result);
1535
1536        if (result == this) {
1537            SimpleBlas.axpy(1.0, other, result);
1538        } else if (result == other) {
1539            SimpleBlas.axpy(1.0, this, result);
1540        } else {
1541            /*SimpleBlas.copy(this, result);
1542            SimpleBlas.axpy(1.0, other, result);*/
1543            JavaBlas.rzgxpy(length, result.data, data, other.data);
1544        }
1545
1546        return result;
1547    }
1548
1549    /** Add a scalar to a matrix (in-place). */
1550    public DoubleMatrix addi(double v, DoubleMatrix result) {
1551        ensureResultLength(null, result);
1552
1553        for (int i = 0; i < length; i++) {
1554            result.put(i, get(i) + v);
1555        }
1556        return result;
1557    }
1558
1559    /** Subtract two matrices (in-place). */
1560    public DoubleMatrix subi(DoubleMatrix other, DoubleMatrix result) {
1561        if (other.isScalar()) {
1562            return subi(other.scalar(), result);
1563        }
1564        if (isScalar()) {
1565            return other.rsubi(scalar(), result);
1566        }
1567
1568        assertSameLength(other);
1569        ensureResultLength(other, result);
1570
1571        if (result == this) {
1572            SimpleBlas.axpy(-1.0, other, result);
1573        } else if (result == other) {
1574            SimpleBlas.scal(-1.0, result);
1575            SimpleBlas.axpy(1.0, this, result);
1576        } else {
1577            SimpleBlas.copy(this, result);
1578            SimpleBlas.axpy(-1.0, other, result);
1579        }
1580        return result;
1581    }
1582
1583    /** Subtract a scalar from a matrix (in-place). */
1584    public DoubleMatrix subi(double v, DoubleMatrix result) {
1585        ensureResultLength(null, result);
1586
1587        for (int i = 0; i < length; i++) {
1588            result.put(i, get(i) - v);
1589        }
1590        return result;
1591    }
1592
1593    /**
1594     * Subtract two matrices, but subtract first from second matrix, that is,
1595     * compute <em>result = other - this</em> (in-place).
1596     * */
1597    public DoubleMatrix rsubi(DoubleMatrix other, DoubleMatrix result) {
1598        return other.subi(this, result);
1599    }
1600
1601    /** Subtract a matrix from a scalar (in-place). */
1602    public DoubleMatrix rsubi(double a, DoubleMatrix result) {
1603        ensureResultLength(null, result);
1604
1605        for (int i = 0; i < length; i++) {
1606            result.put(i, a - get(i));
1607        }
1608        return result;
1609    }
1610
1611    /** Elementwise multiplication (in-place). */
1612    public DoubleMatrix muli(DoubleMatrix other, DoubleMatrix result) {
1613        if (other.isScalar()) {
1614            return muli(other.scalar(), result);
1615        }
1616        if (isScalar()) {
1617            return other.muli(scalar(), result);
1618        }
1619
1620        assertSameLength(other);
1621        ensureResultLength(other, result);
1622
1623        for (int i = 0; i < length; i++) {
1624            result.put(i, get(i) * other.get(i));
1625        }
1626        return result;
1627    }
1628
1629    /** Elementwise multiplication with a scalar (in-place). */
1630    public DoubleMatrix muli(double v, DoubleMatrix result) {
1631        ensureResultLength(null, result);
1632
1633        for (int i = 0; i < length; i++) {
1634            result.put(i, get(i) * v);
1635        }
1636        return result;
1637    }
1638
1639    /** Matrix-matrix multiplication (in-place). */
1640    public DoubleMatrix mmuli(DoubleMatrix other, DoubleMatrix result) {
1641        if (other.isScalar()) {
1642            return muli(other.scalar(), result);
1643        }
1644        if (isScalar()) {
1645            return other.muli(scalar(), result);
1646        }
1647
1648        /* check sizes and resize if necessary */
1649        assertMultipliesWith(other);
1650        if (result.rows != rows || result.columns != other.columns) {
1651            if (result != this && result != other) {
1652                result.resize(rows, other.columns);
1653            } else {
1654                throw new SizeException("Cannot resize result matrix because it is used in-place.");
1655            }
1656        }
1657
1658        if (result == this || result == other) {
1659            /* actually, blas cannot do multiplications in-place. Therefore, we will fake by
1660             * allocating a temporary object on the side and copy the result later.
1661             */
1662            DoubleMatrix temp = new DoubleMatrix(result.rows, result.columns);
1663            if (other.columns == 1) {
1664                SimpleBlas.gemv(1.0, this, other, 0.0, temp);
1665            } else {
1666                SimpleBlas.gemm(1.0, this, other, 0.0, temp);
1667            }
1668            SimpleBlas.copy(temp, result);
1669        } else {
1670            if (other.columns == 1) {
1671                SimpleBlas.gemv(1.0, this, other, 0.0, result);
1672            } else {
1673                SimpleBlas.gemm(1.0, this, other, 0.0, result);
1674            }
1675        }
1676        return result;
1677    }
1678
1679    /** Matrix-matrix multiplication with a scalar (for symmetry, does the
1680     * same as <code>muli(scalar)</code> (in-place).
1681     */
1682    public DoubleMatrix mmuli(double v, DoubleMatrix result) {
1683        return muli(v, result);
1684    }
1685
1686    /** Elementwise division (in-place). */
1687    public DoubleMatrix divi(DoubleMatrix other, DoubleMatrix result) {
1688        if (other.isScalar()) {
1689            return divi(other.scalar(), result);
1690        }
1691        if (isScalar()) {
1692            return other.rdivi(scalar(), result);
1693        }
1694
1695        assertSameLength(other);
1696        ensureResultLength(other, result);
1697
1698        for (int i = 0; i < length; i++) {
1699            result.put(i, get(i) / other.get(i));
1700        }
1701        return result;
1702    }
1703
1704    /** Elementwise division with a scalar (in-place). */
1705    public DoubleMatrix divi(double a, DoubleMatrix result) {
1706        ensureResultLength(null, result);
1707
1708        for (int i = 0; i < length; i++) {
1709            result.put(i, get(i) / a);
1710        }
1711        return result;
1712    }
1713
1714    /**
1715     * Elementwise division, with operands switched. Computes
1716     * <code>result = other / this</code> (in-place). */
1717    public DoubleMatrix rdivi(DoubleMatrix other, DoubleMatrix result) {
1718        return other.divi(this, result);
1719    }
1720
1721    /** (Elementwise) division with a scalar, with operands switched. Computes
1722     * <code>result = a / this</code> (in-place). */
1723    public DoubleMatrix rdivi(double a, DoubleMatrix result) {
1724        ensureResultLength(null, result);
1725
1726        for (int i = 0; i < length; i++) {
1727            result.put(i, a / get(i));
1728        }
1729        return result;
1730    }
1731
1732    /** Negate each element (in-place). */
1733    public DoubleMatrix negi() {
1734        for (int i = 0; i < length; i++) {
1735            put(i, -get(i));
1736        }
1737        return this;
1738    }
1739
1740    /** Negate each element. */
1741    public DoubleMatrix neg() {
1742        return dup().negi();
1743    }
1744
1745    /** Maps zero to 1.0 and all non-zero values to 0.0 (in-place). */
1746    public DoubleMatrix noti() {
1747        for (int i = 0; i < length; i++) {
1748            put(i, get(i) == 0.0 ? 1.0 : 0.0);
1749        }
1750        return this;
1751    }
1752
1753    /** Maps zero to 1.0 and all non-zero values to 0.0. */
1754    public DoubleMatrix not() {
1755        return dup().noti();
1756    }
1757
1758    /** Maps zero to 0.0 and all non-zero values to 1.0 (in-place). */
1759    public DoubleMatrix truthi() {
1760        for (int i = 0; i < length; i++) {
1761            put(i, get(i) == 0.0 ? 0.0 : 1.0);
1762        }
1763        return this;
1764    }
1765
1766    /** Maps zero to 0.0 and all non-zero values to 1.0. */
1767    public DoubleMatrix truth() {
1768        return dup().truthi();
1769    }
1770
1771    public DoubleMatrix isNaNi() {
1772        for (int i = 0; i < length; i++) {
1773            put(i, Double.isNaN(get(i)) ? 1.0 : 0.0);
1774        }
1775        return this;
1776    }
1777
1778    public DoubleMatrix isNaN() {
1779        return dup().isNaNi();
1780    }
1781
1782    public DoubleMatrix isInfinitei() {
1783        for (int i = 0; i < length; i++) {
1784            put(i, Double.isInfinite(get(i)) ? 1.0 : 0.0);
1785        }
1786        return this;
1787    }
1788
1789    public DoubleMatrix isInfinite() {
1790        return dup().isInfinitei();
1791    }
1792
1793    public DoubleMatrix selecti(DoubleMatrix where) {
1794        checkLength(where.length);
1795        for (int i = 0; i < length; i++) {
1796            if (where.get(i) == 0.0) {
1797                put(i, 0.0);
1798            }
1799        }
1800        return this;
1801    }
1802
1803    public DoubleMatrix select(DoubleMatrix where) {
1804        return dup().selecti(where);
1805    }
1806
1807    /****************************************************************
1808     * Rank one-updates
1809     */
1810    /** Computes a rank-1-update A = A + alpha * x * y'. */
1811    public DoubleMatrix rankOneUpdate(double alpha, DoubleMatrix x, DoubleMatrix y) {
1812        if (rows != x.length) {
1813            throw new SizeException("Vector x has wrong length (" + x.length + " != " + rows + ").");
1814        }
1815        if (columns != y.length) {
1816            throw new SizeException("Vector y has wrong length (" + x.length + " != " + columns + ").");
1817        }
1818
1819        SimpleBlas.ger(alpha, x, y, this);
1820        return this;
1821    }
1822
1823    /** Computes a rank-1-update A = A + alpha * x * x'. */
1824    public DoubleMatrix rankOneUpdate(double alpha, DoubleMatrix x) {
1825        return rankOneUpdate(alpha, x, x);
1826    }
1827
1828    /** Computes a rank-1-update A = A + x * x'. */
1829    public DoubleMatrix rankOneUpdate(DoubleMatrix x) {
1830        return rankOneUpdate(1.0, x, x);
1831    }
1832
1833    /** Computes a rank-1-update A = A + x * y'. */
1834    public DoubleMatrix rankOneUpdate(DoubleMatrix x, DoubleMatrix y) {
1835        return rankOneUpdate(1.0, x, y);
1836    }
1837
1838    /****************************************************************
1839     * Logical operations
1840     */
1841    /** Returns the minimal element of the matrix. */
1842    public double min() {
1843        if (isEmpty()) {
1844            return Double.POSITIVE_INFINITY;
1845        }
1846        double v = Double.POSITIVE_INFINITY;
1847        for (int i = 0; i < length; i++) {
1848            if (!Double.isNaN(get(i)) && get(i) < v) {
1849                v = get(i);
1850            }
1851        }
1852
1853        return v;
1854    }
1855
1856    /**
1857     * Returns the linear index of the minimal element. If there are
1858     * more than one elements with this value, the first one is returned.
1859     */
1860    public int argmin() {
1861        if (isEmpty()) {
1862            return -1;
1863        }
1864        double v = Double.POSITIVE_INFINITY;
1865        int a = -1;
1866        for (int i = 0; i < length; i++) {
1867            if (!Double.isNaN(get(i)) && get(i) < v) {
1868                v = get(i);
1869                a = i;
1870            }
1871        }
1872
1873        return a;
1874    }
1875
1876    /**
1877     * Computes the minimum between two matrices. Returns the smaller of the
1878     * corresponding elements in the matrix (in-place).
1879     */
1880    public DoubleMatrix mini(DoubleMatrix other, DoubleMatrix result) {
1881        if (result == this) {
1882            for (int i = 0; i < length; i++) {
1883                if (get(i) > other.get(i)) {
1884                    put(i, other.get(i));
1885                }
1886            }
1887        } else {
1888            for (int i = 0; i < length; i++) {
1889                if (get(i) > other.get(i)) {
1890                    result.put(i, other.get(i));
1891                } else {
1892                    result.put(i, get(i));
1893                }
1894            }
1895        }
1896        return this;
1897    }
1898
1899    /**
1900     * Computes the minimum between two matrices. Returns the smaller of the
1901     * corresponding elements in the matrix (in-place on this).
1902     */
1903    public DoubleMatrix mini(DoubleMatrix other) {
1904        return mini(other, this);
1905    }
1906
1907    /**
1908     * Computes the minimum between two matrices. Returns the smaller of the
1909     * corresponding elements in the matrix (in-place on this).
1910     */
1911    public DoubleMatrix min(DoubleMatrix other) {
1912        return mini(other, new DoubleMatrix(rows, columns));
1913    }
1914
1915    public DoubleMatrix mini(double v, DoubleMatrix result) {
1916        if (result == this) {
1917            for (int i = 0; i < length; i++) {
1918                if (get(i) > v) {
1919                    result.put(i, v);
1920                }
1921            }
1922        } else {
1923            for (int i = 0; i < length; i++) {
1924                if (get(i) > v) {
1925                    result.put(i, v);
1926                } else {
1927                    result.put(i, get(i));
1928                }
1929            }
1930
1931        }
1932        return this;
1933    }
1934
1935    public DoubleMatrix mini(double v) {
1936        return mini(v, this);
1937    }
1938
1939    public DoubleMatrix min(double v) {
1940        return mini(v, new DoubleMatrix(rows, columns));
1941    }
1942
1943    /** Returns the maximal element of the matrix. */
1944    public double max() {
1945        if (isEmpty()) {
1946            return Double.NEGATIVE_INFINITY;
1947        }
1948        double v = Double.NEGATIVE_INFINITY;
1949        for (int i = 0; i < length; i++) {
1950            if (!Double.isNaN(get(i)) && get(i) > v) {
1951                v = get(i);
1952            }
1953        }
1954        return v;
1955    }
1956
1957    /**
1958     * Returns the linear index of the maximal element of the matrix. If
1959     * there are more than one elements with this value, the first one
1960     * is returned.
1961     */
1962    public int argmax() {
1963        if (isEmpty()) {
1964            return -1;
1965        }
1966        double v = Double.NEGATIVE_INFINITY;
1967        int a = -1;
1968        for (int i = 0; i < length; i++) {
1969            if (!Double.isNaN(get(i)) && get(i) > v) {
1970                v = get(i);
1971                a = i;
1972            }
1973        }
1974
1975        return a;
1976    }
1977
1978    /**
1979     * Computes the maximum between two matrices. Returns the larger of the
1980     * corresponding elements in the matrix (in-place).
1981     */
1982    public DoubleMatrix maxi(DoubleMatrix other, DoubleMatrix result) {
1983        if (result == this) {
1984            for (int i = 0; i < length; i++) {
1985                if (get(i) < other.get(i)) {
1986                    put(i, other.get(i));
1987                }
1988            }
1989        } else {
1990            for (int i = 0; i < length; i++) {
1991                if (get(i) < other.get(i)) {
1992                    result.put(i, other.get(i));
1993                } else {
1994                    result.put(i, get(i));
1995                }
1996            }
1997        }
1998        return this;
1999    }
2000
2001    /**
2002     * Computes the maximum between two matrices. Returns the smaller of the
2003     * corresponding elements in the matrix (in-place on this).
2004     */
2005    public DoubleMatrix maxi(DoubleMatrix other) {
2006        return maxi(other, this);
2007    }
2008
2009    /**
2010     * Computes the maximum between two matrices. Returns the larger of the
2011     * corresponding elements in the matrix (in-place on this).
2012     */
2013    public DoubleMatrix max(DoubleMatrix other) {
2014        return maxi(other, new DoubleMatrix(rows, columns));
2015    }
2016
2017    public DoubleMatrix maxi(double v, DoubleMatrix result) {
2018        if (result == this) {
2019            for (int i = 0; i < length; i++) {
2020                if (get(i) < v) {
2021                    result.put(i, v);
2022                }
2023            }
2024        } else {
2025            for (int i = 0; i < length; i++) {
2026                if (get(i) < v) {
2027                    result.put(i, v);
2028                } else {
2029                    result.put(i, get(i));
2030                }
2031            }
2032
2033        }
2034        return this;
2035    }
2036
2037    public DoubleMatrix maxi(double v) {
2038        return maxi(v, this);
2039    }
2040
2041    public DoubleMatrix max(double v) {
2042        return maxi(v, new DoubleMatrix(rows, columns));
2043    }
2044
2045    /** Computes the sum of all elements of the matrix. */
2046    public double sum() {
2047        double s = 0.0;
2048        for (int i = 0; i < length; i++) {
2049            s += get(i);
2050        }
2051        return s;
2052    }
2053
2054    /** Computes the product of all elements of the matrix */
2055    public double prod() {
2056        double p = 1.0;
2057        for (int i = 0; i < length; i++) {
2058            p *= get(i);
2059        }
2060        return p;
2061    }
2062
2063    /**
2064     * Computes the mean value of all elements in the matrix,
2065     * that is, <code>x.sum() / x.length</code>.
2066     */
2067    public double mean() {
2068        return sum() / length;
2069    }
2070
2071    /**
2072     * Computes the cumulative sum, that is, the sum of all elements
2073     * of the matrix up to a given index in linear addressing (in-place).
2074     */
2075    public DoubleMatrix cumulativeSumi() {
2076        double s = 0.0;
2077        for (int i = 0; i < length; i++) {
2078            s += get(i);
2079            put(i, s);
2080        }
2081        return this;
2082    }
2083
2084    /**
2085     * Computes the cumulative sum, that is, the sum of all elements
2086     * of the matrix up to a given index in linear addressing.
2087     */
2088    public DoubleMatrix cumulativeSum() {
2089        return dup().cumulativeSumi();
2090    }
2091
2092    /** The scalar product of this with other. */
2093    public double dot(DoubleMatrix other) {
2094        return SimpleBlas.dot(this, other);
2095    }
2096
2097    /** 
2098     * Computes the projection coefficient of other on this.
2099     *
2100     * The returned scalar times <tt>this</tt> is the orthogonal projection
2101     * of <tt>other</tt> on <tt>this</tt>.
2102     */
2103    public double project(DoubleMatrix other) {
2104        other.checkLength(length);
2105        double norm = 0, dot = 0;
2106        for (int i = 0; i < this.length; i++) {
2107            double x = get(i);
2108            norm += x * x;
2109            dot += x * other.get(i);
2110        }
2111        return dot / norm;
2112    }
2113
2114    /**
2115     * The Euclidean norm of the matrix as vector, also the Frobenius
2116     * norm of the matrix.
2117     */
2118    public double norm2() {
2119        double norm = 0.0;
2120        for (int i = 0; i < length; i++) {
2121            norm += get(i) * get(i);
2122        }
2123        return (double) Math.sqrt(norm);
2124    }
2125
2126    /**
2127     * The maximum norm of the matrix (maximal absolute value of the elements).
2128     */
2129    public double normmax() {
2130        double max = 0.0;
2131        for (int i = 0; i < length; i++) {
2132            double a = Math.abs(get(i));
2133            if (a > max) {
2134                max = a;
2135            }
2136        }
2137        return max;
2138    }
2139
2140    /**
2141     * The 1-norm of the matrix as vector (sum of absolute values of elements).
2142     */
2143    public double norm1() {
2144        double norm = 0.0;
2145        for (int i = 0; i < length; i++) {
2146            norm += Math.abs(get(i));
2147        }
2148        return norm;
2149    }
2150
2151    /**
2152     * Returns the squared (Euclidean) distance.
2153     */
2154    public double squaredDistance(DoubleMatrix other) {
2155        other.checkLength(length);
2156        double sd = 0.0;
2157        for (int i = 0; i < length; i++) {
2158            double d = get(i) - other.get(i);
2159            sd += d * d;
2160        }
2161        return sd;
2162    }
2163
2164    /**
2165     * Returns the (euclidean) distance.
2166     */
2167    public double distance2(DoubleMatrix other) {
2168        return (double) Math.sqrt(squaredDistance(other));
2169    }
2170
2171    /**
2172     * Returns the (1-norm) distance.
2173     */
2174    public double distance1(DoubleMatrix other) {
2175        other.checkLength(length);
2176        double d = 0.0;
2177        for (int i = 0; i < length; i++) {
2178            d += Math.abs(get(i) - other.get(i));
2179        }
2180        return d;
2181    }
2182
2183    /**
2184     * Return a new matrix with all elements sorted.
2185     */
2186    public DoubleMatrix sort() {
2187        double array[] = toArray();
2188        java.util.Arrays.sort(array);
2189        return new DoubleMatrix(rows, columns, array);
2190    }
2191
2192    /**
2193     * Sort elements in-place.
2194     */
2195    public DoubleMatrix sorti() {
2196        Arrays.sort(data);
2197        return this;
2198    }
2199
2200    /**
2201     * Get the sorting permutation.
2202     *
2203     * @return an int[] array such that which indexes the elements in sorted
2204     * order.
2205     */
2206    public int[] sortingPermutation() {
2207        Integer[] indices = new Integer[length];
2208
2209        for (int i = 0; i < length; i++) {
2210            indices[i] = i;
2211        }
2212
2213        final double[] array = data;
2214
2215        Arrays.sort(indices, new Comparator() {
2216
2217            public int compare(Object o1, Object o2) {
2218                int i = (Integer) o1;
2219                int j = (Integer) o2;
2220                if (array[i] < array[j]) {
2221                    return -1;
2222                } else if (array[i] == array[j]) {
2223                    return 0;
2224                } else {
2225                    return 1;
2226                }
2227            }
2228        });
2229
2230        int[] result = new int[length];
2231
2232        for (int i = 0; i < length; i++) {
2233            result[i] = indices[i];
2234        }
2235
2236        return result;
2237    }
2238
2239    /**
2240     * Sort columns (in-place).
2241     */
2242    public DoubleMatrix sortColumnsi() {
2243        for (int i = 0; i < length; i += rows) {
2244            Arrays.sort(data, i, i + rows);
2245        }
2246        return this;
2247    }
2248
2249    /** Sort columns. */
2250    public DoubleMatrix sortColumns() {
2251        return dup().sortColumnsi();
2252    }
2253
2254    /** Return matrix of indices which sort all columns. */
2255    public int[][] columnSortingPermutations() {
2256        int[][] result = new int[columns][];
2257
2258        DoubleMatrix temp = new DoubleMatrix(rows);
2259        for (int c = 0; c < columns; c++) {
2260            result[c] = getColumn(c, temp).sortingPermutation();
2261        }
2262
2263        return result;
2264    }
2265
2266    /** Sort rows (in-place). */
2267    public DoubleMatrix sortRowsi() {
2268        // actually, this is much harder because the data is not consecutive
2269        // in memory...
2270        DoubleMatrix temp = new DoubleMatrix(columns);
2271        for (int r = 0; r < rows; r++) {
2272            putRow(r, getRow(r, temp).sorti());
2273        }
2274        return this;
2275    }
2276
2277    /** Sort rows. */
2278    public DoubleMatrix sortRows() {
2279        return dup().sortRowsi();
2280    }
2281
2282    /** Return matrix of indices which sort all columns. */
2283    public int[][] rowSortingPermutations() {
2284        int[][] result = new int[rows][];
2285
2286        DoubleMatrix temp = new DoubleMatrix(columns);
2287        for (int r = 0; r < rows; r++) {
2288            result[r] = getRow(r, temp).sortingPermutation();
2289        }
2290
2291        return result;
2292    }
2293
2294    /** Return a vector containing the sums of the columns (having number of columns many entries) */
2295    public DoubleMatrix columnSums() {
2296        if (rows == 1) {
2297            return dup();
2298        } else {
2299            DoubleMatrix v = new DoubleMatrix(1, columns);
2300
2301            for (int c = 0; c < columns; c++) {
2302                for (int r = 0; r < rows; r++) {
2303                    v.put(c, v.get(c) + get(r, c));
2304                }
2305            }
2306
2307            return v;
2308        }
2309    }
2310
2311    /** Return a vector containing the means of all columns. */
2312    public DoubleMatrix columnMeans() {
2313        return columnSums().divi(rows);
2314    }
2315
2316    /** Return a vector containing the sum of the rows. */
2317    public DoubleMatrix rowSums() {
2318        if (columns == 1) {
2319            return dup();
2320        } else {
2321            DoubleMatrix v = new DoubleMatrix(rows);
2322
2323            for (int c = 0; c < columns; c++) {
2324                for (int r = 0; r < rows; r++) {
2325                    v.put(r, v.get(r) + get(r, c));
2326                }
2327            }
2328
2329            return v;
2330        }
2331    }
2332
2333    /** Return a vector containing the means of the rows. */
2334    public DoubleMatrix rowMeans() {
2335        return rowSums().divi(columns);
2336    }
2337
2338    /************************************************************************
2339     * Column and rows access.
2340     */
2341
2342    /** Get a copy of a column. */
2343    public DoubleMatrix getColumn(int c) {
2344        return getColumn(c, new DoubleMatrix(rows, 1));
2345    }
2346
2347    /** Copy a column to the given vector. */
2348    public DoubleMatrix getColumn(int c, DoubleMatrix result) {
2349        result.checkLength(rows);
2350        JavaBlas.rcopy(rows, data, index(0, c), 1, result.data, 0, 1);
2351        return result;
2352    }
2353
2354    /** Copy a column back into the matrix. */
2355    public void putColumn(int c, DoubleMatrix v) {
2356        JavaBlas.rcopy(rows, v.data, 0, 1, data, index(0, c), 1);
2357    }
2358
2359    /** Get a copy of a row. */
2360    public DoubleMatrix getRow(int r) {
2361        return getRow(r, new DoubleMatrix(1, columns));
2362    }
2363
2364    /** Copy a row to a given vector. */
2365    public DoubleMatrix getRow(int r, DoubleMatrix result) {
2366        result.checkLength(columns);
2367        JavaBlas.rcopy(columns, data, index(r, 0), rows, result.data, 0, 1);
2368        return result;
2369    }
2370
2371    /** Copy a row back into the matrix. */
2372    public void putRow(int r, DoubleMatrix v) {
2373        JavaBlas.rcopy(columns, v.data, 0, 1, data, index(r, 0), rows);
2374    }
2375
2376    /** Return column-wise minimums. */
2377    public DoubleMatrix columnMins() {
2378        DoubleMatrix mins = new DoubleMatrix(1, columns);
2379        for (int c = 0; c < columns; c++) {
2380            mins.put(c, getColumn(c).min());
2381        }
2382        return mins;
2383    }
2384
2385    /** Return index of minimal element per column. */
2386    public int[] columnArgmins() {
2387        int[] argmins = new int[columns];
2388        for (int c = 0; c < columns; c++) {
2389            argmins[c] = getColumn(c).argmin();
2390        }
2391        return argmins;
2392    }
2393
2394    /** Return column-wise maximums. */
2395    public DoubleMatrix columnMaxs() {
2396        DoubleMatrix maxs = new DoubleMatrix(1, columns);
2397        for (int c = 0; c < columns; c++) {
2398            maxs.put(c, getColumn(c).max());
2399        }
2400        return maxs;
2401    }
2402
2403    /** Return index of minimal element per column. */
2404    public int[] columnArgmaxs() {
2405        int[] argmaxs = new int[columns];
2406        for (int c = 0; c < columns; c++) {
2407            argmaxs[c] = getColumn(c).argmax();
2408        }
2409        return argmaxs;
2410    }
2411
2412    /** Return row-wise minimums. */
2413    public DoubleMatrix rowMins() {
2414        DoubleMatrix mins = new DoubleMatrix(rows);
2415        for (int c = 0; c < rows; c++) {
2416            mins.put(c, getRow(c).min());
2417        }
2418        return mins;
2419    }
2420
2421    /** Return index of minimal element per row. */
2422    public int[] rowArgmins() {
2423        int[] argmins = new int[rows];
2424        for (int c = 0; c < rows; c++) {
2425            argmins[c] = getRow(c).argmin();
2426        }
2427        return argmins;
2428    }
2429
2430    /** Return row-wise maximums. */
2431    public DoubleMatrix rowMaxs() {
2432        DoubleMatrix maxs = new DoubleMatrix(rows);
2433        for (int c = 0; c < rows; c++) {
2434            maxs.put(c, getRow(c).max());
2435        }
2436        return maxs;
2437    }
2438
2439    /** Return index of minimal element per row. */
2440    public int[] rowArgmaxs() {
2441        int[] argmaxs = new int[rows];
2442        for (int c = 0; c < rows; c++) {
2443            argmaxs[c] = getRow(c).argmax();
2444        }
2445        return argmaxs;
2446    }
2447
2448    /**************************************************************************
2449     * Elementwise Functions
2450     */
2451    /** Add a row vector to all rows of the matrix (in place). */
2452    public DoubleMatrix addiRowVector(DoubleMatrix x) {
2453        x.checkLength(columns);
2454        for (int c = 0; c < columns; c++) {
2455            for (int r = 0; r < rows; r++) {
2456                put(r, c, get(r, c) + x.get(c));
2457            }
2458        }
2459        return this;
2460    }
2461
2462    /** Add a row to all rows of the matrix. */
2463    public DoubleMatrix addRowVector(DoubleMatrix x) {
2464        return dup().addiRowVector(x);
2465    }
2466
2467    /** Add a vector to all columns of the matrix (in-place). */
2468    public DoubleMatrix addiColumnVector(DoubleMatrix x) {
2469        x.checkLength(rows);
2470        for (int c = 0; c < columns; c++) {
2471            for (int r = 0; r < rows; r++) {
2472                put(r, c, get(r, c) + x.get(r));
2473            }
2474        }
2475        return this;
2476    }
2477
2478    /** Add a vector to all columns of the matrix. */
2479    public DoubleMatrix addColumnVector(DoubleMatrix x) {
2480        return dup().addiColumnVector(x);
2481    }
2482
2483    /** Subtract a row vector from all rows of the matrix (in-place). */
2484    public DoubleMatrix subiRowVector(DoubleMatrix x) {
2485        // This is a bit crazy, but a row vector must have as length as the columns of the matrix.
2486        x.checkLength(columns);
2487        for (int c = 0; c < columns; c++) {
2488            for (int r = 0; r < rows; r++) {
2489                put(r, c, get(r, c) - x.get(c));
2490            }
2491        }
2492        return this;
2493    }
2494
2495    /** Subtract a row vector from all rows of the matrix. */
2496    public DoubleMatrix subRowVector(DoubleMatrix x) {
2497        return dup().subiRowVector(x);
2498    }
2499
2500    /** Subtract a column vector from all columns of the matrix (in-place). */
2501    public DoubleMatrix subiColumnVector(DoubleMatrix x) {
2502        x.checkLength(rows);
2503        for (int c = 0; c < columns; c++) {
2504            for (int r = 0; r < rows; r++) {
2505                put(r, c, get(r, c) - x.get(r));
2506            }
2507        }
2508        return this;
2509    }
2510
2511    /** Subtract a vector from all columns of the matrix. */
2512    public DoubleMatrix subColumnVector(DoubleMatrix x) {
2513        return dup().subiColumnVector(x);
2514    }
2515
2516    /** Multiply a row by a scalar. */
2517    public DoubleMatrix mulRow(int r, double scale) {
2518        NativeBlas.dscal(columns, scale, data, index(r, 0), rows);
2519        return this;
2520    }
2521
2522    /** Multiply a column by a scalar. */
2523    public DoubleMatrix mulColumn(int c, double scale) {
2524        NativeBlas.dscal(rows, scale, data, index(0, c), 1);
2525        return this;
2526    }
2527
2528    /** Multiply all columns with a column vector (in-place). */
2529    public DoubleMatrix muliColumnVector(DoubleMatrix x) {
2530        x.checkLength(rows);
2531        for (int c = 0; c < columns; c++) {
2532            for (int r = 0; r < rows; r++) {
2533                put(r, c, get(r, c) * x.get(r));
2534            }
2535        }
2536        return this;
2537    }
2538
2539    /** Multiply all columns with a column vector. */
2540    public DoubleMatrix mulColumnVector(DoubleMatrix x) {
2541        return dup().muliColumnVector(x);
2542    }
2543
2544    /** Multiply all rows with a row vector (in-place). */
2545    public DoubleMatrix muliRowVector(DoubleMatrix x) {
2546        x.checkLength(columns);
2547        for (int c = 0; c < columns; c++) {
2548            for (int r = 0; r < rows; r++) {
2549                put(r, c, get(r, c) * x.get(c));
2550            }
2551        }
2552        return this;
2553    }
2554
2555    /** Multiply all rows with a row vector. */
2556    public DoubleMatrix mulRowVector(DoubleMatrix x) {
2557        return dup().muliRowVector(x);
2558    }
2559
2560    public DoubleMatrix diviRowVector(DoubleMatrix x) {
2561        x.checkLength(columns);
2562        for (int c = 0; c < columns; c++) {
2563            for (int r = 0; r < rows; r++) {
2564                put(r, c, get(r, c) / x.get(c));
2565            }
2566        }
2567        return this;
2568    }
2569
2570    public DoubleMatrix divRowVector(DoubleMatrix x) {
2571        return dup().diviRowVector(x);
2572    }
2573
2574    public DoubleMatrix diviColumnVector(DoubleMatrix x) {
2575        x.checkLength(rows);
2576        for (int c = 0; c < columns; c++) {
2577            for (int r = 0; r < rows; r++) {
2578                put(r, c, get(r, c) / x.get(r));
2579            }
2580        }
2581        return this;
2582    }
2583
2584    public DoubleMatrix divColumnVector(DoubleMatrix x) {
2585        return dup().diviColumnVector(x);
2586    }
2587
2588    /**
2589     * Writes out this matrix to the given data stream.
2590     * @param dos the data output stream to write to.
2591     * @throws IOException
2592     */
2593    public void out(DataOutputStream dos) throws IOException {
2594        dos.writeUTF("double");
2595        dos.writeInt(columns);
2596        dos.writeInt(rows);
2597
2598        dos.writeInt(data.length);
2599        for (int i = 0; i < data.length; i++) {
2600            dos.writeDouble(data[i]);
2601        }
2602    }
2603
2604    /**
2605     * Reads in a matrix from the given data stream. Note
2606     * that the old data of this matrix will be discarded.
2607     * @param dis the data input stream to read from.
2608     * @throws IOException
2609     */
2610    public void in(DataInputStream dis) throws IOException {
2611        if (!dis.readUTF().equals("double")) {
2612            throw new IllegalStateException("The matrix in the specified file is not of the correct type!");
2613        }
2614
2615        this.columns = dis.readInt();
2616        this.rows = dis.readInt();
2617
2618        final int MAX = dis.readInt();
2619        data = new double[MAX];
2620        for (int i = 0; i < MAX; i++) {
2621            data[i] = dis.readDouble();
2622        }
2623    }
2624
2625    /**
2626     * Saves this matrix to the specified file.
2627     * @param filename the file to write the matrix in.
2628     * @throws IOException thrown on errors while writing the matrix to the file
2629     */
2630    public void save(String filename) throws IOException {
2631        DataOutputStream dos = new DataOutputStream(new FileOutputStream(filename, false));
2632        this.out(dos);
2633    }
2634
2635    /**
2636     * Loads a matrix from a file into this matrix. Note that the old data
2637     * of this matrix will be discarded.
2638     * @param filename the file to read the matrix from
2639     * @throws IOException thrown on errors while reading the matrix
2640     */
2641    public void load(String filename) throws IOException {
2642        DataInputStream dis = new DataInputStream(new FileInputStream(filename));
2643        this.in(dis);
2644    }
2645
2646    public static DoubleMatrix loadAsciiFile(String filename) throws IOException {
2647        BufferedReader is = new BufferedReader(new InputStreamReader(new FileInputStream(filename)));
2648
2649        // Go through file and count columns and rows. What makes this endeavour a bit difficult is
2650        // that files can have leading or trailing spaces leading to spurious fields
2651        // after String.split().
2652        String line;
2653        int rows = 0;
2654        int columns = -1;
2655        while ((line = is.readLine()) != null) {
2656            String[] elements = line.split("\\s+");
2657            int numElements = elements.length;
2658            if (elements[0].length() == 0) {
2659                numElements--;
2660            }
2661            if (elements[elements.length - 1].length() == 0) {
2662                numElements--;
2663            }
2664
2665            if (columns == -1) {
2666                columns = numElements;
2667            } else {
2668                if (columns != numElements) {
2669                    throw new IOException("Number of elements changes in line " + line + ".");
2670                }
2671            }
2672
2673            rows++;
2674        }
2675        is.close();
2676
2677        // Go through file a second time process the actual data.
2678        is = new BufferedReader(new InputStreamReader(new FileInputStream(filename)));
2679        DoubleMatrix result = new DoubleMatrix(rows, columns);
2680        int r = 0;
2681        while ((line = is.readLine()) != null) {
2682            String[] elements = line.split("\\s+");
2683            int firstElement = (elements[0].length() == 0) ? 1 : 0;
2684            for (int c = 0, cc = firstElement; c < columns; c++, cc++) {
2685                result.put(r, c, Double.valueOf(elements[cc]));
2686            }
2687            r++;
2688        }
2689        return result;
2690    }
2691
2692    public static DoubleMatrix loadCSVFile(String filename) throws IOException {
2693        BufferedReader is = new BufferedReader(new InputStreamReader(new FileInputStream(filename)));
2694
2695        List<DoubleMatrix> rows = new LinkedList<DoubleMatrix>();
2696        String line;
2697        int columns = -1;
2698        while ((line = is.readLine()) != null) {
2699            String[] elements = line.split(",");
2700            int numElements = elements.length;
2701            if (elements[0].length() == 0) {
2702                numElements--;
2703            }
2704            if (elements[elements.length - 1].length() == 0) {
2705                numElements--;
2706            }
2707
2708            if (columns == -1) {
2709                columns = numElements;
2710            } else {
2711                if (columns != numElements) {
2712                    throw new IOException("Number of elements changes in line " + line + ".");
2713                }
2714            }
2715
2716            DoubleMatrix row = new DoubleMatrix(columns);
2717            for (int c = 0; c < columns; c++) {
2718                row.put(c, Double.valueOf(elements[c]));
2719            }
2720            rows.add(row);
2721        }
2722        is.close();
2723
2724        System.out.println("Done reading file");
2725
2726        DoubleMatrix result = new DoubleMatrix(rows.size(), columns);
2727        int r = 0;
2728        Iterator<DoubleMatrix> ri = rows.iterator();
2729        while (ri.hasNext()) {
2730            result.putRow(r, ri.next());
2731            r++;
2732        }
2733        return result;
2734    }
2735
2736    /****************************************************************
2737     * Autogenerated code
2738     */
2739    /***** Code for operators ***************************************/
2740
2741    /* Overloads for the usual arithmetic operations */
2742    /*#
2743    def gen_overloads(base, result_rows, result_cols, verb=''); <<-EOS
2744    #{doc verb.capitalize + " a matrix (in place)."}
2745    public DoubleMatrix #{base}i(DoubleMatrix other) {
2746    return #{base}i(other, this);
2747    }
2748
2749    #{doc verb.capitalize + " a matrix (in place)."}
2750    public DoubleMatrix #{base}(DoubleMatrix other) {
2751    return #{base}i(other, new DoubleMatrix(#{result_rows}, #{result_cols}));
2752    }
2753
2754    #{doc verb.capitalize + " a scalar (in place)."}
2755    public DoubleMatrix #{base}i(double v) {
2756    return #{base}i(v, this);
2757    }
2758
2759    #{doc verb.capitalize + " a scalar."}
2760    public DoubleMatrix #{base}(double v) {
2761    return #{base}i(v, new DoubleMatrix(rows, columns));
2762    }
2763    EOS
2764    end
2765    #*/
2766
2767    /* Generating code for logical operators. This not only generates the stubs
2768     * but really all of the code.
2769     */
2770    /*#
2771    def gen_compare(name, op, cmp); <<-EOS
2772    #{doc 'Test for ' + cmp + ' (in-place).'}
2773    public DoubleMatrix #{name}i(DoubleMatrix other, DoubleMatrix result) {
2774    if (other.isScalar())
2775    return #{name}i(other.scalar(), result);
2776
2777    assertSameLength(other);
2778    ensureResultLength(other, result);
2779
2780    for (int i = 0; i < length; i++)
2781    result.put(i, get(i) #{op} other.get(i) ? 1.0 : 0.0);
2782    return result;
2783    }
2784
2785    #{doc 'Test for ' + cmp + ' (in-place).'}
2786    public DoubleMatrix #{name}i(DoubleMatrix other) {
2787    return #{name}i(other, this);
2788    }
2789
2790    #{doc 'Test for ' + cmp + '.'}
2791    public DoubleMatrix #{name}(DoubleMatrix other) {
2792    return #{name}i(other, new DoubleMatrix(rows, columns));
2793    }
2794
2795    #{doc 'Test for ' + cmp + ' against a scalar (in-place).'}
2796    public DoubleMatrix #{name}i(double value, DoubleMatrix result) {
2797    ensureResultLength(null, result);
2798    for (int i = 0; i < length; i++)
2799    result.put(i, get(i) #{op} value ? 1.0 : 0.0);
2800    return result;
2801    }
2802
2803    #{doc 'Test for ' + cmp + ' against a scalar (in-place).'}
2804    public DoubleMatrix #{name}i(double value) {
2805    return #{name}i(value, this);
2806    }
2807
2808    #{doc 'test for ' + cmp + ' against a scalar.'}
2809    public DoubleMatrix #{name}(double value) {
2810    return #{name}i(value, new DoubleMatrix(rows, columns));
2811    }
2812    EOS
2813    end
2814    #*/
2815    /*#
2816    def gen_logical(name, op, cmp); <<-EOS
2817    #{doc 'Compute elementwise ' + cmp + ' (in-place).'}
2818    public DoubleMatrix #{name}i(DoubleMatrix other, DoubleMatrix result) {
2819    assertSameLength(other);
2820    ensureResultLength(other, result);
2821
2822    for (int i = 0; i < length; i++)
2823    result.put(i, (get(i) != 0.0) #{op} (other.get(i) != 0.0) ? 1.0 : 0.0);
2824    return result;
2825    }
2826
2827    #{doc 'Compute elementwise ' + cmp + ' (in-place).'}
2828    public DoubleMatrix #{name}i(DoubleMatrix other) {
2829    return #{name}i(other, this);
2830    }
2831
2832    #{doc 'Compute elementwise ' + cmp + '.'}
2833    public DoubleMatrix #{name}(DoubleMatrix other) {
2834    return #{name}i(other, new DoubleMatrix(rows, columns));
2835    }
2836
2837    #{doc 'Compute elementwise ' + cmp + ' against a scalar (in-place).'}
2838    public DoubleMatrix #{name}i(double value, DoubleMatrix result) {
2839    ensureResultLength(null, result);
2840    boolean val = (value != 0.0);
2841    for (int i = 0; i < length; i++)
2842    result.put(i, (get(i) != 0.0) #{op} val ? 1.0 : 0.0);
2843    return result;
2844    }
2845
2846    #{doc 'Compute elementwise ' + cmp + ' against a scalar (in-place).'}
2847    public DoubleMatrix #{name}i(double value) {
2848    return #{name}i(value, this);
2849    }
2850
2851    #{doc 'Compute elementwise ' + cmp + ' against a scalar.'}
2852    public DoubleMatrix #{name}(double value) {
2853    return #{name}i(value, new DoubleMatrix(rows, columns));
2854    }
2855    EOS
2856    end
2857    #*/
2858
2859    /*# collect(gen_overloads('add', 'rows', 'columns', 'add'),
2860    gen_overloads('sub', 'rows', 'columns', 'subtract'),
2861    gen_overloads('rsub', 'rows', 'columns', '(right-)subtract'),
2862    gen_overloads('div', 'rows', 'columns', 'elementwise divide by'),
2863    gen_overloads('rdiv', 'rows', 'columns', '(right-)elementwise divide by'),
2864    gen_overloads('mul', 'rows', 'columns', 'elementwise multiply by'),
2865    gen_overloads('mmul', 'rows', 'other.columns', 'matrix-multiply by'),
2866    gen_compare('lt', '<', '"less than"'),
2867    gen_compare('gt', '>', '"greater than"'),
2868    gen_compare('le', '<=', '"less than or equal"'),
2869    gen_compare('ge', '>=', '"greater than or equal"'),
2870    gen_compare('eq', '==', 'equality'),
2871    gen_compare('ne', '!=', 'inequality'),
2872    gen_logical('and', '&', 'logical and'),
2873    gen_logical('or', '|', 'logical or'),
2874    gen_logical('xor', '^', 'logical xor'))
2875    #*/
2876//RJPP-BEGIN------------------------------------------------------------
2877    /** Add a matrix (in place). */
2878    public DoubleMatrix addi(DoubleMatrix other) {
2879    return addi(other, this);
2880    }
2881
2882    /** Add a matrix (in place). */
2883    public DoubleMatrix add(DoubleMatrix other) {
2884    return addi(other, new DoubleMatrix(rows, columns));
2885    }
2886
2887    /** Add a scalar (in place). */
2888    public DoubleMatrix addi(double v) {
2889    return addi(v, this);
2890    }
2891
2892    /** Add a scalar. */
2893    public DoubleMatrix add(double v) {
2894    return addi(v, new DoubleMatrix(rows, columns));
2895    }
2896
2897    /** Subtract a matrix (in place). */
2898    public DoubleMatrix subi(DoubleMatrix other) {
2899    return subi(other, this);
2900    }
2901
2902    /** Subtract a matrix (in place). */
2903    public DoubleMatrix sub(DoubleMatrix other) {
2904    return subi(other, new DoubleMatrix(rows, columns));
2905    }
2906
2907    /** Subtract a scalar (in place). */
2908    public DoubleMatrix subi(double v) {
2909    return subi(v, this);
2910    }
2911
2912    /** Subtract a scalar. */
2913    public DoubleMatrix sub(double v) {
2914    return subi(v, new DoubleMatrix(rows, columns));
2915    }
2916
2917    /** (right-)subtract a matrix (in place). */
2918    public DoubleMatrix rsubi(DoubleMatrix other) {
2919    return rsubi(other, this);
2920    }
2921
2922    /** (right-)subtract a matrix (in place). */
2923    public DoubleMatrix rsub(DoubleMatrix other) {
2924    return rsubi(other, new DoubleMatrix(rows, columns));
2925    }
2926
2927    /** (right-)subtract a scalar (in place). */
2928    public DoubleMatrix rsubi(double v) {
2929    return rsubi(v, this);
2930    }
2931
2932    /** (right-)subtract a scalar. */
2933    public DoubleMatrix rsub(double v) {
2934    return rsubi(v, new DoubleMatrix(rows, columns));
2935    }
2936
2937    /** Elementwise divide by a matrix (in place). */
2938    public DoubleMatrix divi(DoubleMatrix other) {
2939    return divi(other, this);
2940    }
2941
2942    /** Elementwise divide by a matrix (in place). */
2943    public DoubleMatrix div(DoubleMatrix other) {
2944    return divi(other, new DoubleMatrix(rows, columns));
2945    }
2946
2947    /** Elementwise divide by a scalar (in place). */
2948    public DoubleMatrix divi(double v) {
2949    return divi(v, this);
2950    }
2951
2952    /** Elementwise divide by a scalar. */
2953    public DoubleMatrix div(double v) {
2954    return divi(v, new DoubleMatrix(rows, columns));
2955    }
2956
2957    /** (right-)elementwise divide by a matrix (in place). */
2958    public DoubleMatrix rdivi(DoubleMatrix other) {
2959    return rdivi(other, this);
2960    }
2961
2962    /** (right-)elementwise divide by a matrix (in place). */
2963    public DoubleMatrix rdiv(DoubleMatrix other) {
2964    return rdivi(other, new DoubleMatrix(rows, columns));
2965    }
2966
2967    /** (right-)elementwise divide by a scalar (in place). */
2968    public DoubleMatrix rdivi(double v) {
2969    return rdivi(v, this);
2970    }
2971
2972    /** (right-)elementwise divide by a scalar. */
2973    public DoubleMatrix rdiv(double v) {
2974    return rdivi(v, new DoubleMatrix(rows, columns));
2975    }
2976
2977    /** Elementwise multiply by a matrix (in place). */
2978    public DoubleMatrix muli(DoubleMatrix other) {
2979    return muli(other, this);
2980    }
2981
2982    /** Elementwise multiply by a matrix (in place). */
2983    public DoubleMatrix mul(DoubleMatrix other) {
2984    return muli(other, new DoubleMatrix(rows, columns));
2985    }
2986
2987    /** Elementwise multiply by a scalar (in place). */
2988    public DoubleMatrix muli(double v) {
2989    return muli(v, this);
2990    }
2991
2992    /** Elementwise multiply by a scalar. */
2993    public DoubleMatrix mul(double v) {
2994    return muli(v, new DoubleMatrix(rows, columns));
2995    }
2996
2997    /** Matrix-multiply by a matrix (in place). */
2998    public DoubleMatrix mmuli(DoubleMatrix other) {
2999    return mmuli(other, this);
3000    }
3001
3002    /** Matrix-multiply by a matrix (in place). */
3003    public DoubleMatrix mmul(DoubleMatrix other) {
3004    return mmuli(other, new DoubleMatrix(rows, other.columns));
3005    }
3006
3007    /** Matrix-multiply by a scalar (in place). */
3008    public DoubleMatrix mmuli(double v) {
3009    return mmuli(v, this);
3010    }
3011
3012    /** Matrix-multiply by a scalar. */
3013    public DoubleMatrix mmul(double v) {
3014    return mmuli(v, new DoubleMatrix(rows, columns));
3015    }
3016
3017    /** Test for "less than" (in-place). */
3018    public DoubleMatrix lti(DoubleMatrix other, DoubleMatrix result) {
3019    if (other.isScalar())
3020    return lti(other.scalar(), result);
3021
3022    assertSameLength(other);
3023    ensureResultLength(other, result);
3024
3025    for (int i = 0; i < length; i++)
3026    result.put(i, get(i) < other.get(i) ? 1.0 : 0.0);
3027    return result;
3028    }
3029
3030    /** Test for "less than" (in-place). */
3031    public DoubleMatrix lti(DoubleMatrix other) {
3032    return lti(other, this);
3033    }
3034
3035    /** Test for "less than". */
3036    public DoubleMatrix lt(DoubleMatrix other) {
3037    return lti(other, new DoubleMatrix(rows, columns));
3038    }
3039
3040    /** Test for "less than" against a scalar (in-place). */
3041    public DoubleMatrix lti(double value, DoubleMatrix result) {
3042    ensureResultLength(null, result);
3043    for (int i = 0; i < length; i++)
3044    result.put(i, get(i) < value ? 1.0 : 0.0);
3045    return result;
3046    }
3047
3048    /** Test for "less than" against a scalar (in-place). */
3049    public DoubleMatrix lti(double value) {
3050    return lti(value, this);
3051    }
3052
3053    /** test for "less than" against a scalar. */
3054    public DoubleMatrix lt(double value) {
3055    return lti(value, new DoubleMatrix(rows, columns));
3056    }
3057
3058    /** Test for "greater than" (in-place). */
3059    public DoubleMatrix gti(DoubleMatrix other, DoubleMatrix result) {
3060    if (other.isScalar())
3061    return gti(other.scalar(), result);
3062
3063    assertSameLength(other);
3064    ensureResultLength(other, result);
3065
3066    for (int i = 0; i < length; i++)
3067    result.put(i, get(i) > other.get(i) ? 1.0 : 0.0);
3068    return result;
3069    }
3070
3071    /** Test for "greater than" (in-place). */
3072    public DoubleMatrix gti(DoubleMatrix other) {
3073    return gti(other, this);
3074    }
3075
3076    /** Test for "greater than". */
3077    public DoubleMatrix gt(DoubleMatrix other) {
3078    return gti(other, new DoubleMatrix(rows, columns));
3079    }
3080
3081    /** Test for "greater than" against a scalar (in-place). */
3082    public DoubleMatrix gti(double value, DoubleMatrix result) {
3083    ensureResultLength(null, result);
3084    for (int i = 0; i < length; i++)
3085    result.put(i, get(i) > value ? 1.0 : 0.0);
3086    return result;
3087    }
3088
3089    /** Test for "greater than" against a scalar (in-place). */
3090    public DoubleMatrix gti(double value) {
3091    return gti(value, this);
3092    }
3093
3094    /** test for "greater than" against a scalar. */
3095    public DoubleMatrix gt(double value) {
3096    return gti(value, new DoubleMatrix(rows, columns));
3097    }
3098
3099    /** Test for "less than or equal" (in-place). */
3100    public DoubleMatrix lei(DoubleMatrix other, DoubleMatrix result) {
3101    if (other.isScalar())
3102    return lei(other.scalar(), result);
3103
3104    assertSameLength(other);
3105    ensureResultLength(other, result);
3106
3107    for (int i = 0; i < length; i++)
3108    result.put(i, get(i) <= other.get(i) ? 1.0 : 0.0);
3109    return result;
3110    }
3111
3112    /** Test for "less than or equal" (in-place). */
3113    public DoubleMatrix lei(DoubleMatrix other) {
3114    return lei(other, this);
3115    }
3116
3117    /** Test for "less than or equal". */
3118    public DoubleMatrix le(DoubleMatrix other) {
3119    return lei(other, new DoubleMatrix(rows, columns));
3120    }
3121
3122    /** Test for "less than or equal" against a scalar (in-place). */
3123    public DoubleMatrix lei(double value, DoubleMatrix result) {
3124    ensureResultLength(null, result);
3125    for (int i = 0; i < length; i++)
3126    result.put(i, get(i) <= value ? 1.0 : 0.0);
3127    return result;
3128    }
3129
3130    /** Test for "less than or equal" against a scalar (in-place). */
3131    public DoubleMatrix lei(double value) {
3132    return lei(value, this);
3133    }
3134
3135    /** test for "less than or equal" against a scalar. */
3136    public DoubleMatrix le(double value) {
3137    return lei(value, new DoubleMatrix(rows, columns));
3138    }
3139
3140    /** Test for "greater than or equal" (in-place). */
3141    public DoubleMatrix gei(DoubleMatrix other, DoubleMatrix result) {
3142    if (other.isScalar())
3143    return gei(other.scalar(), result);
3144
3145    assertSameLength(other);
3146    ensureResultLength(other, result);
3147
3148    for (int i = 0; i < length; i++)
3149    result.put(i, get(i) >= other.get(i) ? 1.0 : 0.0);
3150    return result;
3151    }
3152
3153    /** Test for "greater than or equal" (in-place). */
3154    public DoubleMatrix gei(DoubleMatrix other) {
3155    return gei(other, this);
3156    }
3157
3158    /** Test for "greater than or equal". */
3159    public DoubleMatrix ge(DoubleMatrix other) {
3160    return gei(other, new DoubleMatrix(rows, columns));
3161    }
3162
3163    /** Test for "greater than or equal" against a scalar (in-place). */
3164    public DoubleMatrix gei(double value, DoubleMatrix result) {
3165    ensureResultLength(null, result);
3166    for (int i = 0; i < length; i++)
3167    result.put(i, get(i) >= value ? 1.0 : 0.0);
3168    return result;
3169    }
3170
3171    /** Test for "greater than or equal" against a scalar (in-place). */
3172    public DoubleMatrix gei(double value) {
3173    return gei(value, this);
3174    }
3175
3176    /** test for "greater than or equal" against a scalar. */
3177    public DoubleMatrix ge(double value) {
3178    return gei(value, new DoubleMatrix(rows, columns));
3179    }
3180
3181    /** Test for equality (in-place). */
3182    public DoubleMatrix eqi(DoubleMatrix other, DoubleMatrix result) {
3183    if (other.isScalar())
3184    return eqi(other.scalar(), result);
3185
3186    assertSameLength(other);
3187    ensureResultLength(other, result);
3188
3189    for (int i = 0; i < length; i++)
3190    result.put(i, get(i) == other.get(i) ? 1.0 : 0.0);
3191    return result;
3192    }
3193
3194    /** Test for equality (in-place). */
3195    public DoubleMatrix eqi(DoubleMatrix other) {
3196    return eqi(other, this);
3197    }
3198
3199    /** Test for equality. */
3200    public DoubleMatrix eq(DoubleMatrix other) {
3201    return eqi(other, new DoubleMatrix(rows, columns));
3202    }
3203
3204    /** Test for equality against a scalar (in-place). */
3205    public DoubleMatrix eqi(double value, DoubleMatrix result) {
3206    ensureResultLength(null, result);
3207    for (int i = 0; i < length; i++)
3208    result.put(i, get(i) == value ? 1.0 : 0.0);
3209    return result;
3210    }
3211
3212    /** Test for equality against a scalar (in-place). */
3213    public DoubleMatrix eqi(double value) {
3214    return eqi(value, this);
3215    }
3216
3217    /** test for equality against a scalar. */
3218    public DoubleMatrix eq(double value) {
3219    return eqi(value, new DoubleMatrix(rows, columns));
3220    }
3221
3222    /** Test for inequality (in-place). */
3223    public DoubleMatrix nei(DoubleMatrix other, DoubleMatrix result) {
3224    if (other.isScalar())
3225    return nei(other.scalar(), result);
3226
3227    assertSameLength(other);
3228    ensureResultLength(other, result);
3229
3230    for (int i = 0; i < length; i++)
3231    result.put(i, get(i) != other.get(i) ? 1.0 : 0.0);
3232    return result;
3233    }
3234
3235    /** Test for inequality (in-place). */
3236    public DoubleMatrix nei(DoubleMatrix other) {
3237    return nei(other, this);
3238    }
3239
3240    /** Test for inequality. */
3241    public DoubleMatrix ne(DoubleMatrix other) {
3242    return nei(other, new DoubleMatrix(rows, columns));
3243    }
3244
3245    /** Test for inequality against a scalar (in-place). */
3246    public DoubleMatrix nei(double value, DoubleMatrix result) {
3247    ensureResultLength(null, result);
3248    for (int i = 0; i < length; i++)
3249    result.put(i, get(i) != value ? 1.0 : 0.0);
3250    return result;
3251    }
3252
3253    /** Test for inequality against a scalar (in-place). */
3254    public DoubleMatrix nei(double value) {
3255    return nei(value, this);
3256    }
3257
3258    /** test for inequality against a scalar. */
3259    public DoubleMatrix ne(double value) {
3260    return nei(value, new DoubleMatrix(rows, columns));
3261    }
3262
3263    /** Compute elementwise logical and (in-place). */
3264    public DoubleMatrix andi(DoubleMatrix other, DoubleMatrix result) {
3265    assertSameLength(other);
3266    ensureResultLength(other, result);
3267
3268    for (int i = 0; i < length; i++)
3269    result.put(i, (get(i) != 0.0) & (other.get(i) != 0.0) ? 1.0 : 0.0);
3270    return result;
3271    }
3272
3273    /** Compute elementwise logical and (in-place). */
3274    public DoubleMatrix andi(DoubleMatrix other) {
3275    return andi(other, this);
3276    }
3277
3278    /** Compute elementwise logical and. */
3279    public DoubleMatrix and(DoubleMatrix other) {
3280    return andi(other, new DoubleMatrix(rows, columns));
3281    }
3282
3283    /** Compute elementwise logical and against a scalar (in-place). */
3284    public DoubleMatrix andi(double value, DoubleMatrix result) {
3285    ensureResultLength(null, result);
3286    boolean val = (value != 0.0);
3287    for (int i = 0; i < length; i++)
3288    result.put(i, (get(i) != 0.0) & val ? 1.0 : 0.0);
3289    return result;
3290    }
3291
3292    /** Compute elementwise logical and against a scalar (in-place). */
3293    public DoubleMatrix andi(double value) {
3294    return andi(value, this);
3295    }
3296
3297    /** Compute elementwise logical and against a scalar. */
3298    public DoubleMatrix and(double value) {
3299    return andi(value, new DoubleMatrix(rows, columns));
3300    }
3301
3302    /** Compute elementwise logical or (in-place). */
3303    public DoubleMatrix ori(DoubleMatrix other, DoubleMatrix result) {
3304    assertSameLength(other);
3305    ensureResultLength(other, result);
3306
3307    for (int i = 0; i < length; i++)
3308    result.put(i, (get(i) != 0.0) | (other.get(i) != 0.0) ? 1.0 : 0.0);
3309    return result;
3310    }
3311
3312    /** Compute elementwise logical or (in-place). */
3313    public DoubleMatrix ori(DoubleMatrix other) {
3314    return ori(other, this);
3315    }
3316
3317    /** Compute elementwise logical or. */
3318    public DoubleMatrix or(DoubleMatrix other) {
3319    return ori(other, new DoubleMatrix(rows, columns));
3320    }
3321
3322    /** Compute elementwise logical or against a scalar (in-place). */
3323    public DoubleMatrix ori(double value, DoubleMatrix result) {
3324    ensureResultLength(null, result);
3325    boolean val = (value != 0.0);
3326    for (int i = 0; i < length; i++)
3327    result.put(i, (get(i) != 0.0) | val ? 1.0 : 0.0);
3328    return result;
3329    }
3330
3331    /** Compute elementwise logical or against a scalar (in-place). */
3332    public DoubleMatrix ori(double value) {
3333    return ori(value, this);
3334    }
3335
3336    /** Compute elementwise logical or against a scalar. */
3337    public DoubleMatrix or(double value) {
3338    return ori(value, new DoubleMatrix(rows, columns));
3339    }
3340
3341    /** Compute elementwise logical xor (in-place). */
3342    public DoubleMatrix xori(DoubleMatrix other, DoubleMatrix result) {
3343    assertSameLength(other);
3344    ensureResultLength(other, result);
3345
3346    for (int i = 0; i < length; i++)
3347    result.put(i, (get(i) != 0.0) ^ (other.get(i) != 0.0) ? 1.0 : 0.0);
3348    return result;
3349    }
3350
3351    /** Compute elementwise logical xor (in-place). */
3352    public DoubleMatrix xori(DoubleMatrix other) {
3353    return xori(other, this);
3354    }
3355
3356    /** Compute elementwise logical xor. */
3357    public DoubleMatrix xor(DoubleMatrix other) {
3358    return xori(other, new DoubleMatrix(rows, columns));
3359    }
3360
3361    /** Compute elementwise logical xor against a scalar (in-place). */
3362    public DoubleMatrix xori(double value, DoubleMatrix result) {
3363    ensureResultLength(null, result);
3364    boolean val = (value != 0.0);
3365    for (int i = 0; i < length; i++)
3366    result.put(i, (get(i) != 0.0) ^ val ? 1.0 : 0.0);
3367    return result;
3368    }
3369
3370    /** Compute elementwise logical xor against a scalar (in-place). */
3371    public DoubleMatrix xori(double value) {
3372    return xori(value, this);
3373    }
3374
3375    /** Compute elementwise logical xor against a scalar. */
3376    public DoubleMatrix xor(double value) {
3377    return xori(value, new DoubleMatrix(rows, columns));
3378    }
3379//RJPP-END--------------------------------------------------------------
3380}