package com.srbenoit.math.linear;

import java.io.BufferedReader;
import java.io.File;
import java.io.FileOutputStream;
import java.io.FileReader;
import java.io.IOException;
import java.io.PrintStream;
import java.text.DecimalFormat;
import java.util.Arrays;
import java.util.logging.Level;
import com.srbenoit.log.LoggedObject;
import com.srbenoit.util.EvaluationException;

/**
 * Manages a matrix of doubles, and includes ways for multiple processes to distribute the work of
 * populating the matrix. A process can request the next uncomputed array element, which marks the
 * element as in-progress. Includes methods to compute the L-U decomposition of the matrix, and
 * perform back-substitution using the result to solve systems of the form A x = b. This also
 * includes methods to save and reload the matrix state, so a lengthy process can be stopped and
 * restarted.
 */
public class DMatrix extends LoggedObject {

    /** characters used to read and write in hex */
    private static final char[] HEX = "0123456789ABCDEF".toCharArray();

    /** the matrix data */
    private final transient double[][] data;

    /** a set of flags for each data element */
    private final transient byte[][] flags;

    /** the row of the first uncompleted element */
    private transient int firstUncRow = 0;

    /** the column of the first uncompleted element */
    private transient int firstUncCol = 0;

    /** the L-U decomposition of the matrix */
    private transient double[][] lud;

    /** the permuted row indexes for the decomposed matrix */
    private transient int[] ludIndexes;

    /** <code>true</code> if an even number of row interchanges; <code>false</code> if odd */
    private transient boolean ludExchangesEven = true;

    /**
     * Constructs a new <code>DMatrix</code> whose values are initialized to zero.
     *
     * @param  numRows     the number of rows
     * @param  numColumns  the number of columns
     */
    public DMatrix(final int numRows, final int numColumns) {

        if (numRows < 1) {
            throw new IllegalArgumentException("Invalid number of rows: " + numRows);
        }

        if (numColumns < 1) {
            throw new IllegalArgumentException("Invalid number of columns: " + numColumns);
        }

        this.data = new double[numRows][numColumns];
        this.flags = new byte[numRows][numColumns];
    }

    /**
     * Constructs a new <code>DMatrix</code> whose values are initialized to provided values
     *
     * @param  numRows     the number of rows
     * @param  numColumns  the number of columns
     * @param  values      the initial values for the matrix, in row-major order
     */
    public DMatrix(final int numRows, final int numColumns, final double... values) {

        this(numRows, numColumns);

        int index;
        byte mask;

        mask = (byte) (0x01 << MatrixElementFlag.COMPLETED.flagBit());

        if (values.length != (numRows * numColumns)) {
            throw new IllegalArgumentException(
                "Number of supplied values does not match matrix size");
        }

        index = 0;

        for (int i = 0; i < numRows; i++) {

            for (int j = 0; j < numColumns; j++) {
                this.data[i][j] = values[index];
                this.flags[i][j] |= mask;
                index++;
            }
        }

        this.firstUncRow = this.data.length;
        this.firstUncCol = 0;
    }

    /**
     * Constructs a new <code>DMatrix</code> whose values are a copy of the values in another
     * <code>DMatrix</code>.
     *
     * @param  matrix  the matrix to copy
     */
    public DMatrix(final DMatrix matrix) {

        // Safe outside synch because matrix size is immutable
        this(matrix.data.length, matrix.data[0].length);

        synchronized (matrix) {

            for (int i = 0; i < this.data.length; i++) {
                System.arraycopy(this.data[i], 0, matrix.data[i], 0, this.data[i].length);
                System.arraycopy(this.flags[i], 0, matrix.flags[i], 0, this.flags[i].length);
            }

            this.firstUncRow = matrix.firstUncRow;
            this.firstUncCol = matrix.firstUncCol;
        }
    }

    /**
     * Gets the number of rows in this matrix.
     *
     * @return  the number of rows
     */
    public int numRows() {

        // Not synchronized because matrix size is immutable
        return this.data.length;
    }

    /**
     * Gets the number of columns in this matrix.
     *
     * @return  the number of columns
     */
    public int numColumns() {

        // Not synchronized because matrix size is immutable
        return this.data[0].length;
    }

    /**
     * Gets the row and column of the next uncompleted (and not currently in progress) element, and
     * mark it as in-progress.
     *
     * @return  a two-integer array containing the row and column of the element (<code>null</code>
     *          if there are no uncompleted elements in the matrix)
     */
    public int[] getUncompletedElement() {

        int[] rowCol = null;

        synchronized (this) {

            if (this.firstUncRow < this.data.length) {
                setFlag(this.firstUncRow, this.firstUncCol, MatrixElementFlag.IN_PROGRESS, true);
                rowCol = new int[] { this.firstUncRow, this.firstUncCol };
                advanceFirstUncompleted();
            }
        }

        return rowCol;
    }

