/*
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
package org.apache.sis.referencing.operation.matrix;

import java.util.Arrays;
import java.util.Objects;
import org.opengis.geometry.Envelope;
import org.opengis.geometry.DirectPosition;
import org.opengis.referencing.cs.AxisDirection;
import org.opengis.referencing.cs.CoordinateSystem;             // For javadoc
import org.opengis.referencing.operation.Matrix;
import org.opengis.referencing.operation.MathTransform;
import org.opengis.geometry.MismatchedDimensionException;
import org.apache.sis.util.Static;
import org.apache.sis.util.CharSequences;
import org.apache.sis.util.ComparisonMode;
import org.apache.sis.util.ArgumentChecks;
import org.apache.sis.util.StringBuilders;
import org.apache.sis.util.resources.Errors;
import org.apache.sis.math.DecimalFunctions;
import org.apache.sis.internal.util.Numerics;
import org.apache.sis.internal.util.DoubleDouble;
import org.apache.sis.internal.referencing.AxisDirections;
import org.apache.sis.internal.referencing.Resources;
import org.apache.sis.internal.referencing.ExtendedPrecisionMatrix;


/**
 * {@link Matrix} factory methods and utilities.
 * This class provides the following methods:
 *
 * <ul>
 *   <li>Creating new matrices of arbitrary size:
 *       {@link #createIdentity createIdentity},
 *       {@link #createDiagonal createDiagonal},
 *       {@link #createZero     createZero},
 *       {@link #create         create},
 *       {@link #copy           copy}.
 *   </li>
 *   <li>Creating new matrices for coordinate operation steps:
 *       {@link #createTransform(Envelope, AxisDirection[], Envelope, AxisDirection[]) createTransform},
 *       {@link #createDimensionSelect createDimensionSelect},
 *       {@link #createPassThrough     createPassThrough}.
 *   </li>
 *   <li>Information:
 *       {@link #isAffine   isAffine},
 *       {@link #isIdentity isIdentity},
 *       {@link #equals(Matrix, Matrix, double, boolean) equals},
 *       {@link #toString(Matrix) toString}.
 *   </li>
 *   <li>Miscellaneous:
 *       {@link #multiply     multiply},
 *       {@link #inverse      inverse},
 *       {@link #unmodifiable unmodifiable},
 *   </li>
 * </ul>
 *
 * @author  Martin Desruisseaux (IRD, Geomatys)
 * @version 1.0
 *
 * @see org.apache.sis.parameter.TensorParameters
 *
 * @since 0.4
 * @module
 */
public final class Matrices extends Static {
    /**
     * Number of spaces to put between columns formatted by {@link #toString(Matrix)}.
     */
    private static final int SPACING = 2;

    /**
     * Do not allows instantiation of this class.
     */
    private Matrices() {
    }

    /**
     * Creates a square identity matrix of size {@code size} × {@code size}.
     * Elements on the diagonal (<var>j</var> == <var>i</var>) are set to 1.
     *
     * <div class="note"><b>Implementation note:</b>
     * For sizes between {@value org.apache.sis.referencing.operation.matrix.Matrix1#SIZE} and
     * {@value org.apache.sis.referencing.operation.matrix.Matrix4#SIZE} inclusive, the matrix
     * is guaranteed to be an instance of one of {@link Matrix1} … {@link Matrix4} subtypes.</div>
     *
     * @param  size  numbers of row and columns. For an affine transform matrix, this is the number of
     *         {@linkplain MathTransform#getSourceDimensions() source} and
     *         {@linkplain MathTransform#getTargetDimensions() target} dimensions + 1.
     * @return an identity matrix of the given size.
     */
    public static MatrixSIS createIdentity(final int size) {
        switch (size) {
            case 1:  return new Matrix1();
            case 2:  return new Matrix2();
            case 3:  return new Matrix3();
            case 4:  return new Matrix4();
            default: return new GeneralMatrix(size, size, true, 1);
        }
    }

    /**
     * Creates a matrix of size {@code numRow} × {@code numCol}.
     * Elements on the diagonal (<var>j</var> == <var>i</var>) are set to 1.
     * The result is an identity matrix if {@code numRow} = {@code numCol}.
     *
     * <div class="note"><b>Implementation note:</b>
     * For {@code numRow} == {@code numCol} with a value between
     * {@value org.apache.sis.referencing.operation.matrix.Matrix1#SIZE} and
     * {@value org.apache.sis.referencing.operation.matrix.Matrix4#SIZE} inclusive, the matrix
     * is guaranteed to be an instance of one of {@link Matrix1} … {@link Matrix4} subtypes.</div>
     *
     * @param  numRow  for a math transform, this is the number of {@linkplain MathTransform#getTargetDimensions() target dimensions} + 1.
     * @param  numCol  for a math transform, this is the number of {@linkplain MathTransform#getSourceDimensions() source dimensions} + 1.
     * @return an identity matrix of the given size.
     */
    public static MatrixSIS createDiagonal(final int numRow, final int numCol) {
        if (numRow == numCol) {
            return createIdentity(numRow);
        } else {
            return new NonSquareMatrix(numRow, numCol, true, 1);
        }
    }

    /**
     * Creates a matrix of size {@code numRow} × {@code numCol} filled with zero values.
     * This constructor is convenient when the caller wants to initialize the matrix elements himself.
     *
     * <div class="note"><b>Implementation note:</b>
     * For {@code numRow} == {@code numCol} with a value between
     * {@value org.apache.sis.referencing.operation.matrix.Matrix1#SIZE} and
     * {@value org.apache.sis.referencing.operation.matrix.Matrix4#SIZE} inclusive, the matrix
     * is guaranteed to be an instance of one of {@link Matrix1} … {@link Matrix4} subtypes.</div>
     *
     * @param  numRow  for a math transform, this is the number of {@linkplain MathTransform#getTargetDimensions() target dimensions} + 1.
     * @param  numCol  for a math transform, this is the number of {@linkplain MathTransform#getSourceDimensions() source dimensions} + 1.
     * @return a matrix of the given size with only zero values.
     */
    public static MatrixSIS createZero(final int numRow, final int numCol) {
        if (numRow == numCol) switch (numRow) {
            case 1:  return new Matrix1(false);
            case 2:  return new Matrix2(false);
            case 3:  return new Matrix3(false);
            case 4:  return new Matrix4(false);
            default: return new GeneralMatrix(numRow, numCol, false, 1);
        }
        return new NonSquareMatrix(numRow, numCol, false, 1);
    }

    /**
     * Creates a matrix filled with zero values, using double-double precision if {@code precision} is {@code true}.
     */
    static MatrixSIS createZero(final int numRow, final int numCol, final boolean precision) {
        return precision ? GeneralMatrix.createExtendedPrecision(numRow, numCol, false) : createZero(numRow, numCol);
    }

