001: /*
002: * JScience - Java(TM) Tools and Libraries for the Advancement of Sciences.
003: * Copyright (C) 2006 - JScience (http://jscience.org/)
004: * All rights reserved.
005: *
006: * Permission to use, copy, modify, and distribute this software is
007: * freely granted, provided that this notice is preserved.
008: */
009: package org.jscience.mathematics.vector;
010:
011: import java.util.Iterator;
012: import java.util.List;
013:
014: import javolution.context.ObjectFactory;
015: import javolution.lang.MathLib;
016: import javolution.util.FastComparator;
017: import javolution.util.FastMap;
018: import javolution.util.FastTable;
019: import javolution.util.Index;
020:
021: import org.jscience.mathematics.structure.Field;
022:
023: /**
024: * <p> This class represents a matrix made of {@link SparseVector sparse
025: * vectors} (as rows). To create a sparse matrix made of column vectors the
026: * {@link #transpose} method can be used.
027: * For example:[code]
028: * SparseVector<Rational> column0 = SparseVector.valueOf(...);
029: * SparseVector<Rational> column1 = SparseVector.valueOf(...);
030: * SparseMatrix<Rational> M = SparseMatrix.valueOf(column0, column1).transpose();
031: * [/code]</p>
032: * <p> As for any concrete {@link org.jscience.mathematics.structure.Structure
033: * structure}, this class is declared <code>final</code> (otherwise most
034: * operations would have to be overridden to return the appropriate type).
035: * Specialized dense matrix should sub-class {@link Matrix} directly.
036: * For example:[code]
037: * // Extension through composition.
038: * final class BandMatrix <F extends Field<F>> extends Matrix<F> {
039: * private SparseMatrix<F> _value;
040: * ...
041: * public BandMatrix opposite() { // Returns the right type.
042: * return BandMatrix.valueOf(_value.opposite());
043: * }
044: * ...
045: * }[/code]
046: * </p>
047: * @author <a href="mailto:jean-marie@dautelle.com">Jean-Marie Dautelle</a>
048: * @version 3.3, January 2, 2007
049: */
050: public final class SparseMatrix<F extends Field<F>> extends Matrix<F> {
051:
052: /**
053: * Holds the number of columns n or the number of rows m if transposed.
054: */
055: int _n;
056:
057: /**
058: * Holds the zero.
059: */
060: F _zero;;
061:
062: /**
063: * Indicates if this matrix is transposed (the rows are then the columns).
064: */
065: boolean _transposed;
066:
067: /**
068: * Holds this matrix rows (or columns when transposed).
069: */
070: final FastTable<SparseVector<F>> _rows = new FastTable<SparseVector<F>>();
071:
072: /**
073: * Returns the sparse square matrix having the specified diagonal
074: * vector. This method is typically used to create an identity matrix.
075: * For example:[code]
076: * SparseMatrix<Real> IDENTITY = Matrix.valueOf(
077: * DenseVector.valueOf({Real.ONE, Real.ONE, Real.ONE}), Real.ZERO);
078: * [/code]
079: *
080: * @param diagonal the diagonal vector.
081: * @param zero value of non-diagonal elements.
082: * @return a square matrix with <code>diagonal</code> on the diagonal and
083: * <code>zero</code> elsewhere.
084: */
085: public static <F extends Field<F>> SparseMatrix<F> valueOf(
086: Vector<F> diagonal, F zero) {
087: int n = diagonal.getDimension();
088: SparseMatrix<F> M = SparseMatrix.newInstance(n, zero, false);
089: for (int i = 0; i < n; i++) {
090: SparseVector<F> row = SparseVector.valueOf(n, zero, i,
091: diagonal.get(i));
092: M._rows.add(row);
093: }
094: return M;
095: }
096:
097: /**
098: * Returns a sparse matrix holding the specified row vectors
099: * (column vectors if {@link #transpose transposed}).
100: *
101: * @param rows the row vectors.
102: * @return the matrix having the specified rows.
103: * @throws DimensionException if the rows do not have the same dimension.
104: */
105: public static <F extends Field<F>> SparseMatrix<F> valueOf(
106: SparseVector<F>... rows) {
107: final int n = rows[0]._dimension;
108: final F zero = rows[0]._zero;
109: SparseMatrix<F> M = SparseMatrix.newInstance(n, zero, false);
110: for (int i = 0, m = rows.length; i < m; i++) {
111: SparseVector<F> rowi = rows[i];
112: if (rowi._dimension != n)
113: throw new DimensionException(
114: "All vectors must have the same dimension.");
115: if (!zero.equals(rowi._zero))
116: throw new DimensionException(
117: "All vectors must have the same zero element.");
118: M._rows.add(rowi);
119: }
120: return M;
121: }
122:
123: /**
124: * Returns a sparse matrix holding the row vectors from the specified
125: * collection (column vectors if {@link #transpose transposed}).
126: *
127: * @param rows the list of row vectors.
128: * @return the matrix having the specified rows.
129: * @throws DimensionException if the rows do not have the same dimension.
130: */
131: public static <F extends Field<F>> SparseMatrix<F> valueOf(
132: List<SparseVector<F>> rows) {
133: final int n = rows.get(0)._dimension;
134: final F zero = rows.get(0)._zero;
135: SparseMatrix<F> M = SparseMatrix.newInstance(n, zero, false);
136: Iterator<SparseVector<F>> iterator = rows.iterator();
137: for (int i = 0, m = rows.size(); i < m; i++) {
138: SparseVector<F> rowi = iterator.next();
139: if (rowi.getDimension() != n)
140: throw new DimensionException(
141: "All vectors must have the same dimension.");
142: if (!zero.equals(rowi._zero))
143: throw new DimensionException(
144: "All vectors must have the same zero element.");
145: M._rows.add(rowi);
146: }
147: return M;
148: }
149:
150: /**
151: * Returns a sparse matrix equivalent to the specified matrix but with
152: * the zero elements removed using the default object equality comparator.
153: *
154: * @param that the matrix to convert.
155: * @param zero the zero element for the sparse vector to return.
156: * @return <code>SparseMatrix.valueOf(that, zero, FastComparator.DEFAULT)</code> or a dense matrix holding the same elements
157: */
158: public static <F extends Field<F>> SparseMatrix<F> valueOf(
159: Matrix<F> that, F zero) {
160: return SparseMatrix.valueOf(that, zero, FastComparator.DEFAULT);
161: }
162:
163: /**
164: * Returns a sparse matrix equivalent to the specified matrix but with
165: * the zero elements removed using the specified object equality comparator.
166: *
167: * @param that the matrix to convert.
168: * @param zero the zero element for the sparse vector to return.
169: * @param comparator the comparator used to determinate zero equality.
170: * @return <code>that</code> or a dense matrix holding the same elements
171: * as the specified matrix.
172: */
173: public static <F extends Field<F>> SparseMatrix<F> valueOf(
174: Matrix<F> that, F zero, FastComparator<? super F> comparator) {
175: if (that instanceof SparseMatrix)
176: return (SparseMatrix<F>) that;
177: int n = that.getNumberOfColumns();
178: int m = that.getNumberOfRows();
179: SparseMatrix<F> M = SparseMatrix.newInstance(n, zero, false);
180: for (int i = 0; i < m; i++) {
181: SparseVector<F> rowi = SparseVector.valueOf(that.getRow(i),
182: zero, comparator);
183: M._rows.add(rowi);
184: }
185: return M;
186: }
187:
188: /**
189: * Returns the value of the non-set elements for this sparse matrix.
190: *
191: * @return the element corresponding to zero.
192: */
193: public F getZero() {
194: return _zero;
195: }
196:
197: @Override
198: public int getNumberOfRows() {
199: return _transposed ? _n : _rows.size();
200: }
201:
202: @Override
203: public int getNumberOfColumns() {
204: return _transposed ? _rows.size() : _n;
205: }
206:
207: @Override
208: public F get(int i, int j) {
209: return _transposed ? _rows.get(j).get(i) : _rows.get(i).get(j);
210: }
211:
212: @Override
213: public SparseVector<F> getRow(int i) {
214: if (!_transposed)
215: return _rows.get(i);
216: // Else transposed.
217: int n = _rows.size();
218: int m = _n;
219: if ((i < 0) || (i >= m))
220: throw new DimensionException();
221: SparseVector<F> V = SparseVector.newInstance(n, _zero);
222: for (int j = 0; j < n; j++) {
223: SparseVector<F> row = _rows.get(j);
224: F e = row._elements.get(Index.valueOf(i));
225: if (e != null) {
226: V._elements.put(Index.valueOf(j), e);
227: }
228: }
229: return V;
230: }
231:
232: @Override
233: public SparseVector<F> getColumn(int j) {
234: if (_transposed)
235: return _rows.get(j);
236: int m = _rows.size();
237: if ((j < 0) || (j >= _n))
238: throw new DimensionException();
239: SparseVector<F> V = SparseVector.newInstance(_n, _zero);
240: for (int i = 0; i < m; i++) {
241: SparseVector<F> row = _rows.get(i);
242: F e = row._elements.get(Index.valueOf(j));
243: if (e != null) {
244: V._elements.put(Index.valueOf(i), e);
245: }
246: }
247: return V;
248: }
249:
250: @Override
251: public SparseVector<F> getDiagonal() {
252: int m = this .getNumberOfRows();
253: int n = this .getNumberOfColumns();
254: int dimension = MathLib.min(m, n);
255: SparseVector<F> V = SparseVector.newInstance(_n, _zero);
256: for (int i = 0; i < dimension; i++) {
257: SparseVector<F> row = _rows.get(i);
258: F e = row._elements.get(Index.valueOf(i));
259: if (e != null) {
260: V._elements.put(Index.valueOf(i), e);
261: }
262: }
263: return V;
264: }
265:
266: @Override
267: public SparseMatrix<F> opposite() {
268: SparseMatrix<F> M = SparseMatrix.newInstance(_n, _zero,
269: _transposed);
270: for (int i = 0, p = _rows.size(); i < p; i++) {
271: M._rows.add(_rows.get(i).opposite());
272: }
273: return M;
274: }
275:
276: @Override
277: public SparseMatrix<F> plus(Matrix<F> that) {
278: if (this .getNumberOfRows() != that.getNumberOfRows())
279: throw new DimensionException();
280: SparseMatrix<F> M = SparseMatrix.newInstance(_n, _zero,
281: _transposed);
282: for (int i = 0, p = _rows.size(); i < p; i++) {
283: M._rows.add(_rows.get(i).plus(
284: _transposed ? that.getColumn(i) : that.getRow(i)));
285: }
286: return M;
287: }
288:
289: @Override
290: public SparseMatrix<F> minus(Matrix<F> that) { // Returns more specialized type.
291: return this .plus(that.opposite());
292: }
293:
294: @Override
295: public SparseMatrix<F> times(F k) {
296: SparseMatrix<F> M = SparseMatrix.newInstance(_n, _zero,
297: _transposed);
298: for (int i = 0, p = _rows.size(); i < p; i++) {
299: M._rows.add(_rows.get(i).times(k));
300: }
301: return M;
302: }
303:
304: @Override
305: public SparseVector<F> times(Vector<F> v) {
306: if (v.getDimension() != this .getNumberOfColumns())
307: throw new DimensionException();
308: final int m = this .getNumberOfRows();
309: SparseVector<F> V = SparseVector.newInstance(m, _zero);
310: for (int i = 0; i < m; i++) {
311: F e = this .getRow(i).times(v);
312: if (!_zero.equals(e)) {
313: V._elements.put(Index.valueOf(i), e);
314: }
315: }
316: return V;
317: }
318:
319: @Override
320: public SparseMatrix<F> times(Matrix<F> that) {
321: final int m = this .getNumberOfRows();
322: final int n = this .getNumberOfColumns();
323: final int p = that.getNumberOfColumns();
324: if (that.getNumberOfRows() != n)
325: throw new DimensionException();
326: // Creates a mxp matrix in transposed form (p columns vectors of size m)
327: FastTable<SparseVector<F>> rows = this .getRows();
328: SparseMatrix<F> M = SparseMatrix.newInstance(m, _zero, true);
329: for (int j = 0; j < p; j++) {
330: Vector<F> thatColj = that.getColumn(j);
331: SparseVector<F> column = SparseVector.newInstance(m, _zero);
332: M._rows.add(column); // M is transposed.
333: for (int i = 0; i < m; i++) {
334: F e = rows.get(i).times(thatColj);
335: if (!_zero.equals(e)) {
336: column._elements.put(Index.valueOf(i), e);
337: }
338: }
339: }
340: return M;
341: }
342:
343: private FastTable<SparseVector<F>> getRows() {
344: if (!_transposed)
345: return _rows;
346: FastTable<SparseVector<F>> rows = FastTable.newInstance();
347: for (int i = 0; i < _n; i++) {
348: rows.add(this .getRow(i));
349: }
350: return rows;
351: }
352:
353: @Override
354: public SparseMatrix<F> inverse() {
355: if (!isSquare())
356: throw new DimensionException("Matrix not square");
357: F detInv = this .determinant().inverse();
358: SparseMatrix<F> A = this .adjoint();
359: // Multiply adjoint elements with 1 / determinant.
360: for (int i = 0, m = A._rows.size(); i < m; i++) {
361: SparseVector<F> row = A._rows.get(i);
362: for (FastMap.Entry<Index, F> e = row._elements.head(), end = row._elements
363: .tail(); (e = e.getNext()) != end;) {
364: F element = e.getValue();
365: e.setValue(detInv.times(element));
366: }
367: }
368: return A;
369: }
370:
371: @Override
372: public F determinant() {
373: if (!isSquare())
374: throw new DimensionException("Matrix not square");
375: if (_n == 1)
376: return this .get(0, 0);
377: // Expansion by minors (also known as Laplacian)
378: // This algorithm is division free but too slow for dense matrix.
379: SparseVector<F> row0 = this .getRow(0);
380: F det = null;
381: for (FastMap.Entry<Index, F> e = row0._elements.head(), end = row0._elements
382: .tail(); (e = e.getNext()) != end;) {
383: int i = e.getKey().intValue();
384: F d = e.getValue().times(cofactor(0, i));
385: if (i % 2 != 0) {
386: d = d.opposite();
387: }
388: det = (det == null) ? d : det.plus(d);
389: }
390: return det == null ? _zero : det;
391: }
392:
393: @Override
394: public Matrix<F> solve(Matrix<F> y) {
395: return this .inverse().times(y);
396: }
397:
398: @Override
399: public SparseMatrix<F> transpose() {
400: SparseMatrix<F> M = SparseMatrix.newInstance(_n, _zero,
401: !_transposed);
402: M._rows.addAll(this ._rows);
403: return M;
404: }
405:
406: @Override
407: public F cofactor(int i, int j) {
408: if (_transposed) {
409: int k = i;
410: i = j;
411: j = k; // Swaps i,j
412: }
413: int m = _rows.size();
414: SparseMatrix<F> M = SparseMatrix.newInstance(m - 1, _zero,
415: _transposed);
416: for (int k1 = 0; k1 < m; k1++) {
417: if (k1 == i)
418: continue;
419: SparseVector<F> row = _rows.get(k1);
420: SparseVector<F> V = SparseVector.newInstance(_n - 1, _zero);
421: M._rows.add(V);
422: for (FastMap.Entry<Index, F> e = row._elements.head(), end = row._elements
423: .tail(); (e = e.getNext()) != end;) {
424: int index = e.getKey().intValue();
425: if (index < j) {
426: V._elements.put(e.getKey(), e.getValue());
427: } else if (index > j) { // Position shifted (index minus one).
428: V._elements.put(Index.valueOf(index - 1), e
429: .getValue());
430: } // Else don't copy element at j.
431: }
432: }
433: return M.determinant();
434: }
435:
436: @Override
437: public SparseMatrix<F> adjoint() {
438: SparseMatrix<F> M = SparseMatrix.newInstance(_n, _zero,
439: _transposed);
440: int m = _rows.size();
441: for (int i = 0; i < m; i++) {
442: SparseVector<F> row = SparseVector.newInstance(_n, _zero);
443: M._rows.add(row);
444: for (int j = 0; j < _n; j++) {
445: F cofactor = _transposed ? cofactor(j, i) : cofactor(i,
446: j);
447: if (!_zero.equals(cofactor)) {
448: row._elements.put(Index.valueOf(j),
449: ((i + j) % 2 == 0) ? cofactor : cofactor
450: .opposite());
451: }
452: }
453: }
454: return M.transpose();
455: }
456:
457: @Override
458: public SparseMatrix<F> tensor(Matrix<F> that) {
459: final int this m = this .getNumberOfRows();
460: final int this n = this .getNumberOfColumns();
461: final int thatm = that.getNumberOfRows();
462: final int thatn = that.getNumberOfColumns();
463: int n = this n * thatn; // Number of columns,
464: int m = this m * thatm; // Number of rows.
465: SparseMatrix<F> M = SparseMatrix.newInstance(n, _zero, false);
466: for (int i = 0; i < m; i++) { // Row index.
467: final int i_rem_thatm = i % thatm;
468: final int i_div_thatm = i / thatm;
469: SparseVector<F> row = SparseVector.newInstance(n, _zero);
470: M._rows.add(row);
471: SparseVector<F> this Row = this .getRow(i_div_thatm);
472: for (FastMap.Entry<Index, F> e = this Row._elements.head(), end = this Row._elements
473: .tail(); (e = e.getNext()) != end;) {
474: F a = e.getValue();
475: int j = e.getKey().intValue();
476: for (int k = 0; k < thatn; k++) {
477: F b = that.get(i_rem_thatm, k);
478: if (!b.equals(_zero)) {
479: row._elements.put(Index.valueOf(j * thatn + k),
480: a.times(b));
481: }
482: }
483: }
484: }
485: return M;
486: }
487:
488: @Override
489: public SparseVector<F> vectorization() {
490: SparseVector<F> V = SparseVector.newInstance(_n
491: * this .getNumberOfRows(), _zero);
492: int offset = 0;
493: for (int j = 0, n = this .getNumberOfColumns(); j < n; j++) {
494: SparseVector<F> column = this .getColumn(j);
495: for (FastMap.Entry<Index, F> e = column._elements.head(), end = column._elements
496: .tail(); (e = e.getNext()) != end;) {
497: V._elements.put(Index.valueOf(e.getKey().intValue()
498: + offset), e.getValue());
499: }
500: offset += this .getNumberOfRows();
501: }
502: return V;
503: }
504:
505: @SuppressWarnings("unchecked")
506: @Override
507: public SparseMatrix<F> copy() {
508: SparseMatrix<F> M = newInstance(_n, (F) _zero.copy(),
509: _transposed);
510: for (SparseVector<F> row : _rows) {
511: M._rows.add(row.copy());
512: }
513: return M;
514: }
515:
516: ///////////////////////
517: // Factory creation. //
518: ///////////////////////
519:
520: @SuppressWarnings("unchecked")
521: static <F extends Field<F>> SparseMatrix<F> newInstance(int n,
522: F zero, boolean transposed) {
523: SparseMatrix<F> M = FACTORY.object();
524: M._n = n;
525: M._zero = zero;
526: M._transposed = transposed;
527: return M;
528: }
529:
530: private static final ObjectFactory<SparseMatrix> FACTORY = new ObjectFactory<SparseMatrix>() {
531: @Override
532: protected SparseMatrix create() {
533: return new SparseMatrix();
534: }
535:
536: @Override
537: protected void cleanup(SparseMatrix matrix) {
538: matrix._rows.reset();
539: }
540: };
541:
542: private SparseMatrix() {
543: }
544:
545: private static final long serialVersionUID = 1L;
546:
547: }
|