    /**
     * Moves the index of the first uncompleted element forward to the next uncompleted (and not
     * in-progress) element.
     */
    private void advanceFirstUncompleted() {

        // NOTE: called only from within synchronized block.

        while (this.firstUncRow < this.data.length) {

            if (isCompleted(this.firstUncRow, this.firstUncCol)
                    || isInProgress(this.firstUncRow, this.firstUncCol)) {

                this.firstUncCol++;

                if (this.firstUncCol >= this.data[0].length) {
                    this.firstUncRow++;
                    this.firstUncCol = 0;
                }
            } else {
                break;
            }
        }
    }

    /**
     * Clears the value in a matrix cell and mark it as not completed, not in progress, and not out
     * of range.
     *
     * @param  row  the index of the data value's row
     * @param  col  the index of the data value's column
     */
    public void clearElementValue(final int row, final int col) {

        synchronized (this) {

            if (isCompleted(row, col)) {
                this.data[row][col] = 0;
                setFlag(row, col, MatrixElementFlag.COMPLETED, false);
                setFlag(row, col, MatrixElementFlag.IN_PROGRESS, false);
                setFlag(row, col, MatrixElementFlag.OUT_OF_RANGE, false);
                clearLud();

                // If this row/column is before the first uncompleted
                // row/column, move the first uncompleted one back to here.
                if (this.firstUncRow > row) {
                    this.firstUncRow = row;
                    this.firstUncCol = col;
                } else if ((this.firstUncRow == row) && (this.firstUncCol > col)) {
                    this.firstUncCol = col;
                }
            }
        }
    }

    /**
     * Copies values from another matrix.
     *
     * @param  src  a matrix from which to copy values
     */
    public void set(final DMatrix src) {

        synchronized (this) {

            if ((src.numRows() == numRows()) && (src.numColumns() == numColumns())) {

                for (int i = 0; i < this.data.length; i++) {
                    System.arraycopy(src.data[i], 0, this.data[i], 0, this.data[i].length);
                    System.arraycopy(src.flags[i], 0, this.flags[i], 0, this.flags[i].length);
                }

                this.firstUncRow = src.firstUncRow;
                this.firstUncCol = src.firstUncCol;
            } else {
                throw new IllegalArgumentException("Matrix size mismatch");
            }
        }
    }

    /**
     * Stores a value in a matrix cell and marks it as completed and no longer in progress.
     *
     * @param  row    the index of the data value's row
     * @param  col    the index of the data value's column
     * @param  value  the data value
     */
    public void set(final int row, final int col, final double value) {

        synchronized (this) {

            if (this.data[row][col] != value) {
                this.data[row][col] = value;
                clearLud();
            }

            setFlag(row, col, MatrixElementFlag.COMPLETED, true);
            setFlag(row, col, MatrixElementFlag.IN_PROGRESS, false);

            if ((row == this.firstUncRow) && (col == this.firstUncCol)) {
                advanceFirstUncompleted();
            }
        }
    }

    /**
     * Gets the data value of a matrix element.
     *
     * @param   row  the row (indexed from 0)
     * @param   col  the column (indexed from 0)
     * @return  the data value stored at that row/column
     */
    public double get(final int row, final int col) {

        synchronized (this) {
            return this.data[row][col];
        }
    }

    /**
     * Gets a particular flag value.
     *
     * @param   row        the row (indexed from 0)
     * @param   col        the column (indexed from 0)
     * @param   whichFlag  the flag to test
     * @return  the flag value
     */
    private boolean getFlag(final int row, final int col, final MatrixElementFlag whichFlag) {

        byte mask;

        mask = (byte) (0x01 << whichFlag.flagBit());

        synchronized (this) {
            return (this.flags[row][col] & mask) == mask;
        }
    }

    /**
     * Tests whether a matrix element is marked as completed.
     *
     * @param   row  the row (indexed from 0)
     * @param   col  the column (indexed from 0)
     * @return  <code>true</code> if completed; <code>false</code> if not
     */
    public boolean isCompleted(final int row, final int col) {

        return getFlag(row, col, MatrixElementFlag.COMPLETED);
    }

    /**
     * Tests whether a matrix element is marked as in progress.
     *
     * @param   row  the row (indexed from 0)
     * @param   col  the column (indexed from 0)
     * @return  <code>true</code> if in progress; <code>false</code> if not
     */
    public boolean isInProgress(final int row, final int col) {

        return getFlag(row, col, MatrixElementFlag.IN_PROGRESS);
    }