    /**
     * Creates a matrix of size {@code numRow} × {@code numCol} initialized to the given elements.
     * The elements array size must be equals to {@code numRow*numCol}. Column indices vary fastest.
     *
     * <div class="note"><b>Implementation note:</b>
     * For {@code numRow} == {@code numCol} with a value between
     * {@value org.apache.sis.referencing.operation.matrix.Matrix1#SIZE} and
     * {@value org.apache.sis.referencing.operation.matrix.Matrix4#SIZE} inclusive, the matrix
     * is guaranteed to be an instance of one of {@link Matrix1} … {@link Matrix4} subtypes.</div>
     *
     * @param  numRow    number of rows.
     * @param  numCol    number of columns.
     * @param  elements  the matrix elements in a row-major array. Column indices vary fastest.
     * @return a matrix initialized to the given elements.
     *
     * @see MatrixSIS#setElements(double[])
     */
    public static MatrixSIS create(final int numRow, final int numCol, final double[] elements) {
        if (numRow == numCol) switch (numRow) {
            case 1:  return new Matrix1(elements);
            case 2:  return new Matrix2(elements);
            case 3:  return new Matrix3(elements);
            case 4:  return new Matrix4(elements);
            default: return new GeneralMatrix(numRow, numCol, elements);
        }
        return new NonSquareMatrix(numRow, numCol, elements);
    }

    /**
     * Creates a matrix of size {@code numRow} × {@code numCol} initialized to the given numbers.
     * The elements array size must be equals to {@code numRow*numCol}. Column indices vary fastest.
     *
     * @param  numRow    number of rows.
     * @param  numCol    number of columns.
     * @param  elements  the matrix elements in a row-major array. Column indices vary fastest.
     * @return a matrix initialized to the given elements.
     */
    public static MatrixSIS create(final int numRow, final int numCol, final Number[] elements) {
        ArgumentChecks.ensureNonNull("elements", elements);
        /*
         * Below is an intantionally undocumented feature. We use those sentinel values as a way to create
         * matrices with extended precision without exposing our double-double arithmetic in public API.
         */
        if (elements == ExtendedPrecisionMatrix.IDENTITY || elements == ExtendedPrecisionMatrix.ZERO) {
            return GeneralMatrix.createExtendedPrecision(numRow, numCol, elements == ExtendedPrecisionMatrix.IDENTITY);
        }
        final GeneralMatrix matrix = GeneralMatrix.createExtendedPrecision(numRow, numCol, false);
        if (matrix.setElements(elements)) {
            /*
             * At least one org.apache.sis.internal.util.DoubleDouble instance has been found,
             * in which case the matrix uses double-double arithmetic.  This case is the main
             * purpose of this method.
             */
            return matrix;
        }
        return create(numRow, numCol, matrix.getElements());
    }

    /**
     * Implementation of {@code createTransform(…)} public methods expecting envelopes and/or axis directions.
     * Argument validity shall be verified by the caller.
     *
     * @param useEnvelopes {@code true} if source and destination envelopes shall be taken in account.
     *        If {@code false}, then source and destination envelopes will be ignored and can be null.
     */
    @SuppressWarnings("null")
    private static MatrixSIS createTransform(final Envelope srcEnvelope, final AxisDirection[] srcAxes,
                                             final Envelope dstEnvelope, final AxisDirection[] dstAxes,
                                             final boolean useEnvelopes)
    {
        final DirectPosition dstCorner, srcCorner, srcOppositeCorner;
        if (useEnvelopes) {
            dstCorner         = dstEnvelope.getLowerCorner();
            srcCorner         = srcEnvelope.getLowerCorner();
            srcOppositeCorner = srcEnvelope.getUpperCorner();
        } else {
            dstCorner = srcCorner = srcOppositeCorner = null;
        }
        /*
         * Unconditionally create extended precision matrix even if standard precision would be
         * enough because callers in other package may perform additional arithmetic operations
         * on it (for example org.apache.sis.referencing.cs.CoordinateSystems.swapAndScaleAxes).
         */
        final MatrixSIS matrix = new GeneralMatrix(dstAxes.length+1, srcAxes.length+1, false, 2);
        /*
         * Maps source axes to destination axes. If no axis is moved (for example if the user
         * want to transform (NORTH,EAST) to (SOUTH,EAST)), then source and destination index
         * will be equal. If some axes are moved (for example if the user want to transform
         * (NORTH,EAST) to (EAST,NORTH)), then coordinates at index {@code srcIndex} will have
         * to be moved at index {@code dstIndex}.
         */
        for (int dstIndex = 0; dstIndex < dstAxes.length; dstIndex++) {
            boolean hasFound = false;
            final AxisDirection dstDir = dstAxes[dstIndex];
            final AxisDirection search = AxisDirections.absolute(dstDir);
            for (int srcIndex = 0; srcIndex < srcAxes.length; srcIndex++) {
                final AxisDirection srcDir = srcAxes[srcIndex];
                if (search.equals(AxisDirections.absolute(srcDir))) {
                    if (hasFound) {
                        throw new IllegalArgumentException(Resources.format(
                                Resources.Keys.ColinearAxisDirections_2, srcDir, dstDir));
                    }
                    hasFound = true;
                    /*
                     * Set the matrix elements. Some matrix elements will never be set.
                     * They will be left to zero, which is their desired value.
                     */
                    final boolean same = srcDir.equals(dstDir);
                    if (useEnvelopes) {
                        /*
                         * See the comment in transform(Envelope, Envelope) for an explanation about why
                         * we use the lower/upper corners instead than getMinimum()/getMaximum() methods.
                         */
                        final DoubleDouble scale = new DoubleDouble(same ? +1d : -1d);
                        scale.multiplyGuessError(dstEnvelope.getSpan(dstIndex));
                        scale.divideGuessError  (srcEnvelope.getSpan(srcIndex));

                        final DoubleDouble translate = new DoubleDouble(scale);
                        translate.multiplyGuessError((same ? srcCorner : srcOppositeCorner).getOrdinate(srcIndex));
                        translate.negate();
                        translate.addGuessError(dstCorner.getOrdinate(dstIndex));

                        matrix.setNumber(dstIndex, srcIndex,       scale);
                        matrix.setNumber(dstIndex, srcAxes.length, translate);
                    } else {
                        matrix.setElement(dstIndex, srcIndex, same ? +1 : -1);
                    }
                }
            }
            if (!hasFound) {
                throw new IllegalArgumentException(Resources.format(
                        Resources.Keys.CanNotMapAxisToDirection_1, dstAxes[dstIndex]));
            }
        }
        matrix.setElement(dstAxes.length, srcAxes.length, 1);
        return matrix;
    }

