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.commonsmath;
025    
026    import org.apache.commons.math.linear.CholeskyDecomposition;
027    import org.apache.commons.math.linear.CholeskyDecompositionImpl;
028    import org.apache.commons.math.linear.EigenDecomposition;
029    import org.apache.commons.math.linear.EigenDecompositionImpl;
030    import org.apache.commons.math.linear.LUDecomposition;
031    import org.apache.commons.math.linear.LUDecompositionImpl;
032    import org.apache.commons.math.linear.QRDecomposition;
033    import org.apache.commons.math.linear.QRDecompositionImpl;
034    import org.apache.commons.math.linear.RealMatrix;
035    import org.apache.commons.math.linear.SingularValueDecomposition;
036    import org.apache.commons.math.linear.SingularValueDecompositionImpl;
037    import org.apache.commons.math.util.MathUtils;
038    import org.ujmp.core.Coordinates;
039    import org.ujmp.core.Matrix;
040    import org.ujmp.core.doublematrix.stub.AbstractDenseDoubleMatrix2D;
041    import org.ujmp.core.exceptions.MatrixException;
042    import org.ujmp.core.interfaces.Wrapper;
043    
044    public abstract class AbstractCommonsMathDenseDoubleMatrix2D extends
045                    AbstractDenseDoubleMatrix2D implements Wrapper<RealMatrix> {
046            private static final long serialVersionUID = -1161807620507675926L;
047    
048            private RealMatrix matrix = null;
049    
050            public AbstractCommonsMathDenseDoubleMatrix2D(RealMatrix matrix) {
051                    this.matrix = matrix;
052            }
053    
054            public RealMatrix getWrappedObject() {
055                    return matrix;
056            }
057    
058            public void setWrappedObject(RealMatrix object) {
059                    this.matrix = object;
060            }
061    
062            public double getDouble(long row, long column) throws MatrixException {
063                    return matrix.getEntry((int) row, (int) column);
064            }
065    
066            public double getDouble(int row, int column) throws MatrixException {
067                    return matrix.getEntry(row, column);
068            }
069    
070            public void setDouble(double value, long row, long column)
071                            throws MatrixException {
072                    matrix.setEntry((int) row, (int) column, value);
073            }
074    
075            public void setDouble(double value, int row, int column)
076                            throws MatrixException {
077                    matrix.setEntry((int) row, (int) column, value);
078            }
079    
080            public long[] getSize() {
081                    return matrix == null ? Coordinates.ZERO2D : new long[] {
082                                    matrix.getRowDimension(), matrix.getColumnDimension() };
083            }
084    
085            public Matrix transpose() {
086                    return CommonsMathDenseDoubleMatrix2DFactory.INSTANCE.dense(matrix
087                                    .transpose());
088            }
089    
090            public Matrix inv() {
091                    return CommonsMathDenseDoubleMatrix2DFactory.INSTANCE
092                                    .dense(new LUDecompositionImpl(matrix).getSolver().getInverse());
093            }
094    
095            public Matrix invSPD() {
096                    try {
097                            return CommonsMathDenseDoubleMatrix2DFactory.INSTANCE
098                                            .dense(new CholeskyDecompositionImpl(matrix).getSolver()
099                                                            .getInverse());
100                    } catch (Exception e) {
101                            throw new MatrixException(e);
102                    }
103            }
104    
105            public Matrix[] lu() {
106                    LUDecomposition lu = new LUDecompositionImpl(matrix);
107                    Matrix l = CommonsMathDenseDoubleMatrix2DFactory.INSTANCE.dense(lu
108                                    .getL());
109                    Matrix u = CommonsMathDenseDoubleMatrix2DFactory.INSTANCE.dense(lu
110                                    .getU());
111                    Matrix p = CommonsMathDenseDoubleMatrix2DFactory.INSTANCE.dense(lu
112                                    .getP());
113                    return new Matrix[] { l, u, p };
114            }
115    
116            public Matrix[] qr() {
117                    QRDecomposition qr = new QRDecompositionImpl(matrix);
118                    Matrix q = CommonsMathDenseDoubleMatrix2DFactory.INSTANCE.dense(qr
119                                    .getQ());
120                    Matrix r = CommonsMathDenseDoubleMatrix2DFactory.INSTANCE.dense(qr
121                                    .getR());
122                    return new Matrix[] { q, r };
123            }
124    
125            public Matrix[] svd() {
126                    SingularValueDecomposition svd = new SingularValueDecompositionImpl(
127                                    matrix);
128                    Matrix u = CommonsMathDenseDoubleMatrix2DFactory.INSTANCE.dense(svd
129                                    .getU());
130                    Matrix s = CommonsMathDenseDoubleMatrix2DFactory.INSTANCE.dense(svd
131                                    .getS());
132                    Matrix v = CommonsMathDenseDoubleMatrix2DFactory.INSTANCE.dense(svd
133                                    .getV());
134                    return new Matrix[] { u, s, v };
135            }
136    
137            public Matrix[] eig() {
138                    EigenDecomposition evd = new EigenDecompositionImpl(matrix,
139                                    MathUtils.EPSILON);
140                    Matrix v = CommonsMathDenseDoubleMatrix2DFactory.INSTANCE.dense(evd
141                                    .getV());
142                    Matrix d = CommonsMathDenseDoubleMatrix2DFactory.INSTANCE.dense(evd
143                                    .getD());
144                    return new Matrix[] { v, d };
145            }
146    
147            public Matrix chol() {
148                    try {
149                            CholeskyDecomposition chol = new CholeskyDecompositionImpl(matrix);
150                            Matrix l = CommonsMathDenseDoubleMatrix2DFactory.INSTANCE
151                                            .dense(chol.getL());
152                            return l;
153                    } catch (Exception e) {
154                            throw new MatrixException(e);
155                    }
156            }
157    
158            public Matrix mtimes(Matrix m2) {
159                    if (m2 instanceof AbstractCommonsMathDenseDoubleMatrix2D) {
160                            return CommonsMathDenseDoubleMatrix2DFactory.INSTANCE
161                                            .dense(matrix
162                                                            .multiply(((AbstractCommonsMathDenseDoubleMatrix2D) m2).matrix));
163                    } else {
164                            return super.mtimes(m2);
165                    }
166            }
167    
168            public Matrix plus(Matrix m2) {
169                    if (m2 instanceof AbstractCommonsMathDenseDoubleMatrix2D) {
170                            return CommonsMathDenseDoubleMatrix2DFactory.INSTANCE.dense(matrix
171                                            .add(((AbstractCommonsMathDenseDoubleMatrix2D) m2).matrix));
172                    } else {
173                            return super.plus(m2);
174                    }
175            }
176    
177            public Matrix minus(Matrix m2) {
178                    if (m2 instanceof AbstractCommonsMathDenseDoubleMatrix2D) {
179                            return CommonsMathDenseDoubleMatrix2DFactory.INSTANCE
180                                            .dense(matrix
181                                                            .subtract(((AbstractCommonsMathDenseDoubleMatrix2D) m2).matrix));
182                    } else {
183                            return super.minus(m2);
184                    }
185            }
186    
187            public Matrix times(double value) {
188                    return CommonsMathDenseDoubleMatrix2DFactory.INSTANCE.dense(matrix
189                                    .scalarMultiply(value));
190            }
191    
192            public Matrix divide(double value) {
193                    return CommonsMathDenseDoubleMatrix2DFactory.INSTANCE.dense(matrix
194                                    .scalarMultiply(1.0 / value));
195            }
196    
197            public Matrix plus(double value) {
198                    return CommonsMathDenseDoubleMatrix2DFactory.INSTANCE.dense(matrix
199                                    .scalarAdd(value));
200            }
201    
202            public Matrix minus(double value) {
203                    return CommonsMathDenseDoubleMatrix2DFactory.INSTANCE.dense(matrix
204                                    .scalarAdd(-value));
205            }
206    
207            public Matrix solve(Matrix b) {
208                    if (b instanceof AbstractCommonsMathDenseDoubleMatrix2D) {
209                            AbstractCommonsMathDenseDoubleMatrix2D b2 = (AbstractCommonsMathDenseDoubleMatrix2D) b;
210                            if (isSquare()) {
211                                    RealMatrix ret = new LUDecompositionImpl(matrix).getSolver()
212                                                    .solve(b2.matrix);
213                                    return CommonsMathDenseDoubleMatrix2DFactory.INSTANCE
214                                                    .dense(ret);
215                            } else {
216                                    RealMatrix ret = new QRDecompositionImpl(matrix).getSolver()
217                                                    .solve(b2.matrix);
218                                    return CommonsMathDenseDoubleMatrix2DFactory.INSTANCE
219                                                    .dense(ret);
220                            }
221                    } else {
222                            return super.solve(b);
223                    }
224            }
225    
226            public Matrix solveSPD(Matrix b) {
227                    try {
228                            if (b instanceof AbstractCommonsMathDenseDoubleMatrix2D) {
229                                    AbstractCommonsMathDenseDoubleMatrix2D b2 = (AbstractCommonsMathDenseDoubleMatrix2D) b;
230                                    RealMatrix ret = new CholeskyDecompositionImpl(matrix)
231                                                    .getSolver().solve(b2.matrix);
232                                    return CommonsMathDenseDoubleMatrix2DFactory.INSTANCE
233                                                    .dense(ret);
234                            } else {
235                                    return super.solve(b);
236                            }
237                    } catch (Exception e) {
238                            throw new MatrixException(e);
239                    }
240            }
241    }