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