    /**
     * Tests whether a matrix element is marked as out of range.
     *
     * @param   row  the row (indexed from 0)
     * @param   col  the column (indexed from 0)
     * @return  <code>true</code> if out of range; <code>false</code> if not
     */
    public boolean isOutOfRange(final int row, final int col) {

        return getFlag(row, col, MatrixElementFlag.OUT_OF_RANGE);
    }

    /**
     * Tests whether a matrix element is marked as tagged.
     *
     * @param   row  the row (indexed from 0)
     * @param   col  the column (indexed from 0)
     * @return  <code>true</code> if tagged; <code>false</code> if not
     */
    public boolean isTagged(final int row, final int col) {

        return getFlag(row, col, MatrixElementFlag.TAGGED);
    }

    /**
     * Sets a particular flag value.
     *
     * @param  row        the row (indexed from 0)
     * @param  col        the column (indexed from 0)
     * @param  whichFlag  the flag to set
     * @param  flagValue  <code>true</code> to flag the value; <code>false</code> to clear the flag
     */
    public void setFlag(final int row, final int col, final MatrixElementFlag whichFlag,
        final boolean flagValue) {

        byte mask;

        mask = (byte) (0x01 << whichFlag.flagBit());

        synchronized (this) {

            if (flagValue) {
                this.flags[row][col] |= mask;
            } else {
                this.flags[row][col] &= (~mask);
            }
        }
    }

    /**
     * Tests whether the entire matrix has been computed.
     *
     * @return  <code>true</code> if matrix is complete; <code>false</code> otherwise
     */
    private boolean isMatrixComplete() {

        byte mask;
        boolean complete = true;

        mask = (byte) (0x01 << MatrixElementFlag.COMPLETED.flagBit());

        synchronized (this) {

            for (int row = 0; row < this.flags.length; row++) {

                for (int col = 0; col < this.flags[0].length; col++) {

                    if ((this.flags[row][col] & mask) != mask) {
                        complete = false;

                        break;
                    }
                }
            }
        }

        return complete;
    }

    /**
     * Generates the <code>String</code> representation of the matrix.
     *
     * @return  the string representation
     */
    @Override public String toString() {

        DecimalFormat format;
        StringBuilder str;
        int numRows;
        int numCols;
        String[][] values;
        int[] widths;
        int delta;
        String crlf;

        format = new DecimalFormat("#.######");
        crlf = System.getProperty("line.separator");

        numRows = this.data.length;
        numCols = this.data[0].length;
        widths = new int[numCols];

        synchronized (this) {

            // Generate String representations of each value and get the
            // maximum width of each
            values = new String[numRows][numCols];

            for (int r = 0; r < numRows; r++) {

                for (int c = 0; c < numCols; c++) {
                    values[r][c] = format.format(this.data[r][c]);

                    if (values[r][c].length() > widths[c]) {
                        widths[c] = values[r][c].length();
                    }
                }
            }

            str = new StringBuilder(50);

            for (int r = 0; r < numRows; r++) {
                str.append('[');
                str.append(' ');

                for (int c = 0; c < numCols; c++) {
                    str.append(values[r][c]);
                    delta = widths[c] - values[r][c].length();

                    while (delta > 0) {
                        str.append(' ');
                        delta--;
                    }

                    str.append(' ');
                }

                str.append(']');
                str.append(crlf);
            }
        }

        return str.toString();
    }

    /**
     * Computes the result of multiplying this matrix on the right by a given matrix.
     *
     * @param   mat  the matrix by which to multiply
     * @return  the product [this] x [mat]
     * @throws  EvaluationException  if the matrix dimensions are not compatible
     */
    public DMatrix mul(final DMatrix mat) throws EvaluationException {

        int thisNumRows;
        int thisNumCols;
        int matNumRows;
        int matNumCols;
        DMatrix product;
        double sum;
        String err;

        // Check the sizes this.numCols == mat.numRows
        thisNumRows = numRows();
        thisNumCols = numColumns();
        matNumRows = mat.numRows();
        matNumCols = mat.numColumns();

        if (thisNumCols != matNumRows) {
            err = "Matrix dimensions not compatible for multiplication (" + thisNumRows + "x"
                + thisNumCols + " matrix multiplied by " + matNumRows + "x" + matNumCols
                + " matrix)";
            LOG.warning(err);
            throw new EvaluationException(err);
        }

        product = new DMatrix(thisNumRows, matNumCols);

        synchronized (this) {

            synchronized (mat) {

                for (int r = 0; r < thisNumRows; r++) {

                    for (int c = 0; c < matNumCols; c++) {
                        sum = 0;

                        for (int k = 0; k < thisNumCols; k++) {
                            sum += get(r, k) * mat.get(k, c);
                        }

                        product.set(r, c, sum);
                    }
                }
            }
        }

        return product;
    }

