001    /*
002     * Copyright (C) 2008-2010 by Holger Arndt
003     *
004     * This file is part of the Universal Java Matrix Package (UJMP).
005     * See the NOTICE file distributed with this work for additional
006     * information regarding copyright ownership and licensing.
007     *
008     * UJMP is free software; you can redistribute it and/or modify
009     * it under the terms of the GNU Lesser General Public License as
010     * published by the Free Software Foundation; either version 2
011     * of the License, or (at your option) any later version.
012     *
013     * UJMP is distributed in the hope that it will be useful,
014     * but WITHOUT ANY WARRANTY; without even the implied warranty of
015     * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
016     * GNU Lesser General Public License for more details.
017     *
018     * You should have received a copy of the GNU Lesser General Public
019     * License along with UJMP; if not, write to the
020     * Free Software Foundation, Inc., 51 Franklin St, Fifth Floor,
021     * Boston, MA  02110-1301  USA
022     */
023    
024    package org.ujmp.jmatharray;
025    
026    import org.math.array.DoubleArray;
027    import org.math.array.LinearAlgebra;
028    import org.ujmp.core.Matrix;
029    import org.ujmp.core.doublematrix.stub.AbstractDenseDoubleMatrix2D;
030    import org.ujmp.core.exceptions.MatrixException;
031    import org.ujmp.core.interfaces.Wrapper;
032    
033    import Jama.CholeskyDecomposition;
034    import Jama.EigenvalueDecomposition;
035    import Jama.LUDecomposition;
036    import Jama.QRDecomposition;
037    import Jama.SingularValueDecomposition;
038    
039    public class JMathArrayDenseDoubleMatrix2D extends AbstractDenseDoubleMatrix2D
040                    implements Wrapper<double[][]> {
041            private static final long serialVersionUID = -3223474248020842822L;
042    
043            private double[][] matrix = null;
044    
045            public JMathArrayDenseDoubleMatrix2D(long... size) {
046                    this.matrix = new double[(int) size[ROW]][(int) size[COLUMN]];
047            }
048    
049            public JMathArrayDenseDoubleMatrix2D(double[][] m) {
050                    this.matrix = m;
051            }
052    
053            public JMathArrayDenseDoubleMatrix2D(Matrix source) throws MatrixException {
054                    this(source.getSize());
055                    for (long[] c : source.availableCoordinates()) {
056                            setDouble(source.getAsDouble(c), c);
057                    }
058            }
059    
060            public double getDouble(long row, long column) {
061                    return matrix[(int) row][(int) column];
062            }
063    
064            public double getDouble(int row, int column) {
065                    return matrix[row][column];
066            }
067    
068            public long[] getSize() {
069                    return new long[] { matrix.length,
070                                    matrix.length == 0 ? 0 : matrix[0].length };
071            }
072    
073            public void setDouble(double value, long row, long column) {
074                    matrix[(int) row][(int) column] = value;
075            }
076    
077            public void setDouble(double value, int row, int column) {
078                    matrix[row][column] = value;
079            }
080    
081            public Matrix[] svd() throws MatrixException {
082                    if (getColumnCount() > getRowCount()) {
083                            SingularValueDecomposition svd = new SingularValueDecomposition(
084                                            new Jama.Matrix(DoubleArray.transpose(matrix)));
085                            Matrix u = new JMathArrayDenseDoubleMatrix2D(svd.getV().getArray());
086                            Matrix s = new JMathArrayDenseDoubleMatrix2D(svd.getS().transpose()
087                                            .getArray());
088                            Matrix v = new JMathArrayDenseDoubleMatrix2D(svd.getU().getArray());
089                            return new Matrix[] { u, s, v };
090                    } else {
091                            SingularValueDecomposition svd = new SingularValueDecomposition(
092                                            new Jama.Matrix(matrix));
093                            Matrix u = new JMathArrayDenseDoubleMatrix2D(svd.getU().getArray());
094                            Matrix s = new JMathArrayDenseDoubleMatrix2D(svd.getS().getArray());
095                            Matrix v = new JMathArrayDenseDoubleMatrix2D(svd.getV().getArray());
096                            return new Matrix[] { u, s, v };
097                    }
098            }
099    
100            public Matrix[] qr() {
101                    if (getRowCount() >= getColumnCount()) {
102                            QRDecomposition qr = new QRDecomposition(new Jama.Matrix(matrix));
103                            Matrix q = new JMathArrayDenseDoubleMatrix2D(qr.getQ().getArray());
104                            Matrix r = new JMathArrayDenseDoubleMatrix2D(qr.getR().getArray());
105                            return new Matrix[] { q, r };
106                    } else {
107                            throw new MatrixException(
108                                            "QR decomposition only works for matrices m>=n");
109                    }
110            }
111    
112            public Matrix[] eig() {
113                    EigenvalueDecomposition eig = new EigenvalueDecomposition(
114                                    new Jama.Matrix(matrix));
115                    Matrix v = new JMathArrayDenseDoubleMatrix2D(eig.getV().getArray());
116                    Matrix d = new JMathArrayDenseDoubleMatrix2D(eig.getD().getArray());
117                    return new Matrix[] { v, d };
118            }
119    
120            public Matrix chol() {
121                    CholeskyDecomposition chol = new CholeskyDecomposition(new Jama.Matrix(
122                                    matrix));
123                    Matrix r = new JMathArrayDenseDoubleMatrix2D(chol.getL().getArray());
124                    return r;
125            }
126    
127            public Matrix[] lu() {
128                    if (getRowCount() >= getColumnCount()) {
129                            LUDecomposition lu = new LUDecomposition(new Jama.Matrix(matrix));
130                            Matrix l = new JMathArrayDenseDoubleMatrix2D(lu.getL().getArray());
131                            Matrix u = new JMathArrayDenseDoubleMatrix2D(lu.getU().getArray());
132                            int m = (int) getRowCount();
133                            int[] piv = lu.getPivot();
134                            Matrix p = new JMathArrayDenseDoubleMatrix2D(m, m);
135                            for (int i = 0; i < m; i++) {
136                                    p.setAsDouble(1, i, piv[i]);
137                            }
138                            return new Matrix[] { l, u, p };
139                    } else {
140                            throw new MatrixException(
141                                            "LU decomposition only works for matrices m>=n");
142                    }
143            }
144    
145            public double[][] getWrappedObject() {
146                    return matrix;
147            }
148    
149            public void setWrappedObject(double[][] object) {
150                    this.matrix = object;
151            }
152    
153            public Matrix transpose() {
154                    return new JMathArrayDenseDoubleMatrix2D(DoubleArray.transpose(matrix));
155            }
156    
157            public Matrix inv() {
158                    return new JMathArrayDenseDoubleMatrix2D(LinearAlgebra.inverse(matrix));
159            }
160    
161            public Matrix invSPD() {
162                    CholeskyDecomposition chol = new CholeskyDecomposition(new Jama.Matrix(
163                                    matrix));
164                    return new JMathArrayDenseDoubleMatrix2D(chol.solve(
165                                    Jama.Matrix.identity(matrix.length, matrix.length)).getArray());
166            }
167    
168            public Matrix plus(double value) {
169                    return new JMathArrayDenseDoubleMatrix2D(DoubleArray.add(DoubleArray
170                                    .copy(matrix), value));
171            }
172    
173            public Matrix minus(double value) {
174                    return new JMathArrayDenseDoubleMatrix2D(DoubleArray.add(DoubleArray
175                                    .copy(matrix), -value));
176            }
177    
178            public Matrix times(double value) {
179                    return new JMathArrayDenseDoubleMatrix2D(LinearAlgebra.times(matrix,
180                                    value));
181            }
182    
183            public Matrix divide(double value) {
184                    return new JMathArrayDenseDoubleMatrix2D(LinearAlgebra.divide(matrix,
185                                    value));
186            }
187    
188            public Matrix mtimes(Matrix m) {
189                    if (m instanceof JMathArrayDenseDoubleMatrix2D) {
190                            return new JMathArrayDenseDoubleMatrix2D(LinearAlgebra.times(
191                                            matrix, ((JMathArrayDenseDoubleMatrix2D) m).matrix));
192                    } else {
193                            return super.mtimes(m);
194                    }
195            }
196    
197            public Matrix copy() {
198                    Matrix m = new JMathArrayDenseDoubleMatrix2D(DoubleArray.copy(matrix));
199                    if (getAnnotation() != null) {
200                            m.setAnnotation(getAnnotation().clone());
201                    }
202                    return m;
203            }
204    
205            public Matrix solve(Matrix b) {
206                    if (b instanceof JMathArrayDenseDoubleMatrix2D) {
207                            JMathArrayDenseDoubleMatrix2D b2 = (JMathArrayDenseDoubleMatrix2D) b;
208                            double[][] x = LinearAlgebra.solve(matrix, b2.matrix);
209                            return new JMathArrayDenseDoubleMatrix2D(x);
210                    } else {
211                            return super.solve(b);
212                    }
213            }
214    
215    }