    /**
     * Creates a transform matrix mapping the given source envelope to the given destination envelope.
     * The given envelopes can have any dimensions, which are handled as below:
     *
     * <ul>
     *   <li>If the two envelopes have the same {@linkplain Envelope#getDimension() dimension},
     *       then the transform is {@linkplain #isAffine(Matrix) affine}.</li>
     *   <li>If the destination envelope has less dimensions than the source envelope,
     *       then trailing dimensions are silently dropped.</li>
     *   <li>If the target envelope has more dimensions than the source envelope,
     *       then the transform will append trailing coordinates with the 0 value.</li>
     * </ul>
     *
     * This method ignores the {@linkplain Envelope#getCoordinateReferenceSystem() envelope CRS}, which may be null.
     * Actually this method is used more often for {@linkplain org.opengis.coverage.grid.GridEnvelope grid envelopes}
     * (which have no CRS) than geodetic envelopes.
     *
     * <h4>Spanning the anti-meridian of a Geographic CRS</h4>
     * If the given envelopes cross the date line, then this method requires their {@code getSpan(int)} method
     * to behave as documented in the {@link org.apache.sis.geometry.AbstractEnvelope#getSpan(int)} javadoc.
     * Furthermore the matrix created by this method will produce expected results only for source or destination
     * points before the date line, since the wrap around operation can not be represented by an affine transform.
     *
     * <h4>Example</h4>
     * Given a source envelope of size 100 × 200 (the units do not matter for this method) and a destination
     * envelope of size 300 × 500, and given {@linkplain Envelope#getLowerCorner() lower corner} translation
     * from (-20, -40) to (-10, -25), then the following method call:
     *
     * {@preformat java
     *   matrix = Matrices.createTransform(
     *           new Envelope2D(null, -20, -40, 100, 200),
     *           new Envelope2D(null, -10, -25, 300, 500));
     * }
     *
     * will return the following square matrix. The transform of the lower corner is given as an example:
     *
     * {@preformat math
     *   ┌     ┐   ┌              ┐   ┌     ┐
     *   │ -10 │   │ 3.0  0    50 │   │ -20 │       // 3.0 is the scale factor from width of 100 to 300
     *   │ -25 │ = │ 0    2.5  75 │ × │ -40 │       // 2.5 is the scale factor from height of 200 to 500
     *   │   1 │   │ 0    0     1 │   │   1 │
     *   └     ┘   └              ┘   └     ┘
     * }
     *
     * @param  srcEnvelope  the source envelope.
     * @param  dstEnvelope  the destination envelope.
     * @return the transform from the given source envelope to target envelope.
     *
     * @see #createTransform(AxisDirection[], AxisDirection[])
     * @see #createTransform(Envelope, AxisDirection[], Envelope, AxisDirection[])
     * @see org.apache.sis.referencing.cs.CoordinateSystems#swapAndScaleAxes(CoordinateSystem, CoordinateSystem)
     */
    public static MatrixSIS createTransform(final Envelope srcEnvelope, final Envelope dstEnvelope) {
        ArgumentChecks.ensureNonNull("srcEnvelope", srcEnvelope);
        ArgumentChecks.ensureNonNull("dstEnvelope", dstEnvelope);
        /*
         * Following code is a simplified version of above createTransform(Envelope, AxisDirection[], ...) method.
         * We need to make sure that those two methods are consistent and compute the matrix values in the same way.
         */
        final int srcDim = srcEnvelope.getDimension();
        final int dstDim = dstEnvelope.getDimension();
        final DirectPosition srcCorner = srcEnvelope.getLowerCorner();
        final DirectPosition dstCorner = dstEnvelope.getLowerCorner();
        final MatrixSIS matrix = createZero(dstDim+1, srcDim+1);
        for (int i = Math.min(srcDim, dstDim); --i >= 0;) {
            /*
             * Note on envelope spanning over the anti-meridian: the GeoAPI javadoc does not mandate the
             * precise behavior of getSpan(int) in such situation.  In the particular case of Apache SIS
             * implementations, the envelope will compute the span correctly (taking in account the wrap
             * around behavior). For non-SIS implementations, we can not know.
             *
             * For the translation term, we really need the lower corner, NOT envelope.getMinimum(i),
             * because we need the starting point, which is not the minimal value when spanning over
             * the anti-meridian.
             */
            final double scale     = dstEnvelope.getSpan(i)   / srcEnvelope.getSpan(i);
            final double translate = dstCorner.getOrdinate(i) - srcCorner.getOrdinate(i)*scale;
            matrix.setElement(i, i,      scale);
            matrix.setElement(i, srcDim, translate);
        }
        matrix.setElement(dstDim, srcDim, 1);
        return matrix;
    }

    /**
     * Creates a transform matrix changing axis order and/or direction. For example the transform may convert
     * (<i>northing</i>, <i>westing</i>) coordinates into (<i>easting</i>, <i>northing</i>) coordinates.
     * This method tries to associate each {@code dstAxes} direction to either an equals {@code srcAxis}
     * direction, or to an opposite {@code srcAxis} direction.
     *
     * <ul>
     *   <li>If some {@code srcAxes} directions can not be mapped to {@code dstAxes} directions, then the transform
     *       will silently drops the coordinates associated to those extra source axis directions.</li>
     *   <li>If some {@code dstAxes} directions can not be mapped to {@code srcAxes} directions,
     *       then an exception will be thrown.</li>
     * </ul>
     *
     * <div class="note"><b>Example:</b>
     * it is legal to transform from (<i>easting</i>, <i>northing</i>, <i>up</i>) to
     * (<i>easting</i>, <i>northing</i>) — this is the first above case — but illegal
     * to transform (<i>easting</i>, <i>northing</i>) to (<i>easting</i>, <i>up</i>).</div>
     *
     * <h4>Example</h4>
     * The following method call:
     *
     * {@preformat java
     *   matrix = Matrices.createTransform(
     *           new AxisDirection[] {AxisDirection.NORTH, AxisDirection.WEST},
     *           new AxisDirection[] {AxisDirection.EAST, AxisDirection.NORTH});
     * }
     *
     * will return the following square matrix, which can be used in coordinate conversions as below:
     *
     * {@preformat math
     *   ┌    ┐   ┌         ┐   ┌    ┐
     *   │ +x │   │ 0 -1  0 │   │  y │
     *   │  y │ = │ 1  0  0 │ × │ -x │
     *   │  1 │   │ 0  0  1 │   │  1 │
     *   └    ┘   └         ┘   └    ┘
     * }
     *
     * @param  srcAxes  the ordered sequence of axis directions for source coordinate system.
     * @param  dstAxes  the ordered sequence of axis directions for destination coordinate system.
     * @return the transform from the given source axis directions to the given target axis directions.
     * @throws IllegalArgumentException if {@code dstAxes} contains at least one axis not found in {@code srcAxes},
     *         or if some colinear axes were found.
     *
     * @see #createTransform(Envelope, Envelope)
     * @see #createTransform(Envelope, AxisDirection[], Envelope, AxisDirection[])
     * @see org.apache.sis.referencing.cs.CoordinateSystems#swapAndScaleAxes(CoordinateSystem, CoordinateSystem)
     */
    public static MatrixSIS createTransform(final AxisDirection[] srcAxes, final AxisDirection[] dstAxes) {
        ArgumentChecks.ensureNonNull("srcAxes", srcAxes);
        ArgumentChecks.ensureNonNull("dstAxes", dstAxes);
        if (Arrays.equals(srcAxes, dstAxes)) {
            /*
             * createTransform(…) may fail if the arrays contain two axes with the same direction, for example
             * AxisDirection.OTHER. This check prevents that failure for the common case of an identity transform.
             */
            return Matrices.createIdentity(srcAxes.length + 1);
        }
        return createTransform(null, srcAxes, null, dstAxes, false);
    }