    /**
     * Computes the result of multiplying this matrix on the right by a given vector.
     *
     * @param   vector  the vector by which to multiply
     * @return  the product [this] x [vector]
     * @throws  EvaluationException  if the vector and matrix sizes are incompatible
     */
    public TupleN transform(final TupleN vector) throws EvaluationException {

        int thisNumRows;
        int thisNumCols;
        TupleN product;
        double sum;
        String err;

        // Check the sizes this.numCols == vector.length
        thisNumCols = numColumns();

        if (thisNumCols != vector.getDimension()) {
            err = "Vector must have the same number of rows as the matrix has columns";
            LOG.warning(err);
            throw new EvaluationException(err);
        }

        thisNumRows = numRows();

        product = new TupleN(vector.getDimension(), false);

        synchronized (this) {

            synchronized (vector) {

                for (int r = 0; r < thisNumRows; r++) {
                    sum = 0;

                    for (int c = 0; c < thisNumCols; c++) {
                        sum += vector.get(c) * get(r, c);
                    }

                    product.set(r, sum);
                }
            }
        }

        return product;
    }

    /**
     * Computes the L-U decomposition of the matrix.
     *
     * @throws  EvaluationException         if the matrix is not square
     * @throws  SingularMatrixException     if the matrix is singular
     * @throws  MatrixNotCompleteException  if the matrix has not been completely computed
     */
    public void luDecompose() throws EvaluationException, SingularMatrixException,
        MatrixNotCompleteException {

        int size;
        double big;
        double temp;
        double[] scaling;
        String err;

        if (!isMatrixComplete()) {
            throw new MatrixNotCompleteException(
                "Cannot compute L-U decomposition of an incomplete matrix");
        }

        size = numRows();

        // L-U Decomposition works only for square matrices!
        if (size != numColumns()) {
            err = "Attempt to perform L-U decomposition on a non-square (" + size + "x"
                + this.data[0].length + ") matrix";
            LOG.warning(err);
            throw new EvaluationException(err);
        }

        synchronized (this) {

            // Prepare the data areas and copy the matrix
            this.lud = new double[size][];

            for (int row = 0; row < size; row++) {
                this.lud[row] = this.data[row].clone();
            }

            this.ludIndexes = new int[size];
            this.ludExchangesEven = true;

            // Loop over rows to get implicit scaling information
            scaling = new double[size];

            // Loop over rows to get implicit scaling information
            for (int row = 0; row < size; row++) {
                big = 0.0;

                for (int col = 0; col < size; col++) {
                    temp = Math.abs(this.lud[row][col]);

                    if (temp > big) {
                        big = temp;
                    }
                }

                if (big == 0.0) {
                    err = "Attempt to perform L-U decomposition on a singular matrix";
                    LOG.warning(err);
                    throw new SingularMatrixException(err);
                }

                scaling[row] = 1.0 / big;
            }

            // This is the loop over columns in Crout's method
            croutsMethod(size, scaling);
        }
    }

    /**
     * Clears the LUD data structures.
     */
    private void clearLud() {

        // NOTE: called only from within synchronized block.

        this.lud = null;
        this.ludIndexes = null;
        this.ludExchangesEven = true;
    }

    /**
     * Implementation of Crout's Method of performing L-U decomposition.
     *
     * @param   size     the size of the (square) matrix
     * @param   scaling  the scaling factors for the rows
     * @throws  SingularMatrixException  if the matrix is singular
     */
    private void croutsMethod(final int size, final double[] scaling)
        throws SingularMatrixException {

        // NOTE: called only from within synchronized block.

        int row;
        int col;
        int rmax;
        double big;
        double dum;
        double sum;
        String err;

        for (col = 0; col < size; col++) {

            for (row = 0; row < col; row++) {
                sum = this.lud[row][col];

                for (int k = 0; k < row; k++) {
                    sum -= this.lud[row][k] * this.lud[k][col];
                }

                this.lud[row][col] = sum;
            }

            // Search for the largest pivot element.
            big = 0.0;
            rmax = 0;

            for (row = col; row < size; row++) {
                sum = this.lud[row][col];

                for (int k = 0; k < col; k++) {
                    sum -= this.lud[row][k] * this.lud[k][col];
                }

                this.lud[row][col] = sum;
                dum = scaling[row] * Math.abs(sum);

                if (dum >= big) {
                    big = dum;
                    rmax = row;
                }
            }

            // If we need to interchange rows, do so.
            if (col != rmax) {
                interchange(size, rmax, col);
                scaling[rmax] = scaling[col];
            }

            this.ludIndexes[col] = rmax;

            // If pivot element is zero, matrix is singular.
            if (this.lud[col][col] == 0) {
                clearLud();

                err = "Attempt to perform L-U decomposition on a singular matrix";
                LOG.warning(err);
                throw new SingularMatrixException(err);
            }

            // Divide by the pivot element.
            if (col != (size - 1)) {
                dum = 1.0 / this.lud[col][col];

                for (row = col + 1; row < size; row++) {
                    this.lud[row][col] *= dum;
                }
            }
        }
    }

