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.jmatrices;
025    
026    import org.jmatrices.dbl.MatrixFactory;
027    import org.jmatrices.dbl.decomposition.CholeskyDecomposition;
028    import org.jmatrices.dbl.decomposition.EigenvalueDecomposition;
029    import org.jmatrices.dbl.decomposition.LUDecomposition;
030    import org.jmatrices.dbl.decomposition.QRDecomposition;
031    import org.jmatrices.dbl.decomposition.SingularValueDecomposition;
032    import org.jmatrices.dbl.operator.MatrixOperator;
033    import org.jmatrices.dbl.transformer.MatrixTransformer;
034    import org.ujmp.core.Coordinates;
035    import org.ujmp.core.Matrix;
036    import org.ujmp.core.doublematrix.stub.AbstractDenseDoubleMatrix2D;
037    import org.ujmp.core.exceptions.MatrixException;
038    import org.ujmp.core.interfaces.Wrapper;
039    
040    public class JMatricesDenseDoubleMatrix2D extends AbstractDenseDoubleMatrix2D
041                    implements Wrapper<org.jmatrices.dbl.Matrix> {
042            private static final long serialVersionUID = 513251881654621L;
043    
044            private org.jmatrices.dbl.Matrix matrix = null;
045    
046            public JMatricesDenseDoubleMatrix2D(org.jmatrices.dbl.Matrix matrix) {
047                    this.matrix = matrix;
048            }
049    
050            public JMatricesDenseDoubleMatrix2D(long... size) {
051                    if (size[ROW] > 0 && size[COLUMN] > 0) {
052                            this.matrix = MatrixFactory.getMatrix((int) size[ROW],
053                                            (int) size[COLUMN], null);
054                    }
055            }
056    
057            public JMatricesDenseDoubleMatrix2D(org.ujmp.core.Matrix source)
058                            throws MatrixException {
059                    this(source.getSize());
060                    for (long[] c : source.availableCoordinates()) {
061                            setDouble(source.getAsDouble(c), c);
062                    }
063            }
064    
065            public double getDouble(long row, long column) {
066                    return matrix.get((int) row + 1, (int) column + 1);
067            }
068    
069            public double getDouble(int row, int column) {
070                    return matrix.get(row + 1, column + 1);
071            }
072    
073            public long[] getSize() {
074                    return matrix == null ? Coordinates.ZERO2D : new long[] {
075                                    matrix.rows(), matrix.cols() };
076            }
077    
078            public void setDouble(double value, long row, long column) {
079                    matrix.set((int) row + 1, (int) column + 1, value);
080            }
081    
082            public void setDouble(double value, int row, int column) {
083                    matrix.set(row + 1, column + 1, value);
084            }
085    
086            public org.jmatrices.dbl.Matrix getWrappedObject() {
087                    return matrix;
088            }
089    
090            public void setWrappedObject(org.jmatrices.dbl.Matrix object) {
091                    this.matrix = object;
092            }
093    
094            public Matrix transpose() {
095                    return new JMatricesDenseDoubleMatrix2D(MatrixTransformer
096                                    .transpose(matrix));
097            }
098    
099            public Matrix inv() {
100                    return new JMatricesDenseDoubleMatrix2D(MatrixTransformer
101                                    .inverse(matrix));
102            }
103    
104            public Matrix[] eig() {
105                    EigenvalueDecomposition evd = new EigenvalueDecomposition(matrix);
106                    Matrix v = new JMatricesDenseDoubleMatrix2D(evd.getV());
107                    Matrix d = new JMatricesDenseDoubleMatrix2D(evd.getD());
108                    return new Matrix[] { v, d };
109            }
110    
111            public Matrix[] qr() {
112                    if (getRowCount() >= getColumnCount()) {
113                            QRDecomposition qr = new QRDecomposition(matrix);
114                            Matrix q = new JMatricesDenseDoubleMatrix2D(qr.getQ());
115                            Matrix r = new JMatricesDenseDoubleMatrix2D(qr.getR());
116                            return new Matrix[] { q, r };
117                    } else {
118                            throw new MatrixException("only allowed for matrices m>=n");
119                    }
120            }
121    
122            public Matrix[] svd() {
123                    if (isSquare()) {
124                            SingularValueDecomposition qr = new SingularValueDecomposition(
125                                            matrix);
126                            Matrix u = new JMatricesDenseDoubleMatrix2D(qr.getU());
127                            Matrix s = new JMatricesDenseDoubleMatrix2D(qr.getS());
128                            Matrix v = new JMatricesDenseDoubleMatrix2D(qr.getV());
129                            return new Matrix[] { u, s, v };
130                    } else {
131                            throw new MatrixException("only allowed for square matrices");
132                    }
133            }
134    
135            public Matrix chol() {
136                    CholeskyDecomposition chol = new CholeskyDecomposition(matrix);
137                    Matrix r = new JMatricesDenseDoubleMatrix2D(chol.getL());
138                    return r;
139            }
140    
141            // some error
142            // public Matrix invSPD() throws MatrixException {
143            // CholeskyDecomposition chol = new CholeskyDecomposition(matrix);
144            // org.jmatrices.dbl.Matrix eye = MatrixFactory.getIdentityMatrix(matrix
145            // .rows(), null);
146            // return new JMatricesDenseDoubleMatrix2D(chol.solve(eye));
147            // }
148    
149            public Matrix[] lu() {
150                    if (getRowCount() >= getColumnCount()) {
151                            LUDecomposition lu = new LUDecomposition(matrix);
152                            Matrix l = new JMatricesDenseDoubleMatrix2D(lu.getL());
153                            Matrix u = new JMatricesDenseDoubleMatrix2D(lu.getU().getSubMatrix(
154                                            1, 1, (int) getColumnCount(), (int) getColumnCount()));
155                            int m = (int) getRowCount();
156                            int[] piv = lu.getPivot();
157                            Matrix p = new JMatricesDenseDoubleMatrix2D(m, m);
158                            for (int i = 0; i < m; i++) {
159                                    p.setAsDouble(1, i, piv[i] - 1);
160                            }
161                            return new Matrix[] { l, u, p };
162                    } else {
163                            throw new MatrixException("only allowed for matrices m>=n");
164                    }
165            }
166    
167            public Matrix mtimes(Matrix m2) {
168                    if (m2 instanceof JMatricesDenseDoubleMatrix2D) {
169                            return new JMatricesDenseDoubleMatrix2D(MatrixOperator.multiply(
170                                            matrix, ((JMatricesDenseDoubleMatrix2D) m2).matrix));
171                    } else {
172                            return super.mtimes(m2);
173                    }
174            }
175    
176            public Matrix solve(Matrix b) {
177                    if (b instanceof JMatricesDenseDoubleMatrix2D) {
178                            JMatricesDenseDoubleMatrix2D b2 = (JMatricesDenseDoubleMatrix2D) b;
179                            if (isSquare()) {
180                                    org.jmatrices.dbl.Matrix x = new LUDecomposition(matrix)
181                                                    .solve(b2.matrix);
182                                    return new JMatricesDenseDoubleMatrix2D(x);
183                            } else {
184                                    org.jmatrices.dbl.Matrix x = new QRDecomposition(matrix)
185                                                    .solve(b2.matrix);
186                                    return new JMatricesDenseDoubleMatrix2D(x);
187                            }
188                    } else {
189                            return super.solve(b);
190                    }
191            }
192    }