    /**
     * Creates a transform matrix mapping the given source envelope to the given destination envelope,
     * combined with changes in axis order and/or direction.
     * Invoking this method is equivalent to concatenating the following matrix transforms:
     *
     * <ul>
     *   <li><code>{@linkplain #createTransform(Envelope, Envelope) createTransform}(srcEnvelope, dstEnvelope)</code></li>
     *   <li><code>{@linkplain #createTransform(AxisDirection[], AxisDirection[]) createTransform}(srcAxes, dstAxes)</code></li>
     * </ul>
     *
     * This method ignores the {@linkplain Envelope#getCoordinateReferenceSystem() envelope CRS}, which may be null.
     * Actually this method is used more often for {@linkplain org.opengis.coverage.grid.GridEnvelope grid envelopes}
     * (which have no CRS) than geodetic envelopes.
     *
     * <h4>Spanning the anti-meridian of a Geographic CRS</h4>
     * If the given envelopes cross the date line, then this method requires their {@code getSpan(int)} method
     * to behave as documented in the {@link org.apache.sis.geometry.AbstractEnvelope#getSpan(int)} javadoc.
     * Furthermore the matrix created by this method will produce expected results only for source or destination
     * points on one side of the date line (depending on whether axis direction is reversed), since the wrap around
     * operation can not be represented by an affine transform.
     *
     * <div class="note"><b>Example:</b>
     * combining the examples documented in the above {@code createTransform(…)} methods, the following method call:
     *
     * {@preformat java
     *   matrix = Matrices.createTransform(
     *           new Envelope2D(null, -40, +20, 200, 100), new AxisDirection[] {AxisDirection.NORTH, AxisDirection.WEST},
     *           new Envelope2D(null, -10, -25, 300, 500), new AxisDirection[] {AxisDirection.EAST, AxisDirection.NORTH});
     * }
     *
     * will return the following square matrix. The transform of a corner is given as an example.
     * Note that the input coordinate values are swapped because of the (<i>North</i>, <i>West</i>) axis directions,
     * and the lower-left corner of the destination envelope is the lower-<em>right</em> corner of the source envelope
     * because of the opposite axis direction.
     *
     * {@preformat math
     *   ┌     ┐   ┌               ┐   ┌     ┐
     *   │ -10 │   │ 0   -3.0  350 │   │ -40 │
     *   │ -25 │ = │ 2.5  0     75 │ × │ 120 │       // 120 is the westernmost source coordinate: (x=20) + (width=100)
     *   │   1 │   │ 0    0      1 │   │   1 │
     *   └     ┘   └               ┘   └     ┘
     * }
     * </div>
     *
     * @param  srcEnvelope  the source envelope.
     * @param  srcAxes      the ordered sequence of axis directions for source coordinate system.
     * @param  dstEnvelope  the destination envelope.
     * @param  dstAxes      the ordered sequence of axis directions for destination coordinate system.
     * @return the transform from the given source envelope and axis directions
     *         to the given envelope and target axis directions.
     * @throws MismatchedDimensionException if an envelope {@linkplain Envelope#getDimension() dimension} does not
     *         match the length of the axis directions sequence.
     * @throws IllegalArgumentException if {@code dstAxes} contains at least one axis not found in {@code srcAxes},
     *         or if some colinear axes were found.
     *
     * @see #createTransform(Envelope, Envelope)
     * @see #createTransform(AxisDirection[], AxisDirection[])
     * @see org.apache.sis.referencing.cs.CoordinateSystems#swapAndScaleAxes(CoordinateSystem, CoordinateSystem)
     */
    public static MatrixSIS createTransform(final Envelope srcEnvelope, final AxisDirection[] srcAxes,
                                            final Envelope dstEnvelope, final AxisDirection[] dstAxes)
    {
        ArgumentChecks.ensureNonNull("srcEnvelope", srcEnvelope);
        ArgumentChecks.ensureNonNull("dstEnvelope", dstEnvelope);
        ArgumentChecks.ensureDimensionMatches("srcEnvelope", srcAxes.length, srcEnvelope);
        ArgumentChecks.ensureDimensionMatches("dstEnvelope", dstAxes.length, dstEnvelope);
        return createTransform(srcEnvelope, srcAxes, dstEnvelope, dstAxes, true);
    }

    /**
     * Creates a matrix for a transform that keep only a subset of source coordinate values.
     * The matrix size will be ({@code selectedDimensions.length} + 1) × ({@code sourceDimensions} + 1).
     * The matrix will contain only zero elements, except for the following cells which will contain 1:
     *
     * <ul>
     *   <li>The last column in the last row.</li>
     *   <li>For any row <var>j</var> other than the last row, the column {@code selectedDimensions[j]}.</li>
     * </ul>
     *
     * <div class="note"><b>Example:</b>
     * given (<var>x</var>,<var>y</var>,<var>z</var>,<var>t</var>) coordinate values, if one wants to keep
     * (<var>y</var>,<var>x</var>,<var>t</var>) coordinates (note the <var>x</var> ↔ <var>y</var> swapping)
     * and discard the <var>z</var> values, then the indices of source coordinates to select are 1 for <var>y</var>,
     * 0 for <var>x</var> and 3 for <var>t</var>. One can use the following method call:
     *
     * {@preformat java
     *   matrix = Matrices.createDimensionSelect(4, new int[] {1, 0, 3});
     * }
     *
     * The above method call will create the following 4×5 matrix,
     * which can be used for converting coordinates as below:
     *
     * {@preformat math
     *   ┌   ┐   ┌           ┐   ┌   ┐
     *   │ y │   │ 0 1 0 0 0 │   │ x │
     *   │ x │   │ 1 0 0 0 0 │   │ y │
     *   │ t │ = │ 0 0 0 1 0 │ × │ z │
     *   │ 1 │   │ 0 0 0 0 1 │   │ t │
     *   └   ┘   └           ┘   │ 1 │
     *                           └   ┘
     * }
     * </div>
     *
     * The inverse of the matrix created by this method will put {@link Double#NaN} values in the extra dimensions.
     * Other dimensions will work as expected.
     *
     * @param  sourceDimensions    the number of dimensions in source coordinates.
     * @param  selectedDimensions  the 0-based indices of source coordinate values to keep.
     *         The length of this array will be the number of dimensions in target coordinates.
     * @return an affine transform matrix keeping only the given source dimensions, and discarding all others.
     * @throws IllegalArgumentException if a value of {@code selectedDimensions} is lower than 0
     *         or not smaller than {@code sourceDimensions}.
     *
     * @see org.apache.sis.referencing.operation.transform.TransformSeparator
     */
    public static MatrixSIS createDimensionSelect(final int sourceDimensions, final int[] selectedDimensions) {
        final int numTargetDim = selectedDimensions.length;
        final MatrixSIS matrix = createZero(numTargetDim+1, sourceDimensions+1);
        for (int j=0; j<numTargetDim; j++) {
            final int i = selectedDimensions[j];
            ArgumentChecks.ensureValidIndex(sourceDimensions, i);
            matrix.setElement(j, i, 1);
        }
        matrix.setElement(numTargetDim, sourceDimensions, 1);
        return matrix;
    }