    /**
     * Performs a row interchange.
     *
     * @param  size    the size of the (square) matrix
     * @param  first   the index of the first row
     * @param  second  the index of the second row
     */
    private void interchange(final int size, final int first, final int second) {

        // NOTE: called only from within synchronized block.

        double dum;

        for (int k = 0; k < size; k++) {
            dum = this.lud[first][k];
            this.lud[first][k] = this.lud[second][k];
            this.lud[second][k] = dum;
        }

        this.ludExchangesEven ^= true;
    }

    /**
     * Solves Ax = b using the L-U decomposition of A computed previously by the <code>
     * luDecompose</code> method.
     *
     * @param   rhs  the right-hand side
     * @return  the x vector that solves the system
     * @throws  EvaluationException         if the matrix is not square
     * @throws  SingularMatrixException     if the matrix is singular
     * @throws  MatrixNotCompleteException  if the matrix has not been completely computed
     */
    public TupleN luSolve(final TupleN rhs) throws EvaluationException, SingularMatrixException,
        MatrixNotCompleteException {

        int size;
        TupleN lhs;
        double sum;
        int indexes;
        int pivot = -1;
        String err;

        synchronized (this) {

            if (this.lud == null) {
                luDecompose();
            }

            size = numColumns();

            if (rhs.getDimension() != size) {
                err = "Right-hand side vector must match number of columns in matrix";
                LOG.warning(err);
                throw new EvaluationException(err);
            }

            lhs = rhs.copy();

            for (int i = 0; i < size; i++) {
                indexes = this.ludIndexes[i];
                sum = lhs.get(indexes);
                lhs.set(indexes, lhs.get(i));

                if (pivot > -1) {

                    for (int j = pivot; j < i; j++) {
                        sum -= this.lud[i][j] * lhs.get(j);
                    }
                } else if (sum != 0) {
                    pivot = i;
                }

                lhs.set(i, sum);
            }

            for (int i = size - 1; i >= 0; i--) {
                sum = lhs.get(i);

                for (int j = i + 1; j < size; j++) {
                    sum -= this.lud[i][j] * lhs.get(j);
                }

                lhs.set(i, sum / this.lud[i][i]);
            }
        }

        return lhs;
    }

    /**
     * Computes the inverse of the matrix.
     *
     * @return  the inverse
     * @throws  EvaluationException         if the matrix is not square
     * @throws  SingularMatrixException     if the matrix is singular
     * @throws  MatrixNotCompleteException  if the matrix has not been completely computed
     */
    public DMatrix inverse() throws EvaluationException, SingularMatrixException,
        MatrixNotCompleteException {

        int size;
        DMatrix inverse;
        TupleN col;
        TupleN vec;

        if (!isMatrixComplete()) {
            throw new MatrixNotCompleteException("Cannot invert an incomplete matrix");
        }

        size = numRows();

        if (size != numColumns()) {
            LOG.log(Level.WARNING, "Attempt to invert  a non-square ({0}x{1}) matrix",
                new Object[] { size, this.data[0].length });
            throw new EvaluationException("Attempt to invert  a non-square (" + size + "x"
                + this.data[0].length + ") matrix");
        }

        synchronized (this) {

            // We use the L-U decomposition to take the inverse.
            if (this.lud == null) {
                luDecompose();
            }

            inverse = new DMatrix(size, size);
            col = new TupleN(size, false);

            for (int c = 0; c < size; c++) {
                col.set(c, 1);
                vec = luSolve(col);
                col.set(c, 0);

                for (int r = 0; r < size; r++) {
                    inverse.set(r, c, vec.get(r));
                }
            }
        }

        return inverse;
    }

    /**
     * Computes the transpose of the matrix.
     *
     * @return  the transpose
     * @throws  MatrixNotCompleteException  if the matrix has not been completely computed
     */
    public DMatrix transpose() throws MatrixNotCompleteException {

        DMatrix transpose;

        if (!isMatrixComplete()) {
            throw new MatrixNotCompleteException("Cannot transpose an incomplete matrix");
        }

        transpose = new DMatrix(this.data[0].length, this.data.length);

        synchronized (this) {

            for (int r = 0; r < this.data.length; r++) {

                for (int c = 0; c < this.data[0].length; c++) {

                    if (isCompleted(r, c)) {
                        transpose.set(c, r, this.data[r][c]);
                    }
                }
            }
        }

        return transpose;
    }

