001    /*
002     * Licensed to the Apache Software Foundation (ASF) under one or more
003     * contributor license agreements.  See the NOTICE file distributed with
004     * this work for additional information regarding copyright ownership.
005     * The ASF licenses this file to You under the Apache License, Version 2.0
006     * (the "License"); you may not use this file except in compliance with
007     * the License.  You may obtain a copy of the License at
008     *
009     *      http://www.apache.org/licenses/LICENSE-2.0
010     *
011     * Unless required by applicable law or agreed to in writing, software
012     * distributed under the License is distributed on an "AS IS" BASIS,
013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014     * See the License for the specific language governing permissions and
015     * limitations under the License.
016     */
017    
018    package org.apache.commons.math.linear;
019    
020    import java.io.IOException;
021    import java.io.ObjectInputStream;
022    import java.io.ObjectOutputStream;
023    import java.lang.reflect.Array;
024    import java.util.Arrays;
025    
026    import org.apache.commons.math.Field;
027    import org.apache.commons.math.FieldElement;
028    import org.apache.commons.math.exception.MathArithmeticException;
029    import org.apache.commons.math.exception.MathIllegalArgumentException;
030    import org.apache.commons.math.exception.OutOfRangeException;
031    import org.apache.commons.math.exception.NoDataException;
032    import org.apache.commons.math.exception.NumberIsTooSmallException;
033    import org.apache.commons.math.exception.NullArgumentException;
034    import org.apache.commons.math.exception.DimensionMismatchException;
035    import org.apache.commons.math.exception.ZeroException;
036    import org.apache.commons.math.exception.util.LocalizedFormats;
037    import org.apache.commons.math.fraction.BigFraction;
038    import org.apache.commons.math.fraction.Fraction;
039    import org.apache.commons.math.util.FastMath;
040    import org.apache.commons.math.util.MathUtils;
041    
042    /**
043     * A collection of static methods that operate on or return matrices.
044     *
045     * @version $Id: MatrixUtils.java 1178009 2011-10-01 15:15:00Z luc $
046     */
047    public class MatrixUtils {
048    
049        /**
050         * Private constructor.
051         */
052        private MatrixUtils() {
053            super();
054        }
055    
056        /**
057         * Returns a {@link RealMatrix} with specified dimensions.
058         * <p>The type of matrix returned depends on the dimension. Below
059         * 2<sup>12</sup> elements (i.e. 4096 elements or 64&times;64 for a
060         * square matrix) which can be stored in a 32kB array, a {@link
061         * Array2DRowRealMatrix} instance is built. Above this threshold a {@link
062         * BlockRealMatrix} instance is built.</p>
063         * <p>The matrix elements are all set to 0.0.</p>
064         * @param rows number of rows of the matrix
065         * @param columns number of columns of the matrix
066         * @return  RealMatrix with specified dimensions
067         * @see #createRealMatrix(double[][])
068         */
069        public static RealMatrix createRealMatrix(final int rows, final int columns) {
070            return (rows * columns <= 4096) ?
071                    new Array2DRowRealMatrix(rows, columns) : new BlockRealMatrix(rows, columns);
072        }
073    
074        /**
075         * Returns a {@link FieldMatrix} with specified dimensions.
076         * <p>The type of matrix returned depends on the dimension. Below
077         * 2<sup>12</sup> elements (i.e. 4096 elements or 64&times;64 for a
078         * square matrix), a {@link FieldMatrix} instance is built. Above
079         * this threshold a {@link BlockFieldMatrix} instance is built.</p>
080         * <p>The matrix elements are all set to field.getZero().</p>
081         * @param <T> the type of the field elements
082         * @param field field to which the matrix elements belong
083         * @param rows number of rows of the matrix
084         * @param columns number of columns of the matrix
085         * @return  FieldMatrix with specified dimensions
086         * @see #createFieldMatrix(FieldElement[][])
087         * @since 2.0
088         */
089        public static <T extends FieldElement<T>> FieldMatrix<T> createFieldMatrix(final Field<T> field,
090                                                                                   final int rows,
091                                                                                   final int columns) {
092            return (rows * columns <= 4096) ?
093                    new Array2DRowFieldMatrix<T>(field, rows, columns) : new BlockFieldMatrix<T>(field, rows, columns);
094        }
095    
096        /**
097         * Returns a {@link RealMatrix} whose entries are the the values in the
098         * the input array.
099         * <p>The type of matrix returned depends on the dimension. Below
100         * 2<sup>12</sup> elements (i.e. 4096 elements or 64&times;64 for a
101         * square matrix) which can be stored in a 32kB array, a {@link
102         * Array2DRowRealMatrix} instance is built. Above this threshold a {@link
103         * BlockRealMatrix} instance is built.</p>
104         * <p>The input array is copied, not referenced.</p>
105         *
106         * @param data input array
107         * @return  RealMatrix containing the values of the array
108         * @throws org.apache.commons.math.exception.DimensionMismatchException
109         * if {@code data} is not rectangular (not all rows have the same length).
110         * @throws NoDataException if a row or column is empty.
111         * @throws NullArgumentException if either {@code data} or {@code data[0]}
112         * is {@code null}.
113         * @see #createRealMatrix(int, int)
114         */
115        public static RealMatrix createRealMatrix(double[][] data) {
116            if (data == null ||
117                data[0] == null) {
118                throw new NullArgumentException();
119            }
120            return (data.length * data[0].length <= 4096) ?
121                    new Array2DRowRealMatrix(data) : new BlockRealMatrix(data);
122        }
123    
124        /**
125         * Returns a {@link FieldMatrix} whose entries are the the values in the
126         * the input array.
127         * <p>The type of matrix returned depends on the dimension. Below
128         * 2<sup>12</sup> elements (i.e. 4096 elements or 64&times;64 for a
129         * square matrix), a {@link FieldMatrix} instance is built. Above
130         * this threshold a {@link BlockFieldMatrix} instance is built.</p>
131         * <p>The input array is copied, not referenced.</p>
132         * @param <T> the type of the field elements
133         * @param data input array
134         * @return a matrix containing the values of the array.
135         * @throws org.apache.commons.math.exception.DimensionMismatchException
136         * if {@code data} is not rectangular (not all rows have the same length).
137         * @throws NoDataException if a row or column is empty.
138         * @throws NullArgumentException if either {@code data} or {@code data[0]}
139         * is {@code null}.
140         * @see #createFieldMatrix(Field, int, int)
141         * @since 2.0
142         */
143        public static <T extends FieldElement<T>> FieldMatrix<T> createFieldMatrix(T[][] data) {
144            if (data == null ||
145                data[0] == null) {
146                throw new NullArgumentException();
147            }
148            return (data.length * data[0].length <= 4096) ?
149                    new Array2DRowFieldMatrix<T>(data) : new BlockFieldMatrix<T>(data);
150        }
151    
152        /**
153         * Returns <code>dimension x dimension</code> identity matrix.
154         *
155         * @param dimension dimension of identity matrix to generate
156         * @return identity matrix
157         * @throws IllegalArgumentException if dimension is not positive
158         * @since 1.1
159         */
160        public static RealMatrix createRealIdentityMatrix(int dimension) {
161            final RealMatrix m = createRealMatrix(dimension, dimension);
162            for (int i = 0; i < dimension; ++i) {
163                m.setEntry(i, i, 1.0);
164            }
165            return m;
166        }
167    
168        /**
169         * Returns <code>dimension x dimension</code> identity matrix.
170         *
171         * @param <T> the type of the field elements
172         * @param field field to which the elements belong
173         * @param dimension dimension of identity matrix to generate
174         * @return identity matrix
175         * @throws IllegalArgumentException if dimension is not positive
176         * @since 2.0
177         */
178        public static <T extends FieldElement<T>> FieldMatrix<T>
179            createFieldIdentityMatrix(final Field<T> field, final int dimension) {
180            final T zero = field.getZero();
181            final T one  = field.getOne();
182            @SuppressWarnings("unchecked")
183            final T[][] d = (T[][]) Array.newInstance(field.getRuntimeClass(), new int[] { dimension, dimension });
184            for (int row = 0; row < dimension; row++) {
185                final T[] dRow = d[row];
186                Arrays.fill(dRow, zero);
187                dRow[row] = one;
188            }
189            return new Array2DRowFieldMatrix<T>(field, d, false);
190        }
191    
192        /**
193         * Returns a diagonal matrix with specified elements.
194         *
195         * @param diagonal diagonal elements of the matrix (the array elements
196         * will be copied)
197         * @return diagonal matrix
198         * @since 2.0
199         */
200        public static RealMatrix createRealDiagonalMatrix(final double[] diagonal) {
201            final RealMatrix m = createRealMatrix(diagonal.length, diagonal.length);
202            for (int i = 0; i < diagonal.length; ++i) {
203                m.setEntry(i, i, diagonal[i]);
204            }
205            return m;
206        }
207    
208        /**
209         * Returns a diagonal matrix with specified elements.
210         *
211         * @param <T> the type of the field elements
212         * @param diagonal diagonal elements of the matrix (the array elements
213         * will be copied)
214         * @return diagonal matrix
215         * @since 2.0
216         */
217        public static <T extends FieldElement<T>> FieldMatrix<T>
218            createFieldDiagonalMatrix(final T[] diagonal) {
219            final FieldMatrix<T> m =
220                createFieldMatrix(diagonal[0].getField(), diagonal.length, diagonal.length);
221            for (int i = 0; i < diagonal.length; ++i) {
222                m.setEntry(i, i, diagonal[i]);
223            }
224            return m;
225        }
226    
227        /**
228         * Creates a {@link RealVector} using the data from the input array.
229         *
230         * @param data the input data
231         * @return a data.length RealVector
232         * @throws NoDataException if {@code data} is empty.
233         * @throws NullArgumentException if {@code data} is {@code null}.
234         */
235        public static RealVector createRealVector(double[] data) {
236            if (data == null) {
237                throw new NullArgumentException();
238            }
239            return new ArrayRealVector(data, true);
240        }
241    
242        /**
243         * Creates a {@link FieldVector} using the data from the input array.
244         *
245         * @param <T> the type of the field elements
246         * @param data the input data
247         * @return a data.length FieldVector
248         * @throws NoDataException if {@code data} is empty.
249         * @throws NullArgumentException if {@code data} is {@code null}.
250         * @throws ZeroException if {@code data} has 0 elements
251         */
252        public static <T extends FieldElement<T>> FieldVector<T> createFieldVector(final T[] data) {
253            if (data == null) {
254                throw new NullArgumentException();
255            }
256            if (data.length == 0) {
257                throw new ZeroException(LocalizedFormats.VECTOR_MUST_HAVE_AT_LEAST_ONE_ELEMENT);
258            }
259            return new ArrayFieldVector<T>(data[0].getField(), data, true);
260        }
261    
262        /**
263         * Create a row {@link RealMatrix} using the data from the input
264         * array.
265         *
266         * @param rowData the input row data
267         * @return a 1 x rowData.length RealMatrix
268         * @throws NoDataException if {@code rowData} is empty.
269         * @throws NullArgumentException if {@code rowData} is {@code null}.
270         */
271        public static RealMatrix createRowRealMatrix(double[] rowData) {
272            if (rowData == null) {
273                throw new NullArgumentException();
274            }
275            final int nCols = rowData.length;
276            final RealMatrix m = createRealMatrix(1, nCols);
277            for (int i = 0; i < nCols; ++i) {
278                m.setEntry(0, i, rowData[i]);
279            }
280            return m;
281        }
282    
283        /**
284         * Create a row {@link FieldMatrix} using the data from the input
285         * array.
286         *
287         * @param <T> the type of the field elements
288         * @param rowData the input row data
289         * @return a 1 x rowData.length FieldMatrix
290         * @throws NoDataException if {@code rowData} is empty.
291         * @throws NullArgumentException if {@code rowData} is {@code null}.
292         */
293        public static <T extends FieldElement<T>> FieldMatrix<T>
294            createRowFieldMatrix(final T[] rowData) {
295            if (rowData == null) {
296                throw new NullArgumentException();
297            }
298            final int nCols = rowData.length;
299            if (nCols == 0) {
300                throw new NoDataException(LocalizedFormats.AT_LEAST_ONE_COLUMN);
301            }
302            final FieldMatrix<T> m = createFieldMatrix(rowData[0].getField(), 1, nCols);
303            for (int i = 0; i < nCols; ++i) {
304                m.setEntry(0, i, rowData[i]);
305            }
306            return m;
307        }
308    
309        /**
310         * Creates a column {@link RealMatrix} using the data from the input
311         * array.
312         *
313         * @param columnData  the input column data
314         * @return a columnData x 1 RealMatrix
315         * @throws NoDataException if {@code columnData} is empty.
316         * @throws NullArgumentException if {@code columnData} is {@code null}.
317         */
318        public static RealMatrix createColumnRealMatrix(double[] columnData) {
319            if (columnData == null) {
320                throw new NullArgumentException();
321            }
322            final int nRows = columnData.length;
323            final RealMatrix m = createRealMatrix(nRows, 1);
324            for (int i = 0; i < nRows; ++i) {
325                m.setEntry(i, 0, columnData[i]);
326            }
327            return m;
328        }
329    
330        /**
331         * Creates a column {@link FieldMatrix} using the data from the input
332         * array.
333         *
334         * @param <T> the type of the field elements
335         * @param columnData  the input column data
336         * @return a columnData x 1 FieldMatrix
337         * @throws NoDataException if {@code data} is empty.
338         * @throws NullArgumentException if {@code columnData} is {@code null}.
339         */
340        public static <T extends FieldElement<T>> FieldMatrix<T>
341            createColumnFieldMatrix(final T[] columnData) {
342            if (columnData == null) {
343                throw new NullArgumentException();
344            }
345            final int nRows = columnData.length;
346            if (nRows == 0) {
347                throw new NoDataException(LocalizedFormats.AT_LEAST_ONE_ROW);
348            }
349            final FieldMatrix<T> m = createFieldMatrix(columnData[0].getField(), nRows, 1);
350            for (int i = 0; i < nRows; ++i) {
351                m.setEntry(i, 0, columnData[i]);
352            }
353            return m;
354        }
355    
356        /**
357         * Check if matrix indices are valid.
358         *
359         * @param m Matrix.
360         * @param row Row index to check.
361         * @param column Column index to check.
362         * @throws OutOfRangeException if {@code row} or {@code column} is not
363         * a valid index.
364         */
365        public static void checkMatrixIndex(final AnyMatrix m,
366                                            final int row, final int column) {
367            checkRowIndex(m, row);
368            checkColumnIndex(m, column);
369        }
370    
371        /**
372         * Check if a row index is valid.
373         *
374         * @param m Matrix.
375         * @param row Row index to check.
376         * @throws OutOfRangeException if {@code row} is not a valid index.
377         */
378        public static void checkRowIndex(final AnyMatrix m, final int row) {
379            if (row < 0 ||
380                row >= m.getRowDimension()) {
381                throw new OutOfRangeException(LocalizedFormats.ROW_INDEX,
382                                              row, 0, m.getRowDimension() - 1);
383            }
384        }
385    
386        /**
387         * Check if a column index is valid.
388         *
389         * @param m Matrix.
390         * @param column Column index to check.
391         * @throws OutOfRangeException if {@code column} is not a valid index.
392         */
393        public static void checkColumnIndex(final AnyMatrix m, final int column) {
394            if (column < 0 || column >= m.getColumnDimension()) {
395                throw new OutOfRangeException(LocalizedFormats.COLUMN_INDEX,
396                                               column, 0, m.getColumnDimension() - 1);
397            }
398        }
399    
400        /**
401         * Check if submatrix ranges indices are valid.
402         * Rows and columns are indicated counting from 0 to {@code n - 1}.
403         *
404         * @param m Matrix.
405         * @param startRow Initial row index.
406         * @param endRow Final row index.
407         * @param startColumn Initial column index.
408         * @param endColumn Final column index.
409         * @throws OutOfRangeException if the indices are invalid.
410         * @throws NumberIsTooSmallException if {@code endRow < startRow} or
411         * {@code endColumn < startColumn}.
412         */
413        public static void checkSubMatrixIndex(final AnyMatrix m,
414                                               final int startRow, final int endRow,
415                                               final int startColumn, final int endColumn) {
416            checkRowIndex(m, startRow);
417            checkRowIndex(m, endRow);
418            if (endRow < startRow) {
419                throw new NumberIsTooSmallException(LocalizedFormats.INITIAL_ROW_AFTER_FINAL_ROW,
420                                                    endRow, startRow, false);
421            }
422    
423            checkColumnIndex(m, startColumn);
424            checkColumnIndex(m, endColumn);
425            if (endColumn < startColumn) {
426                throw new NumberIsTooSmallException(LocalizedFormats.INITIAL_COLUMN_AFTER_FINAL_COLUMN,
427                                                    endColumn, startColumn, false);
428            }
429    
430    
431        }
432    
433        /**
434         * Check if submatrix ranges indices are valid.
435         * Rows and columns are indicated counting from 0 to n-1.
436         *
437         * @param m Matrix.
438         * @param selectedRows Array of row indices.
439         * @param selectedColumns Array of column indices.
440         * @throws NullArgumentException if {@code selectedRows} or
441         * {@code selectedColumns} are {@code null}.
442         * @throws NoDataException if the row or column selections are empty (zero
443         * length).
444         * @throws OutOfRangeException if row or column selections are not valid.
445         */
446        public static void checkSubMatrixIndex(final AnyMatrix m,
447                                               final int[] selectedRows,
448                                               final int[] selectedColumns) {
449            if (selectedRows == null) {
450                throw new NullArgumentException();
451            }
452            if (selectedColumns == null) {
453                throw new NullArgumentException();
454            }
455            if (selectedRows.length == 0) {
456                throw new NoDataException(LocalizedFormats.EMPTY_SELECTED_ROW_INDEX_ARRAY);
457            }
458            if (selectedColumns.length == 0) {
459                throw new NoDataException(LocalizedFormats.EMPTY_SELECTED_COLUMN_INDEX_ARRAY);
460            }
461    
462            for (final int row : selectedRows) {
463                checkRowIndex(m, row);
464            }
465            for (final int column : selectedColumns) {
466                checkColumnIndex(m, column);
467            }
468        }
469    
470        /**
471         * Check if matrices are addition compatible.
472         *
473         * @param left Left hand side matrix.
474         * @param right Right hand side matrix.
475         * @throws MatrixDimensionMismatchException if the matrices are not addition compatible.
476         */
477        public static void checkAdditionCompatible(final AnyMatrix left, final AnyMatrix right) {
478            if ((left.getRowDimension()    != right.getRowDimension()) ||
479                (left.getColumnDimension() != right.getColumnDimension())) {
480                throw new MatrixDimensionMismatchException(left.getRowDimension(), left.getColumnDimension(),
481                                                           right.getRowDimension(), right.getColumnDimension());
482            }
483        }
484    
485        /**
486         * Check if matrices are subtraction compatible
487         *
488         * @param left Left hand side matrix.
489         * @param right Right hand side matrix.
490         * @throws MatrixDimensionMismatchException if the matrices are not addition compatible.
491         */
492        public static void checkSubtractionCompatible(final AnyMatrix left, final AnyMatrix right) {
493            if ((left.getRowDimension()    != right.getRowDimension()) ||
494                (left.getColumnDimension() != right.getColumnDimension())) {
495                throw new MatrixDimensionMismatchException(left.getRowDimension(), left.getColumnDimension(),
496                                                           right.getRowDimension(), right.getColumnDimension());
497            }
498        }
499    
500        /**
501         * Check if matrices are multiplication compatible
502         *
503         * @param left Left hand side matrix.
504         * @param right Right hand side matrix.
505         * @throws DimensionMismatchException if matrices are not multiplication compatible.
506         */
507        public static void checkMultiplicationCompatible(final AnyMatrix left, final AnyMatrix right) {
508            if (left.getColumnDimension() != right.getRowDimension()) {
509                throw new DimensionMismatchException(left.getColumnDimension(),
510                                                     right.getRowDimension());
511            }
512        }
513    
514        /**
515         * Convert a {@link FieldMatrix}/{@link Fraction} matrix to a {@link RealMatrix}.
516         * @param m Matrix to convert.
517         * @return the converted matrix.
518         */
519        public static Array2DRowRealMatrix fractionMatrixToRealMatrix(final FieldMatrix<Fraction> m) {
520            final FractionMatrixConverter converter = new FractionMatrixConverter();
521            m.walkInOptimizedOrder(converter);
522            return converter.getConvertedMatrix();
523        }
524    
525        /** Converter for {@link FieldMatrix}/{@link Fraction}. */
526        private static class FractionMatrixConverter extends DefaultFieldMatrixPreservingVisitor<Fraction> {
527            /** Converted array. */
528            private double[][] data;
529            /** Simple constructor. */
530            public FractionMatrixConverter() {
531                super(Fraction.ZERO);
532            }
533    
534            /** {@inheritDoc} */
535            @Override
536            public void start(int rows, int columns,
537                              int startRow, int endRow, int startColumn, int endColumn) {
538                data = new double[rows][columns];
539            }
540    
541            /** {@inheritDoc} */
542            @Override
543            public void visit(int row, int column, Fraction value) {
544                data[row][column] = value.doubleValue();
545            }
546    
547            /**
548             * Get the converted matrix.
549             *
550             * @return the converted matrix.
551             */
552            Array2DRowRealMatrix getConvertedMatrix() {
553                return new Array2DRowRealMatrix(data, false);
554            }
555    
556        }
557    
558        /**
559         * Convert a {@link FieldMatrix}/{@link BigFraction} matrix to a {@link RealMatrix}.
560         *
561         * @param m Matrix to convert.
562         * @return the converted matrix.
563         */
564        public static Array2DRowRealMatrix bigFractionMatrixToRealMatrix(final FieldMatrix<BigFraction> m) {
565            final BigFractionMatrixConverter converter = new BigFractionMatrixConverter();
566            m.walkInOptimizedOrder(converter);
567            return converter.getConvertedMatrix();
568        }
569    
570        /** Converter for {@link FieldMatrix}/{@link BigFraction}. */
571        private static class BigFractionMatrixConverter extends DefaultFieldMatrixPreservingVisitor<BigFraction> {
572            /** Converted array. */
573            private double[][] data;
574            /** Simple constructor. */
575            public BigFractionMatrixConverter() {
576                super(BigFraction.ZERO);
577            }
578    
579            /** {@inheritDoc} */
580            @Override
581            public void start(int rows, int columns,
582                              int startRow, int endRow, int startColumn, int endColumn) {
583                data = new double[rows][columns];
584            }
585    
586            /** {@inheritDoc} */
587            @Override
588            public void visit(int row, int column, BigFraction value) {
589                data[row][column] = value.doubleValue();
590            }
591    
592            /**
593             * Get the converted matrix.
594             *
595             * @return the converted matrix.
596             */
597            Array2DRowRealMatrix getConvertedMatrix() {
598                return new Array2DRowRealMatrix(data, false);
599            }
600        }
601    
602        /** Serialize a {@link RealVector}.
603         * <p>
604         * This method is intended to be called from within a private
605         * <code>writeObject</code> method (after a call to
606         * <code>oos.defaultWriteObject()</code>) in a class that has a
607         * {@link RealVector} field, which should be declared <code>transient</code>.
608         * This way, the default handling does not serialize the vector (the {@link
609         * RealVector} interface is not serializable by default) but this method does
610         * serialize it specifically.
611         * </p>
612         * <p>
613         * The following example shows how a simple class with a name and a real vector
614         * should be written:
615         * <pre><code>
616         * public class NamedVector implements Serializable {
617         *
618         *     private final String name;
619         *     private final transient RealVector coefficients;
620         *
621         *     // omitted constructors, getters ...
622         *
623         *     private void writeObject(ObjectOutputStream oos) throws IOException {
624         *         oos.defaultWriteObject();  // takes care of name field
625         *         MatrixUtils.serializeRealVector(coefficients, oos);
626         *     }
627         *
628         *     private void readObject(ObjectInputStream ois) throws ClassNotFoundException, IOException {
629         *         ois.defaultReadObject();  // takes care of name field
630         *         MatrixUtils.deserializeRealVector(this, "coefficients", ois);
631         *     }
632         *
633         * }
634         * </code></pre>
635         * </p>
636         *
637         * @param vector real vector to serialize
638         * @param oos stream where the real vector should be written
639         * @exception IOException if object cannot be written to stream
640         * @see #deserializeRealVector(Object, String, ObjectInputStream)
641         */
642        public static void serializeRealVector(final RealVector vector,
643                                               final ObjectOutputStream oos)
644            throws IOException {
645            final int n = vector.getDimension();
646            oos.writeInt(n);
647            for (int i = 0; i < n; ++i) {
648                oos.writeDouble(vector.getEntry(i));
649            }
650        }
651    
652        /** Deserialize  a {@link RealVector} field in a class.
653         * <p>
654         * This method is intended to be called from within a private
655         * <code>readObject</code> method (after a call to
656         * <code>ois.defaultReadObject()</code>) in a class that has a
657         * {@link RealVector} field, which should be declared <code>transient</code>.
658         * This way, the default handling does not deserialize the vector (the {@link
659         * RealVector} interface is not serializable by default) but this method does
660         * deserialize it specifically.
661         * </p>
662         * @param instance instance in which the field must be set up
663         * @param fieldName name of the field within the class (may be private and final)
664         * @param ois stream from which the real vector should be read
665         * @exception ClassNotFoundException if a class in the stream cannot be found
666         * @exception IOException if object cannot be read from the stream
667         * @see #serializeRealVector(RealVector, ObjectOutputStream)
668         */
669        public static void deserializeRealVector(final Object instance,
670                                                 final String fieldName,
671                                                 final ObjectInputStream ois)
672          throws ClassNotFoundException, IOException {
673            try {
674    
675                // read the vector data
676                final int n = ois.readInt();
677                final double[] data = new double[n];
678                for (int i = 0; i < n; ++i) {
679                    data[i] = ois.readDouble();
680                }
681    
682                // create the instance
683                final RealVector vector = new ArrayRealVector(data, false);
684    
685                // set up the field
686                final java.lang.reflect.Field f =
687                    instance.getClass().getDeclaredField(fieldName);
688                f.setAccessible(true);
689                f.set(instance, vector);
690    
691            } catch (NoSuchFieldException nsfe) {
692                IOException ioe = new IOException();
693                ioe.initCause(nsfe);
694                throw ioe;
695            } catch (IllegalAccessException iae) {
696                IOException ioe = new IOException();
697                ioe.initCause(iae);
698                throw ioe;
699            }
700    
701        }
702    
703        /** Serialize a {@link RealMatrix}.
704         * <p>
705         * This method is intended to be called from within a private
706         * <code>writeObject</code> method (after a call to
707         * <code>oos.defaultWriteObject()</code>) in a class that has a
708         * {@link RealMatrix} field, which should be declared <code>transient</code>.
709         * This way, the default handling does not serialize the matrix (the {@link
710         * RealMatrix} interface is not serializable by default) but this method does
711         * serialize it specifically.
712         * </p>
713         * <p>
714         * The following example shows how a simple class with a name and a real matrix
715         * should be written:
716         * <pre><code>
717         * public class NamedMatrix implements Serializable {
718         *
719         *     private final String name;
720         *     private final transient RealMatrix coefficients;
721         *
722         *     // omitted constructors, getters ...
723         *
724         *     private void writeObject(ObjectOutputStream oos) throws IOException {
725         *         oos.defaultWriteObject();  // takes care of name field
726         *         MatrixUtils.serializeRealMatrix(coefficients, oos);
727         *     }
728         *
729         *     private void readObject(ObjectInputStream ois) throws ClassNotFoundException, IOException {
730         *         ois.defaultReadObject();  // takes care of name field
731         *         MatrixUtils.deserializeRealMatrix(this, "coefficients", ois);
732         *     }
733         *
734         * }
735         * </code></pre>
736         * </p>
737         *
738         * @param matrix real matrix to serialize
739         * @param oos stream where the real matrix should be written
740         * @exception IOException if object cannot be written to stream
741         * @see #deserializeRealMatrix(Object, String, ObjectInputStream)
742         */
743        public static void serializeRealMatrix(final RealMatrix matrix,
744                                               final ObjectOutputStream oos)
745            throws IOException {
746            final int n = matrix.getRowDimension();
747            final int m = matrix.getColumnDimension();
748            oos.writeInt(n);
749            oos.writeInt(m);
750            for (int i = 0; i < n; ++i) {
751                for (int j = 0; j < m; ++j) {
752                    oos.writeDouble(matrix.getEntry(i, j));
753                }
754            }
755        }
756    
757        /** Deserialize  a {@link RealMatrix} field in a class.
758         * <p>
759         * This method is intended to be called from within a private
760         * <code>readObject</code> method (after a call to
761         * <code>ois.defaultReadObject()</code>) in a class that has a
762         * {@link RealMatrix} field, which should be declared <code>transient</code>.
763         * This way, the default handling does not deserialize the matrix (the {@link
764         * RealMatrix} interface is not serializable by default) but this method does
765         * deserialize it specifically.
766         * </p>
767         * @param instance instance in which the field must be set up
768         * @param fieldName name of the field within the class (may be private and final)
769         * @param ois stream from which the real matrix should be read
770         * @exception ClassNotFoundException if a class in the stream cannot be found
771         * @exception IOException if object cannot be read from the stream
772         * @see #serializeRealMatrix(RealMatrix, ObjectOutputStream)
773         */
774        public static void deserializeRealMatrix(final Object instance,
775                                                 final String fieldName,
776                                                 final ObjectInputStream ois)
777          throws ClassNotFoundException, IOException {
778            try {
779    
780                // read the matrix data
781                final int n = ois.readInt();
782                final int m = ois.readInt();
783                final double[][] data = new double[n][m];
784                for (int i = 0; i < n; ++i) {
785                    final double[] dataI = data[i];
786                    for (int j = 0; j < m; ++j) {
787                        dataI[j] = ois.readDouble();
788                    }
789                }
790    
791                // create the instance
792                final RealMatrix matrix = new Array2DRowRealMatrix(data, false);
793    
794                // set up the field
795                final java.lang.reflect.Field f =
796                    instance.getClass().getDeclaredField(fieldName);
797                f.setAccessible(true);
798                f.set(instance, matrix);
799    
800            } catch (NoSuchFieldException nsfe) {
801                IOException ioe = new IOException();
802                ioe.initCause(nsfe);
803                throw ioe;
804            } catch (IllegalAccessException iae) {
805                IOException ioe = new IOException();
806                ioe.initCause(iae);
807                throw ioe;
808            }
809        }
810    
811        /**Solve  a  system of composed of a Lower Triangular Matrix
812         * {@link RealMatrix}.
813         * <p>
814         * This method is called to solve systems of equations which are
815         * of the lower triangular form. The matrix {@link RealMatrix}
816         * is assumed, though not checked, to be in lower triangular form.
817         * The vector {@link RealVector} is overwritten with the solution.
818         * The matrix is checked that it is square and its dimensions match
819         * the length of the vector.
820         * </p>
821         * @param rm RealMatrix which is lower triangular
822         * @param b  RealVector this is overwritten
823         * @exception IllegalArgumentException if the matrix and vector are not conformable
824         * @exception ArithmeticException there is a zero or near zero on the diagonal of rm
825         */
826        public static void solveLowerTriangularSystem( RealMatrix rm, RealVector b){
827            if ((rm == null) || (b == null) || ( rm.getRowDimension() != b.getDimension())) {
828                throw new MathIllegalArgumentException(LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE,
829                        (rm == null) ? 0 : rm.getRowDimension(),
830                        (b == null) ? 0 : b.getDimension());
831            }
832            if( rm.getColumnDimension() != rm.getRowDimension() ){
833                throw new MathIllegalArgumentException(LocalizedFormats.DIMENSIONS_MISMATCH_2x2,
834                        rm.getRowDimension(),rm.getRowDimension(),
835                        rm.getRowDimension(),rm.getColumnDimension());
836            }
837            int rows = rm.getRowDimension();
838            for( int i = 0 ; i < rows ; i++ ){
839                double diag = rm.getEntry(i, i);
840                if( FastMath.abs(diag) < MathUtils.SAFE_MIN ){
841                    throw new MathArithmeticException(LocalizedFormats.ZERO_DENOMINATOR);
842                }
843                double bi = b.getEntry(i)/diag;
844                b.setEntry(i,  bi );
845                for( int j = i+1; j< rows; j++ ){
846                    b.setEntry(j, b.getEntry(j)-bi*rm.getEntry(j,i)  );
847                }
848            }
849        }
850    
851        /** Solver a  system composed  of an Upper Triangular Matrix
852         * {@link RealMatrix}.
853         * <p>
854         * This method is called to solve systems of equations which are
855         * of the lower triangular form. The matrix {@link RealMatrix}
856         * is assumed, though not checked, to be in upper triangular form.
857         * The vector {@link RealVector} is overwritten with the solution.
858         * The matrix is checked that it is square and its dimensions match
859         * the length of the vector.
860         * </p>
861         * @param rm RealMatrix which is upper triangular
862         * @param b  RealVector this is overwritten
863         * @exception IllegalArgumentException if the matrix and vector are not conformable
864         * @exception ArithmeticException there is a zero or near zero on the diagonal of rm
865         */
866        public static void solveUpperTriangularSystem( RealMatrix rm, RealVector b){
867            if ((rm == null) || (b == null) || ( rm.getRowDimension() != b.getDimension())) {
868                throw new MathIllegalArgumentException(LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE,
869                        (rm == null) ? 0 : rm.getRowDimension(),
870                        (b == null) ? 0 : b.getDimension());
871            }
872            if( rm.getColumnDimension() != rm.getRowDimension() ){
873                throw new MathIllegalArgumentException(LocalizedFormats.DIMENSIONS_MISMATCH_2x2,
874                        rm.getRowDimension(),rm.getRowDimension(),
875                        rm.getRowDimension(),rm.getColumnDimension());
876            }
877            int rows = rm.getRowDimension();
878            for( int i = rows-1 ; i >-1 ; i-- ){
879                double diag = rm.getEntry(i, i);
880                if( FastMath.abs(diag) < MathUtils.SAFE_MIN ){
881                    throw new MathArithmeticException(LocalizedFormats.ZERO_DENOMINATOR);
882                }
883                double bi = b.getEntry(i)/diag;
884                b.setEntry(i,  bi );
885                for( int j = i-1; j>-1; j-- ){
886                    b.setEntry(j, b.getEntry(j)-bi*rm.getEntry(j,i)  );
887                }
888            }
889        }
890    }