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