    /**
     * Creates a matrix which converts a subset of coordinates using the transform given by another matrix.
     * For example giving (<var>latitude</var>, <var>longitude</var>, <var>height</var>) coordinates,
     * a pass through operation can convert the height values from feet to metres without affecting
     * the (<var>latitude</var>, <var>longitude</var>) values.
     *
     * <p>The given sub-matrix shall have the following properties:</p>
     * <ul>
     *   <li>The last row often (but not necessarily) contains 0 values everywhere except in the last column.</li>
     *   <li>Values in the last column are translation terms, except in the last row.</li>
     *   <li>All other values are scale or shear terms.</li>
     * </ul>
     *
     * A square matrix complying with the above conditions is often {@linkplain #isAffine(Matrix) affine},
     * but this is not mandatory
     * (for example a <cite>perspective transform</cite> may contain non-zero values in the last row).
     *
     * <p>This method builds a new matrix with the following content:</p>
     * <ul>
     *   <li>An amount of {@code firstAffectedCoordinate} rows and columns are inserted before the first
     *       row and columns of the sub-matrix. The elements for the new rows and columns are set to 1
     *       on the diagonal, and 0 elsewhere.</li>
     *   <li>The sub-matrix - except for its last row and column - is copied in the new matrix starting
     *       at index ({@code firstAffectedCoordinate}, {@code firstAffectedCoordinate}).</li>
     *   <li>An amount of {@code numTrailingCoordinates} rows and columns are appended after the above sub-matrix.
     *       Their elements are set to 1 on the pseudo-diagonal ending in the lower-right corner, and 0 elsewhere.</li>
     *   <li>The last sub-matrix row is copied in the last row of the new matrix, and the last sub-matrix column
     *       is copied in the last column of the sub-matrix.</li>
     * </ul>
     *
     * <div class="note"><b>Example:</b>
     * given the following sub-matrix which converts height values from feet to metres before to subtracts 25 metres:
     *
     * {@preformat math
     *   ┌    ┐   ┌             ┐   ┌   ┐
     *   │ z' │ = │ 0.3048  -25 │ × │ z │
     *   │ 1  │   │ 0         1 │   │ 1 │
     *   └    ┘   └             ┘   └   ┘
     * }
     *
     * Then a call to {@code Matrices.createPassThrough(2, subMatrix, 1)} will return the following matrix,
     * which can be used for converting the height (<var>z</var>) without affecting the other coordinate values
     * (<var>x</var>,<var>y</var>,<var>t</var>):
     *
     * {@preformat math
     *   ┌    ┐   ┌                      ┐   ┌   ┐
     *   │ x  │   │ 1  0  0       0    0 │   │ x │
     *   │ y  │   │ 0  1  0       0    0 │   │ y │
     *   │ z' │ = │ 0  0  0.3048  0  -25 │ × │ z │
     *   │ t  │   │ 0  0  0       1    0 │   │ t │
     *   │ 1  │   │ 0  0  0       0    1 │   │ 1 │
     *   └    ┘   └                      ┘   └   ┘
     * }
     * </div>
     *
     * @param  firstAffectedCoordinate  the lowest index of the affected coordinates.
     * @param  subMatrix                the matrix to use for affected coordinates.
     * @param  numTrailingCoordinates   number of trailing coordinates to pass through.
     * @return a matrix for the same transform than the given matrix,
     *         augmented with leading and trailing pass-through coordinates.
     *
     * @see org.apache.sis.referencing.operation.transform.DefaultMathTransformFactory#createPassThroughTransform(int, MathTransform, int)
     */
    public static MatrixSIS createPassThrough(final int firstAffectedCoordinate,
            final Matrix subMatrix, final int numTrailingCoordinates)
    {
        ArgumentChecks.ensureNonNull ("subMatrix",               subMatrix);
        ArgumentChecks.ensurePositive("firstAffectedCoordinate", firstAffectedCoordinate);
        ArgumentChecks.ensurePositive("numTrailingCoordinates",  numTrailingCoordinates);
        final int  expansion = firstAffectedCoordinate + numTrailingCoordinates;
        int sourceDimensions = subMatrix.getNumCol();           // Will become the number of dimensions later.
        int targetDimensions = subMatrix.getNumRow();
        /*
         * Get data from the source matrix, together with the error terms if present.
         * The 'stride' and 'length' values will be used for computing indices in that array.
         * The DoubleDouble temporary object is used only if the array contains error terms.
         */
        final int      stride  = sourceDimensions;
        final int      length  = sourceDimensions * targetDimensions;
        final double[] sources = getExtendedElements(subMatrix);
        final DoubleDouble transfer = (sources.length > length) ? new DoubleDouble() : null;
        final MatrixSIS matrix = createZero(targetDimensions-- + expansion,
                                            sourceDimensions-- + expansion,
                                            transfer != null);
        /*
         * Following code processes from upper row to lower row.
         * First, set the diagonal elements on leading new dimensions.
         */
        for (int j=0; j<firstAffectedCoordinate; j++) {
            matrix.setElement(j, j, 1);
        }
        /*
         * Copy the sub-matrix, with special case for the translation terms
         * which are unconditionally stored in the last column.
         */
        final int lastColumn = sourceDimensions + expansion;
        matrix.setElements(sources, length, stride, transfer,
                0,                     0,                               // Source (row, colum)
                firstAffectedCoordinate, firstAffectedCoordinate,       // Target (row, column)
                targetDimensions,        sourceDimensions);             // Number of rows and columns to copy.
        matrix.setElements(sources, length, stride, transfer,
                0,                       sourceDimensions,              // Source (row, colum):  last column
                firstAffectedCoordinate, lastColumn,                    // Target (row, column): part of last column
                targetDimensions,        1);                            // Copy some rows of only 1 column.
        /*
         * Set the pseudo-diagonal elements on the trailing new dimensions.
         * 'diff' is zero for a square matrix and non-zero for rectangular matrix.
         */
        final int diff = targetDimensions - sourceDimensions;
        for (int i=lastColumn - numTrailingCoordinates; i<lastColumn; i++) {
            matrix.setElement(diff + i, i, 1);
        }
        /*
         * Copy the last row from the sub-matrix. In the usual affine transform,
         * this row contains only 0 element except for the last one, which is 1.
         */
        final int lastRow = targetDimensions + expansion;
        matrix.setElements(sources, length, stride, transfer,
                targetDimensions, 0,                                // Source (row, colum):  last row
                lastRow,          firstAffectedCoordinate,          // Target (row, column): part of last row
                1,                sourceDimensions);                // Copy some columns of only 1 row.
        matrix.setElements(sources, length, stride, transfer,
                targetDimensions, sourceDimensions,
                lastRow,          lastColumn,
                1,                1);
        return matrix;
    }

