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.jama;
025    
026    import org.ujmp.core.Matrix;
027    import org.ujmp.core.doublematrix.stub.AbstractDenseDoubleMatrix2D;
028    import org.ujmp.core.exceptions.MatrixException;
029    import org.ujmp.core.interfaces.Wrapper;
030    
031    import Jama.CholeskyDecomposition;
032    import Jama.EigenvalueDecomposition;
033    import Jama.LUDecomposition;
034    import Jama.QRDecomposition;
035    import Jama.SingularValueDecomposition;
036    
037    public class JamaDenseDoubleMatrix2D extends AbstractDenseDoubleMatrix2D
038                    implements Wrapper<Jama.Matrix> {
039            private static final long serialVersionUID = -6065454603299978242L;
040    
041            private Jama.Matrix matrix = null;
042    
043            public JamaDenseDoubleMatrix2D(long... size) {
044                    this.matrix = new Jama.Matrix((int) size[ROW], (int) size[COLUMN]);
045            }
046    
047            public JamaDenseDoubleMatrix2D(Jama.Matrix matrix) {
048                    this.matrix = matrix;
049            }
050    
051            public JamaDenseDoubleMatrix2D(Matrix source) throws MatrixException {
052                    this(source.getSize());
053                    for (long[] c : source.availableCoordinates()) {
054                            setDouble(source.getAsDouble(c), c);
055                    }
056            }
057    
058            public static Jama.Matrix identity(int m, int n) {
059                    Jama.Matrix A = new Jama.Matrix(m, n);
060                    double[][] X = A.getArray();
061                    for (int i = 0; i < m; i++) {
062                            for (int j = 0; j < n; j++) {
063                                    X[i][j] = (i == j ? 1.0 : 0.0);
064                            }
065                    }
066                    return A;
067            }
068    
069            public Matrix inv() throws MatrixException {
070                    return new JamaDenseDoubleMatrix2D(matrix.inverse());
071            }
072    
073            public Matrix invSPD() throws MatrixException {
074                    CholeskyDecomposition chol = new CholeskyDecomposition(matrix);
075                    return new JamaDenseDoubleMatrix2D(chol.solve(Jama.Matrix.identity(
076                                    matrix.getRowDimension(), matrix.getRowDimension())));
077            }
078    
079            public Matrix[] svd() throws MatrixException {
080                    if (getColumnCount() > getRowCount()) {
081                            SingularValueDecomposition svd = new SingularValueDecomposition(
082                                            matrix.transpose());
083                            Matrix u = new JamaDenseDoubleMatrix2D(svd.getV());
084                            Matrix s = new JamaDenseDoubleMatrix2D(svd.getS().transpose());
085                            Matrix v = new JamaDenseDoubleMatrix2D(svd.getU());
086                            return new Matrix[] { u, s, v };
087                    } else {
088                            SingularValueDecomposition svd = new SingularValueDecomposition(
089                                            matrix);
090                            Matrix u = new JamaDenseDoubleMatrix2D(svd.getU());
091                            Matrix s = new JamaDenseDoubleMatrix2D(svd.getS());
092                            Matrix v = new JamaDenseDoubleMatrix2D(svd.getV());
093                            return new Matrix[] { u, s, v };
094                    }
095            }
096    
097            public double getDouble(long row, long column) {
098                    return matrix.get((int) row, (int) column);
099            }
100    
101            public double getDouble(int row, int column) {
102                    return matrix.get(row, column);
103            }
104    
105            public long[] getSize() {
106                    return new long[] { matrix.getRowDimension(),
107                                    matrix.getColumnDimension() };
108            }
109    
110            public void setDouble(double value, long row, long column) {
111                    matrix.set((int) row, (int) column, value);
112            }
113    
114            public void setDouble(double value, int row, int column) {
115                    matrix.set(row, column, value);
116            }
117    
118            public Jama.Matrix getWrappedObject() {
119                    return matrix;
120            }
121    
122            public void setWrappedObject(Jama.Matrix object) {
123                    this.matrix = object;
124            }
125    
126            public final Matrix copy() throws MatrixException {
127                    Matrix m = new JamaDenseDoubleMatrix2D(matrix.copy());
128                    if (getAnnotation() != null) {
129                            m.setAnnotation(getAnnotation().clone());
130                    }
131                    return m;
132            }
133    
134            public Matrix transpose() {
135                    return new JamaDenseDoubleMatrix2D(matrix.transpose());
136            }
137    
138            public Matrix[] qr() {
139                    if (getRowCount() >= getColumnCount()) {
140                            QRDecomposition qr = new QRDecomposition(matrix);
141                            Matrix q = new JamaDenseDoubleMatrix2D(qr.getQ());
142                            Matrix r = new JamaDenseDoubleMatrix2D(qr.getR());
143                            return new Matrix[] { q, r };
144                    } else {
145                            throw new MatrixException(
146                                            "QR decomposition only works for matrices m>=n");
147                    }
148            }
149    
150            public Matrix[] lu() {
151                    LUDecomposition lu = new LUDecomposition(matrix);
152                    Matrix l = new JamaDenseDoubleMatrix2D(lu.getL());
153                    Matrix u = new JamaDenseDoubleMatrix2D(lu.getU());
154                    int m = (int) getRowCount();
155                    int[] piv = lu.getPivot();
156                    Matrix p = new JamaDenseDoubleMatrix2D(m, m);
157                    for (int i = 0; i < m; i++) {
158                            p.setAsDouble(1, i, piv[i]);
159                    }
160                    return new Matrix[] { l, u, p };
161            }
162    
163            public Matrix[] eig() {
164                    EigenvalueDecomposition eig = new EigenvalueDecomposition(matrix);
165                    Matrix v = new JamaDenseDoubleMatrix2D(eig.getV());
166                    Matrix d = new JamaDenseDoubleMatrix2D(eig.getD());
167                    return new Matrix[] { v, d };
168            }
169    
170            public Matrix chol() {
171                    CholeskyDecomposition chol = new CholeskyDecomposition(matrix);
172                    Matrix r = new JamaDenseDoubleMatrix2D(chol.getL());
173                    return r;
174            }
175    
176            public Matrix mtimes(Matrix m) {
177                    if (m instanceof JamaDenseDoubleMatrix2D) {
178                            return new JamaDenseDoubleMatrix2D(matrix
179                                            .times(((JamaDenseDoubleMatrix2D) m).matrix));
180                    } else {
181                            return super.mtimes(m);
182                    }
183            }
184    
185            public Matrix times(double value) {
186                    return new JamaDenseDoubleMatrix2D(matrix.times(value));
187            }
188    
189            public Matrix divide(double value) {
190                    return new JamaDenseDoubleMatrix2D(matrix.times(1.0 / value));
191            }
192    
193            public double det() {
194                    return matrix.det();
195            }
196    
197            public Matrix plus(Matrix m) {
198                    if (m instanceof JamaDenseDoubleMatrix2D) {
199                            return new JamaDenseDoubleMatrix2D(matrix
200                                            .plus(((JamaDenseDoubleMatrix2D) m).matrix));
201                    } else {
202                            return super.plus(m);
203                    }
204            }
205    
206            public Matrix minus(Matrix m) {
207                    if (m instanceof JamaDenseDoubleMatrix2D) {
208                            return new JamaDenseDoubleMatrix2D(matrix
209                                            .minus(((JamaDenseDoubleMatrix2D) m).matrix));
210                    } else {
211                            return super.minus(m);
212                    }
213            }
214    
215            public Matrix solve(Matrix b) {
216                    if (b instanceof JamaDenseDoubleMatrix2D) {
217                            JamaDenseDoubleMatrix2D b2 = (JamaDenseDoubleMatrix2D) b;
218                            Jama.Matrix x = matrix.solve(b2.matrix);
219                            return new JamaDenseDoubleMatrix2D(x);
220                    } else {
221                            return super.solve(b);
222                    }
223            }
224    
225            public Matrix solveSPD(Matrix b) {
226                    if (b instanceof JamaDenseDoubleMatrix2D) {
227                            JamaDenseDoubleMatrix2D b2 = (JamaDenseDoubleMatrix2D) b;
228                            CholeskyDecomposition chol = new CholeskyDecomposition(matrix);
229                            return new JamaDenseDoubleMatrix2D(chol.solve(b2.matrix));
230                    } else {
231                            return super.solve(b);
232                    }
233            }
234    
235    }