    /**
     * Tests whether all entries in the matrix are within \epsilon of the entries in a test matrix.
     * The test performed on each matrix element is "if (Math.abs(thisMatrix[r][c] -
     * refMatrix[r][c]) < epsilon)". If this statement is true for all entries, this method returns
     * true.
     *
     * @param   ref      the reference matrix to compare to
     * @param   epsilon  the tolerance
     * @return  <code>true</code> if all entries are within tolerance, <code>false</code> otherwise
     * @throws  MatrixNotCompleteException  if the matrix has not been completely computed
     */
    public boolean epsilonEquals(final DMatrix ref, final double epsilon)
        throws MatrixNotCompleteException {

        boolean equals;

        if (!isMatrixComplete() || !ref.isMatrixComplete()) {
            throw new MatrixNotCompleteException("Cannot test equality of an incomplete matrix");
        }

        synchronized (this) {

            synchronized (ref) {

                if ((this.data.length == ref.data.length)
                        && (this.data[0].length == ref.data[0].length)) {
                    equals = true;

outer:
                    for (int r = 0; r < this.data.length; r++) {

                        for (int c = 0; c < this.data[0].length; c++) {

                            if (Math.abs(this.data[r][c] - ref.data[r][c]) > epsilon) {
                                equals = false;

                                break outer;
                            }
                        }
                    }
                } else {
                    equals = false;
                }
            }
        }

        return equals;
    }

    /**
     * Tests whether all entries in the matrix are equal to the entries in a test matrix.
     *
     * @param   ref  the matrix to compare to
     * @return  <code>true</code> if matrix sizes and all entries are equal, <code>false</code>
     *          otherwise
     * @throws  MatrixNotCompleteException  if the matrix has not been completely computed
     */
    public boolean equalsMatrix(final DMatrix ref) throws MatrixNotCompleteException {

        boolean equals;

        if (!isMatrixComplete() || !ref.isMatrixComplete()) {
            throw new MatrixNotCompleteException("Cannot test equality of an incomplete matrix");
        }

        synchronized (this) {

            synchronized (ref) {

                if ((this.data.length == ref.data.length)
                        && (this.data[0].length == ref.data[0].length)) {
                    equals = true;

outer:
                    for (int r = 0; r < this.data.length; r++) {

                        for (int c = 0; c < this.data[0].length; c++) {

                            if (this.data[r][c] != ref.data[r][c]) {
                                equals = false;

                                break outer;
                            }
                        }
                    }
                } else {
                    equals = false;
                }
            }
        }

        return equals;
    }

    /**
     * Tests whether an object is a <code>DMatrix</code> and is equal to this matrix.
     *
     * @param   obj  the object to compare to
     * @return  <code>true</code> if the object is a <code>DMatrix</code> and matrix sizes and all
     *          entries are equal, <code>false</code> otherwise
     */
    @Override public boolean equals(final Object obj) {

        boolean equals;

        if (obj instanceof DMatrix) {

            try {
                equals = equalsMatrix((DMatrix) obj);
            } catch (MatrixNotCompleteException e) {
                equals = false;
            }
        } else {
            equals = false;
        }

        return equals;
    }

    /**
     * Computes the hash code of the object.
     *
     * @return  the hash code
     */
    @Override public int hashCode() {

        synchronized (this) {
            return this.data.hashCode();
        }
    }

    /**
     * Saves the matrix, including the completion status of all values, to a file. The in-progress
     * status of values, and any LU decomposition data is not saved with the matrix. The matrix is
     * saved in a format that will allow GnuPlot to display it as a two-dimensional data set using
     * the 'SPLOT' command. Completion status of values is stored in the file as GnuPlot comments.
     *
     * @param   file  the file to which to write
     * @throws  IOException  if there is an error writing to the file
     */
    public void save(final File file) throws IOException {

        File tmp;
        FileOutputStream fos;
        PrintStream out;
        int row;
        int col;

        // If there is an existing file, we write this out to a new temporary
        // file, then replace the existing file only on successful write.
        if (file.exists()) {
            tmp = File.createTempFile("DMatrix", "dat", file.getParentFile());
        } else {
            tmp = file;
        }

        fos = new FileOutputStream(tmp);
        out = new PrintStream(fos);

        out.println("# DMatrix v1.0 storage file - do not alter comments");

        out.print("# numRows=");
        out.println(this.data.length);

        out.print("# numCols=");
        out.println(this.data[0].length);

        // Output the completed status of the matrix
        synchronized (this) {

            for (row = 0; row < this.data.length; row++) {
                out.print("# ");
                out.print(row);
                out.print(":");

                for (col = 0; col < this.data[0].length; col++) {
                    out.print(HEX[this.flags[row][col] & 0x0F]);
                }

                out.println("");
            }

            // Now output the rows
            for (row = 0; row < this.data.length; row++) {
                out.println("");

                for (col = 0; col < this.data[0].length; col++) {
                    out.print(row);
                    out.print(' ');
                    out.print(col);
                    out.print(' ');
                    out.println(this.data[row][col]);
                }
            }
        }

        out.println("# END");

        out.close();
        fos.close();

        if (!tmp.equals(file)) {

            if (file.exists()) {
                LOG.log(Level.FINE, "Deleting {0}", file.getAbsolutePath());

                if (!file.delete()) {
                    throw new IOException("Unable to overwrite '" + file.getAbsolutePath() + "'");
                }
            }

            LOG.log(Level.FINE, "Renaming to {0}", file.getAbsolutePath());

            if (!tmp.renameTo(file)) {
                throw new IOException("Unable to rename temp file to '" + file.getAbsolutePath()
                    + "'");
            }
        }
    }