    /**
     * Returns a matrix with the same content than the given matrix but a different size, assuming an affine transform.
     * This method can be invoked for adding or removing the <strong>last</strong> dimensions of an affine transform.
     * More specifically:
     *
     * <ul class="verbose">
     *   <li>If the given {@code numCol} is <var>n</var> less than the number of columns in the given matrix,
     *       then the <var>n</var> columns <em>before the last column</em> are removed.
     *       The last column is left unchanged because it is assumed to contain the translation terms.</li>
     *   <li>If the given {@code numCol} is <var>n</var> more than the number of columns in the given matrix,
     *       then <var>n</var> columns are inserted <em>before the last column</em>.
     *       All values in the new columns will be zero.</li>
     *   <li>If the given {@code numRow} is <var>n</var> less than the number of rows in the given matrix,
     *       then the <var>n</var> rows <em>before the last row</em> are removed.
     *       The last row is left unchanged because it is assumed to contain the usual [0 0 0 … 1] terms.</li>
     *   <li>If the given {@code numRow} is <var>n</var> more than the number of rows in the given matrix,
     *       then <var>n</var> rows are inserted <em>before the last row</em>.
     *       The corresponding offset and scale factors will be 0 and 1 respectively.
     *       In other words, new dimensions are propagated unchanged.</li>
     * </ul>
     *
     * @param  matrix  the matrix to resize. This matrix will never be changed.
     * @param  numRow  the new number of rows. This is equal to the desired number of target dimensions plus 1.
     * @param  numCol  the new number of columns. This is equal to the desired number of source dimensions plus 1.
     * @return a new matrix of the given size, or the given {@code matrix} if no resizing was needed.
     */
    public static Matrix resizeAffine(final Matrix matrix, int numRow, int numCol) {
        ArgumentChecks.ensureNonNull         ("matrix", matrix);
        ArgumentChecks.ensureStrictlyPositive("numRow", numRow);
        ArgumentChecks.ensureStrictlyPositive("numCol", numCol);
        int srcRow = matrix.getNumRow();
        int srcCol = matrix.getNumCol();
        if (numRow == srcRow && numCol == srcCol) {
            return matrix;
        }
        final int      stride  = srcCol;
        final int      length  = srcCol * srcRow;
        final double[] sources = getExtendedElements(matrix);
        final DoubleDouble transfer = (sources.length > length) ? new DoubleDouble() : null;
        final MatrixSIS resized = createZero(numRow, numCol, transfer != null);
        final int copyRow = Math.min(--numRow, --srcRow);
        final int copyCol = Math.min(--numCol, --srcCol);
        for (int j=copyRow; j<numRow; j++) {
            resized.setElement(j, j, 1);
        }
        resized.setElements(sources, length, stride, transfer, 0,      0,       0,      0,       copyRow, copyCol);    // Shear and scale terms.
        resized.setElements(sources, length, stride, transfer, 0,      srcCol,  0,      numCol,  copyRow, 1);          // Translation column.
        resized.setElements(sources, length, stride, transfer, srcRow, 0,       numRow, 0,       1,       copyCol);    // Last row.
        resized.setElements(sources, length, stride, transfer, srcRow, srcCol,  numRow, numCol,  1,       1);          // Last row.
        return resized;
    }

    /**
     * Returns {@code true} if the given matrix is likely to use extended precision.
     * A value of {@code true} is not a guarantee that the matrix uses extended precision,
     * but a value of {@code false} is a guarantee that it does not.
     */
    private static boolean isExtendedPrecision(final Matrix matrix) {
        return (matrix instanceof MatrixSIS) ? ((MatrixSIS) matrix).isExtendedPrecision() :
               (matrix instanceof ExtendedPrecisionMatrix);  // Not guarantee that the matrix really uses double-double.
    }

    /**
     * Returns the elements of the given matrix, together with error terms if available.
     */
    private static double[] getExtendedElements(final Matrix matrix) {
        if (matrix instanceof ExtendedPrecisionMatrix) {
            return ((ExtendedPrecisionMatrix) matrix).getExtendedElements();
        } else {
            return MatrixSIS.castOrCopy(matrix).getElements();
        }
    }

    /**
     * Creates a new matrix which is a copy of the given matrix.
     *
     * <div class="note"><b>Implementation note:</b>
     * For square matrix with a size between {@value org.apache.sis.referencing.operation.matrix.Matrix1#SIZE}
     * and {@value org.apache.sis.referencing.operation.matrix.Matrix4#SIZE} inclusive, the returned matrix is
     * usually an instance of one of {@link Matrix1} … {@link Matrix4} subtypes.</div>
     *
     * @param  matrix  the matrix to copy, or {@code null}.
     * @return a copy of the given matrix, or {@code null} if the given matrix was null.
     *
     * @see MatrixSIS#clone()
     * @see MatrixSIS#castOrCopy(Matrix)
     */
    public static MatrixSIS copy(final Matrix matrix) {
        if (matrix == null) {
            return null;
        }
        final int size = matrix.getNumRow();
        if (size != matrix.getNumCol()) {
            return new NonSquareMatrix(matrix);
        }
        if (!isExtendedPrecision(matrix)) {
            switch (size) {
                case 1: return new Matrix1(matrix);
                case 2: return new Matrix2(matrix);
                case 3: return new Matrix3(matrix);
                case 4: return new Matrix4(matrix);
            }
        }
        return new GeneralMatrix(matrix);
    }

    /**
     * Returns an unmodifiable view of the given matrix. The returned matrix is immutable
     * only if the given {@code matrix} is not modified anymore after this method call.
     *
     * @param  matrix  the matrix for which to get an unmodifiable view, or {@code null}.
     * @return a unmodifiable view of the given matrix, or {@code null} if the given matrix was null.
     *
     * @since 0.6
     */
    public static MatrixSIS unmodifiable(final Matrix matrix) {
        if (matrix == null || matrix instanceof UnmodifiableMatrix) {
            return (MatrixSIS) matrix;
        } else {
            return new UnmodifiableMatrix(matrix);
        }
    }

    /**
     * Returns a new matrix which is the result of multiplying the first matrix with the second one.
     * In other words, returns {@code m1} × {@code m2}.
     *
     * @param  m1  the first matrix to multiply.
     * @param  m2  the second matrix to multiply.
     * @return the result of {@code m1} × {@code m2}.
     * @throws MismatchedMatrixSizeException if the number of columns in {@code m1} is not equals to the
     *         number of rows in {@code m2}.
     *
     * @see MatrixSIS#multiply(Matrix)
     *
     * @since 0.6
     */
    public static MatrixSIS multiply(final Matrix m1, final Matrix m2) throws MismatchedMatrixSizeException {
        if (m1 instanceof MatrixSIS) {
            return ((MatrixSIS) m1).multiply(m2);           // Maybe the subclass overrides that method.
        }
        final int nc = m2.getNumCol();
        MatrixSIS.ensureNumRowMatch(m1.getNumCol(), m2.getNumRow(), nc);
        final GeneralMatrix result = GeneralMatrix.createExtendedPrecision(m1.getNumRow(), nc, false);
        result.setToProduct(m1, m2);
        return result;
    }

    /**
     * Returns the inverse of the given matrix.
     *
     * @param  matrix  the matrix to inverse, or {@code null}.
     * @return the inverse of this matrix, or {@code null} if the given matrix was null.
     * @throws NoninvertibleMatrixException if the given matrix is not invertible.
     *
     * @see MatrixSIS#inverse()
     *
     * @since 0.6
     */
    public static MatrixSIS inverse(final Matrix matrix) throws NoninvertibleMatrixException {
        if (matrix == null) {
            return null;
        } else if (matrix instanceof MatrixSIS) {
            return ((MatrixSIS) matrix).inverse();  // Maybe the subclass override that method.
        } else if (matrix.getNumRow() != matrix.getNumCol()) {
            return new NonSquareMatrix(matrix).inverse();
        } else {
            return Solver.inverse(matrix, true);
        }
    }

