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