    /**
     * Loads the matrix, including the completion status of all values, from a file.
     *
     * @param   file  the file from which to load
     * @return  the loaded matrix
     * @throws  IOException  if there is an error reading from the file
     */
    public static DMatrix load(final File file) throws IOException {

        FileReader reader;
        BufferedReader buf;
        String line;
        int numRows;
        int numCols;
        DMatrix matrix;
        int row;
        int col;
        String expect;
        char[] chars;
        String[] fields;

        reader = new FileReader(file);
        buf = new BufferedReader(reader);

        // Read the first-line comment
        line = buf.readLine();

        if (!line.startsWith("# DMatrix v1.0")) {
            throw new IOException(
                "Invlaid matrix file - first line does not begin with '# DMatrix v1.0'.");
        }

        // Read the number of rows and columns
        line = buf.readLine();

        if (!line.startsWith("# numRows=")) {
            throw new IOException(
                "Invlaid matrix file - second line does not contain '# numRows=...'");
        }

        numRows = Integer.parseInt(line.substring(10));

        line = buf.readLine();

        if (!line.startsWith("# numCols=")) {
            throw new IOException(
                "Invlaid matrix file - third line does not contain '# numCols=...'");
        }

        numCols = Integer.parseInt(line.substring(10));

        matrix = new DMatrix(numRows, numCols);

        // Read in the completed status of the matrix
        for (row = 0; row < matrix.data.length; row++) {
            line = buf.readLine();
            expect = "# " + row + ":";

            if (!line.startsWith(expect)) {
                throw new IOException("Invlaid matrix file - line " + (3 + row)
                    + " does not contain '" + expect + "'");
            }

            chars = line.substring(expect.length()).toCharArray();

            if (chars.length != matrix.data[0].length) {
                throw new IOException("Invlaid matrix file - line " + (3 + row) + " does not have "
                    + matrix.data[0].length + " flag values");
            }

            for (col = 0; col < matrix.data[0].length; col++) {
                matrix.flags[row][col] = -1;

                for (int i = 0; i < HEX.length; i++) {

                    if (HEX[i] == chars[col]) {
                        matrix.flags[row][col] = (byte) i;

                        break;
                    }
                }

                if (matrix.flags[row][col] == -1) {
                    throw new IOException("Invlaid matrix file - line " + (3 + row)
                        + " has invalid flag value: '" + chars[col] + "'");
                }
            }
        }

        // Read in the matrix data values.
        for (row = 0; row < matrix.data.length; row++) {
            line = buf.readLine();

            if (line.length() > 0) {
                throw new IOException(
                    "Invlaid matrix file - expected blank line before matrix row " + row);
            }

            for (col = 0; col < matrix.data[0].length; col++) {
                line = buf.readLine();
                fields = line.split(" ");

                if (fields.length != 3) {
                    throw new IOException("Invlaid matrix file - matrix row " + row + " column "
                        + col + " does not have 3 fields");
                }

                if (Integer.parseInt(fields[0]) != row) {
                    throw new IOException("Invlaid matrix file - matrix row " + row + " column "
                        + col + " has bad row number");
                }

                if (Integer.parseInt(fields[1]) != col) {
                    throw new IOException("Invlaid matrix file - matrix row " + row + " column "
                        + col + " has bad column number");
                }

                matrix.data[row][col] = Double.parseDouble(fields[2]);
            }
        }

        // Check for proper termination.
        line = buf.readLine();

        if (!line.startsWith("# END")) {
            throw new IOException(
                "Invlaid matrix file - line after matrix does not begin with '# END'.");
        }

        buf.close();
        reader.close();

        // Finally, compute the first uncompleted cell
        matrix.firstUncRow = 0;
        matrix.firstUncCol = 0;
        matrix.advanceFirstUncompleted();

        return matrix;
    }