    /**
     * Returns {@code true} if the given matrix represents an affine transform.
     * A transform is affine if the matrix is square and its last row contains
     * only zeros, except in the last column which contains 1.
     *
     * @param  matrix  the matrix to test.
     * @return {@code true} if the matrix represents an affine transform.
     *
     * @see MatrixSIS#isAffine()
     */
    public static boolean isAffine(final Matrix matrix) {
        if (matrix instanceof MatrixSIS) {
            return ((MatrixSIS) matrix).isAffine();
        } else {
            return MatrixSIS.isAffine(matrix);
        }
    }

    /**
     * Returns {@code true} if the given matrix represents a translation.
     * This method returns {@code true} if the given matrix {@linkplain #isAffine(Matrix) is affine}
     * and differs from the identity matrix only in the last column.
     *
     * @param  matrix  the matrix to test.
     * @return {@code true} if the matrix represents a translation.
     *
     * @since 0.7
     */
    public static boolean isTranslation(final Matrix matrix) {
        if (!isAffine(matrix)) {
            return false;
        }
        final int numRow = matrix.getNumRow() - 1;      // Excluding translation column.
        final int numCol = matrix.getNumCol() - 1;      // Excluding last row in affine transform.
        for (int j=0; j<numRow; j++) {
            for (int i=0; i<numCol; i++) {
                if (matrix.getElement(j,i) != ((i == j) ? 1 : 0)) {
                    return false;
                }
            }
        }
        return true;
    }

    /**
     * Returns {@code true} if the given matrix is close to an identity matrix, given a tolerance threshold.
     * This method is equivalent to computing the difference between the given matrix and an identity matrix
     * of identical size, and returning {@code true} if and only if all differences are smaller than or equal
     * to {@code tolerance}.
     *
     * <p><b>Caution:</b> {@linkplain org.apache.sis.referencing.datum.BursaWolfParameters Bursa-Wolf parameters},
     * when represented as a matrix, are close to an identity transform and could easily be confused with rounding
     * errors. In case of doubt, it is often safer to use the strict {@link MatrixSIS#isIdentity()} method instead
     * than this one.</p>
     *
     * @param  matrix     the matrix to test for identity.
     * @param  tolerance  the tolerance value, or 0 for a strict comparison.
     * @return {@code true} if this matrix is close to the identity matrix given the tolerance threshold.
     *
     * @see MatrixSIS#isIdentity()
     */
    public static boolean isIdentity(final Matrix matrix, final double tolerance) {
        final int numRow = matrix.getNumRow();
        final int numCol = matrix.getNumCol();
        if (numRow != numCol) {
            return false;
        }
        for (int j=0; j<numRow; j++) {
            for (int i=0; i<numCol; i++) {
                double e = matrix.getElement(j,i);
                if (i == j) {
                    e--;
                }
                if (!(Math.abs(e) <= tolerance)) {              // Uses '!' in order to catch NaN values.
                    return false;
                }
            }
        }
        return true;
    }

    /**
     * Compares the given matrices for equality, using the given relative or absolute tolerance threshold.
     * The matrix elements are compared as below:
     *
     * <ul>
     *   <li>{@link Double#NaN} values are considered equals to all other NaN values</li>
     *   <li>Infinite values are considered equal to other infinite values of the same sign</li>
     *   <li>All other values are considered equal if the absolute value of their difference is
     *       smaller than or equals to the threshold described below.</li>
     * </ul>
     *
     * If {@code relative} is {@code true}, then for any pair of values <var>v1</var><sub>j,i</sub>
     * and <var>v2</var><sub>j,i</sub> to compare, the tolerance threshold is scaled by
     * {@code max(abs(v1), abs(v2))}. Otherwise the threshold is used as-is.
     *
     * @param  m1        the first matrix to compare, or {@code null}.
     * @param  m2        the second matrix to compare, or {@code null}.
     * @param  epsilon   the tolerance value.
     * @param  relative  if {@code true}, then the tolerance value is relative to the magnitude
     *                   of the matrix elements being compared.
     * @return {@code true} if the values of the two matrix do not differ by a quantity greater
     *         than the given tolerance threshold.
     *
     * @see MatrixSIS#equals(Matrix, double)
     */
    public static boolean equals(final Matrix m1, final Matrix m2, final double epsilon, final boolean relative) {
        if (m1 != m2) {
            if (m1 == null || m2 == null) {
                return false;
            }
            final int numRow = m1.getNumRow();
            if (numRow != m2.getNumRow()) {
                return false;
            }
            final int numCol = m1.getNumCol();
            if (numCol != m2.getNumCol()) {
                return false;
            }
            for (int j=0; j<numRow; j++) {
                for (int i=0; i<numCol; i++) {
                    final double v1 = m1.getElement(j, i);
                    final double v2 = m2.getElement(j, i);
                    double tolerance = epsilon;
                    if (relative) {
                        tolerance *= Math.max(Math.abs(v1), Math.abs(v2));
                    }
                    if (!(Math.abs(v1 - v2) <= tolerance)) {
                        if (Numerics.equals(v1, v2)) {
                            // Special case for NaN and infinite values.
                            continue;
                        }
                        return false;
                    }
                }
            }
        }
        return true;
    }

    /**
     * Compares the given matrices for equality, using the given comparison strictness level.
     * To be considered equal, the two matrices must met the following conditions, which depend
     * on the {@code mode} argument:
     *
     * <ul>
     *   <li>{@link ComparisonMode#STRICT STRICT}:
     *       the two matrices must be of the same class, have the same size and the same element values.</li>
     *   <li>{@link ComparisonMode#BY_CONTRACT BY_CONTRACT}:
     *       the two matrices must have the same size and the same element values,
     *       but are not required to be the same implementation class (any {@link Matrix} is okay).</li>
     *   <li>{@link ComparisonMode#IGNORE_METADATA IGNORE_METADATA}: same as {@code BY_CONTRACT},
     *       since matrices have no metadata.</li>
     *   <li>{@link ComparisonMode#APPROXIMATE APPROXIMATE}:
     *       the two matrices must have the same size, but the element values can differ up to some threshold.
     *       The threshold value is determined empirically and may change in any future SIS versions.
     *       For more control, use {@link #equals(Matrix, Matrix, double, boolean)} instead.</li>
     * </ul>
     *
     * @param  m1    the first matrix to compare, or {@code null}.
     * @param  m2    the second matrix to compare, or {@code null}.
     * @param  mode  the strictness level of the comparison.
     * @return {@code true} if both matrices are equal.
     *
     * @see MatrixSIS#equals(Object, ComparisonMode)
     */
    public static boolean equals(final Matrix m1, final Matrix m2, final ComparisonMode mode) {
        switch (mode) {
            case STRICT:          return Objects.equals(m1, m2);
            case BY_CONTRACT:     // Fall through
            case IGNORE_METADATA: return equals(m1, m2, 0, false);
            case DEBUG:           // Fall through
            case ALLOW_VARIANT:   // Fall through
            case APPROXIMATE:     return equals(m1, m2, Numerics.COMPARISON_THRESHOLD, true);
            default: throw new IllegalArgumentException(Errors.format(
                    Errors.Keys.UnknownEnumValue_2, ComparisonMode.class, mode));
        }
    }

