You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by er...@apache.org on 2021/08/11 11:31:33 UTC
[commons-math] branch master updated: Delete spurious files.
This is an automated email from the ASF dual-hosted git repository.
erans pushed a commit to branch master
in repository https://gitbox.apache.org/repos/asf/commons-math.git
The following commit(s) were added to refs/heads/master by this push:
new 9b7a2c8 Delete spurious files.
9b7a2c8 is described below
commit 9b7a2c8edca579929625f46f55b67eaffb3244ad
Author: Gilles Sadowski <gi...@gmail.com>
AuthorDate: Wed Aug 11 13:30:44 2021 +0200
Delete spurious files.
Files were committed by mistake.
---
.../legacy/linear/EigenDecomposition.java_PRINT | 945 ---------------------
.../scalar/noderiv/SimplexOptimizer.java.DEBUG | 339 --------
2 files changed, 1284 deletions(-)
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/linear/EigenDecomposition.java_PRINT b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/linear/EigenDecomposition.java_PRINT
deleted file mode 100644
index b74e3e0..0000000
--- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/linear/EigenDecomposition.java_PRINT
+++ /dev/null
@@ -1,945 +0,0 @@
-/*
- * 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.commons.math4.linear;
-
-import org.apache.commons.numbers.complex.Complex;
-import org.apache.commons.numbers.core.Precision;
-import org.apache.commons.math4.exception.DimensionMismatchException;
-import org.apache.commons.math4.exception.MathArithmeticException;
-import org.apache.commons.math4.exception.MathUnsupportedOperationException;
-import org.apache.commons.math4.exception.MaxCountExceededException;
-import org.apache.commons.math4.exception.util.LocalizedFormats;
-import org.apache.commons.math4.util.FastMath;
-
-/**
- * Calculates the eigen decomposition of a real matrix.
- * <p>
- * The eigen decomposition of matrix A is a set of two matrices:
- * V and D such that A = V × D × V<sup>T</sup>.
- * A, V and D are all m × m matrices.
- * <p>
- * This class is similar in spirit to the {@code EigenvalueDecomposition}
- * class from the <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a>
- * library, with the following changes:
- * <ul>
- * <li>a {@link #getVT() getVt} method has been added,</li>
- * <li>two {@link #getRealEigenvalue(int) getRealEigenvalue} and
- * {@link #getImagEigenvalue(int) getImagEigenvalue} methods to pick up a
- * single eigenvalue have been added,</li>
- * <li>a {@link #getEigenvector(int) getEigenvector} method to pick up a
- * single eigenvector has been added,</li>
- * <li>a {@link #getDeterminant() getDeterminant} method has been added.</li>
- * <li>a {@link #getSolver() getSolver} method has been added.</li>
- * </ul>
- * <p>
- * As of 3.1, this class supports general real matrices (both symmetric and non-symmetric):
- * <p>
- * If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is diagonal
- * and the eigenvector matrix V is orthogonal, i.e.
- * {@code A = V.multiply(D.multiply(V.transpose()))} and
- * {@code V.multiply(V.transpose())} equals the identity matrix.
- * </p>
- * <p>
- * If A is not symmetric, then the eigenvalue matrix D is block diagonal with the real
- * eigenvalues in 1-by-1 blocks and any complex eigenvalues, lambda + i*mu, in 2-by-2
- * blocks:
- * <pre>
- * [lambda, mu ]
- * [ -mu, lambda]
- * </pre>
- * The columns of V represent the eigenvectors in the sense that {@code A*V = V*D},
- * i.e. A.multiply(V) equals V.multiply(D).
- * The matrix V may be badly conditioned, or even singular, so the validity of the
- * equation {@code A = V*D*inverse(V)} depends upon the condition of V.
- * <p>
- * This implementation is based on the paper by A. Drubrulle, R.S. Martin and
- * J.H. Wilkinson "The Implicit QL Algorithm" in Wilksinson and Reinsch (1971)
- * Handbook for automatic computation, vol. 2, Linear algebra, Springer-Verlag,
- * New-York.
- *
- * @see <a href="http://mathworld.wolfram.com/EigenDecomposition.html">MathWorld</a>
- * @see <a href="http://en.wikipedia.org/wiki/Eigendecomposition_of_a_matrix">Wikipedia</a>
- * @since 2.0 (changed to concrete class in 3.0)
- */
-public class EigenDecomposition {
- /** Internally used epsilon criteria. */
- private static final double EPSILON = 1e-12;
- /** Maximum number of iterations accepted in the implicit QL transformation */
- private static final byte MAX_ITER = 30;
- /** Main diagonal of the tridiagonal matrix. */
- private double[] main;
- /** Secondary diagonal of the tridiagonal matrix. */
- private double[] secondary;
- /**
- * Transformer to tridiagonal (may be null if matrix is already
- * tridiagonal).
- */
- private TriDiagonalTransformer transformer;
- /** Real part of the realEigenvalues. */
- private double[] realEigenvalues;
- /** Imaginary part of the realEigenvalues. */
- private double[] imagEigenvalues;
- /** Eigenvectors. */
- private ArrayRealVector[] eigenvectors;
- /** Cached value of V. */
- private RealMatrix cachedV;
- /** Cached value of D. */
- private RealMatrix cachedD;
- /** Cached value of Vt. */
- private RealMatrix cachedVt;
- /** Whether the matrix is symmetric. */
- private final boolean isSymmetric;
-
- /**
- * Calculates the eigen decomposition of the given real matrix.
- * <p>
- * Supports decomposition of a general matrix since 3.1.
- *
- * @param matrix Matrix to decompose.
- * @throws MaxCountExceededException if the algorithm fails to converge.
- * @throws MathArithmeticException if the decomposition of a general matrix
- * results in a matrix with zero norm
- * @since 3.1
- */
- public EigenDecomposition(final RealMatrix matrix)
- throws MathArithmeticException {
- final double symTol = 10 * matrix.getRowDimension() * matrix.getColumnDimension() * Precision.EPSILON;
- isSymmetric = MatrixUtils.isSymmetric(matrix, symTol);
- if (isSymmetric) {
- transformToTridiagonal(matrix);
- findEigenVectors(transformer.getQ().getData());
- } else {
- final SchurTransformer t = transformToSchur(matrix);
- findEigenVectorsFromSchur(t);
- }
- }
-
- /**
- * Calculates the eigen decomposition of the symmetric tridiagonal
- * matrix. The Householder matrix is assumed to be the identity matrix.
- *
- * @param main Main diagonal of the symmetric tridiagonal form.
- * @param secondary Secondary of the tridiagonal form.
- * @throws MaxCountExceededException if the algorithm fails to converge.
- * @since 3.1
- */
- public EigenDecomposition(final double[] main, final double[] secondary) {
- isSymmetric = true;
- this.main = main.clone();
- this.secondary = secondary.clone();
- transformer = null;
- final int size = main.length;
- final double[][] z = new double[size][size];
- for (int i = 0; i < size; i++) {
- z[i][i] = 1.0;
- }
- findEigenVectors(z);
- }
-
- /**
- * Gets the matrix V of the decomposition.
- * V is an orthogonal matrix, i.e. its transpose is also its inverse.
- * The columns of V are the eigenvectors of the original matrix.
- * No assumption is made about the orientation of the system axes formed
- * by the columns of V (e.g. in a 3-dimension space, V can form a left-
- * or right-handed system).
- *
- * @return the V matrix.
- */
- public RealMatrix getV() {
-
- if (cachedV == null) {
- final int m = eigenvectors.length;
- cachedV = MatrixUtils.createRealMatrix(m, m);
- for (int k = 0; k < m; ++k) {
- cachedV.setColumnVector(k, eigenvectors[k]);
- }
- }
- // return the cached matrix
- return cachedV;
- }
-
- /**
- * Gets the block diagonal matrix D of the decomposition.
- * D is a block diagonal matrix.
- * Real eigenvalues are on the diagonal while complex values are on
- * 2x2 blocks { {real +imaginary}, {-imaginary, real} }.
- *
- * @return the D matrix.
- *
- * @see #getRealEigenvalues()
- * @see #getImagEigenvalues()
- */
- public RealMatrix getD() {
-
- if (cachedD == null) {
- // cache the matrix for subsequent calls
- cachedD = MatrixUtils.createRealMatrixWithDiagonal(realEigenvalues);
-
- for (int i = 0; i < imagEigenvalues.length; i++) {
- if (Precision.compareTo(imagEigenvalues[i], 0.0, EPSILON) > 0) {
- cachedD.setEntry(i, i+1, imagEigenvalues[i]);
- } else if (Precision.compareTo(imagEigenvalues[i], 0.0, EPSILON) < 0) {
- cachedD.setEntry(i, i-1, imagEigenvalues[i]);
- }
- }
- }
- return cachedD;
- }
-
- /**
- * Gets the transpose of the matrix V of the decomposition.
- * V is an orthogonal matrix, i.e. its transpose is also its inverse.
- * The columns of V are the eigenvectors of the original matrix.
- * No assumption is made about the orientation of the system axes formed
- * by the columns of V (e.g. in a 3-dimension space, V can form a left-
- * or right-handed system).
- *
- * @return the transpose of the V matrix.
- */
- public RealMatrix getVT() {
-
- if (cachedVt == null) {
- final int m = eigenvectors.length;
- cachedVt = MatrixUtils.createRealMatrix(m, m);
- for (int k = 0; k < m; ++k) {
- cachedVt.setRowVector(k, eigenvectors[k]);
- }
- }
-
- // return the cached matrix
- return cachedVt;
- }
-
- /**
- * Returns whether the calculated eigen values are complex or real.
- * <p>The method performs a zero check for each element of the
- * {@link #getImagEigenvalues()} array and returns {@code true} if any
- * element is not equal to zero.
- *
- * @return {@code true} if the eigen values are complex, {@code false} otherwise
- * @since 3.1
- */
- public boolean hasComplexEigenvalues() {
- for (int i = 0; i < imagEigenvalues.length; i++) {
- if (!Precision.equals(imagEigenvalues[i], 0.0, EPSILON)) {
- return true;
- }
- }
- return false;
- }
-
- /**
- * Gets a copy of the real parts of the eigenvalues of the original matrix.
- *
- * @return a copy of the real parts of the eigenvalues of the original matrix.
- *
- * @see #getD()
- * @see #getRealEigenvalue(int)
- * @see #getImagEigenvalues()
- */
- public double[] getRealEigenvalues() {
- return realEigenvalues.clone();
- }
-
- /**
- * Returns the real part of the i<sup>th</sup> eigenvalue of the original
- * matrix.
- *
- * @param i index of the eigenvalue (counting from 0)
- * @return real part of the i<sup>th</sup> eigenvalue of the original
- * matrix.
- *
- * @see #getD()
- * @see #getRealEigenvalues()
- * @see #getImagEigenvalue(int)
- */
- public double getRealEigenvalue(final int i) {
- return realEigenvalues[i];
- }
-
- /**
- * Gets a copy of the imaginary parts of the eigenvalues of the original
- * matrix.
- *
- * @return a copy of the imaginary parts of the eigenvalues of the original
- * matrix.
- *
- * @see #getD()
- * @see #getImagEigenvalue(int)
- * @see #getRealEigenvalues()
- */
- public double[] getImagEigenvalues() {
- return imagEigenvalues.clone();
- }
-
- /**
- * Gets the imaginary part of the i<sup>th</sup> eigenvalue of the original
- * matrix.
- *
- * @param i Index of the eigenvalue (counting from 0).
- * @return the imaginary part of the i<sup>th</sup> eigenvalue of the original
- * matrix.
- *
- * @see #getD()
- * @see #getImagEigenvalues()
- * @see #getRealEigenvalue(int)
- */
- public double getImagEigenvalue(final int i) {
- return imagEigenvalues[i];
- }
-
- /**
- * Gets a copy of the i<sup>th</sup> eigenvector of the original matrix.
- *
- * @param i Index of the eigenvector (counting from 0).
- * @return a copy of the i<sup>th</sup> eigenvector of the original matrix.
- * @see #getD()
- */
- public RealVector getEigenvector(final int i) {
- return eigenvectors[i].copy();
- }
-
- /**
- * Computes the determinant of the matrix.
- *
- * @return the determinant of the matrix.
- */
- public double getDeterminant() {
- double determinant = 1;
- for (double lambda : realEigenvalues) {
- determinant *= lambda;
- }
- return determinant;
- }
-
- /**
- * Computes the square-root of the matrix.
- * This implementation assumes that the matrix is symmetric and positive
- * definite.
- *
- * @return the square-root of the matrix.
- * @throws MathUnsupportedOperationException if the matrix is not
- * symmetric or not positive definite.
- * @since 3.1
- */
- public RealMatrix getSquareRoot() {
- if (!isSymmetric) {
- throw new MathUnsupportedOperationException();
- }
-
- final double[] sqrtEigenValues = new double[realEigenvalues.length];
- for (int i = 0; i < realEigenvalues.length; i++) {
- final double eigen = realEigenvalues[i];
- if (eigen <= 0) {
- throw new MathUnsupportedOperationException();
- }
- sqrtEigenValues[i] = FastMath.sqrt(eigen);
- }
- final RealMatrix sqrtEigen = MatrixUtils.createRealDiagonalMatrix(sqrtEigenValues);
- final RealMatrix v = getV();
- final RealMatrix vT = getVT();
-
- return v.multiply(sqrtEigen).multiply(vT);
- }
-
- /**
- * Gets a solver for finding the A × X = B solution in exact
- * linear sense.
- * <p>
- * Since 3.1, eigen decomposition of a general matrix is supported,
- * but the {@link DecompositionSolver} only supports real eigenvalues.
- *
- * @return a solver
- * @throws MathUnsupportedOperationException if the decomposition resulted in
- * complex eigenvalues
- */
- public DecompositionSolver getSolver() {
- if (hasComplexEigenvalues()) {
- throw new MathUnsupportedOperationException();
- }
- return new Solver(realEigenvalues, imagEigenvalues, eigenvectors);
- }
-
- /** Specialized solver. */
- private static class Solver implements DecompositionSolver {
- /** Real part of the realEigenvalues. */
- private final double[] realEigenvalues;
- /** Imaginary part of the realEigenvalues. */
- private final double[] imagEigenvalues;
- /** Eigenvectors. */
- private final ArrayRealVector[] eigenvectors;
-
- /**
- * Builds a solver from decomposed matrix.
- *
- * @param realEigenvalues Real parts of the eigenvalues.
- * @param imagEigenvalues Imaginary parts of the eigenvalues.
- * @param eigenvectors Eigenvectors.
- */
- private Solver(final double[] realEigenvalues,
- final double[] imagEigenvalues,
- final ArrayRealVector[] eigenvectors) {
- this.realEigenvalues = realEigenvalues;
- this.imagEigenvalues = imagEigenvalues;
- this.eigenvectors = eigenvectors;
- }
-
- /**
- * Solves the linear equation A × X = B for symmetric matrices A.
- * <p>
- * This method only finds exact linear solutions, i.e. solutions for
- * which ||A × X - B|| is exactly 0.
- * </p>
- *
- * @param b Right-hand side of the equation A × X = B.
- * @return a Vector X that minimizes the two norm of A × X - B.
- *
- * @throws DimensionMismatchException if the matrices dimensions do not match.
- * @throws SingularMatrixException if the decomposed matrix is singular.
- */
- @Override
- public RealVector solve(final RealVector b) {
- if (!isNonSingular()) {
- throw new SingularMatrixException();
- }
-
- final int m = realEigenvalues.length;
- if (b.getDimension() != m) {
- throw new DimensionMismatchException(b.getDimension(), m);
- }
-
- final double[] bp = new double[m];
- for (int i = 0; i < m; ++i) {
- final ArrayRealVector v = eigenvectors[i];
- final double[] vData = v.getDataRef();
- final double s = v.dotProduct(b) / realEigenvalues[i];
- for (int j = 0; j < m; ++j) {
- bp[j] += s * vData[j];
- }
- }
-
- return new ArrayRealVector(bp, false);
- }
-
- /** {@inheritDoc} */
- @Override
- public RealMatrix solve(RealMatrix b) {
-
- if (!isNonSingular()) {
- throw new SingularMatrixException();
- }
-
- final int m = realEigenvalues.length;
- if (b.getRowDimension() != m) {
- throw new DimensionMismatchException(b.getRowDimension(), m);
- }
-
- final int nColB = b.getColumnDimension();
- final double[][] bp = new double[m][nColB];
- final double[] tmpCol = new double[m];
- for (int k = 0; k < nColB; ++k) {
- for (int i = 0; i < m; ++i) {
- tmpCol[i] = b.getEntry(i, k);
- bp[i][k] = 0;
- }
- for (int i = 0; i < m; ++i) {
- final ArrayRealVector v = eigenvectors[i];
- final double[] vData = v.getDataRef();
- double s = 0;
- for (int j = 0; j < m; ++j) {
- s += v.getEntry(j) * tmpCol[j];
- }
- s /= realEigenvalues[i];
- for (int j = 0; j < m; ++j) {
- bp[j][k] += s * vData[j];
- }
- }
- }
-
- return new Array2DRowRealMatrix(bp, false);
-
- }
-
- /**
- * Checks whether the decomposed matrix is non-singular.
- *
- * @return true if the decomposed matrix is non-singular.
- */
- @Override
- public boolean isNonSingular() {
- double largestEigenvalueNorm = 0.0;
- // Looping over all values (in case they are not sorted in decreasing
- // order of their norm).
- for (int i = 0; i < realEigenvalues.length; ++i) {
- largestEigenvalueNorm = FastMath.max(largestEigenvalueNorm, eigenvalueNorm(i));
- }
- // Corner case: zero matrix, all exactly 0 eigenvalues
- if (largestEigenvalueNorm == 0.0) {
- return false;
- }
- for (int i = 0; i < realEigenvalues.length; ++i) {
- // Looking for eigenvalues that are 0, where we consider anything much much smaller
- // than the largest eigenvalue to be effectively 0.
- if (Precision.equals(eigenvalueNorm(i) / largestEigenvalueNorm, 0, EPSILON)) {
- return false;
- }
- }
- return true;
- }
-
- /**
- * @param i which eigenvalue to find the norm of
- * @return the norm of ith (complex) eigenvalue.
- */
- private double eigenvalueNorm(int i) {
- final double re = realEigenvalues[i];
- final double im = imagEigenvalues[i];
- return FastMath.sqrt(re * re + im * im);
- }
-
- /**
- * Get the inverse of the decomposed matrix.
- *
- * @return the inverse matrix.
- * @throws SingularMatrixException if the decomposed matrix is singular.
- */
- @Override
- public RealMatrix getInverse() {
- if (!isNonSingular()) {
- throw new SingularMatrixException();
- }
-
- final int m = realEigenvalues.length;
- final double[][] invData = new double[m][m];
-
- for (int i = 0; i < m; ++i) {
- final double[] invI = invData[i];
- for (int j = 0; j < m; ++j) {
- double invIJ = 0;
- for (int k = 0; k < m; ++k) {
- final double[] vK = eigenvectors[k].getDataRef();
- invIJ += vK[i] * vK[j] / realEigenvalues[k];
- }
- invI[j] = invIJ;
- }
- }
- return MatrixUtils.createRealMatrix(invData);
- }
- }
-
- /**
- * Transforms the matrix to tridiagonal form.
- *
- * @param matrix Matrix to transform.
- */
- private void transformToTridiagonal(final RealMatrix matrix) {
- // transform the matrix to tridiagonal
- transformer = new TriDiagonalTransformer(matrix);
- main = transformer.getMainDiagonalRef();
- secondary = transformer.getSecondaryDiagonalRef();
- }
-
- /**
- * Find eigenvalues and eigenvectors (Dubrulle et al., 1971)
- *
- * @param householderMatrix Householder matrix of the transformation
- * to tridiagonal form.
- */
- private void findEigenVectors(final double[][] householderMatrix) {
- final double[][]z = householderMatrix.clone();
- final int n = main.length;
- realEigenvalues = new double[n];
- imagEigenvalues = new double[n];
- final double[] e = new double[n];
- for (int i = 0; i < n - 1; i++) {
- realEigenvalues[i] = main[i];
- e[i] = secondary[i];
- }
- realEigenvalues[n - 1] = main[n - 1];
- e[n - 1] = 0;
-
- // Determine the largest main and secondary value in absolute term.
- double maxAbsoluteValue = 0;
- for (int i = 0; i < n; i++) {
- if (FastMath.abs(realEigenvalues[i]) > maxAbsoluteValue) {
- maxAbsoluteValue = FastMath.abs(realEigenvalues[i]);
- }
- if (FastMath.abs(e[i]) > maxAbsoluteValue) {
- maxAbsoluteValue = FastMath.abs(e[i]);
- }
- }
- // Make null any main and secondary value too small to be significant
- if (maxAbsoluteValue != 0) {
- for (int i=0; i < n; i++) {
- if (FastMath.abs(realEigenvalues[i]) <= Precision.EPSILON * maxAbsoluteValue) {
- realEigenvalues[i] = 0;
- }
- if (FastMath.abs(e[i]) <= Precision.EPSILON * maxAbsoluteValue) {
- e[i]=0;
- }
- }
- }
-
- for (int j = 0; j < n; j++) {
- int its = 0;
- int m;
- do {
- for (m = j; m < n - 1; m++) {
- double delta = FastMath.abs(realEigenvalues[m]) +
- FastMath.abs(realEigenvalues[m + 1]);
- if (FastMath.abs(e[m]) + delta == delta) {
- break;
- }
- }
- if (m != j) {
- if (its == MAX_ITER) {
- throw new MaxCountExceededException(LocalizedFormats.CONVERGENCE_FAILED,
- MAX_ITER);
- }
- its++;
- double q = (realEigenvalues[j + 1] - realEigenvalues[j]) / (2 * e[j]);
- double t = FastMath.sqrt(1 + q * q);
- if (q < 0.0) {
- q = realEigenvalues[m] - realEigenvalues[j] + e[j] / (q - t);
- } else {
- q = realEigenvalues[m] - realEigenvalues[j] + e[j] / (q + t);
- }
- double u = 0.0;
- double s = 1.0;
- double c = 1.0;
- int i;
- for (i = m - 1; i >= j; i--) {
- double p = s * e[i];
- double h = c * e[i];
- if (FastMath.abs(p) >= FastMath.abs(q)) {
- c = q / p;
- t = FastMath.sqrt(c * c + 1.0);
- e[i + 1] = p * t;
- s = 1.0 / t;
- c *= s;
- } else {
- s = p / q;
- t = FastMath.sqrt(s * s + 1.0);
- e[i + 1] = q * t;
- c = 1.0 / t;
- s *= c;
- }
- if (e[i + 1] == 0.0) {
- realEigenvalues[i + 1] -= u;
- e[m] = 0.0;
- break;
- }
- q = realEigenvalues[i + 1] - u;
- t = (realEigenvalues[i] - q) * s + 2.0 * c * h;
- u = s * t;
- realEigenvalues[i + 1] = q + u;
- q = c * t - h;
- for (int ia = 0; ia < n; ia++) {
- p = z[ia][i + 1];
- z[ia][i + 1] = s * z[ia][i] + c * p;
- z[ia][i] = c * z[ia][i] - s * p;
- }
- }
- if (t == 0.0 && i >= j) {
- continue;
- }
- realEigenvalues[j] -= u;
- e[j] = q;
- e[m] = 0.0;
- }
- } while (m != j);
- }
-
- //Sort the eigen values (and vectors) in increase order
- for (int i = 0; i < n; i++) {
- int k = i;
- double p = realEigenvalues[i];
- for (int j = i + 1; j < n; j++) {
- if (realEigenvalues[j] > p) {
- k = j;
- p = realEigenvalues[j];
- }
- }
- if (k != i) {
- realEigenvalues[k] = realEigenvalues[i];
- realEigenvalues[i] = p;
- for (int j = 0; j < n; j++) {
- p = z[j][i];
- z[j][i] = z[j][k];
- z[j][k] = p;
- }
- }
- }
-
- // Determine the largest eigen value in absolute term.
- maxAbsoluteValue = 0;
- for (int i = 0; i < n; i++) {
- if (FastMath.abs(realEigenvalues[i]) > maxAbsoluteValue) {
- maxAbsoluteValue=FastMath.abs(realEigenvalues[i]);
- }
- }
- // Make null any eigen value too small to be significant
- if (maxAbsoluteValue != 0.0) {
- for (int i=0; i < n; i++) {
- if (FastMath.abs(realEigenvalues[i]) < Precision.EPSILON * maxAbsoluteValue) {
- realEigenvalues[i] = 0;
- }
- }
- }
- eigenvectors = new ArrayRealVector[n];
- final double[] tmp = new double[n];
- for (int i = 0; i < n; i++) {
- for (int j = 0; j < n; j++) {
- tmp[j] = z[j][i];
- }
- eigenvectors[i] = new ArrayRealVector(tmp);
- }
- }
-
- /**
- * Transforms the matrix to Schur form and calculates the eigenvalues.
- *
- * @param matrix Matrix to transform.
- * @return the {@link SchurTransformer Shur transform} for this matrix
- */
- private SchurTransformer transformToSchur(final RealMatrix matrix) {
- final SchurTransformer schurTransform = new SchurTransformer(matrix);
- final double[][] matT = schurTransform.getT().getData();
-
- realEigenvalues = new double[matT.length];
- imagEigenvalues = new double[matT.length];
-
- for (int i = 0; i < realEigenvalues.length; i++) {
- if (i == (realEigenvalues.length - 1) ||
- Precision.equals(matT[i + 1][i], 0.0, EPSILON)) {
- realEigenvalues[i] = matT[i][i];
- } else {
- final double x = matT[i + 1][i + 1];
- final double p = 0.5 * (matT[i][i] - x);
- final double z = FastMath.sqrt(FastMath.abs(p * p + matT[i + 1][i] * matT[i][i + 1]));
- realEigenvalues[i] = x + p;
- imagEigenvalues[i] = z;
- realEigenvalues[i + 1] = x + p;
- imagEigenvalues[i + 1] = -z;
- i++;
- }
- }
- return schurTransform;
- }
-
- /**
- * Performs a division of two complex numbers.
- *
- * @param xr real part of the first number
- * @param xi imaginary part of the first number
- * @param yr real part of the second number
- * @param yi imaginary part of the second number
- * @return result of the complex division
- */
- private Complex cdiv(final double xr, final double xi,
- final double yr, final double yi) {
- return Complex.ofCartesian(xr, xi).divide(Complex.ofCartesian(yr, yi));
- }
-
- /**
- * Find eigenvectors from a matrix transformed to Schur form.
- *
- * @param schur the schur transformation of the matrix
- * @throws MathArithmeticException if the Schur form has a norm of zero
- */
- private void findEigenVectorsFromSchur(final SchurTransformer schur)
- throws MathArithmeticException {
- final double[][] matrixT = schur.getT().getData();
- final double[][] matrixP = schur.getP().getData();
-
- final int n = matrixT.length;
-
- // compute matrix norm
- double norm = 0.0;
- for (int i = 0; i < n; i++) {
- for (int j = FastMath.max(i - 1, 0); j < n; j++) {
- norm += FastMath.abs(matrixT[i][j]);
- }
- }
-
- // we can not handle a matrix with zero norm
- if (Precision.equals(norm, 0.0, EPSILON)) {
- throw new MathArithmeticException(LocalizedFormats.ZERO_NORM);
- }
-
- // Backsubstitute to find vectors of upper triangular form
-
- double r = 0.0;
- double s = 0.0;
- double z = 0.0;
-
- for (int idx = n - 1; idx >= 0; idx--) {
- double p = realEigenvalues[idx];
- double q = imagEigenvalues[idx];
-
- if (Precision.equals(q, 0.0)) {
- // Real vector
- int l = idx;
- matrixT[idx][idx] = 1.0;
- for (int i = idx - 1; i >= 0; i--) {
- double w = matrixT[i][i] - p;
- r = 0.0;
- for (int j = l; j <= idx; j++) {
- r += matrixT[i][j] * matrixT[j][idx];
- }
- if (Precision.compareTo(imagEigenvalues[i], 0.0, EPSILON) < 0) {
- z = w;
- s = r;
- } else {
- l = i;
- if (Precision.equals(imagEigenvalues[i], 0.0)) {
- if (w != 0.0) {
- matrixT[i][idx] = -r / w;
- } else {
- matrixT[i][idx] = -r / (Precision.EPSILON * norm);
- }
- } else {
- // Solve real equations
- double x = matrixT[i][i + 1];
- double y = matrixT[i + 1][i];
- q = (realEigenvalues[i] - p) * (realEigenvalues[i] - p) +
- imagEigenvalues[i] * imagEigenvalues[i];
- double t = (x * s - z * r) / q;
- matrixT[i][idx] = t;
- if (FastMath.abs(x) > FastMath.abs(z)) {
- matrixT[i + 1][idx] = (-r - w * t) / x;
- } else {
- matrixT[i + 1][idx] = (-s - y * t) / z;
- }
- }
-
- // Overflow control
- double t = FastMath.abs(matrixT[i][idx]);
- if ((Precision.EPSILON * t) * t > 1) {
- for (int j = i; j <= idx; j++) {
- matrixT[j][idx] /= t;
- }
- }
- }
- }
- } else if (q < 0.0) {
- // Complex vector
- int l = idx - 1;
-
- // Last vector component imaginary so matrix is triangular
- if (FastMath.abs(matrixT[idx][idx - 1]) > FastMath.abs(matrixT[idx - 1][idx])) {
- matrixT[idx - 1][idx - 1] = q / matrixT[idx][idx - 1];
- matrixT[idx - 1][idx] = -(matrixT[idx][idx] - p) / matrixT[idx][idx - 1];
- } else {
- final Complex result = cdiv(0.0, -matrixT[idx - 1][idx],
- matrixT[idx - 1][idx - 1] - p, q);
- matrixT[idx - 1][idx - 1] = result.getReal();
- matrixT[idx - 1][idx] = result.getImaginary();
- }
-
- matrixT[idx][idx - 1] = 0.0;
- matrixT[idx][idx] = 1.0;
-
- for (int i = idx - 2; i >= 0; i--) {
- double ra = 0.0;
- double sa = 0.0;
- for (int j = l; j <= idx; j++) {
- ra += matrixT[i][j] * matrixT[j][idx - 1];
- sa += matrixT[i][j] * matrixT[j][idx];
- }
- double w = matrixT[i][i] - p;
-
- if (Precision.compareTo(imagEigenvalues[i], 0.0, EPSILON) < 0) {
- z = w;
- r = ra;
- s = sa;
- } else {
- l = i;
- if (Precision.equals(imagEigenvalues[i], 0.0)) {
- final Complex c = cdiv(-ra, -sa, w, q);
- matrixT[i][idx - 1] = c.getReal();
- matrixT[i][idx] = c.getImaginary();
- } else {
- // Solve complex equations
- double x = matrixT[i][i + 1];
- double y = matrixT[i + 1][i];
- double vr = (realEigenvalues[i] - p) * (realEigenvalues[i] - p) +
- imagEigenvalues[i] * imagEigenvalues[i] - q * q;
- final double vi = (realEigenvalues[i] - p) * 2.0 * q;
- if (Precision.equals(vr, 0.0) && Precision.equals(vi, 0.0)) {
- vr = Precision.EPSILON * norm *
- (FastMath.abs(w) + FastMath.abs(q) + FastMath.abs(x) +
- FastMath.abs(y) + FastMath.abs(z));
- }
- final Complex c = cdiv(x * r - z * ra + q * sa,
- x * s - z * sa - q * ra, vr, vi);
- matrixT[i][idx - 1] = c.getReal();
- matrixT[i][idx] = c.getImaginary();
-
- if (FastMath.abs(x) > (FastMath.abs(z) + FastMath.abs(q))) {
- matrixT[i + 1][idx - 1] = (-ra - w * matrixT[i][idx - 1] +
- q * matrixT[i][idx]) / x;
- matrixT[i + 1][idx] = (-sa - w * matrixT[i][idx] -
- q * matrixT[i][idx - 1]) / x;
- } else {
- final Complex c2 = cdiv(-r - y * matrixT[i][idx - 1],
- -s - y * matrixT[i][idx], z, q);
- matrixT[i + 1][idx - 1] = c2.getReal();
- matrixT[i + 1][idx] = c2.getImaginary();
- }
- }
-
- // Overflow control
- double t = FastMath.max(FastMath.abs(matrixT[i][idx - 1]),
- FastMath.abs(matrixT[i][idx]));
- if ((Precision.EPSILON * t) * t > 1) {
- for (int j = i; j <= idx; j++) {
- matrixT[j][idx - 1] /= t;
- matrixT[j][idx] /= t;
- }
- }
- }
- }
- }
- }
-
- // Back transformation to get eigenvectors of original matrix
- for (int j = n - 1; j >= 0; j--) {
- for (int i = 0; i <= n - 1; i++) {
- z = 0.0;
- for (int k = 0; k <= FastMath.min(j, n - 1); k++) {
- z += matrixP[i][k] * matrixT[k][j];
- }
- matrixP[i][j] = z;
- }
- }
-
- eigenvectors = new ArrayRealVector[n];
- final double[] tmp = new double[n];
- for (int i = 0; i < n; i++) {
- System.out.println("Eigenvector " + i + ": "); // XXX
- for (int j = 0; j < n; j++) {
- tmp[j] = matrixP[j][i];
- System.out.print(tmp[j] + "\t"); // XXX
- }
- System.out.println(); // XXX
- eigenvectors[i] = new ArrayRealVector(tmp);
- }
- }
-}
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/SimplexOptimizer.java.DEBUG b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/SimplexOptimizer.java.DEBUG
deleted file mode 100644
index c6b4420..0000000
--- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/SimplexOptimizer.java.DEBUG
+++ /dev/null
@@ -1,339 +0,0 @@
-/*
- * 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.commons.math4.optim.nonlinear.scalar.noderiv;
-
-import java.util.Comparator;
-
-import org.apache.commons.rng.UniformRandomProvider;
-import org.apache.commons.rng.simple.RandomSource;
-import org.apache.commons.math4.analysis.MultivariateFunction;
-import org.apache.commons.math4.exception.MathUnsupportedOperationException;
-import org.apache.commons.math4.exception.NullArgumentException;
-import org.apache.commons.math4.exception.util.LocalizedFormats;
-import org.apache.commons.math4.optim.ConvergenceChecker;
-import org.apache.commons.math4.optim.OptimizationData;
-import org.apache.commons.math4.optim.PointValuePair;
-import org.apache.commons.math4.optim.SimpleValueChecker;
-import org.apache.commons.math4.optim.nonlinear.scalar.GoalType;
-import org.apache.commons.math4.optim.nonlinear.scalar.SimulatedAnnealing;
-import org.apache.commons.math4.optim.nonlinear.scalar.MultivariateOptimizer;
-
-/**
- * This class implements simplex-based direct search optimization.
- *
- * <p>
- * Direct search methods only use objective function values, they do
- * not need derivatives and don't either try to compute approximation
- * of the derivatives. According to a 1996 paper by Margaret H. Wright
- * (<a href="http://cm.bell-labs.com/cm/cs/doc/96/4-02.ps.gz">Direct
- * Search Methods: Once Scorned, Now Respectable</a>), they are used
- * when either the computation of the derivative is impossible (noisy
- * functions, unpredictable discontinuities) or difficult (complexity,
- * computation cost). In the first cases, rather than an optimum, a
- * <em>not too bad</em> point is desired. In the latter cases, an
- * optimum is desired but cannot be reasonably found. In all cases
- * direct search methods can be useful.
- * </p>
- * <p>
- * Simplex-based direct search methods are based on comparison of
- * the objective function values at the vertices of a simplex (which is a
- * set of n+1 points in dimension n) that is updated by the algorithms
- * steps.
- * </p>
- * <p>
- * The simplex update procedure ({@link NelderMeadSimplex} or
- * {@link MultiDirectionalSimplex}) must be passed to the
- * {@code optimize} method.
- * </p>
- * <p>
- * Each call to {@code optimize} will re-use the start configuration of
- * the current simplex and move it such that its first vertex is at the
- * provided start point of the optimization.
- * If the {@code optimize} method is called to solve a different problem
- * and the number of parameters change, the simplex must be re-initialized
- * to one with the appropriate dimensions.
- * </p>
- * <p>
- * Convergence is checked by providing the <em>worst</em> points of
- * previous and current simplex to the convergence checker, not the best
- * ones.
- * </p>
- * <p>
- * This simplex optimizer implementation does not directly support constrained
- * optimization with simple bounds; so, for such optimizations, either a more
- * dedicated algorithm must be used like
- * {@link CMAESOptimizer} or {@link BOBYQAOptimizer}, or the objective
- * function must be wrapped in an adapter like
- * {@link org.apache.commons.math4.optim.nonlinear.scalar.MultivariateFunctionMappingAdapter
- * MultivariateFunctionMappingAdapter} or
- * {@link org.apache.commons.math4.optim.nonlinear.scalar.MultivariateFunctionPenaltyAdapter
- * MultivariateFunctionPenaltyAdapter}.
- * <br>
- * The call to {@link #optimize(OptimizationData[]) optimize} will throw
- * {@link MathUnsupportedOperationException} if bounds are passed to it.
- * </p>
- *
- * @since 3.0
- */
-public class SimplexOptimizer extends MultivariateOptimizer {
- /** Simplex update rule. */
- private AbstractSimplex simplex;
- /** Simulated annealing setup. */
- private SimulatedAnnealing annealing;
- /** Overall best. */
- private PointValuePair best;
-
- /**
- * @param checker Convergence checker.
- */
- public SimplexOptimizer(ConvergenceChecker<PointValuePair> checker) {
- super(checker);
- }
-
- /**
- * @param rel Relative threshold.
- * @param abs Absolute threshold.
- */
- public SimplexOptimizer(double rel, double abs) {
- this(new SimpleValueChecker(rel, abs));
- }
-
- /**
- * {@inheritDoc}
- *
- * @param optData Optimization data. In addition to those documented in
- * {@link MultivariateOptimizer#parseOptimizationData(OptimizationData[])
- * MultivariateOptimizer}, this method will register the following data:
- * <ul>
- * <li>{@link AbstractSimplex}</li>
- * <li>{@link SimulatedAnnealing}</li>
- * </ul>
- * @return {@inheritDoc}
- */
- @Override
- public PointValuePair optimize(OptimizationData... optData) {
- // Set up base class and perform computation.
- return super.optimize(optData);
- }
-
- /** {@inheritDoc} */
- @Override
- protected PointValuePair doOptimize() {
- checkParameters();
-
- // Indirect call to "computeObjectiveValue" in order to update the
- // evaluations counter.
- final MultivariateFunction evalFunc
- = new MultivariateFunction() {
- /** {@inheritDoc} */
- @Override
- public double value(double[] point) {
- return computeObjectiveValue(point);
- }
- };
-
- final boolean isMinim = getGoalType() == GoalType.MINIMIZE;
- final Comparator<PointValuePair> comparator
- = new Comparator<PointValuePair>() {
- /** {@inheritDoc} */
- @Override
- public int compare(final PointValuePair o1,
- final PointValuePair o2) {
- final double v1 = o1.getValue();
- final double v2 = o2.getValue();
- return isMinim ? Double.compare(v1, v2) : Double.compare(v2, v1);
- }
- };
-
- // Initialize search.
- simplex.build(getStartPoint());
- simplex.evaluate(evalFunc, comparator);
- final UniformRandomProvider rng = annealing != null ?
- RandomSource.create(RandomSource.KISS) :
- null;
-
- PointValuePair[] previous = null;
- int iteration = 0;
- final ConvergenceChecker<PointValuePair> checker = getConvergenceChecker();
- while (true) {
- iteration = getIterations();
- if (iteration > 0) {
- boolean converged = true;
- for (int i = 0; i < simplex.getSize(); i++) {
- PointValuePair prev = previous[i];
- converged = converged &&
- checker.converged(iteration, prev, simplex.getPoint(i));
-
- if (!converged) {
- // Short circuit, since "converged" will stay "false".
- break;
- }
- }
- if (converged) {
- System.out.println(" saAcceptCount=" + saAcceptCount); // XXX
- // We have found an optimum.
- return best;
- }
- }
-
- // We still need to search.
- previous = simplex.getPoints();
- simplex.iterate(evalFunc, comparator);
-
- // Track best point.
- final int bestIndex = 0; // Index of best point.
- if (best == null ||
- comparator.compare(best, simplex.getPoint(bestIndex)) > 0) {
- best = simplex.getPoint(bestIndex);
- }
-
- if (annealing != null) {
- // Simulated annealing step.
- simulatedAnnealing(iteration,
- evalFunc,
- isMinim,
- rng);
- }
-
- incrementIterationCount();
- }
- }
-
- /**
- * Scans the list of (required and optional) optimization data that
- * characterize the problem.
- *
- * @param optData Optimization data.
- * The following data will be looked for:
- * <ul>
- * <li>{@link AbstractSimplex}</li>
- * <li>{@link SimulatedAnnealing}</li>
- * </ul>
- */
- @Override
- protected void parseOptimizationData(OptimizationData... optData) {
- // Allow base class to register its own data.
- super.parseOptimizationData(optData);
-
- // The existing values (as set by the previous call) are reused if
- // not provided in the argument list.
- for (OptimizationData data : optData) {
- if (data instanceof AbstractSimplex) {
- simplex = (AbstractSimplex) data;
- continue;
- }
- if (data instanceof SimulatedAnnealing) {
- annealing = (SimulatedAnnealing) data;
- continue;
- }
- }
- }
-
- /**
- * @throws MathUnsupportedOperationException if bounds were passed to the
- * {@link #optimize(OptimizationData[]) optimize} method.
- * @throws NullArgumentException if no initial simplex was passed to the
- * {@link #optimize(OptimizationData[]) optimize} method.
- */
- private void checkParameters() {
- if (simplex == null) {
- throw new NullArgumentException();
- }
- if (getLowerBound() != null ||
- getUpperBound() != null) {
- throw new MathUnsupportedOperationException(LocalizedFormats.CONSTRAINT);
- }
- }
-
- private int saAcceptCount = 0; // XXX
- /**
- * Perform the annealing step (possibly replacing simplex's best point).
- *
- * @param iteration Current iteration.
- * @param evalFunc Evaluation function.
- * @param isMinim Whether a minimization is performed.
- * @param rng RNG.
- */
- private void simulatedAnnealing(int iteration,
- MultivariateFunction evalFunc,
- boolean isMinim,
- UniformRandomProvider rng) {
- if (iteration > annealing.getIterations()) {
- return; // Do nothing.
- }
-
- // Construct alternative state.
- final int bestIndex = 0; // Index of best point.
- final PointValuePair alt = alternativeState(simplex,
- bestIndex,
- evalFunc,
- rng);
-
- if (annealing.accept(simplex.getPoint(bestIndex).getValue(),
- alt.getValue(),
- isMinim,
- iteration)) {
- System.out.println("eO=" + simplex.getPoint(bestIndex).getValue()); // XXX
- System.out.println("eN=" + alt.getValue()); // XXX
- ++saAcceptCount; // XXX
-
- // Modify best point of the current simplex.
- simplex.setPoint(bestIndex, alt);
- }
- }
-
- /**
- * Creates a state that could replace the one stored at
- * {@code replaceIndex} in the given {@code simplex}.
- *
- * @param simplex Current simplex.
- * @param replaceIndex Index of the simplex point that will potentially be
- * replace by the new state.
- * @param evalFunc Evaluation function.
- * @param rng RNG.
- * @return a new state.
- */
- private static PointValuePair alternativeState(AbstractSimplex simplex,
- int replaceIndex,
- MultivariateFunction evalFunc,
- UniformRandomProvider rng) {
- final PointValuePair[] points = simplex.getPoints();
- final int numPoints = points.length;
- final int spaceDim = numPoints - 1;
-
- // Compute mean coordinate offsets from the point to replace
- // to all the other points.
- final double[] coord = new double[spaceDim];
- final double[] replaceCoord = points[replaceIndex].getPointRef();
- for (int j = 0; j < numPoints; j++) {
- if (j == replaceIndex) {
- continue;
- }
- final double[] c = points[j].getPointRef();
- for (int i = 0; i < spaceDim; i++) {
- coord[i] += c[i] - replaceCoord[i];
- }
- }
-
- for (int i = 0; i < spaceDim; i++) {
- coord[i] /= numPoints; // Mean coordinate offset.
- coord[i] = replaceCoord[i] + (rng.nextDouble() - 0.5) * coord[i];
- }
-
- return new PointValuePair(coord, evalFunc.value(coord), false);
- }
-}