    /**
     * Generates a list of ranges to use to color map the matrix. Each range will have roughly the
     * same number of matrix values in it, and will be in ascending order by value. If the matrix
     * has many cells with the same data value, it may not be possible to distribute the color
     * ranges evenly. If the matrix is not complete, only the completed data is used. If there are
     * fewer values in the matrix than the number of ranges requested, all ranges after a certain
     * point will have the same value.
     *
     * @param   numRanges  the number of ranges to generate
     * @return  the color range table - an array of doubles whose values are the boundaries of each
     *          color range; the last entry in the array is at least as large as the largest data
     *          value in the matrix
     */
    public double[] colorRanges(final int numRanges) {

        double[] ranges;
        int total;
        int count;
        double[] values;
        double which;

        ranges = new double[numRanges];

        // We scan the matrix copying all distinct data values into a list.
        synchronized (this) {
            total = 0;

            for (int r = 0; r < this.data.length; r++) {

                for (int c = 0; c < this.data[0].length; c++) {

                    if (isCompleted(r, c)) {
                        total++;
                    }
                }
            }

            values = new double[total];
            count = 0;

            for (int r = 0; r < this.data.length; r++) {

                for (int c = 0; c < this.data[0].length; c++) {

                    if (isCompleted(r, c)) {
                        values[count] = this.data[r][c];
                        count++;
                    }
                }
            }

            // Now we sort the data values into ascending order
            Arrays.sort(values);

            // If there are more ranges (colors) than the length of the values
            // array, we still want to use all the colors, so we just divide
            // the range linearly. Otherwise, we distribute ranges evenly among
            // the values.
            if (numRanges > values.length) {

                for (int i = 0; i < values.length; i++) {
                    which = i * numRanges / values.length;
                    Arrays.fill(ranges, (int) which, numRanges - 1, values[i]);
                }
            } else {

                for (double i = 0; i < numRanges; i++) {
                    which = i * values.length / numRanges;

                    ranges[(int) i] = (which < (values.length - 2)) ? values[(int) which + 1]
                                                                    : values[values.length - 1];
                }
            }
        }

        return ranges;
    }

    /**
     * Main method for testing.
     *
     * @param  args  command-line arguments
     */
    public static void main(final String... args) {

        DMatrix mat;
        DMatrix mat2;
        int[] rowCol;

        mat = new DMatrix(4, 4);
        mat.set(0, 0, 1);
        mat.set(1, 1, Math.PI);
        mat.set(2, 2, Math.E);
        mat.set(3, 3, Math.sqrt(2));

        try {
            mat.save(new File("/imp/testMatrix.dat"));

            LOG.fine(mat.toString());

            mat2 = load(new File("/imp/testMatrix.dat"));

            LOG.fine(mat2.toString());

            rowCol = mat2.getUncompletedElement();
            LOG.log(Level.FINE, "A: {0}, {1}", new Object[] { rowCol[0], rowCol[1] });

            rowCol = mat2.getUncompletedElement();
            LOG.log(Level.FINE, "B: {0}, {1}", new Object[] { rowCol[0], rowCol[1] });

            rowCol = mat2.getUncompletedElement();
            LOG.log(Level.FINE, "C: {0}, {1}", new Object[] { rowCol[0], rowCol[1] });

            rowCol = mat2.getUncompletedElement();
            LOG.log(Level.FINE, "D: {0}, {1}", new Object[] { rowCol[0], rowCol[1] });

            rowCol = mat2.getUncompletedElement();
            LOG.log(Level.FINE, "E: {0}, {1}", new Object[] { rowCol[0], rowCol[1] });

            rowCol = mat2.getUncompletedElement();
            LOG.log(Level.FINE, "F: {0}, {1}", new Object[] { rowCol[0], rowCol[1] });

            rowCol = mat2.getUncompletedElement();
            LOG.log(Level.FINE, "G: {0}, {1}", new Object[] { rowCol[0], rowCol[1] });

            rowCol = mat2.getUncompletedElement();
            LOG.log(Level.FINE, "H: {0}, {1}", new Object[] { rowCol[0], rowCol[1] });

            rowCol = mat2.getUncompletedElement();
            LOG.log(Level.FINE, "I: {0}, {1}", new Object[] { rowCol[0], rowCol[1] });

            rowCol = mat2.getUncompletedElement();
            LOG.log(Level.FINE, "J: {0}, {1}", new Object[] { rowCol[0], rowCol[1] });

            rowCol = mat2.getUncompletedElement();
            LOG.log(Level.FINE, "K: {0}, {1}", new Object[] { rowCol[0], rowCol[1] });

            rowCol = mat2.getUncompletedElement();
            LOG.log(Level.FINE, "L: {0}, {1}", new Object[] { rowCol[0], rowCol[1] });

            rowCol = mat2.getUncompletedElement();
            LOG.fine(Boolean.toString(rowCol == null));
        } catch (IOException e) {
            LOG.throwing("DMatrix", "main", e);
        }
    }
}