    /**
     * Returns a unlocalized string representation of the given matrix.
     * For each column, the numbers are aligned on the decimal separator.
     *
     * <p>The current implementation formats ±0 and ±1 without trailing {@code ".0"}, and all other values with
     * a per-column uniform number of fraction digits. The ±0 and ±1 values are treated especially because they
     * usually imply a <cite>"no scale"</cite>, <cite>"no translation"</cite> or <cite>"orthogonal axes"</cite>
     * meaning. A matrix in SIS is often populated mostly by ±0 and ±1 values, with a few "interesting" values.
     * The simpler ±0 and ±1 formatting makes easier to spot the "interesting" values.</p>
     *
     * <p>The following example shows the string representation of an affine transform which swap
     * (<var>latitude</var>, <var>longitude</var>) axes, converts degrees to radians and converts
     * height values from feet to metres:</p>
     *
     * {@preformat math
     *   ┌                                                       ┐
     *   │ 0                     0.017453292519943295  0       0 │
     *   │ 0.017453292519943295  0                     0       0 │
     *   │ 0                     0                     0.3048  0 │
     *   │ 0                     0                     0       1 │
     *   └                                                       ┘
     * }
     *
     * <div class="note"><b>Note:</b>
     * Formatting on a per-column basis is convenient for the kind of matrices used in referencing by coordinates,
     * because each column is typically a displacement vector in a different dimension of the source coordinate
     * reference system. In addition, the last column is often a translation vector having a magnitude very
     * different than the other columns.</div>
     *
     * @param  matrix  the matrix for which to get a string representation.
     * @return a string representation of the given matrix.
     */
    public static String toString(final Matrix matrix) {
        final int numRow = matrix.getNumRow();
        final int numCol = matrix.getNumCol();
        final String[]  elements            = new String [numCol * numRow];     // String representation of matrix values.
        final boolean[] noFractionDigits    = new boolean[numCol * numRow];     // Whether to remove the trailing ".0" for a given number.
        final boolean[] hasDecimalSeparator = new boolean[numCol];              // Whether the column has at least one number where fraction digits are shown.
        final byte[] maximumFractionDigits  = new byte   [numCol];              // The greatest amount of fraction digits found in a column.
        final byte[] maximumPaddingZeros    = new byte   [numCol * numRow];     // Maximal amount of zeros that we can append before to exceed the IEEE 754 accuracy.
        final byte[] widthBeforeFraction    = new byte   [numCol];              // Number of characters before the fraction digits: spacing + ('-') + integerDigits + '.'
        final byte[] columnWidth            = new byte   [numCol];              // Total column width.
        int totalWidth = 1;
        /*
         * Create now the string representation of all matrix elements and measure the width
         * of the integer field and the fraction field, then the total width of each column.
         */
        int spacing = 1;                // Spacing is 1 before the first column only, then SPACING for other columns.
        for (int i=0; i<numCol; i++) {
            for (int j=0; j<numRow; j++) {
                final int flatIndex = j*numCol + i;
                final double value  = matrix.getElement(j,i);
                String element = Double.toString(value);
                final int width;
                /*
                 * Special case for ±0 and ±1 (because those values appear very often and have
                 * a particular meaning): for those values, we will ignore the fraction digits.
                 * For all other values, we will format all fraction digits.
                 */
                if (value == -1 || value == 0 || value == +1) {
                    noFractionDigits[flatIndex] = true;
                    width = spacing + element.length() - 2;           // The -2 is for ignoring the trailing ".0"
                    widthBeforeFraction[i] = (byte) Math.max(widthBeforeFraction[i], width);
                } else {
                    /*
                     * All values other than ±0 and ±1. If the values is NaN or infinity (in which case there is
                     * no decimal separator), give all spaces to the "before fraction" side for right-alignment.
                     */
                    int s = element.lastIndexOf('.');
                    if (s < 0) {
                        element = CharSequences.replace(element, "Infinity", "∞").toString();
                        width = spacing + element.length();
                        widthBeforeFraction[i] = (byte) Math.max(widthBeforeFraction[i], width);
                    } else {
                        /*
                         * All values other than ±0, ±1, NaN and infinity. We store separately the width before
                         * and after the decimal separator. The width before the separator contains the spacing
                         * between cells.
                         */
                        hasDecimalSeparator[i] = true;
                        final int numFractionDigits = element.length() - ++s;
                        width = (widthBeforeFraction  [i] = (byte) Math.max(widthBeforeFraction  [i], spacing + s))
                              + (maximumFractionDigits[i] = (byte) Math.max(maximumFractionDigits[i], numFractionDigits));
                        /*
                         * If the number use exponential notation, we will not be allowed to append any zero.
                         * Otherwise we will append some zeros for right-alignment, but without exceeding the
                         * IEEE 754 'double' accuracy for not giving a false sense of precision.
                         */
                        if (element.indexOf('E') < 0) {
                            final int accuracy = -DecimalFunctions.floorLog10(Math.ulp(value));
                            maximumPaddingZeros[flatIndex] = (byte) (accuracy - numFractionDigits);
                        }
                    }
                }
                columnWidth[i] = (byte) Math.max(columnWidth[i], width);
                elements[flatIndex] = element;
            }
            totalWidth += columnWidth[i];
            spacing = SPACING;                              // Spacing before all columns after the first one.
        }
        /*
         * Now append the formatted elements with the appropriate amount of spaces before each value,
         * and trailling zeros after each value except ±0, ±1, NaN and infinities.
         */
        final String   lineSeparator = System.lineSeparator();
        final CharSequence whiteLine = CharSequences.spaces(totalWidth);
        final StringBuilder   buffer = new StringBuilder((totalWidth + 2 + lineSeparator.length()) * (numRow + 2));
        buffer.append('┌').append(whiteLine).append('┐').append(lineSeparator);
        int flatIndex = 0;
        for (int j=0; j<numRow; j++) {
            buffer.append('│');
            for (int i=0; i<numCol; i++) {
                final String element = elements[flatIndex];
                final int width = element.length();
                int spaces, s = element.lastIndexOf('.');
                if (s >= 0) {
                    if (hasDecimalSeparator[i]) s++;
                    spaces = widthBeforeFraction[i] - s;    // Number of spaces for alignment on the decimal separator
                } else {
                    spaces = columnWidth[i] - width;        // Number of spaces for right alignment (NaN or ∞ cases)
                }
                buffer.append(CharSequences.spaces(spaces)).append(element);
                if (s >= 0) {
                    /*
                     * Append trailing spaces for ±0 and ±1 values,
                     * or trailing zeros for all other real values.
                     */
                    s += maximumFractionDigits[i] - width;
                    if (noFractionDigits[flatIndex]) {
                        buffer.setLength(buffer.length() - 2);      // Erase the trailing ".0"
                        s += 2;
                    } else {
                        int n = Math.min(s, maximumPaddingZeros[flatIndex]);
                        StringBuilders.repeat(buffer, '0', n);
                        s -= n;
                    }
                    buffer.append(CharSequences.spaces(s));
                }
                flatIndex++;
            }
            buffer.append(" │").append(lineSeparator);
        }
        return buffer.append('└').append(whiteLine).append('┘').append(lineSeparator).toString();
    }
}