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