diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/analysis/interpolation/UnivariatePeriodicInterpolator.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/analysis/interpolation/UnivariatePeriodicInterpolator.java
index ef3122a138..942d055563 100644
--- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/analysis/interpolation/UnivariatePeriodicInterpolator.java
+++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/analysis/interpolation/UnivariatePeriodicInterpolator.java
@@ -114,7 +114,7 @@ public UnivariateFunction interpolate(double[] xval,
             y[index] = yval[i];
         }
 
-        SortInPlace.ASCENDING.accept(x, y);
+        SortInPlace.ASCENDING.apply(x, y);
 
         final UnivariateFunction f = interpolator.interpolate(x, y);
         return new UnivariateFunction() {
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/analysis/polynomials/PolynomialFunctionLagrangeForm.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/analysis/polynomials/PolynomialFunctionLagrangeForm.java
index c0fe244b84..14f94b7c69 100644
--- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/analysis/polynomials/PolynomialFunctionLagrangeForm.java
+++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/analysis/polynomials/PolynomialFunctionLagrangeForm.java
@@ -78,7 +78,7 @@ public PolynomialFunctionLagrangeForm(double x[], double y[])
         coefficientsComputed = false;
 
         if (!verifyInterpolationArray(x, y, false)) {
-            SortInPlace.ASCENDING.accept(this.x, this.y);
+            SortInPlace.ASCENDING.apply(this.x, this.y);
             // Second check in case some abscissa is duplicated.
             verifyInterpolationArray(this.x, this.y, true);
         }
@@ -183,7 +183,7 @@ public static double evaluate(double x[], double y[], double z)
         System.arraycopy(x, 0, xNew, 0, x.length);
         System.arraycopy(y, 0, yNew, 0, y.length);
 
-        SortInPlace.ASCENDING.accept(xNew, yNew);
+        SortInPlace.ASCENDING.apply(xNew, yNew);
         // Second check in case some abscissa is duplicated.
         verifyInterpolationArray(xNew, yNew, true);
         return evaluateInternal(xNew, yNew, z);
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/analysis/polynomials/PolynomialsUtils.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/analysis/polynomials/PolynomialsUtils.java
index 44da4c1451..6430d3140f 100644
--- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/analysis/polynomials/PolynomialsUtils.java
+++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/analysis/polynomials/PolynomialsUtils.java
@@ -366,7 +366,9 @@ public static double[] shift(final double[] coefficients,
     private static PolynomialFunction buildPolynomial(final int degree,
                                                       final List<BigFraction> coefficients,
                                                       final RecurrenceCoefficientsGenerator generator) {
-
+        // Synchronizing on a method parameter is not safe; however, in this
+        // case, the lock object is an immutable field that belongs to this
+        // class.
         synchronized (coefficients) {
             final int maxDegree = (int) AccurateMath.floor(AccurateMath.sqrt(2 * coefficients.size())) - 1;
             if (degree > maxDegree) {
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/distribution/EmpiricalDistribution.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/distribution/EmpiricalDistribution.java
index c8f5021533..65e580883c 100644
--- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/distribution/EmpiricalDistribution.java
+++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/distribution/EmpiricalDistribution.java
@@ -46,40 +46,46 @@
 
 /**
  * <p>Represents an <a href="http://en.wikipedia.org/wiki/Empirical_distribution_function">
- * empirical probability distribution</a> -- a probability distribution derived
+ * empirical probability distribution</a>: Probability distribution derived
  * from observed data without making any assumptions about the functional form
  * of the population distribution that the data come from.</p>
  *
- * <p>An <code>EmpiricalDistribution</code> maintains data structures, called
- * <i>distribution digests</i>, that describe empirical distributions and
- * support the following operations: <ul>
- * <li>loading the distribution from a file of observed data values</li>
- * <li>dividing the input data into "bin ranges" and reporting bin frequency
- *     counts (data for histogram)</li>
- * <li>reporting univariate statistics describing the full set of data values
- *     as well as the observations within each bin</li>
- * <li>generating random values from the distribution</li>
+ * <p>An {@code EmpiricalDistribution} maintains data structures called
+ * <i>distribution digests</i> that describe empirical distributions and
+ * support the following operations:
+ * <ul>
+ *  <li>loading the distribution from a file of observed data values</li>
+ *  <li>dividing the input data into "bin ranges" and reporting bin frequency
+ *      counts (data for histogram)</li>
+ *  <li>reporting univariate statistics describing the full set of data values
+ *      as well as the observations within each bin</li>
+ *  <li>generating random values from the distribution</li>
  * </ul>
- * Applications can use <code>EmpiricalDistribution</code> to build grouped
+ *
+ * Applications can use {@code EmpiricalDistribution} to build grouped
  * frequency histograms representing the input data or to generate random values
- * "like" those in the input file -- i.e., the values generated will follow the
+ * "like" those in the input file, i.e. the values generated will follow the
  * distribution of the values in the file.
  *
  * <p>The implementation uses what amounts to the
  * <a href="http://nedwww.ipac.caltech.edu/level5/March02/Silverman/Silver2_6.html">
  * Variable Kernel Method</a> with Gaussian smoothing:<p>
  * <strong>Digesting the input file</strong>
- * <ol><li>Pass the file once to compute min and max.</li>
- * <li>Divide the range from min-max into <code>binCount</code> "bins."</li>
- * <li>Pass the data file again, computing bin counts and univariate
- *     statistics (mean, std dev.) for each of the bins </li>
- * <li>Divide the interval (0,1) into subintervals associated with the bins,
- *     with the length of a bin's subinterval proportional to its count.</li></ol>
- * <strong>Generating random values from the distribution</strong><ol>
- * <li>Generate a uniformly distributed value in (0,1) </li>
- * <li>Select the subinterval to which the value belongs.
- * <li>Generate a random Gaussian value with mean = mean of the associated
- *     bin and std dev = std dev of associated bin.</li></ol>
+ * <ol>
+ *  <li>Pass the file once to compute min and max.</li>
+ *  <li>Divide the range from min to max into {@code binCount} bins.</li>
+ *  <li>Pass the data file again, computing bin counts and univariate
+ *      statistics (mean and std dev.) for each bin.</li>
+ *  <li>Divide the interval (0,1) into subintervals associated with the bins,
+ *      with the length of a bin's subinterval proportional to its count.</li>
+ * </ol>
+ * <strong>Generating random values from the distribution</strong>
+ * <ol>
+ *  <li>Generate a uniformly distributed value in (0,1) </li>
+ *  <li>Select the subinterval to which the value belongs.
+ *  <li>Generate a random Gaussian value with mean = mean of the associated
+ *      bin and std dev = std dev of associated bin.</li>
+ * </ol>
  *
  * <p>EmpiricalDistribution implements the {@link ContinuousDistribution} interface
  * as follows.  Given x within the range of values in the dataset, let B
@@ -91,49 +97,36 @@
  * grouped frequency distribution at the bin endpoints and interpolates within
  * bins using within-bin kernels.</p>
  *
- *<strong>USAGE NOTES:</strong><ul>
- *<li>The <code>binCount</code> is set by default to 1000.  A good rule of thumb
- *    is to set the bin count to approximately the length of the input file divided
- *    by 10. </li>
- *<li>The input file <i>must</i> be a plain text file containing one valid numeric
- *    entry per line.</li>
+ * <strong>USAGE NOTES:</strong>
+ * <ul>
+ * <li>The {@code binCount} is set by default to 1000.  A good rule of thumb
+ *     is to set the bin count to approximately the length of the input file divided
+ *     by 10. </li>
+ * <li>The input file <i>must</i> be a plain text file containing one valid numeric
+ *     entry per line.</li>
  * </ul>
- *
  */
 public class EmpiricalDistribution extends AbstractRealDistribution
     implements ContinuousDistribution {
-
     /** Default bin count. */
     public static final int DEFAULT_BIN_COUNT = 1000;
-
     /** Character set for file input. */
     private static final String FILE_CHARSET = "US-ASCII";
-
-    /** Serializable version identifier. */
-    private static final long serialVersionUID = 5729073523949762654L;
-
-    /** List of SummaryStatistics objects characterizing the bins. */
+    /** Bins' characteristics. */
     private final List<SummaryStatistics> binStats;
-
     /** Sample statistics. */
     private SummaryStatistics sampleStats;
-
     /** Max loaded value. */
     private double max = Double.NEGATIVE_INFINITY;
-
     /** Min loaded value. */
     private double min = Double.POSITIVE_INFINITY;
-
     /** Grid size. */
     private double delta;
-
-    /** number of bins. */
+    /** Number of bins. */
     private final int binCount;
-
-    /** is the distribution loaded? */
+    /** Whether the distribution is loaded. */
     private boolean loaded;
-
-    /** upper bounds of subintervals in (0,1) "belonging" to the bins. */
+    /** Upper bounds of subintervals in (0,1) belonging to the bins. */
     private double[] upperBounds;
 
     /**
@@ -247,11 +240,10 @@ public void load(File file) throws IOException {
     }
 
     /**
-     * Provides methods for computing <code>sampleStats</code> and
-     * <code>beanStats</code> abstracting the source of data.
+     * Provides methods for computing {@code sampleStats} and
+     * {@code beanStats} abstracting the source of data.
      */
-    private abstract class DataAdapter{
-
+    private abstract class DataAdapter {
         /**
          * Compute bin stats.
          *
@@ -265,16 +257,14 @@ private abstract class DataAdapter{
          * @throws IOException if an error occurs computing sample stats
          */
         public abstract void computeStats() throws IOException;
-
     }
 
     /**
-     * <code>DataAdapter</code> for data provided through some input stream.
+     * {@code DataAdapter} for data provided through some input stream.
      */
-    private class StreamDataAdapter extends DataAdapter{
-
+    private class StreamDataAdapter extends DataAdapter {
         /** Input stream providing access to the data. */
-        private BufferedReader inputStream;
+        private final BufferedReader inputStream;
 
         /**
          * Create a StreamDataAdapter from a BufferedReader.
@@ -282,7 +272,6 @@ private class StreamDataAdapter extends DataAdapter{
          * @param in BufferedReader input stream
          */
         StreamDataAdapter(BufferedReader in){
-            super();
             inputStream = in;
         }
 
@@ -298,7 +287,6 @@ public void computeBinStats() throws IOException {
             }
 
             inputStream.close();
-            inputStream = null;
         }
 
         /** {@inheritDoc} */
@@ -312,15 +300,13 @@ public void computeStats() throws IOException {
                 sampleStats.addValue(val);
             }
             inputStream.close();
-            inputStream = null;
         }
     }
 
     /**
-     * <code>DataAdapter</code> for data provided as array of doubles.
+     * {@code DataAdapter} for data provided as array of doubles.
      */
     private class ArrayDataAdapter extends DataAdapter {
-
         /** Array of input  data values. */
         private final double[] inputArray;
 
@@ -331,7 +317,6 @@ private class ArrayDataAdapter extends DataAdapter {
          * @throws NullArgumentException if in is null
          */
         ArrayDataAdapter(double[] in) {
-            super();
             NullArgumentException.check(in);
             inputArray = in;
         }
@@ -349,8 +334,7 @@ public void computeStats() throws IOException {
         @Override
         public void computeBinStats() throws IOException {
             for (int i = 0; i < inputArray.length; i++) {
-                SummaryStatistics stats =
-                    binStats.get(findBin(inputArray[i]));
+                SummaryStatistics stats = binStats.get(findBin(inputArray[i]));
                 stats.addValue(inputArray[i]);
             }
         }
@@ -362,12 +346,11 @@ public void computeBinStats() throws IOException {
      * @param da object providing access to the data
      * @throws IOException  if an IO error occurs
      */
-    private void fillBinStats(final DataAdapter da)
-        throws IOException {
+    private void fillBinStats(final DataAdapter da) throws IOException {
         // Set up grid
         min = sampleStats.getMin();
         max = sampleStats.getMax();
-        delta = (max - min)/binCount;
+        delta = (max - min) / binCount;
 
         // Initialize binStats ArrayList
         if (!binStats.isEmpty()) {
@@ -375,7 +358,7 @@ private void fillBinStats(final DataAdapter da)
         }
         for (int i = 0; i < binCount; i++) {
             SummaryStatistics stats = new SummaryStatistics();
-            binStats.add(i,stats);
+            binStats.add(i, stats);
         }
 
         // Filling data in binStats Array
@@ -383,13 +366,12 @@ private void fillBinStats(final DataAdapter da)
 
         // Assign upperBounds based on bin counts
         upperBounds = new double[binCount];
-        upperBounds[0] =
-        ((double) binStats.get(0).getN()) / (double) sampleStats.getN();
-        for (int i = 1; i < binCount-1; i++) {
-            upperBounds[i] = upperBounds[i-1] +
-            ((double) binStats.get(i).getN()) / (double) sampleStats.getN();
+        upperBounds[0] = binStats.get(0).getN() / (double) sampleStats.getN();
+        for (int i = 1; i < binCount - 1; i++) {
+            upperBounds[i] = upperBounds[i - 1] +
+                binStats.get(i).getN() / (double) sampleStats.getN();
         }
-        upperBounds[binCount-1] = 1.0d;
+        upperBounds[binCount - 1] = 1d;
     }
 
     /**
@@ -399,9 +381,8 @@ private void fillBinStats(final DataAdapter da)
      * @return the index of the bin containing the value
      */
     private int findBin(double value) {
-        return AccurateMath.min(
-                AccurateMath.max((int) AccurateMath.ceil((value - min) / delta) - 1, 0),
-                binCount - 1);
+        return AccurateMath.min(AccurateMath.max((int) AccurateMath.ceil((value - min) / delta) - 1, 0),
+                                binCount - 1);
     }
 
     /**
@@ -490,7 +471,7 @@ public boolean isLoaded() {
         return loaded;
     }
 
-    // Distribution methods ---------------------------
+    // Distribution methods.
 
     /**
      * {@inheritDoc}
@@ -588,21 +569,22 @@ public double cumulativeProbability(double x) {
      */
     @Override
     public double inverseCumulativeProbability(final double p) {
-        if (p < 0.0 || p > 1.0) {
+        if (p < 0 ||
+            p > 1) {
             throw new OutOfRangeException(p, 0, 1);
         }
 
-        if (p == 0.0) {
+        if (p == 0) {
             return getSupportLowerBound();
         }
 
-        if (p == 1.0) {
+        if (p == 1) {
             return getSupportUpperBound();
         }
 
         int i = 0;
         while (cumBinP(i) < p) {
-            i++;
+            ++i;
         }
 
         final ContinuousDistribution kernel = getKernel(binStats.get(i));
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/ml/clustering/Cluster.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/ml/clustering/Cluster.java
index e336b80270..5aedc6b36e 100644
--- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/ml/clustering/Cluster.java
+++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/ml/clustering/Cluster.java
@@ -17,7 +17,6 @@
 
 package org.apache.commons.math4.legacy.ml.clustering;
 
-import java.io.Serializable;
 import java.util.ArrayList;
 import java.util.List;
 
@@ -26,10 +25,7 @@
  * @param <T> the type of points that can be clustered
  * @since 3.2
  */
-public class Cluster<T extends Clusterable> implements Serializable {
-
-    /** Serializable version identifier. */
-    private static final long serialVersionUID = -3442297081515880464L;
+public class Cluster<T extends Clusterable> {
 
     /** The points contained in this cluster. */
     private final List<T> points;
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/ml/distance/DistanceMeasure.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/ml/distance/DistanceMeasure.java
index c8fabc2d73..a5a22f15a8 100644
--- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/ml/distance/DistanceMeasure.java
+++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/ml/distance/DistanceMeasure.java
@@ -16,8 +16,6 @@
  */
 package org.apache.commons.math4.legacy.ml.distance;
 
-import java.io.Serializable;
-
 import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
 
 /**
@@ -25,7 +23,7 @@
  *
  * @since 3.2
  */
-public interface DistanceMeasure extends Serializable {
+public interface DistanceMeasure {
 
     /**
      * Compute the distance between two n-dimensional vectors.
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/AbstractSimplex.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/AbstractSimplex.java
deleted file mode 100644
index bd1d555c9d..0000000000
--- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/AbstractSimplex.java
+++ /dev/null
@@ -1,345 +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.legacy.optim.nonlinear.scalar.noderiv;
-
-import java.util.Arrays;
-import java.util.Comparator;
-
-import org.apache.commons.math4.legacy.analysis.MultivariateFunction;
-import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
-import org.apache.commons.math4.legacy.exception.MathIllegalArgumentException;
-import org.apache.commons.math4.legacy.exception.NotStrictlyPositiveException;
-import org.apache.commons.math4.legacy.exception.NullArgumentException;
-import org.apache.commons.math4.legacy.exception.OutOfRangeException;
-import org.apache.commons.math4.legacy.exception.ZeroException;
-import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
-import org.apache.commons.math4.legacy.optim.OptimizationData;
-import org.apache.commons.math4.legacy.optim.PointValuePair;
-
-/**
- * This class implements the simplex concept.
- * It is intended to be used in conjunction with {@link SimplexOptimizer}.
- * <br>
- * The initial configuration of the simplex is set by the constructors
- * {@link #AbstractSimplex(double[])} or {@link #AbstractSimplex(double[][])}.
- * The other {@link #AbstractSimplex(int) constructor} will set all steps
- * to 1, thus building a default configuration from a unit hypercube.
- * <br>
- * Users <em>must</em> call the {@link #build(double[]) build} method in order
- * to create the data structure that will be acted on by the other methods of
- * this class.
- *
- * @see SimplexOptimizer
- * @since 3.0
- */
-public abstract class AbstractSimplex implements OptimizationData {
-    /** Simplex. */
-    private PointValuePair[] simplex;
-    /** Start simplex configuration. */
-    private double[][] startConfiguration;
-    /** Simplex dimension (must be equal to {@code simplex.length - 1}). */
-    private final int dimension;
-
-    /**
-     * Build a unit hypercube simplex.
-     *
-     * @param n Dimension of the simplex.
-     */
-    protected AbstractSimplex(int n) {
-        this(n, 1d);
-    }
-
-    /**
-     * Build a hypercube simplex with the given side length.
-     *
-     * @param n Dimension of the simplex.
-     * @param sideLength Length of the sides of the hypercube.
-     */
-    protected AbstractSimplex(int n,
-                              double sideLength) {
-        this(createHypercubeSteps(n, sideLength));
-    }
-
-    /**
-     * The start configuration for simplex is built from a box parallel to
-     * the canonical axes of the space. The simplex is the subset of vertices
-     * of a box parallel to the canonical axes. It is built as the path followed
-     * while traveling from one vertex of the box to the diagonally opposite
-     * vertex moving only along the box edges. The first vertex of the box will
-     * be located at the start point of the optimization.
-     * As an example, in dimension 3 a simplex has 4 vertices. Setting the
-     * steps to (1, 10, 2) and the start point to (1, 1, 1) would imply the
-     * start simplex would be: { (1, 1, 1), (2, 1, 1), (2, 11, 1), (2, 11, 3) }.
-     * The first vertex would be set to the start point at (1, 1, 1) and the
-     * last vertex would be set to the diagonally opposite vertex at (2, 11, 3).
-     *
-     * @param steps Steps along the canonical axes representing box edges. They
-     * may be negative but not zero.
-     * @throws NullArgumentException if {@code steps} is {@code null}.
-     * @throws ZeroException if one of the steps is zero.
-     */
-    protected AbstractSimplex(final double[] steps) {
-        if (steps == null) {
-            throw new NullArgumentException();
-        }
-        if (steps.length == 0) {
-            throw new ZeroException();
-        }
-        dimension = steps.length;
-
-        // Only the relative position of the n final vertices with respect
-        // to the first one are stored.
-        startConfiguration = new double[dimension][dimension];
-        for (int i = 0; i < dimension; i++) {
-            final double[] vertexI = startConfiguration[i];
-            for (int j = 0; j < i + 1; j++) {
-                if (steps[j] == 0) {
-                    throw new ZeroException(LocalizedFormats.EQUAL_VERTICES_IN_SIMPLEX);
-                }
-                System.arraycopy(steps, 0, vertexI, 0, j + 1);
-            }
-        }
-    }
-
-    /**
-     * The real initial simplex will be set up by moving the reference
-     * simplex such that its first point is located at the start point of the
-     * optimization.
-     *
-     * @param referenceSimplex Reference simplex.
-     * @throws NotStrictlyPositiveException if the reference simplex does not
-     * contain at least one point.
-     * @throws DimensionMismatchException if there is a dimension mismatch
-     * in the reference simplex.
-     * @throws IllegalArgumentException if one of its vertices is duplicated.
-     */
-    protected AbstractSimplex(final double[][] referenceSimplex) {
-        if (referenceSimplex.length <= 0) {
-            throw new NotStrictlyPositiveException(LocalizedFormats.SIMPLEX_NEED_ONE_POINT,
-                                                   referenceSimplex.length);
-        }
-        dimension = referenceSimplex.length - 1;
-
-        // Only the relative position of the n final vertices with respect
-        // to the first one are stored.
-        startConfiguration = new double[dimension][dimension];
-        final double[] ref0 = referenceSimplex[0];
-
-        // Loop over vertices.
-        for (int i = 0; i < referenceSimplex.length; i++) {
-            final double[] refI = referenceSimplex[i];
-
-            // Safety checks.
-            if (refI.length != dimension) {
-                throw new DimensionMismatchException(refI.length, dimension);
-            }
-            for (int j = 0; j < i; j++) {
-                final double[] refJ = referenceSimplex[j];
-                boolean allEquals = true;
-                for (int k = 0; k < dimension; k++) {
-                    if (refI[k] != refJ[k]) {
-                        allEquals = false;
-                        break;
-                    }
-                }
-                if (allEquals) {
-                    throw new MathIllegalArgumentException(LocalizedFormats.EQUAL_VERTICES_IN_SIMPLEX,
-                                                           i, j);
-                }
-            }
-
-            // Store vertex i position relative to vertex 0 position.
-            if (i > 0) {
-                final double[] confI = startConfiguration[i - 1];
-                for (int k = 0; k < dimension; k++) {
-                    confI[k] = refI[k] - ref0[k];
-                }
-            }
-        }
-    }
-
-    /**
-     * Get simplex dimension.
-     *
-     * @return the dimension of the simplex.
-     */
-    public int getDimension() {
-        return dimension;
-    }
-
-    /**
-     * Get simplex size.
-     * After calling the {@link #build(double[]) build} method, this method will
-     * will be equivalent to {@code getDimension() + 1}.
-     *
-     * @return the size of the simplex.
-     */
-    public int getSize() {
-        return simplex.length;
-    }
-
-    /**
-     * Compute the next simplex of the algorithm.
-     *
-     * @param evaluationFunction Evaluation function.
-     * @param comparator Comparator to use to sort simplex vertices from best
-     * to worst.
-     * @throws org.apache.commons.math4.legacy.exception.TooManyEvaluationsException
-     * if the algorithm fails to converge.
-     */
-    public abstract void iterate(MultivariateFunction evaluationFunction,
-                                 Comparator<PointValuePair> comparator);
-
-    /**
-     * Build an initial simplex.
-     *
-     * @param startPoint First point of the simplex.
-     * @throws DimensionMismatchException if the start point does not match
-     * simplex dimension.
-     */
-    public void build(final double[] startPoint) {
-        if (dimension != startPoint.length) {
-            throw new DimensionMismatchException(dimension, startPoint.length);
-        }
-
-        // Set first vertex.
-        simplex = new PointValuePair[dimension + 1];
-        simplex[0] = new PointValuePair(startPoint, Double.NaN);
-
-        // Set remaining vertices.
-        for (int i = 0; i < dimension; i++) {
-            final double[] confI = startConfiguration[i];
-            final double[] vertexI = new double[dimension];
-            for (int k = 0; k < dimension; k++) {
-                vertexI[k] = startPoint[k] + confI[k];
-            }
-            simplex[i + 1] = new PointValuePair(vertexI, Double.NaN);
-        }
-    }
-
-    /**
-     * Evaluate all the non-evaluated points of the simplex.
-     *
-     * @param evaluationFunction Evaluation function.
-     * @param comparator Comparator to use to sort simplex vertices from best to worst.
-     * @throws org.apache.commons.math4.legacy.exception.TooManyEvaluationsException
-     * if the maximal number of evaluations is exceeded.
-     */
-    public void evaluate(final MultivariateFunction evaluationFunction,
-                         final Comparator<PointValuePair> comparator) {
-        // Evaluate the objective function at all non-evaluated simplex points.
-        for (int i = 0; i < simplex.length; i++) {
-            final PointValuePair vertex = simplex[i];
-            final double[] point = vertex.getPointRef();
-            if (Double.isNaN(vertex.getValue())) {
-                simplex[i] = new PointValuePair(point, evaluationFunction.value(point), false);
-            }
-        }
-
-        // Sort the simplex from best to worst.
-        Arrays.sort(simplex, comparator);
-    }
-
-    /**
-     * Replace the worst point of the simplex by a new point.
-     *
-     * @param pointValuePair Point to insert.
-     * @param comparator Comparator to use for sorting the simplex vertices
-     * from best to worst.
-     */
-    protected void replaceWorstPoint(PointValuePair pointValuePair,
-                                     final Comparator<PointValuePair> comparator) {
-        for (int i = 0; i < dimension; i++) {
-            if (comparator.compare(simplex[i], pointValuePair) > 0) {
-                PointValuePair tmp = simplex[i];
-                simplex[i] = pointValuePair;
-                pointValuePair = tmp;
-            }
-        }
-        simplex[dimension] = pointValuePair;
-    }
-
-    /**
-     * Get the points of the simplex.
-     *
-     * @return all the simplex points.
-     */
-    public PointValuePair[] getPoints() {
-        final PointValuePair[] copy = new PointValuePair[simplex.length];
-        System.arraycopy(simplex, 0, copy, 0, simplex.length);
-        return copy;
-    }
-
-    /**
-     * Get the simplex point stored at the requested {@code index}.
-     *
-     * @param index Location.
-     * @return the point at location {@code index}.
-     */
-    public PointValuePair getPoint(int index) {
-        if (index < 0 ||
-            index >= simplex.length) {
-            throw new OutOfRangeException(index, 0, simplex.length - 1);
-        }
-        return simplex[index];
-    }
-
-    /**
-     * Store a new point at location {@code index}.
-     * Note that no deep-copy of {@code point} is performed.
-     *
-     * @param index Location.
-     * @param point New value.
-     */
-    protected void setPoint(int index, PointValuePair point) {
-        if (index < 0 ||
-            index >= simplex.length) {
-            throw new OutOfRangeException(index, 0, simplex.length - 1);
-        }
-        simplex[index] = point;
-    }
-
-    /**
-     * Replace all points.
-     * Note that no deep-copy of {@code points} is performed.
-     *
-     * @param points New Points.
-     */
-    protected void setPoints(PointValuePair[] points) {
-        if (points.length != simplex.length) {
-            throw new DimensionMismatchException(points.length, simplex.length);
-        }
-        simplex = points;
-    }
-
-    /**
-     * Create steps for a unit hypercube.
-     *
-     * @param n Dimension of the hypercube.
-     * @param sideLength Length of the sides of the hypercube.
-     * @return the steps.
-     */
-    private static double[] createHypercubeSteps(int n,
-                                                 double sideLength) {
-        final double[] steps = new double[n];
-        for (int i = 0; i < n; i++) {
-            steps[i] = sideLength;
-        }
-        return steps;
-    }
-}
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/MultiDirectionalSimplex.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/MultiDirectionalSimplex.java
deleted file mode 100644
index da47eaa543..0000000000
--- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/MultiDirectionalSimplex.java
+++ /dev/null
@@ -1,215 +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.legacy.optim.nonlinear.scalar.noderiv;
-
-import java.util.Comparator;
-
-import org.apache.commons.math4.legacy.analysis.MultivariateFunction;
-import org.apache.commons.math4.legacy.optim.PointValuePair;
-
-/**
- * This class implements the multi-directional direct search method.
- *
- * @since 3.0
- */
-public class MultiDirectionalSimplex extends AbstractSimplex {
-    /** Default value for {@link #khi}: {@value}. */
-    private static final double DEFAULT_KHI = 2;
-    /** Default value for {@link #gamma}: {@value}. */
-    private static final double DEFAULT_GAMMA = 0.5;
-    /** Expansion coefficient. */
-    private final double khi;
-    /** Contraction coefficient. */
-    private final double gamma;
-
-    /**
-     * Build a multi-directional simplex with default coefficients.
-     * The default values are 2.0 for khi and 0.5 for gamma.
-     *
-     * @param n Dimension of the simplex.
-     */
-    public MultiDirectionalSimplex(final int n) {
-        this(n, 1d);
-    }
-
-    /**
-     * Build a multi-directional simplex with default coefficients.
-     * The default values are 2.0 for khi and 0.5 for gamma.
-     *
-     * @param n Dimension of the simplex.
-     * @param sideLength Length of the sides of the default (hypercube)
-     * simplex. See {@link AbstractSimplex#AbstractSimplex(int,double)}.
-     */
-    public MultiDirectionalSimplex(final int n, double sideLength) {
-        this(n, sideLength, DEFAULT_KHI, DEFAULT_GAMMA);
-    }
-
-    /**
-     * Build a multi-directional simplex with specified coefficients.
-     *
-     * @param n Dimension of the simplex. See
-     * {@link AbstractSimplex#AbstractSimplex(int,double)}.
-     * @param khi Expansion coefficient.
-     * @param gamma Contraction coefficient.
-     */
-    public MultiDirectionalSimplex(final int n,
-                                   final double khi, final double gamma) {
-        this(n, 1d, khi, gamma);
-    }
-
-    /**
-     * Build a multi-directional simplex with specified coefficients.
-     *
-     * @param n Dimension of the simplex. See
-     * {@link AbstractSimplex#AbstractSimplex(int,double)}.
-     * @param sideLength Length of the sides of the default (hypercube)
-     * simplex. See {@link AbstractSimplex#AbstractSimplex(int,double)}.
-     * @param khi Expansion coefficient.
-     * @param gamma Contraction coefficient.
-     */
-    public MultiDirectionalSimplex(final int n, double sideLength,
-                                   final double khi, final double gamma) {
-        super(n, sideLength);
-
-        this.khi   = khi;
-        this.gamma = gamma;
-    }
-
-    /**
-     * Build a multi-directional simplex with default coefficients.
-     * The default values are 2.0 for khi and 0.5 for gamma.
-     *
-     * @param steps Steps along the canonical axes representing box edges.
-     * They may be negative but not zero. See
-     */
-    public MultiDirectionalSimplex(final double[] steps) {
-        this(steps, DEFAULT_KHI, DEFAULT_GAMMA);
-    }
-
-    /**
-     * Build a multi-directional simplex with specified coefficients.
-     *
-     * @param steps Steps along the canonical axes representing box edges.
-     * They may be negative but not zero. See
-     * {@link AbstractSimplex#AbstractSimplex(double[])}.
-     * @param khi Expansion coefficient.
-     * @param gamma Contraction coefficient.
-     */
-    public MultiDirectionalSimplex(final double[] steps,
-                                   final double khi, final double gamma) {
-        super(steps);
-
-        this.khi   = khi;
-        this.gamma = gamma;
-    }
-
-    /**
-     * Build a multi-directional simplex with default coefficients.
-     * The default values are 2.0 for khi and 0.5 for gamma.
-     *
-     * @param referenceSimplex Reference simplex. See
-     * {@link AbstractSimplex#AbstractSimplex(double[][])}.
-     */
-    public MultiDirectionalSimplex(final double[][] referenceSimplex) {
-        this(referenceSimplex, DEFAULT_KHI, DEFAULT_GAMMA);
-    }
-
-    /**
-     * Build a multi-directional simplex with specified coefficients.
-     *
-     * @param referenceSimplex Reference simplex. See
-     * {@link AbstractSimplex#AbstractSimplex(double[][])}.
-     * @param khi Expansion coefficient.
-     * @param gamma Contraction coefficient.
-     * @throws org.apache.commons.math4.legacy.exception.NotStrictlyPositiveException
-     * if the reference simplex does not contain at least one point.
-     * @throws org.apache.commons.math4.legacy.exception.DimensionMismatchException
-     * if there is a dimension mismatch in the reference simplex.
-     */
-    public MultiDirectionalSimplex(final double[][] referenceSimplex,
-                                   final double khi, final double gamma) {
-        super(referenceSimplex);
-
-        this.khi   = khi;
-        this.gamma = gamma;
-    }
-
-    /** {@inheritDoc} */
-    @Override
-    public void iterate(final MultivariateFunction evaluationFunction,
-                        final Comparator<PointValuePair> comparator) {
-        // Save the original simplex.
-        final PointValuePair[] original = getPoints();
-        final PointValuePair best = original[0];
-
-        // Perform a reflection step.
-        final PointValuePair reflected = evaluateNewSimplex(evaluationFunction,
-                                                                original, 1, comparator);
-        if (comparator.compare(reflected, best) < 0) {
-            // Compute the expanded simplex.
-            final PointValuePair[] reflectedSimplex = getPoints();
-            final PointValuePair expanded = evaluateNewSimplex(evaluationFunction,
-                                                                   original, khi, comparator);
-            if (comparator.compare(reflected, expanded) <= 0) {
-                // Keep the reflected simplex.
-                setPoints(reflectedSimplex);
-            }
-            // Keep the expanded simplex.
-            return;
-        }
-
-        // Compute the contracted simplex.
-        evaluateNewSimplex(evaluationFunction, original, gamma, comparator);
-
-    }
-
-    /**
-     * Compute and evaluate a new simplex.
-     *
-     * @param evaluationFunction Evaluation function.
-     * @param original Original simplex (to be preserved).
-     * @param coeff Linear coefficient.
-     * @param comparator Comparator to use to sort simplex vertices from best
-     * to poorest.
-     * @return the best point in the transformed simplex.
-     * @throws org.apache.commons.math4.legacy.exception.TooManyEvaluationsException
-     * if the maximal number of evaluations is exceeded.
-     */
-    private PointValuePair evaluateNewSimplex(final MultivariateFunction evaluationFunction,
-                                                  final PointValuePair[] original,
-                                                  final double coeff,
-                                                  final Comparator<PointValuePair> comparator) {
-        final double[] xSmallest = original[0].getPointRef();
-        // Perform a linear transformation on all the simplex points,
-        // except the first one.
-        setPoint(0, original[0]);
-        final int dim = getDimension();
-        for (int i = 1; i < getSize(); i++) {
-            final double[] xOriginal = original[i].getPointRef();
-            final double[] xTransformed = new double[dim];
-            for (int j = 0; j < dim; j++) {
-                xTransformed[j] = xSmallest[j] + coeff * (xSmallest[j] - xOriginal[j]);
-            }
-            setPoint(i, new PointValuePair(xTransformed, Double.NaN, false));
-        }
-
-        // Evaluate the simplex.
-        evaluate(evaluationFunction, comparator);
-
-        return getPoint(0);
-    }
-}
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/MultiDirectionalTransform.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/MultiDirectionalTransform.java
new file mode 100644
index 0000000000..fa1c849159
--- /dev/null
+++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/MultiDirectionalTransform.java
@@ -0,0 +1,119 @@
+/*
+ * 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.legacy.optim.nonlinear.scalar.noderiv;
+
+import java.util.Comparator;
+import java.util.function.UnaryOperator;
+
+import org.apache.commons.math4.legacy.analysis.MultivariateFunction;
+import org.apache.commons.math4.legacy.optim.PointValuePair;
+import org.apache.commons.math4.legacy.optim.OptimizationData;
+
+/**
+ * Multi-directional search method.
+ */
+public class MultiDirectionalTransform
+    implements Simplex.TransformFactory,
+               OptimizationData {
+    /** Default value for {@link #gamma}: {@value}. */
+    private static final double DEFAULT_GAMMA = 2;
+    /** Default value for {@link #rho}: {@value}. */
+    private static final double DEFAULT_RHO = 0.5;
+    /** Expansion coefficient. */
+    private final double gamma;
+    /** Contraction coefficient. */
+    private final double rho;
+
+    /**
+     * @param gamma Expansion coefficient.
+     * @param rho Contraction coefficient.
+     */
+    public MultiDirectionalTransform(double gamma,
+                                     double rho) {
+        this.gamma = gamma;
+        this.rho = rho;
+    }
+
+    /**
+     * Transform with default values.
+     */
+    public MultiDirectionalTransform() {
+        this(DEFAULT_GAMMA,
+             DEFAULT_RHO);
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public UnaryOperator<Simplex> apply(final MultivariateFunction evaluationFunction,
+                                        final Comparator<PointValuePair> comparator) {
+        return original -> {
+            final PointValuePair best = original.get(0);
+
+            // Perform a reflection step.
+            final Simplex reflectedSimplex = transform(original,
+                                                       1);
+            final PointValuePair reflectedBest = reflectedSimplex.get(0);
+
+            if (comparator.compare(reflectedBest, best) < 0) {
+                // Compute the expanded simplex.
+                final Simplex expandedSimplex = transform(original,
+                                                          gamma);
+                final PointValuePair expandedBest = expandedSimplex.get(0);
+
+                return comparator.compare(reflectedBest, expandedBest) <= 0 ?
+                    reflectedSimplex :
+                    expandedSimplex;
+            } else {
+                // Compute the contracted simplex.
+                return transform(original,
+                                 rho);
+            }
+        };
+    }
+
+    /**
+     * Computes and evaluates a new simplex.
+     *
+     * @param original Original simplex.
+     * @param coeff Linear coefficient.
+     * @return the transformed simplex.
+     * @throws org.apache.commons.math4.legacy.exception.TooManyEvaluationsException
+     * if the maximal number of evaluations is exceeded.
+     */
+    private Simplex transform(Simplex original,
+                              double coeff) {
+        // Transformed simplex is the result a linear transformation on all
+        // points except the first one.
+        final int dim = original.getDimension();
+        final int len = original.getSize();
+        final double[][] simplex = new double[len][];
+
+        simplex[0] = original.get(0).getPoint();
+        final double[] xSmallest = simplex[0];
+
+        for (int i = 1; i < len; i++) {
+            final double[] xOriginal = original.get(i).getPoint();
+            final double[] xTransformed = new double[dim];
+            for (int j = 0; j < dim; j++) {
+                xTransformed[j] = xSmallest[j] + coeff * (xSmallest[j] - xOriginal[j]);
+            }
+            simplex[i] = xTransformed;
+        }
+
+        return Simplex.of(simplex);
+    }
+}
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/NelderMeadSimplex.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/NelderMeadSimplex.java
deleted file mode 100644
index bc923fa6c6..0000000000
--- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/NelderMeadSimplex.java
+++ /dev/null
@@ -1,280 +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.legacy.optim.nonlinear.scalar.noderiv;
-
-import java.util.Comparator;
-
-import org.apache.commons.math4.legacy.analysis.MultivariateFunction;
-import org.apache.commons.math4.legacy.optim.PointValuePair;
-
-/**
- * This class implements the Nelder-Mead simplex algorithm.
- *
- * @since 3.0
- */
-public class NelderMeadSimplex extends AbstractSimplex {
-    /** Default value for {@link #rho}: {@value}. */
-    private static final double DEFAULT_RHO = 1;
-    /** Default value for {@link #khi}: {@value}. */
-    private static final double DEFAULT_KHI = 2;
-    /** Default value for {@link #gamma}: {@value}. */
-    private static final double DEFAULT_GAMMA = 0.5;
-    /** Default value for {@link #sigma}: {@value}. */
-    private static final double DEFAULT_SIGMA = 0.5;
-    /** Reflection coefficient. */
-    private final double rho;
-    /** Expansion coefficient. */
-    private final double khi;
-    /** Contraction coefficient. */
-    private final double gamma;
-    /** Shrinkage coefficient. */
-    private final double sigma;
-
-    /**
-     * Build a Nelder-Mead simplex with default coefficients.
-     * The default coefficients are 1.0 for rho, 2.0 for khi and 0.5
-     * for both gamma and sigma.
-     *
-     * @param n Dimension of the simplex.
-     */
-    public NelderMeadSimplex(final int n) {
-        this(n, 1d);
-    }
-
-    /**
-     * Build a Nelder-Mead simplex with default coefficients.
-     * The default coefficients are 1.0 for rho, 2.0 for khi and 0.5
-     * for both gamma and sigma.
-     *
-     * @param n Dimension of the simplex.
-     * @param sideLength Length of the sides of the default (hypercube)
-     * simplex. See {@link AbstractSimplex#AbstractSimplex(int,double)}.
-     */
-    public NelderMeadSimplex(final int n, double sideLength) {
-        this(n, sideLength,
-             DEFAULT_RHO, DEFAULT_KHI, DEFAULT_GAMMA, DEFAULT_SIGMA);
-    }
-
-    /**
-     * Build a Nelder-Mead simplex with specified coefficients.
-     *
-     * @param n Dimension of the simplex. See
-     * {@link AbstractSimplex#AbstractSimplex(int,double)}.
-     * @param sideLength Length of the sides of the default (hypercube)
-     * simplex. See {@link AbstractSimplex#AbstractSimplex(int,double)}.
-     * @param rho Reflection coefficient.
-     * @param khi Expansion coefficient.
-     * @param gamma Contraction coefficient.
-     * @param sigma Shrinkage coefficient.
-     */
-    public NelderMeadSimplex(final int n, double sideLength,
-                             final double rho, final double khi,
-                             final double gamma, final double sigma) {
-        super(n, sideLength);
-
-        this.rho = rho;
-        this.khi = khi;
-        this.gamma = gamma;
-        this.sigma = sigma;
-    }
-
-    /**
-     * Build a Nelder-Mead simplex with specified coefficients.
-     *
-     * @param n Dimension of the simplex. See
-     * {@link AbstractSimplex#AbstractSimplex(int)}.
-     * @param rho Reflection coefficient.
-     * @param khi Expansion coefficient.
-     * @param gamma Contraction coefficient.
-     * @param sigma Shrinkage coefficient.
-     */
-    public NelderMeadSimplex(final int n,
-                             final double rho, final double khi,
-                             final double gamma, final double sigma) {
-        this(n, 1d, rho, khi, gamma, sigma);
-    }
-
-    /**
-     * Build a Nelder-Mead simplex with default coefficients.
-     * The default coefficients are 1.0 for rho, 2.0 for khi and 0.5
-     * for both gamma and sigma.
-     *
-     * @param steps Steps along the canonical axes representing box edges.
-     * They may be negative but not zero. See
-     */
-    public NelderMeadSimplex(final double[] steps) {
-        this(steps, DEFAULT_RHO, DEFAULT_KHI, DEFAULT_GAMMA, DEFAULT_SIGMA);
-    }
-
-    /**
-     * Build a Nelder-Mead simplex with specified coefficients.
-     *
-     * @param steps Steps along the canonical axes representing box edges.
-     * They may be negative but not zero. See
-     * {@link AbstractSimplex#AbstractSimplex(double[])}.
-     * @param rho Reflection coefficient.
-     * @param khi Expansion coefficient.
-     * @param gamma Contraction coefficient.
-     * @param sigma Shrinkage coefficient.
-     * @throws IllegalArgumentException if one of the steps is zero.
-     */
-    public NelderMeadSimplex(final double[] steps,
-                             final double rho, final double khi,
-                             final double gamma, final double sigma) {
-        super(steps);
-
-        this.rho = rho;
-        this.khi = khi;
-        this.gamma = gamma;
-        this.sigma = sigma;
-    }
-
-    /**
-     * Build a Nelder-Mead simplex with default coefficients.
-     * The default coefficients are 1.0 for rho, 2.0 for khi and 0.5
-     * for both gamma and sigma.
-     *
-     * @param referenceSimplex Reference simplex. See
-     * {@link AbstractSimplex#AbstractSimplex(double[][])}.
-     */
-    public NelderMeadSimplex(final double[][] referenceSimplex) {
-        this(referenceSimplex, DEFAULT_RHO, DEFAULT_KHI, DEFAULT_GAMMA, DEFAULT_SIGMA);
-    }
-
-    /**
-     * Build a Nelder-Mead simplex with specified coefficients.
-     *
-     * @param referenceSimplex Reference simplex. See
-     * {@link AbstractSimplex#AbstractSimplex(double[][])}.
-     * @param rho Reflection coefficient.
-     * @param khi Expansion coefficient.
-     * @param gamma Contraction coefficient.
-     * @param sigma Shrinkage coefficient.
-     * @throws org.apache.commons.math4.legacy.exception.NotStrictlyPositiveException
-     * if the reference simplex does not contain at least one point.
-     * @throws org.apache.commons.math4.legacy.exception.DimensionMismatchException
-     * if there is a dimension mismatch in the reference simplex.
-     */
-    public NelderMeadSimplex(final double[][] referenceSimplex,
-                             final double rho, final double khi,
-                             final double gamma, final double sigma) {
-        super(referenceSimplex);
-
-        this.rho = rho;
-        this.khi = khi;
-        this.gamma = gamma;
-        this.sigma = sigma;
-    }
-
-    /** {@inheritDoc} */
-    @Override
-    public void iterate(final MultivariateFunction evaluationFunction,
-                        final Comparator<PointValuePair> comparator) {
-        // The simplex has n + 1 points if dimension is n.
-        final int n = getDimension();
-
-        // Interesting values.
-        final PointValuePair best = getPoint(0);
-        final PointValuePair secondBest = getPoint(n - 1);
-        final PointValuePair worst = getPoint(n);
-        final double[] xWorst = worst.getPointRef();
-
-        // Compute the centroid of the best vertices (dismissing the worst
-        // point at index n).
-        final double[] centroid = new double[n];
-        for (int i = 0; i < n; i++) {
-            final double[] x = getPoint(i).getPointRef();
-            for (int j = 0; j < n; j++) {
-                centroid[j] += x[j];
-            }
-        }
-        final double scaling = 1.0 / n;
-        for (int j = 0; j < n; j++) {
-            centroid[j] *= scaling;
-        }
-
-        // compute the reflection point
-        final double[] xR = new double[n];
-        for (int j = 0; j < n; j++) {
-            xR[j] = centroid[j] + rho * (centroid[j] - xWorst[j]);
-        }
-        final PointValuePair reflected
-            = new PointValuePair(xR, evaluationFunction.value(xR), false);
-
-        if (comparator.compare(best, reflected) <= 0 &&
-            comparator.compare(reflected, secondBest) < 0) {
-            // Accept the reflected point.
-            replaceWorstPoint(reflected, comparator);
-        } else if (comparator.compare(reflected, best) < 0) {
-            // Compute the expansion point.
-            final double[] xE = new double[n];
-            for (int j = 0; j < n; j++) {
-                xE[j] = centroid[j] + khi * (xR[j] - centroid[j]);
-            }
-            final PointValuePair expanded
-                = new PointValuePair(xE, evaluationFunction.value(xE), false);
-
-            if (comparator.compare(expanded, reflected) < 0) {
-                // Accept the expansion point.
-                replaceWorstPoint(expanded, comparator);
-            } else {
-                // Accept the reflected point.
-                replaceWorstPoint(reflected, comparator);
-            }
-        } else {
-            if (comparator.compare(reflected, worst) < 0) {
-                // Perform an outside contraction.
-                final double[] xC = new double[n];
-                for (int j = 0; j < n; j++) {
-                    xC[j] = centroid[j] + gamma * (xR[j] - centroid[j]);
-                }
-                final PointValuePair outContracted
-                    = new PointValuePair(xC, evaluationFunction.value(xC), false);
-                if (comparator.compare(outContracted, reflected) <= 0) {
-                    // Accept the contraction point.
-                    replaceWorstPoint(outContracted, comparator);
-                    return;
-                }
-            } else {
-                // Perform an inside contraction.
-                final double[] xC = new double[n];
-                for (int j = 0; j < n; j++) {
-                    xC[j] = centroid[j] - gamma * (centroid[j] - xWorst[j]);
-                }
-                final PointValuePair inContracted
-                    = new PointValuePair(xC, evaluationFunction.value(xC), false);
-
-                if (comparator.compare(inContracted, worst) < 0) {
-                    // Accept the contraction point.
-                    replaceWorstPoint(inContracted, comparator);
-                    return;
-                }
-            }
-
-            // Perform a shrink.
-            final double[] xSmallest = getPoint(0).getPointRef();
-            for (int i = 1; i <= n; i++) {
-                final double[] x = getPoint(i).getPoint();
-                for (int j = 0; j < n; j++) {
-                    x[j] = xSmallest[j] + sigma * (x[j] - xSmallest[j]);
-                }
-                setPoint(i, new PointValuePair(x, Double.NaN, false));
-            }
-            evaluate(evaluationFunction, comparator);
-        }
-    }
-}
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/NelderMeadTransform.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/NelderMeadTransform.java
new file mode 100644
index 0000000000..1fbb88a63d
--- /dev/null
+++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/NelderMeadTransform.java
@@ -0,0 +1,182 @@
+/*
+ * 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.legacy.optim.nonlinear.scalar.noderiv;
+
+import java.util.Comparator;
+import java.util.function.UnaryOperator;
+
+import org.apache.commons.math4.legacy.analysis.MultivariateFunction;
+import org.apache.commons.math4.legacy.optim.PointValuePair;
+import org.apache.commons.math4.legacy.optim.OptimizationData;
+
+/**
+ * <a href="https://en.wikipedia.org/wiki/Nelder%E2%80%93Mead_method">Nelder-Mead method</a>.
+ */
+public class NelderMeadTransform
+    implements Simplex.TransformFactory,
+               OptimizationData {
+    /** Default value for {@link #alpha}: {@value}. */
+    private static final double DEFAULT_ALPHA = 1;
+    /** Default value for {@link #gamma}: {@value}. */
+    private static final double DEFAULT_GAMMA = 2;
+    /** Default value for {@link #rho}: {@value}. */
+    private static final double DEFAULT_RHO = 0.5;
+    /** Default value for {@link #sigma}: {@value}. */
+    private static final double DEFAULT_SIGMA = 0.5;
+    /** Reflection coefficient. */
+    private final double alpha;
+    /** Expansion coefficient. */
+    private final double gamma;
+    /** Contraction coefficient. */
+    private final double rho;
+    /** Shrinkage coefficient. */
+    private final double sigma;
+
+    /**
+     * @param alpha Reflection coefficient.
+     * @param gamma Expansion coefficient.
+     * @param rho Contraction coefficient.
+     * @param sigma Shrinkage coefficient.
+     */
+    public NelderMeadTransform(double alpha,
+                               double gamma,
+                               double rho,
+                               double sigma) {
+        this.alpha = alpha;
+        this.gamma = gamma;
+        this.rho = rho;
+        this.sigma = sigma;
+    }
+
+    /**
+     * Transform with default values.
+     */
+    public NelderMeadTransform() {
+        this(DEFAULT_ALPHA,
+             DEFAULT_GAMMA,
+             DEFAULT_RHO,
+             DEFAULT_SIGMA);
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public UnaryOperator<Simplex> apply(final MultivariateFunction evaluationFunction,
+                                        final Comparator<PointValuePair> comparator) {
+        return original -> {
+            Simplex newSimplex = original;
+
+            // The simplex has n + 1 points if dimension is n.
+            final int n = newSimplex.getDimension();
+
+            // Interesting values.
+            final PointValuePair best = newSimplex.get(0);
+            final PointValuePair secondWorst = newSimplex.get(n - 1);
+            final PointValuePair worst = newSimplex.get(n);
+            final double[] xWorst = worst.getPoint();
+
+            // Centroid of the best vertices, dismissing the worst point (at index n).
+            final double[] centroid = new double[n];
+            for (int i = 0; i < n; i++) {
+                final double[] x = newSimplex.get(i).getPoint();
+                for (int j = 0; j < n; j++) {
+                    centroid[j] += x[j];
+                }
+            }
+            final double scaling = 1d / n;
+            for (int j = 0; j < n; j++) {
+                centroid[j] *= scaling;
+            }
+
+            // Reflection.
+            final PointValuePair reflected = newPoint(centroid,
+                                                      -alpha,
+                                                      xWorst,
+                                                      evaluationFunction);
+            if (comparator.compare(reflected, secondWorst) < 0 &&
+                comparator.compare(best, reflected) <= 0) {
+                return newSimplex.withReplacement(n, reflected);
+            }
+
+            if (comparator.compare(reflected, best) < 0) {
+                // Expansion.
+                final PointValuePair expanded = newPoint(centroid,
+                                                         -gamma,
+                                                         xWorst,
+                                                         evaluationFunction);
+                if (comparator.compare(expanded, reflected) < 0) {
+                    return newSimplex.withReplacement(n, expanded);
+                } else {
+                    return newSimplex.withReplacement(n, reflected);
+                }
+            }
+
+            if (comparator.compare(reflected, worst) < 0) {
+                // Outside contraction.
+                final PointValuePair contracted = newPoint(centroid,
+                                                           rho,
+                                                           reflected.getPoint(),
+                                                           evaluationFunction);
+                if (comparator.compare(contracted, reflected) < 0) {
+                    return newSimplex.withReplacement(n, contracted); // Accept contracted point.
+                }
+            } else {
+                // Inside contraction.
+                final PointValuePair contracted = newPoint(centroid,
+                                                           rho,
+                                                           xWorst,
+                                                           evaluationFunction);
+                if (comparator.compare(contracted, worst) < 0) {
+                    return newSimplex.withReplacement(n, contracted); // Accept contracted point.
+                }
+            }
+
+            // Shrink.
+            final int len = newSimplex.getSize();
+            final double[] xBest = best.getPoint();
+            for (int i = 1; i < len; i++) {
+                newSimplex = newSimplex.withReplacement(i, newPoint(xBest,
+                                                                    sigma,
+                                                                    newSimplex.get(i).getPoint(),
+                                                                    evaluationFunction));
+            }
+            return newSimplex;
+        };
+    }
+
+    /**
+     * Evaluates point with coordinates \( a_i + s (b_i - a_i) \).
+     *
+     * @param a Cartesian coordinates.
+     * @param s Scaling factor.
+     * @param b Cartesian coordinates.
+     * @param eval Evaluation function.
+     * @return a new point.
+     */
+    private static PointValuePair newPoint(double[] a,
+                                           double s,
+                                           double[] b,
+                                           MultivariateFunction eval) {
+        final int dim = a.length;
+        final double[] r = new double[dim];
+        for (int i = 0; i < dim; i++) {
+            final double m = a[i];
+            r[i] = m + s * (b[i] - m);
+        }
+
+        return new PointValuePair(r, eval.value(r), false);
+    }
+}
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/Simplex.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/Simplex.java
new file mode 100644
index 0000000000..63e60434b8
--- /dev/null
+++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/Simplex.java
@@ -0,0 +1,302 @@
+/*
+ * 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.legacy.optim.nonlinear.scalar.noderiv;
+
+import java.util.List;
+import java.util.ArrayList;
+import java.util.Arrays;
+import java.util.Comparator;
+import java.util.Collections;
+import java.util.function.BiFunction;
+import java.util.function.UnaryOperator;
+
+import org.apache.commons.math4.legacy.analysis.MultivariateFunction;
+import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
+import org.apache.commons.math4.legacy.exception.MathIllegalArgumentException;
+import org.apache.commons.math4.legacy.exception.NotStrictlyPositiveException;
+import org.apache.commons.math4.legacy.exception.ZeroException;
+import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
+import org.apache.commons.math4.legacy.optim.OptimizationData;
+import org.apache.commons.math4.legacy.optim.PointValuePair;
+
+/**
+ * Represents a simplex.
+ *
+ * @see SimplexOptimizer
+ * @see MultiDirectionalTransform
+ * @see NelderMeadTransform
+ */
+public final class Simplex implements OptimizationData {
+    /** Coordinates. */
+    private final List<PointValuePair> points;
+
+    /**
+     * Builds from a given set of coordinates.
+     *
+     * @param referenceSimplex Reference simplex.
+     * @throws NotStrictlyPositiveException if the reference simplex does not
+     * contain at least one point.
+     * @throws DimensionMismatchException if there is a dimension mismatch
+     * in the reference simplex.
+     * @throws IllegalArgumentException if one of its vertices is duplicated.
+     */
+    private Simplex(double[][] referenceSimplex) {
+        if (referenceSimplex.length <= 0) {
+            throw new NotStrictlyPositiveException(LocalizedFormats.SIMPLEX_NEED_ONE_POINT,
+                                                   referenceSimplex.length);
+        }
+        final int len = referenceSimplex.length;
+        points = new ArrayList<>(len);
+
+        final int dim = len - 1;
+
+        // Loop over vertices.
+        for (int i = 0; i < len; i++) {
+            final double[] refI = referenceSimplex[i];
+
+            // Safety checks.
+            if (refI.length != dim) {
+                throw new DimensionMismatchException(refI.length, dim);
+            }
+            for (int j = 1; j < i; j++) {
+                final double[] refJ = referenceSimplex[j];
+                boolean allEquals = true;
+                for (int k = 0; k < dim; k++) {
+                    if (refI[k] != refJ[k]) {
+                        allEquals = false;
+                        break;
+                    }
+                }
+                if (allEquals) {
+                    throw new MathIllegalArgumentException(LocalizedFormats.EQUAL_VERTICES_IN_SIMPLEX,
+                                                           i, j);
+                }
+            }
+
+            points.add(new PointValuePair(refI, Double.NaN));
+        }
+    }
+
+    /**
+     * Builds from an existing simplex.
+     *
+     * @param points Simplex data.  Reference will be stored in the newly
+     * constructed instance.
+     */
+    private Simplex(List<PointValuePair> points) {
+        this.points = points;
+    }
+
+    /**
+     * Builds from a given set of coordinates.
+     *
+     * @param simplex Reference simplex.
+     * @return a new instance.
+     * @throws NotStrictlyPositiveException if the reference simplex does not
+     * contain at least one point.
+     * @throws DimensionMismatchException if there is a dimension mismatch
+     * in the reference simplex.
+     * @throws IllegalArgumentException if one of its vertices is duplicated.
+     */
+    public static Simplex of(double[][] simplex) {
+        return new Simplex(simplex);
+    }
+
+    /**
+     * Builds a unit hypercube simplex.
+     *
+     * @param n Dimension of the simplex.
+     * @return a new instance.
+     */
+    public static Simplex of(int n) {
+        return of(n, 1d);
+    }
+
+    /**
+     * Builds a hypercube simplex with the given side length.
+     *
+     * @param n Dimension of the simplex.
+     * @param sideLength Length of the sides of the hypercube.
+     * @return a new instance.
+     */
+    public static Simplex of(int n,
+                             double sideLength) {
+        final double[] steps = new double[n];
+        Arrays.fill(steps, sideLength);
+        return of(steps);
+    }
+
+    /**
+     * The start configuration for simplex is built from a box parallel to
+     * the canonical axes of the space. The simplex is the subset of vertices
+     * of a box parallel to the canonical axes. It is built as the path followed
+     * while traveling from one vertex of the box to the diagonally opposite
+     * vertex moving only along the box edges. The first vertex of the box will
+     * be located at the start point of the optimization.
+     * As an example, in dimension 3 a simplex has 4 vertices. Setting the
+     * steps to (1, 10, 2) and the start point to (1, 1, 1) would imply the
+     * start simplex would be: { (1, 1, 1), (2, 1, 1), (2, 11, 1), (2, 11, 3) }.
+     * The first vertex would be set to the start point at (1, 1, 1) and the
+     * last vertex would be set to the diagonally opposite vertex at (2, 11, 3).
+     *
+     * @param steps Steps along the canonical axes representing box edges.
+     * They may be negative but not zero.
+     * @throws ZeroException if one of the steps is zero.
+     * @return a new instance.
+     */
+    public static Simplex of(double[] steps) {
+        if (steps.length == 0) {
+            throw new ZeroException();
+        }
+        final int dim = steps.length;
+        final int len = dim + 1;
+
+        // Only the relative position of the n final vertices with respect
+        // to the first one are stored.
+        final double[][] simplex = new double[len][dim];
+        for (int i = 1; i < len; i++) { // First point is the origin (zero).
+            final double[] vertexI = simplex[i];
+            for (int j = 0; j < i; j++) {
+                if (steps[j] == 0) {
+                    throw new ZeroException(LocalizedFormats.EQUAL_VERTICES_IN_SIMPLEX);
+                }
+                System.arraycopy(steps, 0, vertexI, 0, j + 1);
+            }
+        }
+
+        return new Simplex(simplex);
+    }
+
+    /**
+     * Returns the space dimension.
+     *
+     * @return the dimension of the simplex.
+     */
+    public int getDimension() {
+        return points.size() - 1;
+    }
+
+    /**
+     * Returns the number of vertices.
+     *
+     * @return the size of the simplex.
+     */
+    public int getSize() {
+        return points.size();
+    }
+
+    /**
+     * Evaluates the (non-evaluated) simplex points and returns a new instance
+     * with vertices sorted from best to worst.
+     *
+     * @param function Evaluation function.
+     * @param comparator Comparator for sorting vertices, from best to worst.
+     * @return a new instance in which the vertices are sorted according to
+     * the given {@code comparator}.
+     */
+    public Simplex evaluate(MultivariateFunction function,
+                            Comparator<PointValuePair> comparator) {
+        final List<PointValuePair> newPoints = new ArrayList<>(points.size());
+        for (PointValuePair pv : points) {
+            final double[] coord = pv.getPoint();
+            final double value = Double.isNaN(pv.getValue()) ?
+                function.value(coord) :
+                pv.getValue();
+
+            newPoints.add(new PointValuePair(coord, value, false));
+        }
+
+        Collections.sort(newPoints, comparator);
+        return new Simplex(newPoints);
+    }
+
+    /**
+     * Retrieves a copy of the simplex point stored at {@code index}.
+     *
+     * @param index Location.
+     * @return the point at location {@code index}.
+     */
+    public PointValuePair get(int index) {
+        final PointValuePair p = points.get(index);
+        return new PointValuePair(p.getPoint(), p.getValue());
+    }
+
+    /**
+     * Translates the simplex such that its first point is at the given {@code point}.
+     *
+     * @param point Coordinates of the new simplex's first point.
+     * @return a new instance.
+     * @throws DimensionMismatchException if the start point does not match
+     * simplex dimension.
+     */
+    public Simplex translate(final double[] point) {
+        final int dim = getDimension();
+        if (dim != point.length) {
+            throw new DimensionMismatchException(dim, point.length);
+        }
+
+        final int len = points.size();
+        final double[][] coordinates = new double[len][dim];
+        final double[] current0 = points.get(0).getPoint(); // Current first point.
+
+        // Set new vertices.
+        for (int i = 0; i < len; i++) {
+            final double[] currentI = points.get(i).getPoint();
+            final double[] newI = coordinates[i];
+            for (int k = 0; k < dim; k++) {
+                newI[k] = point[k] + currentI[k] - current0[k];
+            }
+        }
+
+        return of(coordinates);
+    }
+
+    /**
+     * Builds a new simplex where the given {@code point} replaces the
+     * one at {@code index} in this instance.
+     *
+     * @param index Index of the point to replace.
+     * @param point Replacement for the point currently at {@code index}.
+     * @return a new instance.
+     * @throws IndexOutOfBoundsException if {@code index} is out of bounds.
+     */
+    public Simplex withReplacement(int index,
+                                   PointValuePair point) {
+        final int len = points.size();
+        if (index < 0 ||
+            index >= len) {
+            throw new IndexOutOfBoundsException("index: " + index);
+        }
+
+        final List<PointValuePair> newPoints = new ArrayList<>(len);
+        for (int i = 0; i < len; i++) {
+            final PointValuePair pv = i == index ?
+                point :
+                points.get(i);
+            newPoints.add(new PointValuePair(pv.getPoint(), pv.getValue(), false));
+        }
+
+        return new Simplex(newPoints);
+    }
+
+    /**
+     * Generator of simplex transform.
+     */
+    public interface TransformFactory
+        extends BiFunction<MultivariateFunction, Comparator<PointValuePair>, UnaryOperator<Simplex>> {}
+}
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/SimplexOptimizer.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/SimplexOptimizer.java
index 48aaafd49a..60af50d08e 100644
--- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/SimplexOptimizer.java
+++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/SimplexOptimizer.java
@@ -17,10 +17,10 @@
 package org.apache.commons.math4.legacy.optim.nonlinear.scalar.noderiv;
 
 import java.util.Comparator;
+import java.util.function.UnaryOperator;
 
 import org.apache.commons.math4.legacy.analysis.MultivariateFunction;
 import org.apache.commons.math4.legacy.exception.MathUnsupportedOperationException;
-import org.apache.commons.math4.legacy.exception.NullArgumentException;
 import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
 import org.apache.commons.math4.legacy.optim.ConvergenceChecker;
 import org.apache.commons.math4.legacy.optim.OptimizationData;
@@ -33,62 +33,57 @@
  * 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>
+ * 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>
- *  The simplex update procedure ({@link NelderMeadSimplex} or
- * {@link MultiDirectionalSimplex})  must be passed to the
- * {@code optimize} method.
- * </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>
- *  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>
+ * In addition to those documented in
+ * {@link MultivariateOptimizer#optimize(OptimizationData[]) MultivariateOptimizer},
+ * an instance of this class will register the following data:
+ * <ul>
+ *  <li>{@link Simplex}</li>
+ *  <li>{@link Simplex.TransformFactory} (either {@link NelderMeadTransform}
+ *   or {@link MultiDirectionalTransform})</li>
+ * </ul>
+ *
  * <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>
+ * 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>
- *  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.legacy.optim.nonlinear.scalar.MultivariateFunctionMappingAdapter
- *  MultivariateFunctionMappingAdapter} or
- *  {@link org.apache.commons.math4.legacy.optim.nonlinear.scalar.MultivariateFunctionPenaltyAdapter
- *  MultivariateFunctionPenaltyAdapter}.
- *  <br>
- *  The call to {@link #optimize(OptimizationData[]) optimize} will throw
- *  {@link MathUnsupportedOperationException} if bounds are passed to it.
- * </p>
+ * Convergence is considered achieved when <em>all</em> the simplex points
+ * have converged.
  *
- * @since 3.0
+ * <p>
+ * This implementation does not directly support constrained optimization
+ * with simple bounds.
+ * The call to {@link #optimize(OptimizationData[]) optimize} will throw
+ * {@link MathUnsupportedOperationException} if bounds are passed to it.
  */
 public class SimplexOptimizer extends MultivariateOptimizer {
-    /** Simplex update rule. */
-    private AbstractSimplex simplex;
+    /** Simplex update function factory. */
+    private Simplex.TransformFactory updateRule;
+    /** Current simplex. */
+    private Simplex simplex;
 
     /**
      * @param checker Convergence checker.
@@ -101,27 +96,11 @@ public SimplexOptimizer(ConvergenceChecker<PointValuePair> checker) {
      * @param rel Relative threshold.
      * @param abs Absolute threshold.
      */
-    public SimplexOptimizer(double rel, double abs) {
+    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>
-     * </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() {
@@ -129,58 +108,42 @@ protected PointValuePair doOptimize() {
 
         // 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 MultivariateFunction evalFunc = this::computeObjectiveValue;
 
         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);
-            }
+        final Comparator<PointValuePair> comparator = (o1, o2) -> {
+            final double v1 = o1.getValue();
+            final double v2 = o2.getValue();
+            return isMinim ? Double.compare(v1, v2) : Double.compare(v2, v1);
         };
 
+        final UnaryOperator<Simplex> update = updateRule.apply(evalFunc, comparator);
+
         // Initialize search.
-        simplex.build(getStartPoint());
-        simplex.evaluate(evalFunc, comparator);
+        simplex = simplex.translate(getStartPoint()).evaluate(evalFunc, comparator);
 
-        PointValuePair[] previous = null;
-        int iteration = 0;
+        Simplex previous = null;
         final ConvergenceChecker<PointValuePair> checker = getConvergenceChecker();
         while (true) {
-            iteration = getIterations();
-            if (iteration > 0) {
+            if (previous != null) { // Skip check at first iteration.
                 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".
+                    if (!checker.converged(getIterations(),
+                                           previous.get(i),
+                                           simplex.get(i))) {
+                        converged = false;
                         break;
                     }
                 }
                 if (converged) {
                     // We have found an optimum.
-                    return simplex.getPoint(0);
+                    return simplex.get(0);
                 }
             }
 
             // We still need to search.
-            previous = simplex.getPoints();
-            simplex.iterate(evalFunc, comparator);
+            previous = simplex;
+            simplex = update.apply(simplex).evaluate(evalFunc, comparator);
 
             incrementIterationCount();
         }
@@ -193,7 +156,8 @@ public int compare(final PointValuePair o1,
      * @param optData Optimization data.
      * The following data will be looked for:
      * <ul>
-     *  <li>{@link AbstractSimplex}</li>
+     *  <li>{@link Simplex}</li>
+     *  <li>{@link Simplex.TransformFactory}</li>
      * </ul>
      */
     @Override
@@ -201,14 +165,13 @@ 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.
+        // 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;
-                // If more data must be parsed, this statement _must_ be
-                // changed to "continue".
-                break;
+            if (data instanceof Simplex) {
+                simplex = (Simplex) data;
+            } else if (data instanceof Simplex.TransformFactory) {
+                updateRule = (Simplex.TransformFactory) data;
             }
         }
     }
@@ -216,12 +179,15 @@ protected void parseOptimizationData(OptimizationData... optData) {
     /**
      * @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.
+     * @throws NullPointerException if no initial simplex or no transform rule
+     * was passed to the {@link #optimize(OptimizationData[]) optimize} method.
      */
     private void checkParameters() {
+        if (updateRule == null) {
+            throw new NullPointerException("No update rule");
+        }
         if (simplex == null) {
-            throw new NullArgumentException();
+            throw new NullPointerException("No initial simplex");
         }
         if (getLowerBound() != null ||
             getUpperBound() != null) {
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/random/HaltonSequenceGenerator.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/quasirandom/HaltonSequenceGenerator.java
similarity index 80%
rename from commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/random/HaltonSequenceGenerator.java
rename to commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/quasirandom/HaltonSequenceGenerator.java
index 8a175618bf..8562714f0f 100644
--- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/random/HaltonSequenceGenerator.java
+++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/quasirandom/HaltonSequenceGenerator.java
@@ -14,14 +14,15 @@
  * See the License for the specific language governing permissions and
  * limitations under the License.
  */
-package org.apache.commons.math4.legacy.random;
-
-import java.util.function.Supplier;
+package org.apache.commons.math4.legacy.quasirandom;
 
 import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
+import org.apache.commons.math4.legacy.exception.NotPositiveException;
 import org.apache.commons.math4.legacy.exception.NullArgumentException;
 import org.apache.commons.math4.legacy.exception.OutOfRangeException;
 
+import java.util.Arrays;
+
 /**
  * Implementation of a Halton sequence.
  * <p>
@@ -43,7 +44,7 @@
  * The generator supports two modes:
  * <ul>
  *   <li>sequential generation of points: {@link #get()}</li>
- *   <li>random access to the i-th point in the sequence: {@link #skipTo(int)}</li>
+ *   <li>random access to the i-th point in the sequence: {@link LowDiscrepancySequence#jump(long)}</li>
  * </ul>
  *
  * @see <a href="http://en.wikipedia.org/wiki/Halton_sequence">Halton sequence (Wikipedia)</a>
@@ -51,7 +52,7 @@
  * On the Halton sequence and its scramblings</a>
  * @since 3.3
  */
-public class HaltonSequenceGenerator implements Supplier<double[]> {
+public class HaltonSequenceGenerator implements LowDiscrepancySequence {
 
     /** The first 40 primes. */
     private static final int[] PRIMES = new int[] {
@@ -71,7 +72,7 @@ public class HaltonSequenceGenerator implements Supplier<double[]> {
     private final int dimension;
 
     /** The current index in the sequence. */
-    private int count;
+    private long count;
 
     /** The base numbers for each component. */
     private final int[] base;
@@ -89,6 +90,13 @@ public HaltonSequenceGenerator(final int dimension) {
         this(dimension, PRIMES, WEIGHTS);
     }
 
+    private HaltonSequenceGenerator(HaltonSequenceGenerator original){
+        this.base = original.base;
+        this.count = original.count;
+        this.dimension = original.dimension;
+        this.weight = Arrays.copyOf(original.weight,original.weight.length);
+    }
+
     /**
      * Construct a new Halton sequence generator with the given base numbers and weights for each dimension.
      * The length of the bases array defines the space dimension and is required to be &gt; 0.
@@ -123,12 +131,12 @@ public HaltonSequenceGenerator(final int dimension, final int[] bases, final int
     public double[] get() {
         final double[] v = new double[dimension];
         for (int i = 0; i < dimension; i++) {
-            int index = count;
+            long index = count;
             double f = 1.0 / base[i];
 
             int j = 0;
             while (index > 0) {
-                final int digit = scramble(i, j, base[i], index % base[i]);
+                final long digit = scramble(i, j, base[i], index % base[i]);
                 v[i] += f * digit;
                 index /= base[i]; // floor( index / base )
                 f /= base[i];
@@ -151,31 +159,43 @@ public double[] get() {
      * @param digit the j-th digit
      * @return the scrambled digit
      */
-    protected int scramble(final int i, final int j, final int b, final int digit) {
+    protected long scramble(final int i, final int j, final int b, final long digit) {
         return weight != null ? (weight[i] * digit) % b : digit;
     }
 
     /**
-     * Skip to the i-th point in the Halton sequence.
+§     * Jump to the i-th point in the Halton sequence.
      * <p>
      * This operation can be performed in O(1).
      *
      * @param index the index in the sequence to skip to
-     * @return the i-th point in the Halton sequence
-     * @throws org.apache.commons.math4.legacy.exception.NotPositiveException NotPositiveException if index &lt; 0
+     * @return the copy of this sequence
+     * @throws org.apache.commons.math4.legacy.exception.NotPositiveException if index &lt; 0
+     */
+    @Override
+    public LowDiscrepancySequence jump(final long index) throws NotPositiveException {
+        if(index < 0) {
+            throw new NotPositiveException(index);
+        }
+        HaltonSequenceGenerator copy = this.copy();
+        copy.performJump(index);
+        return copy;
+    }
+
+    /**
+     *  Do jump at the index.
+     * @param index
      */
-    public double[] skipTo(final int index) {
+    private void performJump(long index) {
         count = index;
-        return get();
     }
 
     /**
-     * Returns the index i of the next point in the Halton sequence that will be returned
-     * by calling {@link #get()}.
-     *
-     * @return the index of the next point
+     * Private constructor avoid side effects.
+     * @return copy of HaltonSequenceGenerator
      */
-    public int getNextIndex() {
-        return count;
+    private HaltonSequenceGenerator copy() {
+        return new HaltonSequenceGenerator(this);
     }
+
 }
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/random/NormalizedRandomGenerator.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/quasirandom/LowDiscrepancySequence.java
similarity index 53%
rename from commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/random/NormalizedRandomGenerator.java
rename to commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/quasirandom/LowDiscrepancySequence.java
index e41a008c44..0d6c834c82 100644
--- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/random/NormalizedRandomGenerator.java
+++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/quasirandom/LowDiscrepancySequence.java
@@ -14,24 +14,19 @@
  * See the License for the specific language governing permissions and
  * limitations under the License.
  */
+package org.apache.commons.math4.legacy.quasirandom;
 
-package org.apache.commons.math4.legacy.random;
+import java.util.function.Supplier;
 
-/**
- * This interface represent a normalized random generator for
- * scalars.
- * Normalized generator provide null mean and unit standard deviation scalars.
- * @since 1.2
+/** Interface to Low Discrepancy Sequence generator and supplier.
+ *  Supplier of a low discrepancy vectors.
+ *  Offers navigation through underlying sequence.
  */
-public interface NormalizedRandomGenerator {
-
-  /** Generate a random scalar with null mean and unit standard deviation.
-   * <p>This method does <strong>not</strong> specify the shape of the
-   * distribution, it is the implementing class that provides it. The
-   * only contract here is to generate numbers with null mean and unit
-   * standard deviation.</p>
-   * @return a random scalar with null mean and unit standard deviation
-   */
-  double nextNormalizedDouble();
-
+public interface LowDiscrepancySequence extends Supplier<double[]> {
+    /**
+     *  Skip to the index in the sequence.
+     * @param index of the element to skip to
+     * @return T element at the index
+     */
+    LowDiscrepancySequence jump(long index);
 }
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/random/SobolSequenceGenerator.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/quasirandom/SobolSequenceGenerator.java
similarity index 88%
rename from commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/random/SobolSequenceGenerator.java
rename to commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/quasirandom/SobolSequenceGenerator.java
index 8269eec95c..e89fdf4b4d 100644
--- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/random/SobolSequenceGenerator.java
+++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/quasirandom/SobolSequenceGenerator.java
@@ -14,7 +14,7 @@
  * See the License for the specific language governing permissions and
  * limitations under the License.
  */
-package org.apache.commons.math4.legacy.random;
+package org.apache.commons.math4.legacy.quasirandom;
 
 import java.io.BufferedReader;
 import java.io.IOException;
@@ -24,10 +24,10 @@
 import java.util.Arrays;
 import java.util.NoSuchElementException;
 import java.util.StringTokenizer;
-import java.util.function.Supplier;
 
 import org.apache.commons.math4.legacy.exception.MathInternalError;
 import org.apache.commons.math4.legacy.exception.MathParseException;
+import org.apache.commons.math4.legacy.exception.NotPositiveException;
 import org.apache.commons.math4.legacy.exception.NotStrictlyPositiveException;
 import org.apache.commons.math4.legacy.exception.OutOfRangeException;
 import org.apache.commons.math4.legacy.core.jdkmath.AccurateMath;
@@ -45,7 +45,7 @@
  * The generator supports two modes:
  * <ul>
  *   <li>sequential generation of points: {@link #get()}</li>
- *   <li>random access to the i-th point in the sequence: {@link #skipTo(int)}</li>
+ *   <li>random access to the i-th point in the sequence: {@link LowDiscrepancySequence#jump(long)}</li>
  * </ul>
  *
  * @see <a href="http://en.wikipedia.org/wiki/Sobol_sequence">Sobol sequence (Wikipedia)</a>
@@ -53,7 +53,8 @@
  *
  * @since 3.3
  */
-public class SobolSequenceGenerator implements Supplier<double[]> {
+public class SobolSequenceGenerator implements LowDiscrepancySequence {
+
     /** The number of bits to use. */
     private static final int BITS = 52;
 
@@ -73,7 +74,7 @@ public class SobolSequenceGenerator implements Supplier<double[]> {
     private final int dimension;
 
     /** The current index in the sequence. */
-    private int count;
+    private long count;
 
     /** The direction vector for each component. */
     private final long[][] direction;
@@ -171,6 +172,13 @@ public SobolSequenceGenerator(final int dimension, final InputStream is)
         }
     }
 
+    private SobolSequenceGenerator(SobolSequenceGenerator sobolSequenceGenerator) {
+        this.dimension = sobolSequenceGenerator.dimension;
+        this.count = sobolSequenceGenerator.count;
+        this.direction = Arrays.copyOf(sobolSequenceGenerator.direction, sobolSequenceGenerator.direction.length);
+        this.x = Arrays.copyOf(sobolSequenceGenerator.x, sobolSequenceGenerator.x.length);
+    }
+
     /**
      * Load the direction vector for each dimension from the given stream.
      * <p>
@@ -261,7 +269,7 @@ public double[] get() {
 
         // find the index c of the rightmost 0
         int c = 1;
-        int value = count - 1;
+        long value = count - 1;
         while ((value & 1) == 1) {
             value >>= 1;
             c++;
@@ -281,15 +289,30 @@ public double[] get() {
      * This operation can be performed in O(1).
      *
      * @param index the index in the sequence to skip to
-     * @return the i-th point in the Sobol sequence
-     * @throws org.apache.commons.math4.legacy.exception.NotPositiveException NotPositiveException if index &lt; 0
+     * @return the sequence it self
+     * @throws NotPositiveException if index &lt; 0
      */
-    public double[] skipTo(final int index) {
+    @Override
+    public LowDiscrepancySequence jump(final long index) throws NotPositiveException {
+        if(index < 0) {
+            throw new NotPositiveException(index);
+        }
+
+        SobolSequenceGenerator copy = this.copy();
+        copy.performJump(index);
+        return copy;
+    }
+
+    /**
+     *  Do jump at the index.
+     * @param index
+     */
+    private void performJump(long index) {
         if (index == 0) {
             // reset x vector
             Arrays.fill(x, 0);
         } else {
-            final int i = index - 1;
+            final long i = index - 1;
             final long grayCode = i ^ (i >> 1); // compute the gray code of i = i XOR floor(i / 2)
             for (int j = 0; j < dimension; j++) {
                 long result = 0;
@@ -307,17 +330,14 @@ public double[] skipTo(final int index) {
             }
         }
         count = index;
-        return get();
     }
 
     /**
-     * Returns the index i of the next point in the Sobol sequence that will be returned
-     * by calling {@link #get()}.
-     *
-     * @return the index of the next point
+     * Private constructor avoid side effects.
+     * @return copy of LowDiscrepancySequence
      */
-    public int getNextIndex() {
-        return count;
+    private SobolSequenceGenerator copy() {
+        return new SobolSequenceGenerator(this);
     }
 
 }
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/quasirandom/package-info.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/quasirandom/package-info.java
new file mode 100644
index 0000000000..ef4d50192e
--- /dev/null
+++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/quasirandom/package-info.java
@@ -0,0 +1,20 @@
+/*
+ * 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.
+ */
+/**
+ * Contains quasi random sequences.
+ */
+package org.apache.commons.math4.legacy.quasirandom;
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/random/StableRandomGenerator.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/random/StableRandomGenerator.java
deleted file mode 100644
index cfb6b4f48c..0000000000
--- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/random/StableRandomGenerator.java
+++ /dev/null
@@ -1,133 +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.legacy.random;
-
-import org.apache.commons.rng.UniformRandomProvider;
-import org.apache.commons.math4.legacy.exception.NullArgumentException;
-import org.apache.commons.math4.legacy.exception.OutOfRangeException;
-import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
-import org.apache.commons.math4.legacy.core.jdkmath.AccurateMath;
-
-/**
- * <p>This class provides a stable normalized random generator. It samples from a stable
- * distribution with location parameter 0 and scale 1.</p>
- *
- * <p>The implementation uses the Chambers-Mallows-Stuck method as described in
- * <i>Handbook of computational statistics: concepts and methods</i> by
- * James E. Gentle, Wolfgang H&auml;rdle, Yuichi Mori.</p>
- *
- * @since 3.0
- */
-public class StableRandomGenerator implements NormalizedRandomGenerator {
-    /** Underlying generator. */
-    private final UniformRandomProvider generator;
-    /** stability parameter. */
-    private final double alpha;
-    /** skewness parameter. */
-    private final double beta;
-    /** cache of expression value used in generation. */
-    private final double zeta;
-
-    /**
-     * Create a new generator.
-     *
-     * @param generator Underlying random generator
-     * @param alpha Stability parameter. Must be in range (0, 2]
-     * @param beta Skewness parameter. Must be in range [-1, 1]
-     * @throws NullArgumentException if generator is null
-     * @throws OutOfRangeException if {@code alpha <= 0} or {@code alpha > 2}
-     * or {@code beta < -1} or {@code beta > 1}
-     */
-    public StableRandomGenerator(final UniformRandomProvider generator,
-                                 final double alpha, final double beta)
-        throws NullArgumentException, OutOfRangeException {
-        if (generator == null) {
-            throw new NullArgumentException();
-        }
-
-        if (!(alpha > 0d && alpha <= 2d)) {
-            throw new OutOfRangeException(LocalizedFormats.OUT_OF_RANGE_LEFT,
-                    alpha, 0, 2);
-        }
-
-        if (!(beta >= -1d && beta <= 1d)) {
-            throw new OutOfRangeException(LocalizedFormats.OUT_OF_RANGE_SIMPLE,
-                    beta, -1, 1);
-        }
-
-        this.generator = generator;
-        this.alpha = alpha;
-        this.beta = beta;
-        if (alpha < 2d && beta != 0d) {
-            zeta = beta * AccurateMath.tan(AccurateMath.PI * alpha / 2);
-        } else {
-            zeta = 0d;
-        }
-    }
-
-    /**
-     * Generate a random scalar with zero location and unit scale.
-     *
-     * @return a random scalar with zero location and unit scale
-     */
-    @Override
-    public double nextNormalizedDouble() {
-        // we need 2 uniform random numbers to calculate omega and phi
-        double omega = -AccurateMath.log(generator.nextDouble());
-        double phi = AccurateMath.PI * (generator.nextDouble() - 0.5);
-
-        // Normal distribution case (Box-Muller algorithm)
-        if (alpha == 2d) {
-            return AccurateMath.sqrt(2d * omega) * AccurateMath.sin(phi);
-        }
-
-        double x;
-        // when beta = 0, zeta is zero as well
-        // Thus we can exclude it from the formula
-        if (beta == 0d) {
-            // Cauchy distribution case
-            if (alpha == 1d) {
-                x = AccurateMath.tan(phi);
-            } else {
-                x = AccurateMath.pow(omega * AccurateMath.cos((1 - alpha) * phi),
-                    1d / alpha - 1d) *
-                    AccurateMath.sin(alpha * phi) /
-                    AccurateMath.pow(AccurateMath.cos(phi), 1d / alpha);
-            }
-        } else {
-            // Generic stable distribution
-            double cosPhi = AccurateMath.cos(phi);
-            // to avoid rounding errors around alpha = 1
-            if (AccurateMath.abs(alpha - 1d) > 1e-8) {
-                double alphaPhi = alpha * phi;
-                double invAlphaPhi = phi - alphaPhi;
-                x = (AccurateMath.sin(alphaPhi) + zeta * AccurateMath.cos(alphaPhi)) / cosPhi *
-                    (AccurateMath.cos(invAlphaPhi) + zeta * AccurateMath.sin(invAlphaPhi)) /
-                     AccurateMath.pow(omega * cosPhi, (1 - alpha) / alpha);
-            } else {
-                double betaPhi = AccurateMath.PI / 2 + beta * phi;
-                x = 2d / AccurateMath.PI * (betaPhi * AccurateMath.tan(phi) - beta *
-                    AccurateMath.log(AccurateMath.PI / 2d * omega * cosPhi / betaPhi));
-
-                if (alpha != 1d) {
-                    x += beta * AccurateMath.tan(AccurateMath.PI * alpha / 2);
-                }
-            }
-        }
-        return x;
-    }
-}
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/descriptive/rank/PSquarePercentile.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/descriptive/rank/PSquarePercentile.java
index 01664fab22..5ffc871a22 100644
--- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/descriptive/rank/PSquarePercentile.java
+++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/descriptive/rank/PSquarePercentile.java
@@ -767,7 +767,7 @@ private double estimate() {
                             new double[] { xval[1], xval[1 + delta] };
                     final double[] yBad =
                             new double[] { yval[1], yval[1 + delta] };
-                    SortInPlace.ASCENDING.accept(xBad, yBad);// since d can be +/- 1
+                    SortInPlace.ASCENDING.apply(xBad, yBad);// since d can be +/- 1
                     univariateFunction = linear.interpolate(xBad, yBad);
                     markerHeight = univariateFunction.value(xD);
                 }
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/descriptive/rank/Percentile.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/descriptive/rank/Percentile.java
index 61c5dc5166..8d2ac8883e 100644
--- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/descriptive/rank/Percentile.java
+++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/descriptive/rank/Percentile.java
@@ -1150,7 +1150,7 @@ protected double index(final double p, final int length) {
             @Override
             public double evaluate(final double[] work, final double[] sampleWeights,
                                    final double p) {
-                SortInPlace.ASCENDING.accept(work, sampleWeights);
+                SortInPlace.ASCENDING.apply(work, sampleWeights);
                 double[] sk = new double[work.length];
                 for(int k = 0; k < work.length; k++) {
                     sk[k] = 0;
diff --git a/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/distribution/EmpiricalDistributionTest.java b/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/distribution/EmpiricalDistributionTest.java
index 70b9f1d8b0..b86b7c348c 100644
--- a/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/distribution/EmpiricalDistributionTest.java
+++ b/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/distribution/EmpiricalDistributionTest.java
@@ -43,7 +43,6 @@
 import org.junit.Assert;
 import org.junit.Before;
 import org.junit.Test;
-import org.junit.Ignore;
 
 /**
  * Test cases for the {@link EmpiricalDistribution} class.
@@ -667,7 +666,6 @@ public void testMath1431() {
         }
     }
 
-    @Ignore
     @Test
     public void testMath1462() {
         final double[] data = {
@@ -686,16 +684,23 @@ public void testMath1462() {
             6461.3944, 6384.1345
         };
 
-        final EmpiricalDistribution ed = new EmpiricalDistribution(data.length);
+        final EmpiricalDistribution ed = new EmpiricalDistribution(data.length / 10);
         ed.load(data);
 
-        final double p50 = ed.inverseCumulativeProbability(0.5);
-        final double p51 = ed.inverseCumulativeProbability(0.51111);
-        final double p49 = ed.inverseCumulativeProbability(0.49999);
+        double v;
+        double p;
 
-        Assert.assertTrue(p51 < 6350);
-        Assert.assertTrue(p49 < 6341);
-        Assert.assertTrue(p50 < 7000);
+        p = 0.49999;
+        v = ed.inverseCumulativeProbability(p);
+        Assert.assertTrue("p=" + p + " => v=" + v, v < 6344);
+
+        p = 0.5;
+        v = ed.inverseCumulativeProbability(p);
+        Assert.assertTrue("p=" + p + " => v=" + v, v < 7000);
+
+        p = 0.51111;
+        v = ed.inverseCumulativeProbability(p);
+        Assert.assertTrue("p=" + p + " => v=" + v, v < 6350);
     }
 
     /**
diff --git a/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/MultiStartMultivariateOptimizerTest.java b/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/MultiStartMultivariateOptimizerTest.java
index f324c3885b..a431635ec9 100644
--- a/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/MultiStartMultivariateOptimizerTest.java
+++ b/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/MultiStartMultivariateOptimizerTest.java
@@ -25,8 +25,9 @@
 import org.apache.commons.math4.legacy.optim.SimpleValueChecker;
 import org.apache.commons.math4.legacy.optim.nonlinear.scalar.gradient.CircleScalar;
 import org.apache.commons.math4.legacy.optim.nonlinear.scalar.gradient.NonLinearConjugateGradientOptimizer;
-import org.apache.commons.math4.legacy.optim.nonlinear.scalar.noderiv.NelderMeadSimplex;
+import org.apache.commons.math4.legacy.optim.nonlinear.scalar.noderiv.NelderMeadTransform;
 import org.apache.commons.math4.legacy.optim.nonlinear.scalar.noderiv.SimplexOptimizer;
+import org.apache.commons.math4.legacy.optim.nonlinear.scalar.noderiv.Simplex;
 import org.apache.commons.rng.UniformRandomProvider;
 import org.apache.commons.rng.simple.RandomSource;
 import org.apache.commons.rng.sampling.distribution.GaussianSampler;
@@ -59,6 +60,7 @@ public void testCircleFitting() {
             = optimizer.optimize(new MaxEval(1000),
                                  circle.getObjectiveFunction(),
                                  circle.getObjectiveFunctionGradient(),
+                                 new NelderMeadTransform(),
                                  GoalType.MINIMIZE,
                                  new InitialGuess(new double[] { 98.680, 47.345 }));
         Assert.assertEquals(1000, optimizer.getMaxEvaluations());
@@ -83,7 +85,7 @@ public void testRosenbrock() {
         Rosenbrock rosenbrock = new Rosenbrock();
         SimplexOptimizer underlying
             = new SimplexOptimizer(new SimpleValueChecker(-1, 1e-3));
-        NelderMeadSimplex simplex = new NelderMeadSimplex(new double[][] {
+        final Simplex simplex = Simplex.of(new double[][] {
                 { -1.2,  1.0 },
                 { 0.9, 1.2 } ,
                 {  3.5, -2.3 }
@@ -99,6 +101,7 @@ public void testRosenbrock() {
                                  new ObjectiveFunction(rosenbrock),
                                  GoalType.MINIMIZE,
                                  simplex,
+                                 new NelderMeadTransform(),
                                  new InitialGuess(new double[] { -1.2, 1.0 }));
         Assert.assertEquals(nbStarts, optimizer.getOptima().length);
 
diff --git a/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/MultivariateFunctionMappingAdapterTest.java b/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/MultivariateFunctionMappingAdapterTest.java
index b79508bd34..ecc0617d8a 100644
--- a/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/MultivariateFunctionMappingAdapterTest.java
+++ b/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/MultivariateFunctionMappingAdapterTest.java
@@ -20,9 +20,9 @@
 import org.apache.commons.math4.legacy.optim.InitialGuess;
 import org.apache.commons.math4.legacy.optim.MaxEval;
 import org.apache.commons.math4.legacy.optim.PointValuePair;
-import org.apache.commons.math4.legacy.optim.nonlinear.scalar.noderiv.AbstractSimplex;
-import org.apache.commons.math4.legacy.optim.nonlinear.scalar.noderiv.NelderMeadSimplex;
 import org.apache.commons.math4.legacy.optim.nonlinear.scalar.noderiv.SimplexOptimizer;
+import org.apache.commons.math4.legacy.optim.nonlinear.scalar.noderiv.Simplex;
+import org.apache.commons.math4.legacy.optim.nonlinear.scalar.noderiv.NelderMeadTransform;
 import org.junit.Assert;
 import org.junit.Test;
 
@@ -36,7 +36,7 @@ public void testStartSimplexInsideRange() {
                                                      biQuadratic.getUpper());
 
         SimplexOptimizer optimizer = new SimplexOptimizer(1e-10, 1e-30);
-        final AbstractSimplex simplex = new NelderMeadSimplex(new double[][] {
+        final Simplex simplex = Simplex.of(new double[][] {
                 wrapped.boundedToUnbounded(new double[] { 1.5, 2.75 }),
                 wrapped.boundedToUnbounded(new double[] { 1.5, 2.95 }),
                 wrapped.boundedToUnbounded(new double[] { 1.7, 2.90 })
@@ -46,6 +46,7 @@ public void testStartSimplexInsideRange() {
             = optimizer.optimize(new MaxEval(300),
                                  new ObjectiveFunction(wrapped),
                                  simplex,
+                                 new NelderMeadTransform(),
                                  GoalType.MINIMIZE,
                                  new InitialGuess(wrapped.boundedToUnbounded(new double[] { 1.5, 2.25 })));
         final double[] bounded = wrapped.unboundedToBounded(optimum.getPoint());
@@ -63,7 +64,7 @@ public void testOptimumOutsideRange() {
                                                      biQuadratic.getUpper());
 
         SimplexOptimizer optimizer = new SimplexOptimizer(1e-10, 1e-30);
-        final AbstractSimplex simplex = new NelderMeadSimplex(new double[][] {
+        final Simplex simplex = Simplex.of(new double[][] {
                 wrapped.boundedToUnbounded(new double[] { 1.5, 2.75 }),
                 wrapped.boundedToUnbounded(new double[] { 1.5, 2.95 }),
                 wrapped.boundedToUnbounded(new double[] { 1.7, 2.90 })
@@ -73,6 +74,7 @@ public void testOptimumOutsideRange() {
             = optimizer.optimize(new MaxEval(100),
                                  new ObjectiveFunction(wrapped),
                                  simplex,
+                                 new NelderMeadTransform(),
                                  GoalType.MINIMIZE,
                                  new InitialGuess(wrapped.boundedToUnbounded(new double[] { 1.5, 2.25 })));
         final double[] bounded = wrapped.unboundedToBounded(optimum.getPoint());
@@ -92,7 +94,7 @@ public void testUnbounded() {
                                                      biQuadratic.getUpper());
 
         SimplexOptimizer optimizer = new SimplexOptimizer(1e-10, 1e-30);
-        final AbstractSimplex simplex = new NelderMeadSimplex(new double[][] {
+        final Simplex simplex = Simplex.of(new double[][] {
                 wrapped.boundedToUnbounded(new double[] { 1.5, 2.75 }),
                 wrapped.boundedToUnbounded(new double[] { 1.5, 2.95 }),
                 wrapped.boundedToUnbounded(new double[] { 1.7, 2.90 })
@@ -102,6 +104,7 @@ public void testUnbounded() {
             = optimizer.optimize(new MaxEval(300),
                                  new ObjectiveFunction(wrapped),
                                  simplex,
+                                 new NelderMeadTransform(),
                                  GoalType.MINIMIZE,
                                  new InitialGuess(wrapped.boundedToUnbounded(new double[] { 1.5, 2.25 })));
         final double[] bounded = wrapped.unboundedToBounded(optimum.getPoint());
@@ -121,7 +124,7 @@ public void testHalfBounded() {
                                                      biQuadratic.getUpper());
 
         SimplexOptimizer optimizer = new SimplexOptimizer(1e-13, 1e-30);
-        final AbstractSimplex simplex = new NelderMeadSimplex(new double[][] {
+        final Simplex simplex = Simplex.of(new double[][] {
                 wrapped.boundedToUnbounded(new double[] { 1.5, 2.75 }),
                 wrapped.boundedToUnbounded(new double[] { 1.5, 2.95 }),
                 wrapped.boundedToUnbounded(new double[] { 1.7, 2.90 })
@@ -131,6 +134,7 @@ public void testHalfBounded() {
             = optimizer.optimize(new MaxEval(200),
                                  new ObjectiveFunction(wrapped),
                                  simplex,
+                                 new NelderMeadTransform(),
                                  GoalType.MINIMIZE,
                                  new InitialGuess(wrapped.boundedToUnbounded(new double[] { 1.5, 2.25 })));
         final double[] bounded = wrapped.unboundedToBounded(optimum.getPoint());
diff --git a/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/MultivariateFunctionPenaltyAdapterTest.java b/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/MultivariateFunctionPenaltyAdapterTest.java
index 26ceb21f60..01ac25676a 100644
--- a/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/MultivariateFunctionPenaltyAdapterTest.java
+++ b/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/MultivariateFunctionPenaltyAdapterTest.java
@@ -21,9 +21,9 @@
 import org.apache.commons.math4.legacy.optim.MaxEval;
 import org.apache.commons.math4.legacy.optim.PointValuePair;
 import org.apache.commons.math4.legacy.optim.SimplePointChecker;
-import org.apache.commons.math4.legacy.optim.nonlinear.scalar.noderiv.AbstractSimplex;
-import org.apache.commons.math4.legacy.optim.nonlinear.scalar.noderiv.NelderMeadSimplex;
 import org.apache.commons.math4.legacy.optim.nonlinear.scalar.noderiv.SimplexOptimizer;
+import org.apache.commons.math4.legacy.optim.nonlinear.scalar.noderiv.Simplex;
+import org.apache.commons.math4.legacy.optim.nonlinear.scalar.noderiv.NelderMeadTransform;
 import org.junit.Assert;
 import org.junit.Test;
 
@@ -38,12 +38,12 @@ public void testStartSimplexInsideRange() {
                                                        1000.0, new double[] { 100.0, 100.0 });
 
         SimplexOptimizer optimizer = new SimplexOptimizer(1e-10, 1e-30);
-        final AbstractSimplex simplex = new NelderMeadSimplex(new double[] { 1.0, 0.5 });
 
         final PointValuePair optimum
             = optimizer.optimize(new MaxEval(300),
                                  new ObjectiveFunction(wrapped),
-                                 simplex,
+                                 Simplex.of(new double[] { 1.0, 0.5 }),
+                                 new NelderMeadTransform(),
                                  GoalType.MINIMIZE,
                                  new InitialGuess(new double[] { 1.5, 2.25 }));
 
@@ -61,12 +61,12 @@ public void testStartSimplexOutsideRange() {
                                                        1000.0, new double[] { 100.0, 100.0 });
 
         SimplexOptimizer optimizer = new SimplexOptimizer(1e-10, 1e-30);
-        final AbstractSimplex simplex = new NelderMeadSimplex(new double[] { 1.0, 0.5 });
 
         final PointValuePair optimum
             = optimizer.optimize(new MaxEval(300),
                                  new ObjectiveFunction(wrapped),
-                                 simplex,
+                                 Simplex.of(new double[] { 1.0, 0.5 }),
+                                 new NelderMeadTransform(),
                                  GoalType.MINIMIZE,
                                  new InitialGuess(new double[] { -1.5, 4.0 }));
 
@@ -84,12 +84,12 @@ public void testOptimumOutsideRange() {
                                                       1000.0, new double[] { 100.0, 100.0 });
 
         SimplexOptimizer optimizer = new SimplexOptimizer(new SimplePointChecker<PointValuePair>(1.0e-11, 1.0e-20));
-        final AbstractSimplex simplex = new NelderMeadSimplex(new double[] { 1.0, 0.5 });
 
         final PointValuePair optimum
             = optimizer.optimize(new MaxEval(600),
                                  new ObjectiveFunction(wrapped),
-                                 simplex,
+                                 Simplex.of(new double[] { 1.0, 0.5 }),
+                                 new NelderMeadTransform(),
                                  GoalType.MINIMIZE,
                                  new InitialGuess(new double[] { -1.5, 4.0 }));
 
@@ -109,12 +109,12 @@ public void testUnbounded() {
                                                      1000.0, new double[] { 100.0, 100.0 });
 
         SimplexOptimizer optimizer = new SimplexOptimizer(1e-10, 1e-30);
-        final AbstractSimplex simplex = new NelderMeadSimplex(new double[] { 1.0, 0.5 });
 
         final PointValuePair optimum
             = optimizer.optimize(new MaxEval(300),
                                  new ObjectiveFunction(wrapped),
-                                 simplex,
+                                 Simplex.of(new double[] { 1.0, 0.5 }),
+                                 new NelderMeadTransform(),
                                  GoalType.MINIMIZE,
                                  new InitialGuess(new double[] { -1.5, 4.0 }));
 
@@ -134,12 +134,12 @@ public void testHalfBounded() {
                                                        1000.0, new double[] { 100.0, 100.0 });
 
         SimplexOptimizer optimizer = new SimplexOptimizer(new SimplePointChecker<PointValuePair>(1.0e-10, 1.0e-20));
-        final AbstractSimplex simplex = new NelderMeadSimplex(new double[] { 1.0, 0.5 });
 
         final PointValuePair optimum
             = optimizer.optimize(new MaxEval(400),
                                  new ObjectiveFunction(wrapped),
-                                 simplex,
+                                 Simplex.of(new double[] { 1.0, 0.5 }),
+                                 new NelderMeadTransform(),
                                  GoalType.MINIMIZE,
                                  new InitialGuess(new double[] { -1.5, 4.0 }));
 
@@ -197,7 +197,5 @@ public double getBoundedXOptimum() {
         public double getBoundedYOptimum() {
             return (yOptimum < yMin) ? yMin : ((yOptimum > yMax) ? yMax : yOptimum);
         }
-
     }
-
 }
diff --git a/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/OptimTestUtils.java b/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/OptimTestUtils.java
index 5739e557ee..2f5d9cbfde 100644
--- a/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/OptimTestUtils.java
+++ b/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/OptimTestUtils.java
@@ -238,32 +238,17 @@ public double value(double[] x) {
         }
     }
 
+    // https://www.sfu.ca/~ssurjano/rastr.html
     static class Rastrigin implements MultivariateFunction {
-        private double axisratio;
-        private double amplitude;
-
-        Rastrigin() {
-            this(1, 10);
-        }
-
-        Rastrigin(double axisratio, double amplitude) {
-            this.axisratio = axisratio;
-            this.amplitude = amplitude;
-        }
-
+        private static final double A = 10;
         @Override
         public double value(double[] x) {
-            double f = 0;
-            double fac;
-            for (int i = 0; i < x.length; ++i) {
-                fac = AccurateMath.pow(axisratio, (i - 1.) / (x.length - 1.));
-                if (i == 0 && x[i] < 0) {
-                    fac *= 1.;
-                }
-                f += fac * fac * x[i] * x[i] + amplitude
-                    * (1. - AccurateMath.cos(2. * AccurateMath.PI * fac * x[i]));
+            double sum = 0;
+            for (int i = 0; i < x.length; i++) {
+                final double xi = x[i];
+                sum += xi * xi - A * AccurateMath.cos(2 * Math.PI * xi);
             }
-            return f;
+            return A * x.length + sum;
         }
     }
 
@@ -347,12 +332,23 @@ static double[] point(int n, double value) {
     static double[] point(int n,
                           double value,
                           double jitter) {
+        final double[] ds = new double[n];
+        Arrays.fill(ds, value);
+        return point(ds, jitter);
+    }
+
+    /**
+     * @param point Point.
+     * @param jitter Random noise to add to each {@code value} element.
+     */
+    static double[] point(double[] value,
+                          double jitter) {
         final ContinuousUniformSampler s = new ContinuousUniformSampler(rng(),
                                                                         -jitter,
                                                                         jitter);
-        final double[] ds = new double[n];
-        for (int i = 0; i < n; i++) {
-            ds[i] = value + s.sample();
+        final double[] ds = new double[value.length];
+        for (int i = 0; i < value.length; i++) {
+            ds[i] = value[i] + s.sample();
         }
         return ds;
     }
diff --git a/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/SimplexOptimizerMultiDirectionalTest.java b/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/SimplexOptimizerMultiDirectionalTest.java
index 4b16492ad7..571d42153f 100644
--- a/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/SimplexOptimizerMultiDirectionalTest.java
+++ b/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/SimplexOptimizerMultiDirectionalTest.java
@@ -17,18 +17,19 @@
 
 package org.apache.commons.math4.legacy.optim.nonlinear.scalar.noderiv;
 
+import java.util.Arrays;
 import org.apache.commons.math4.legacy.analysis.MultivariateFunction;
 import org.apache.commons.math4.legacy.exception.MathUnsupportedOperationException;
 import org.apache.commons.math4.legacy.optim.InitialGuess;
 import org.apache.commons.math4.legacy.optim.MaxEval;
 import org.apache.commons.math4.legacy.optim.PointValuePair;
 import org.apache.commons.math4.legacy.optim.SimpleBounds;
-import org.apache.commons.math4.legacy.optim.SimpleValueChecker;
 import org.apache.commons.math4.legacy.optim.nonlinear.scalar.GoalType;
 import org.apache.commons.math4.legacy.optim.nonlinear.scalar.ObjectiveFunction;
 import org.apache.commons.math4.legacy.core.MathArrays;
 import org.junit.Assert;
 import org.junit.Test;
+import org.junit.Ignore;
 
 public class SimplexOptimizerMultiDirectionalTest {
     private static final int DIM = 13;
@@ -42,139 +43,22 @@ public void testBoundsUnsupported() {
                            new ObjectiveFunction(fourExtrema),
                            GoalType.MINIMIZE,
                            new InitialGuess(new double[] { -3, 0 }),
-                           new NelderMeadSimplex(new double[] { 0.2, 0.2 }),
+                           Simplex.of(new double[] { 0.2, 0.2 }),
+                           new MultiDirectionalTransform(),
                            new SimpleBounds(new double[] { -5, -1 },
                                             new double[] { 5, 1 }));
     }
 
-    @Test
-    public void testMinimize1() {
-        SimplexOptimizer optimizer = new SimplexOptimizer(1e-11, 1e-30);
-        final OptimTestUtils.FourExtrema fourExtrema = new OptimTestUtils.FourExtrema();
-
-        final PointValuePair optimum
-            = optimizer.optimize(new MaxEval(200),
-                                 new ObjectiveFunction(fourExtrema),
-                                 GoalType.MINIMIZE,
-                                 new InitialGuess(new double[] { -3, 0 }),
-                                 new MultiDirectionalSimplex(new double[] { 0.2, 0.2 }));
-        Assert.assertEquals(fourExtrema.xM, optimum.getPoint()[0], 4e-6);
-        Assert.assertEquals(fourExtrema.yP, optimum.getPoint()[1], 3e-6);
-        Assert.assertEquals(fourExtrema.valueXmYp, optimum.getValue(), 8e-13);
-        Assert.assertTrue(optimizer.getEvaluations() > 120);
-        Assert.assertTrue(optimizer.getEvaluations() < 150);
-
-        // Check that the number of iterations is updated (MATH-949).
-        Assert.assertTrue(optimizer.getIterations() > 0);
-    }
-
-    @Test
-    public void testMinimize2() {
-        SimplexOptimizer optimizer = new SimplexOptimizer(1e-11, 1e-30);
-        final OptimTestUtils.FourExtrema fourExtrema = new OptimTestUtils.FourExtrema();
-
-        final PointValuePair optimum
-            = optimizer.optimize(new MaxEval(200),
-                                 new ObjectiveFunction(fourExtrema),
-                                 GoalType.MINIMIZE,
-                                 new InitialGuess(new double[] { 1, 0 }),
-                                 new MultiDirectionalSimplex(new double[] { 0.2, 0.2 }));
-        Assert.assertEquals(fourExtrema.xP, optimum.getPoint()[0], 2e-8);
-        Assert.assertEquals(fourExtrema.yM, optimum.getPoint()[1], 3e-6);
-        Assert.assertEquals(fourExtrema.valueXpYm, optimum.getValue(), 2e-12);
-        Assert.assertTrue(optimizer.getEvaluations() > 120);
-        Assert.assertTrue(optimizer.getEvaluations() < 150);
-
-        // Check that the number of iterations is updated (MATH-949).
-        Assert.assertTrue(optimizer.getIterations() > 0);
-    }
-
-    @Test
-    public void testMaximize1() {
-        SimplexOptimizer optimizer = new SimplexOptimizer(1e-11, 1e-30);
-        final OptimTestUtils.FourExtrema fourExtrema = new OptimTestUtils.FourExtrema();
-
-        final PointValuePair optimum
-            = optimizer.optimize(new MaxEval(200),
-                                 new ObjectiveFunction(fourExtrema),
-                                 GoalType.MAXIMIZE,
-                                 new InitialGuess(new double[] { -3.0, 0.0 }),
-                                 new MultiDirectionalSimplex(new double[] { 0.2, 0.2 }));
-        Assert.assertEquals(fourExtrema.xM, optimum.getPoint()[0], 7e-7);
-        Assert.assertEquals(fourExtrema.yM, optimum.getPoint()[1], 3e-7);
-        Assert.assertEquals(fourExtrema.valueXmYm, optimum.getValue(), 2e-14);
-        Assert.assertTrue(optimizer.getEvaluations() > 120);
-        Assert.assertTrue(optimizer.getEvaluations() < 150);
-
-        // Check that the number of iterations is updated (MATH-949).
-        Assert.assertTrue(optimizer.getIterations() > 0);
-    }
-
-    @Test
-    public void testMaximize2() {
-        SimplexOptimizer optimizer = new SimplexOptimizer(new SimpleValueChecker(1e-15, 1e-30));
-        final OptimTestUtils.FourExtrema fourExtrema = new OptimTestUtils.FourExtrema();
-
-        final PointValuePair optimum
-            = optimizer.optimize(new MaxEval(200),
-                                 new ObjectiveFunction(fourExtrema),
-                                 GoalType.MAXIMIZE,
-                                 new InitialGuess(new double[] { 1, 0 }),
-                                 new MultiDirectionalSimplex(new double[] { 0.2, 0.2 }));
-        Assert.assertEquals(fourExtrema.xP, optimum.getPoint()[0], 2e-8);
-        Assert.assertEquals(fourExtrema.yP, optimum.getPoint()[1], 3e-6);
-        Assert.assertEquals(fourExtrema.valueXpYp, optimum.getValue(), 2e-12);
-        Assert.assertTrue(optimizer.getEvaluations() > 180);
-        Assert.assertTrue(optimizer.getEvaluations() < 220);
-
-        // Check that the number of iterations is updated (MATH-949).
-        Assert.assertTrue(optimizer.getIterations() > 0);
-    }
-
-    @Test
-    public void testRosenbrock() {
-        final OptimTestUtils.Rosenbrock rosenbrock = new OptimTestUtils.Rosenbrock();
-        SimplexOptimizer optimizer = new SimplexOptimizer(-1, 1e-3);
-        PointValuePair optimum
-           = optimizer.optimize(new MaxEval(100),
-                                new ObjectiveFunction(rosenbrock),
-                                GoalType.MINIMIZE,
-                                new InitialGuess(new double[] { -1.2, 1 }),
-                                new MultiDirectionalSimplex(new double[][] {
-                                        { -1.2,  1.0 },
-                                        { 0.9, 1.2 },
-                                        {  3.5, -2.3 } }));
-        Assert.assertTrue(optimizer.getEvaluations() > 50);
-        Assert.assertTrue(optimizer.getEvaluations() < 100);
-        Assert.assertTrue(optimum.getValue() > 1e-2);
-    }
-
-    @Test
-    public void testPowell() {
-        final OptimTestUtils.Powell powell = new OptimTestUtils.Powell();
-        SimplexOptimizer optimizer = new SimplexOptimizer(-1, 1e-3);
-        PointValuePair optimum
-            = optimizer.optimize(new MaxEval(1000),
-                                 new ObjectiveFunction(powell),
-                                 GoalType.MINIMIZE,
-                                 new InitialGuess(new double[] { 3, -1, 0, 1 }),
-                                 new MultiDirectionalSimplex(4));
-        Assert.assertTrue(optimizer.getEvaluations() > 800);
-        Assert.assertTrue(optimizer.getEvaluations() < 900);
-        Assert.assertTrue(optimum.getValue() > 1e-2);
-    }
-
     @Test
     public void testMath283() {
-        // fails because MultiDirectional.iterateSimplex is looping forever
-        // the while(true) should be replaced with a convergence check
         SimplexOptimizer optimizer = new SimplexOptimizer(1e-14, 1e-14);
         final OptimTestUtils.Gaussian2D function = new OptimTestUtils.Gaussian2D(0, 0, 1);
         PointValuePair estimate = optimizer.optimize(new MaxEval(1000),
                                                      new ObjectiveFunction(function),
                                                      GoalType.MAXIMIZE,
                                                      new InitialGuess(function.getMaximumPosition()),
-                                                     new MultiDirectionalSimplex(2));
+                                                     Simplex.of(2),
+                                                     new MultiDirectionalTransform());
         final double EPSILON = 1e-5;
         final double expectedMaximum = function.getMaximum();
         final double actualMaximum = estimate.getValue();
@@ -186,92 +70,157 @@ public void testMath283() {
         Assert.assertEquals(expectedPosition[1], actualPosition[1], EPSILON );
     }
 
+    @Test
+    public void testFourExtremaMinimize1() {
+        final OptimTestUtils.FourExtrema f = new OptimTestUtils.FourExtrema();
+        doTest(f,
+               OptimTestUtils.point(new double[] {-3, 0}, 1e-1),
+               GoalType.MINIMIZE,
+               105,
+               Simplex.of(OptimTestUtils.point(2, 0.2, 1e-2)),
+               new PointValuePair(new double[] {f.xM, f.yP}, f.valueXmYp),
+               1e-6);
+    }
+    @Test
+    public void testFourExtremaMaximize1() {
+        final OptimTestUtils.FourExtrema f = new OptimTestUtils.FourExtrema();
+        doTest(f,
+               OptimTestUtils.point(new double[] {-3, 0}, 1e-1),
+               GoalType.MAXIMIZE,
+               100,
+               Simplex.of(OptimTestUtils.point(2, 0.2, 1e-2)),
+               new PointValuePair(new double[] {f.xM, f.yM}, f.valueXmYm),
+               1e-6);
+    }
+    @Test
+    public void testFourExtremaMinimize2() {
+        final OptimTestUtils.FourExtrema f = new OptimTestUtils.FourExtrema();
+        doTest(f,
+               OptimTestUtils.point(new double[] {1, 0}, 1e-1),
+               GoalType.MINIMIZE,
+               100,
+               Simplex.of(OptimTestUtils.point(2, 0.2, 1e-2)),
+               new PointValuePair(new double[] {f.xP, f.yM}, f.valueXpYm),
+               1e-6);
+    }
+    @Test
+    public void testFourExtremaMaximize2() {
+        final OptimTestUtils.FourExtrema f = new OptimTestUtils.FourExtrema();
+        doTest(f,
+               OptimTestUtils.point(new double[] {1, 0}, 1e-1),
+               GoalType.MAXIMIZE,
+               110,
+               Simplex.of(OptimTestUtils.point(2, 0.2, 1e-2)),
+               new PointValuePair(new double[] {f.xP, f.yP}, f.valueXpYp),
+               1e-6);
+    }
+
+    @Test
+    public void testRosenbrock() {
+        doTest(new OptimTestUtils.Rosenbrock(),
+               OptimTestUtils.point(new double[] { -1.2, 1 }, 1e-1),
+               GoalType.MINIMIZE,
+               190,
+               new PointValuePair(OptimTestUtils.point(2, 1.0), 0.0),
+               1e-4);
+    }
+
+    @Test
+    public void testPowell() {
+        doTest(new OptimTestUtils.Powell(),
+               OptimTestUtils.point(new double[] { 3, -1, 0, 1 }, 1e-1),
+               GoalType.MINIMIZE,
+               420,
+               new PointValuePair(OptimTestUtils.point(4, 0.0), 0.0),
+               2e-4);
+    }
+
     @Test
     public void testRosen() {
         doTest(new OptimTestUtils.Rosen(),
                OptimTestUtils.point(DIM, 0.1),
                GoalType.MINIMIZE,
-               183861,
+               186915,
                new PointValuePair(OptimTestUtils.point(DIM, 1.0), 0.0),
                1e-4);
     }
 
-    @Test
+    @Ignore@Test
     public void testEllipse() {
         doTest(new OptimTestUtils.Elli(),
-               OptimTestUtils.point(DIM, 1.0),
+               OptimTestUtils.point(DIM, 1.0, 1e-1),
                GoalType.MINIMIZE,
-               873,
+               911,
                new PointValuePair(OptimTestUtils.point(DIM, 0.0), 0.0),
-               1e-14);
+               1e-9);
      }
 
-    @Test
+    @Ignore@Test
     public void testElliRotated() {
         doTest(new OptimTestUtils.ElliRotated(),
-               OptimTestUtils.point(DIM, 1.0),
+               OptimTestUtils.point(DIM, 1.0, 1e-1),
                GoalType.MINIMIZE,
-               873,
+               911,
                new PointValuePair(OptimTestUtils.point(DIM, 0.0), 0.0),
-               1e-14);
+               1e-9);
     }
 
     @Test
     public void testCigar() {
         doTest(new OptimTestUtils.Cigar(),
-               OptimTestUtils.point(DIM, 1.0),
+               OptimTestUtils.point(DIM, 1.0, 1e-1),
                GoalType.MINIMIZE,
-               925,
+               7000,
                new PointValuePair(OptimTestUtils.point(DIM,0.0), 0.0),
-               1e-14);
+               1e-6);
     }
 
-    @Test
+    @Ignore@Test
     public void testTwoAxes() {
         doTest(new OptimTestUtils.TwoAxes(),
-               OptimTestUtils.point(DIM, 1.0),
+               OptimTestUtils.point(DIM, 1.0, 1e-1),
                GoalType.MINIMIZE,
-               1159,
+               1219,
                new PointValuePair(OptimTestUtils.point(DIM, 0.0), 0.0),
-               1e-14);
+               1e-12);
      }
 
-    @Test
+    @Ignore@Test
     public void testCigTab() {
         doTest(new OptimTestUtils.CigTab(),
-               OptimTestUtils.point(DIM, 1.0),
+               OptimTestUtils.point(DIM, 1.0, 1e-1),
                GoalType.MINIMIZE,
-               795,
+               827,
                new PointValuePair(OptimTestUtils.point(DIM, 0.0), 0.0),
-               1e-14);
+               1e-8);
      }
 
     @Test
     public void testSphere() {
         doTest(new OptimTestUtils.Sphere(),
-               OptimTestUtils.point(DIM, 1.0),
+               OptimTestUtils.point(DIM, 1.0, 1e-1),
                GoalType.MINIMIZE,
-               665,
+               2580,
                new PointValuePair(OptimTestUtils.point(DIM, 0.0), 0.0),
-               1e-14);
+               1e-6);
     }
 
-    @Test
+    @Ignore@Test
     public void testTablet() {
         doTest(new OptimTestUtils.Tablet(),
-               OptimTestUtils.point(DIM, 1.0),
+               OptimTestUtils.point(DIM, 1.0, 1e-1),
                GoalType.MINIMIZE,
-               873,
+               911,
                new PointValuePair(OptimTestUtils.point(DIM, 0.0), 0.0),
-               1e-14);
+               1e-9);
     }
 
-    @Test
+    @Ignore@Test
     public void testDiffPow() {
         doTest(new OptimTestUtils.DiffPow(),
-               OptimTestUtils.point(DIM, 1.0),
+               OptimTestUtils.point(DIM, 1.0, 1e-1),
                GoalType.MINIMIZE,
-               1000,
+               631,
                new PointValuePair(OptimTestUtils.point(DIM, 0.0), 0.0),
                1e-14);
     }
@@ -281,27 +230,27 @@ public void testSsDiffPow() {
         doTest(new OptimTestUtils.SsDiffPow(),
                OptimTestUtils.point(DIM / 2, 1.0, 1e-1),
                GoalType.MINIMIZE,
-               5000,
+               4000,
                new PointValuePair(OptimTestUtils.point(DIM / 2, 0.0), 0.0),
                1e-3);
     }
 
-    @Test
+    @Ignore@Test
     public void testAckley() {
         doTest(new OptimTestUtils.Ackley(),
                OptimTestUtils.point(DIM, 1.0, 1e-1),
                GoalType.MINIMIZE,
-               7500,
+               7900,
                new PointValuePair(OptimTestUtils.point(DIM, 0.0), 0d),
                1e-11);
     }
 
-    @Test
+    @Ignore@Test
     public void testRastrigin() {
         doTest(new OptimTestUtils.Rastrigin(),
-               OptimTestUtils.point(DIM, 1.0, 1e-1),
+               OptimTestUtils.point(DIM, 1.0, 2e-1),
                GoalType.MINIMIZE,
-               5000,
+               4600,
                new PointValuePair(OptimTestUtils.point(DIM, 0.0), 0d),
                1e-6);
     }
@@ -320,15 +269,53 @@ private void doTest(MultivariateFunction func,
                         int maxEvaluations,
                         PointValuePair expected,
                         double tol) {
-        final int dim = startPoint.length;
-        final SimplexOptimizer optim = new SimplexOptimizer(1e-10, 1e-12);
-        final PointValuePair result = optim.optimize(new MaxEval(maxEvaluations),
+        doTest(func,
+               startPoint,
+               goal,
+               maxEvaluations,
+               Simplex.of(startPoint.length, 1),
+               expected,
+               tol);
+    }
+
+    /**
+     * @param func Function to optimize.
+     * @param startPoint Starting point.
+     * @param goal Minimization or maximization.
+     * @param maxEvaluations Maximum number of evaluations.
+     * @param simplexSteps Initial simplex.
+     * @param expected Expected optimum.
+     * @param tol Tolerance for checking that the optimum is correct.
+     */
+    private void doTest(MultivariateFunction func,
+                        double[] startPoint,
+                        GoalType goal,
+                        int maxEvaluations,
+                        Simplex simplex,
+                        PointValuePair expected,
+                        double tol) {
+        final int maxEval = Math.max(maxEvaluations, 12000);
+        final SimplexOptimizer optim = new SimplexOptimizer(1e-13, 1e-14);
+        final PointValuePair result = optim.optimize(new MaxEval(maxEval),
                                                      new ObjectiveFunction(func),
                                                      goal,
                                                      new InitialGuess(startPoint),
-                                                     new MultiDirectionalSimplex(dim, 1));
+                                                     simplex,
+                                                     new NelderMeadTransform());
+        final String name = func.getClass().getSimpleName();
+
+        final double[] endPoint = result.getPoint();
+        final double funcValue = result.getValue();
+        Assert.assertEquals(name + ": value at " + Arrays.toString(endPoint),
+                            expected.getValue(),
+                            funcValue, 1e-2);
+
         final double dist = MathArrays.distance(expected.getPoint(),
-                                                result.getPoint());
-        Assert.assertEquals(0d, dist, tol);
+                                                endPoint);
+        Assert.assertEquals(name + ": distance to optimum", 0d, dist, tol);
+
+        final int nEval = optim.getEvaluations();
+        Assert.assertTrue(name + ": nEval=" + nEval,
+                          nEval < maxEvaluations);
     }
 }
diff --git a/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/SimplexOptimizerNelderMeadTest.java b/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/SimplexOptimizerNelderMeadTest.java
index 42d3d99bff..e619c32c15 100644
--- a/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/SimplexOptimizerNelderMeadTest.java
+++ b/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/SimplexOptimizerNelderMeadTest.java
@@ -17,7 +17,7 @@
 
 package org.apache.commons.math4.legacy.optim.nonlinear.scalar.noderiv;
 
-
+import java.util.Arrays;
 import org.apache.commons.math4.legacy.analysis.MultivariateFunction;
 import org.apache.commons.math4.legacy.analysis.MultivariateVectorFunction;
 import org.apache.commons.math4.legacy.exception.MathUnsupportedOperationException;
@@ -48,129 +48,12 @@ public void testBoundsUnsupported() {
                            new ObjectiveFunction(fourExtrema),
                            GoalType.MINIMIZE,
                            new InitialGuess(new double[] { -3, 0 }),
-                           new NelderMeadSimplex(new double[] { 0.2, 0.2 }),
+                           Simplex.of(new double[] { 0.2, 0.2 }),
+                           new NelderMeadTransform(),
                            new SimpleBounds(new double[] { -5, -1 },
                                             new double[] { 5, 1 }));
     }
 
-    @Test
-    public void testMinimize1() {
-        SimplexOptimizer optimizer = new SimplexOptimizer(1e-10, 1e-30);
-        final OptimTestUtils.FourExtrema fourExtrema = new OptimTestUtils.FourExtrema();
-
-        final PointValuePair optimum
-            = optimizer.optimize(new MaxEval(100),
-                                 new ObjectiveFunction(fourExtrema),
-                                 GoalType.MINIMIZE,
-                                 new InitialGuess(new double[] { -3, 0 }),
-                                 new NelderMeadSimplex(new double[] { 0.2, 0.2 }));
-        Assert.assertEquals(fourExtrema.xM, optimum.getPoint()[0], 2e-7);
-        Assert.assertEquals(fourExtrema.yP, optimum.getPoint()[1], 2e-5);
-        Assert.assertEquals(fourExtrema.valueXmYp, optimum.getValue(), 6e-12);
-        Assert.assertTrue(optimizer.getEvaluations() > 60);
-        Assert.assertTrue(optimizer.getEvaluations() < 90);
-
-        // Check that the number of iterations is updated (MATH-949).
-        Assert.assertTrue(optimizer.getIterations() > 0);
-    }
-
-    @Test
-    public void testMinimize2() {
-        SimplexOptimizer optimizer = new SimplexOptimizer(1e-10, 1e-30);
-        final OptimTestUtils.FourExtrema fourExtrema = new OptimTestUtils.FourExtrema();
-
-        final PointValuePair optimum
-            = optimizer.optimize(new MaxEval(100),
-                                 new ObjectiveFunction(fourExtrema),
-                                 GoalType.MINIMIZE,
-                                 new InitialGuess(new double[] { 1, 0 }),
-                                 new NelderMeadSimplex(new double[] { 0.2, 0.2 }));
-        Assert.assertEquals(fourExtrema.xP, optimum.getPoint()[0], 5e-6);
-        Assert.assertEquals(fourExtrema.yM, optimum.getPoint()[1], 6e-6);
-        Assert.assertEquals(fourExtrema.valueXpYm, optimum.getValue(), 1e-11);
-        Assert.assertTrue(optimizer.getEvaluations() > 60);
-        Assert.assertTrue(optimizer.getEvaluations() < 90);
-
-        // Check that the number of iterations is updated (MATH-949).
-        Assert.assertTrue(optimizer.getIterations() > 0);
-    }
-
-    @Test
-    public void testMaximize1() {
-        SimplexOptimizer optimizer = new SimplexOptimizer(1e-10, 1e-30);
-        final OptimTestUtils.FourExtrema fourExtrema = new OptimTestUtils.FourExtrema();
-
-        final PointValuePair optimum
-            = optimizer.optimize(new MaxEval(100),
-                                 new ObjectiveFunction(fourExtrema),
-                                 GoalType.MAXIMIZE,
-                                 new InitialGuess(new double[] { -3, 0 }),
-                                 new NelderMeadSimplex(new double[] { 0.2, 0.2 }));
-        Assert.assertEquals(fourExtrema.xM, optimum.getPoint()[0], 1e-5);
-        Assert.assertEquals(fourExtrema.yM, optimum.getPoint()[1], 3e-6);
-        Assert.assertEquals(fourExtrema.valueXmYm, optimum.getValue(), 3e-12);
-        Assert.assertTrue(optimizer.getEvaluations() > 60);
-        Assert.assertTrue(optimizer.getEvaluations() < 90);
-
-        // Check that the number of iterations is updated (MATH-949).
-        Assert.assertTrue(optimizer.getIterations() > 0);
-    }
-
-    @Test
-    public void testMaximize2() {
-        SimplexOptimizer optimizer = new SimplexOptimizer(1e-10, 1e-30);
-        final OptimTestUtils.FourExtrema fourExtrema = new OptimTestUtils.FourExtrema();
-
-        final PointValuePair optimum
-            = optimizer.optimize(new MaxEval(100),
-                                 new ObjectiveFunction(fourExtrema),
-                                 GoalType.MAXIMIZE,
-                                 new InitialGuess(new double[] { 1, 0 }),
-                                 new NelderMeadSimplex(new double[] { 0.2, 0.2 }));
-        Assert.assertEquals(fourExtrema.xP, optimum.getPoint()[0], 4e-6);
-        Assert.assertEquals(fourExtrema.yP, optimum.getPoint()[1], 5e-6);
-        Assert.assertEquals(fourExtrema.valueXpYp, optimum.getValue(), 7e-12);
-        Assert.assertTrue(optimizer.getEvaluations() > 60);
-        Assert.assertTrue(optimizer.getEvaluations() < 90);
-
-        // Check that the number of iterations is updated (MATH-949).
-        Assert.assertTrue(optimizer.getIterations() > 0);
-    }
-
-    @Test
-    public void testRosenbrock() {
-
-        OptimTestUtils.Rosenbrock rosenbrock = new OptimTestUtils.Rosenbrock();
-        SimplexOptimizer optimizer = new SimplexOptimizer(-1, 1e-3);
-        PointValuePair optimum
-        = optimizer.optimize(new MaxEval(100),
-                             new ObjectiveFunction(rosenbrock),
-                             GoalType.MINIMIZE,
-                             new InitialGuess(new double[] { -1.2, 1 }),
-                                new NelderMeadSimplex(new double[][] {
-                                        { -1.2,  1 },
-                                        { 0.9, 1.2 },
-                                        {  3.5, -2.3 } }));
-        Assert.assertTrue(optimizer.getEvaluations() > 40);
-        Assert.assertTrue(optimizer.getEvaluations() < 50);
-        Assert.assertTrue(optimum.getValue() < 8e-4);
-    }
-
-    @Test
-    public void testPowell() {
-        OptimTestUtils.Powell powell = new OptimTestUtils.Powell();
-        SimplexOptimizer optimizer = new SimplexOptimizer(-1, 1e-3);
-        PointValuePair optimum =
-            optimizer.optimize(new MaxEval(200),
-                               new ObjectiveFunction(powell),
-                               GoalType.MINIMIZE,
-                               new InitialGuess(new double[] { 3, -1, 0, 1 }),
-                               new NelderMeadSimplex(4));
-        Assert.assertTrue(optimizer.getEvaluations() > 110);
-        Assert.assertTrue(optimizer.getEvaluations() < 130);
-        Assert.assertTrue(optimum.getValue() < 2e-3);
-    }
-
     @Test
     public void testLeastSquares1() {
         final RealMatrix factors
@@ -190,11 +73,13 @@ public double[] value(double[] variables) {
                                new ObjectiveFunction(ls),
                                GoalType.MINIMIZE,
                                new InitialGuess(new double[] { 10, 10 }),
-                               new NelderMeadSimplex(2));
-        Assert.assertEquals( 2, optimum.getPointRef()[0], 3e-5);
+                               Simplex.of(2),
+                               new NelderMeadTransform());
+        Assert.assertEquals( 2, optimum.getPointRef()[0], 1e-3);
         Assert.assertEquals(-3, optimum.getPointRef()[1], 4e-4);
-        Assert.assertTrue(optimizer.getEvaluations() > 60);
-        Assert.assertTrue(optimizer.getEvaluations() < 80);
+        final int nEval = optimizer.getEvaluations();
+        Assert.assertTrue("nEval=" + nEval, nEval > 60);
+        Assert.assertTrue("nEval=" + nEval, nEval < 80);
         Assert.assertTrue(optimum.getValue() < 1.0e-6);
     }
 
@@ -217,11 +102,13 @@ public double[] value(double[] variables) {
                                new ObjectiveFunction(ls),
                                GoalType.MINIMIZE,
                                new InitialGuess(new double[] { 10, 10 }),
-                               new NelderMeadSimplex(2));
-        Assert.assertEquals( 2, optimum.getPointRef()[0], 5e-5);
+                               Simplex.of(2),
+                               new NelderMeadTransform());
+        Assert.assertEquals( 2, optimum.getPointRef()[0], 1e-4);
         Assert.assertEquals(-3, optimum.getPointRef()[1], 8e-4);
-        Assert.assertTrue(optimizer.getEvaluations() > 60);
-        Assert.assertTrue(optimizer.getEvaluations() < 80);
+        final int nEval = optimizer.getEvaluations();
+        Assert.assertTrue("nEval=" + nEval, nEval > 70);
+        Assert.assertTrue("nEval=" + nEval, nEval < 85);
         Assert.assertTrue(optimum.getValue() < 1e-6);
     }
 
@@ -246,11 +133,13 @@ public double[] value(double[] variables) {
                                  new ObjectiveFunction(ls),
                                  GoalType.MINIMIZE,
                                  new InitialGuess(new double[] { 10, 10 }),
-                                 new NelderMeadSimplex(2));
-        Assert.assertEquals( 2, optimum.getPointRef()[0], 2e-3);
-        Assert.assertEquals(-3, optimum.getPointRef()[1], 8e-4);
-        Assert.assertTrue(optimizer.getEvaluations() > 60);
-        Assert.assertTrue(optimizer.getEvaluations() < 80);
+                                 Simplex.of(2),
+                                 new NelderMeadTransform());
+        Assert.assertEquals( 2, optimum.getPointRef()[0], 1e-2);
+        Assert.assertEquals(-3, optimum.getPointRef()[1], 1e-2);
+        final int nEval = optimizer.getEvaluations();
+        Assert.assertTrue("nEval=" + nEval, nEval > 60);
+        Assert.assertTrue("nEval=" + nEval, nEval < 80);
         Assert.assertTrue(optimum.getValue() < 1e-6);
     }
 
@@ -262,7 +151,73 @@ public void testMaxIterations() {
                            new ObjectiveFunction(powell),
                            GoalType.MINIMIZE,
                            new InitialGuess(new double[] { 3, -1, 0, 1 }),
-                           new NelderMeadSimplex(4));
+                           Simplex.of(4),
+                           new NelderMeadTransform());
+    }
+
+    @Test
+    public void testFourExtremaMinimize1() {
+        final OptimTestUtils.FourExtrema f = new OptimTestUtils.FourExtrema();
+        doTest(f,
+               OptimTestUtils.point(new double[] {-3, 0}, 1e-1),
+               GoalType.MINIMIZE,
+               105,
+               Simplex.of(OptimTestUtils.point(2, 0.2, 1e-2)),
+               new PointValuePair(new double[] {f.xM, f.yP}, f.valueXmYp),
+               1e-6);
+    }
+    @Test
+    public void testFourExtremaMaximize1() {
+        final OptimTestUtils.FourExtrema f = new OptimTestUtils.FourExtrema();
+        doTest(f,
+               OptimTestUtils.point(new double[] {-3, 0}, 1e-1),
+               GoalType.MAXIMIZE,
+               100,
+               Simplex.of(OptimTestUtils.point(2, 0.2, 1e-2)),
+               new PointValuePair(new double[] {f.xM, f.yM}, f.valueXmYm),
+               1e-6);
+    }
+    @Test
+    public void testFourExtremaMinimize2() {
+        final OptimTestUtils.FourExtrema f = new OptimTestUtils.FourExtrema();
+        doTest(f,
+               OptimTestUtils.point(new double[] {1, 0}, 1e-1),
+               GoalType.MINIMIZE,
+               100,
+               Simplex.of(OptimTestUtils.point(2, 0.2, 1e-2)),
+               new PointValuePair(new double[] {f.xP, f.yM}, f.valueXpYm),
+               1e-6);
+    }
+    @Test
+    public void testFourExtremaMaximize2() {
+        final OptimTestUtils.FourExtrema f = new OptimTestUtils.FourExtrema();
+        doTest(f,
+               OptimTestUtils.point(new double[] {1, 0}, 1e-1),
+               GoalType.MAXIMIZE,
+               110,
+               Simplex.of(OptimTestUtils.point(2, 0.2, 1e-2)),
+               new PointValuePair(new double[] {f.xP, f.yP}, f.valueXpYp),
+               1e-6);
+    }
+
+    @Test
+    public void testRosenbrock() {
+        doTest(new OptimTestUtils.Rosenbrock(),
+               OptimTestUtils.point(new double[] { -1.2, 1 }, 1e-1),
+               GoalType.MINIMIZE,
+               180,
+               new PointValuePair(OptimTestUtils.point(2, 1.0), 0.0),
+               1e-4);
+    }
+
+    @Test
+    public void testPowell() {
+        doTest(new OptimTestUtils.Powell(),
+               OptimTestUtils.point(new double[] { 3, -1, 0, 1 }, 1e-1),
+               GoalType.MINIMIZE,
+               420,
+               new PointValuePair(OptimTestUtils.point(4, 0.0), 0.0),
+               2e-4);
     }
 
     @Test
@@ -275,19 +230,17 @@ public void testRosen() {
                1e-6);
     }
 
-    @Ignore
-    @Test
+    @Ignore@Test
     public void testEllipse() {
         doTest(new OptimTestUtils.Elli(),
                OptimTestUtils.point(DIM, 1.0, 2e-1),
                GoalType.MINIMIZE,
-               12000,
+               15000,
                new PointValuePair(OptimTestUtils.point(DIM, 0.0), 0.0),
                1e-6);
-     }
+    }
 
-    @Ignore
-    @Test
+    @Ignore@Test
     public void testElliRotated() {
         doTest(new OptimTestUtils.ElliRotated(),
                OptimTestUtils.point(DIM, 1.0, 1e-1),
@@ -307,8 +260,7 @@ public void testCigar() {
                1e-6);
     }
 
-    @Ignore
-    @Test
+    @Ignore@Test
     public void testTwoAxes() {
         doTest(new OptimTestUtils.TwoAxes(),
                OptimTestUtils.point(DIM, 1.0, 1e-1),
@@ -316,10 +268,9 @@ public void testTwoAxes() {
                3451,
                new PointValuePair(OptimTestUtils.point(DIM, 0.0), 0.0),
                1e-14);
-     }
+    }
 
-    @Ignore
-    @Test
+    @Ignore@Test
     public void testCigTab() {
         doTest(new OptimTestUtils.CigTab(),
                OptimTestUtils.point(DIM, 1.0, 1e-1),
@@ -327,7 +278,7 @@ public void testCigTab() {
                7000,
                new PointValuePair(OptimTestUtils.point(DIM, 0.0), 0.0),
                1e-4);
-     }
+    }
 
     @Test
     public void testSphere() {
@@ -339,8 +290,7 @@ public void testSphere() {
                1e-6);
     }
 
-    @Ignore
-    @Test
+    @Ignore@Test
     public void testTablet() {
         doTest(new OptimTestUtils.Tablet(),
                OptimTestUtils.point(DIM, 1.0, 1e-1),
@@ -350,8 +300,7 @@ public void testTablet() {
                1e-14);
     }
 
-    @Ignore
-    @Test
+    @Ignore@Test
     public void testDiffPow() {
         doTest(new OptimTestUtils.DiffPow(),
                OptimTestUtils.point(DIM, 1.0, 2e-1),
@@ -371,24 +320,22 @@ public void testSsDiffPow() {
                1e-3);
     }
 
-    @Ignore
-    @Test
+    @Ignore@Test
     public void testAckley() {
         doTest(new OptimTestUtils.Ackley(),
                OptimTestUtils.point(DIM, 1.0, 2e-1),
                GoalType.MINIMIZE,
-               10000,
+               7900,
                new PointValuePair(OptimTestUtils.point(DIM, 0.0), 0.0),
                1e-11);
     }
 
-    @Ignore
-    @Test
+    @Ignore@Test
     public void testRastrigin() {
         doTest(new OptimTestUtils.Rastrigin(),
                OptimTestUtils.point(DIM, 1.0, 2e-1),
                GoalType.MINIMIZE,
-               10000,
+               4600,
                new PointValuePair(OptimTestUtils.point(DIM, 0.0), 0.0),
                0);
     }
@@ -407,15 +354,53 @@ private void doTest(MultivariateFunction func,
                         int maxEvaluations,
                         PointValuePair expected,
                         double tol) {
-        final int dim = startPoint.length;
+        doTest(func,
+               startPoint,
+               goal,
+               maxEvaluations,
+               Simplex.of(startPoint.length, 1),
+               expected,
+               tol);
+    }
+
+    /**
+     * @param func Function to optimize.
+     * @param startPoint Starting point.
+     * @param goal Minimization or maximization.
+     * @param maxEvaluations Maximum number of evaluations.
+     * @param simplexSteps Initial simplex.
+     * @param expected Expected optimum.
+     * @param tol Tolerance for checking that the optimum is correct.
+     */
+    private void doTest(MultivariateFunction func,
+                        double[] startPoint,
+                        GoalType goal,
+                        int maxEvaluations,
+                        Simplex simplex,
+                        PointValuePair expected,
+                        double tol) {
+        final int maxEval = Math.max(maxEvaluations, 12000);
         final SimplexOptimizer optim = new SimplexOptimizer(1e-13, 1e-14);
-        final PointValuePair result = optim.optimize(new MaxEval(maxEvaluations),
+        final PointValuePair result = optim.optimize(new MaxEval(maxEval),
                                                      new ObjectiveFunction(func),
                                                      goal,
                                                      new InitialGuess(startPoint),
-                                                     new NelderMeadSimplex(dim, 1));
+                                                     simplex,
+                                                     new NelderMeadTransform());
+        final String name = func.getClass().getSimpleName();
+
+        final double[] endPoint = result.getPoint();
+        final double funcValue = result.getValue();
+        Assert.assertEquals(name + ": value at " + Arrays.toString(endPoint),
+                            expected.getValue(),
+                            funcValue, 1e-2);
+
         final double dist = MathArrays.distance(expected.getPoint(),
-                                                result.getPoint());
-        Assert.assertEquals(0d, dist, tol);
+                                                endPoint);
+        Assert.assertEquals(name + ": distance to optimum", 0d, dist, tol);
+
+        final int nEval = optim.getEvaluations();
+        Assert.assertTrue(name + ": nEval=" + nEval,
+                          nEval < maxEvaluations);
     }
 }
diff --git a/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/random/HaltonSequenceGeneratorTest.java b/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/quasirandom/HaltonSequenceGeneratorTest.java
similarity index 85%
rename from commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/random/HaltonSequenceGeneratorTest.java
rename to commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/quasirandom/HaltonSequenceGeneratorTest.java
index 4dd800b41b..084f3f223f 100644
--- a/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/random/HaltonSequenceGeneratorTest.java
+++ b/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/quasirandom/HaltonSequenceGeneratorTest.java
@@ -14,12 +14,13 @@
  * See the License for the specific language governing permissions and
  * limitations under the License.
  */
-package org.apache.commons.math4.legacy.random;
+package org.apache.commons.math4.legacy.quasirandom;
 
-import org.junit.Assert;
 import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
+import org.apache.commons.math4.legacy.exception.NotPositiveException;
 import org.apache.commons.math4.legacy.exception.NullArgumentException;
 import org.apache.commons.math4.legacy.exception.OutOfRangeException;
+import org.junit.Assert;
 import org.junit.Before;
 import org.junit.Test;
 
@@ -63,7 +64,6 @@ public void test3DReference() {
         for (int i = 0; i < referenceValues.length; i++) {
             double[] result = generator.get();
             Assert.assertArrayEquals(referenceValues[i], result, 1e-3);
-            Assert.assertEquals(i + 1, generator.getNextIndex());
         }
     }
 
@@ -73,7 +73,6 @@ public void test2DUnscrambledReference() {
         for (int i = 0; i < referenceValuesUnscrambled.length; i++) {
             double[] result = generator.get();
             Assert.assertArrayEquals(referenceValuesUnscrambled[i], result, 1e-3);
-            Assert.assertEquals(i + 1, generator.getNextIndex());
         }
     }
 
@@ -119,16 +118,27 @@ public void testConstructor2() throws Exception{
     }
 
     @Test
-    public void testSkip() {
-        double[] result = generator.skipTo(5);
+    public void testJump() {
+        LowDiscrepancySequence copyOfSeq = generator.jump(5);
+        double[] result = copyOfSeq.get();
         Assert.assertArrayEquals(referenceValues[5], result, 1e-3);
-        Assert.assertEquals(6, generator.getNextIndex());
+
 
         for (int i = 6; i < referenceValues.length; i++) {
-            result = generator.get();
+            result = copyOfSeq.get();
             Assert.assertArrayEquals(referenceValues[i], result, 1e-3);
-            Assert.assertEquals(i + 1, generator.getNextIndex());
         }
     }
 
+    @Test(expected = NotPositiveException.class)
+    public void testJumpNegativeIndex() {
+        LowDiscrepancySequence copyOfSeq = generator.jump(-5);
+    }
+
+    @Test
+    public void testFirstSupplying() {
+        LowDiscrepancySequence sequence = new HaltonSequenceGenerator(3);
+        Assert.assertArrayEquals(new double[]{0.0, 0.0, 0.0}, sequence.get(),1e-6);
+    }
+
 }
diff --git a/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/random/SobolSequenceGeneratorTest.java b/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/quasirandom/SobolSequenceGeneratorTest.java
similarity index 82%
rename from commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/random/SobolSequenceGeneratorTest.java
rename to commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/quasirandom/SobolSequenceGeneratorTest.java
index 11c45a30de..57be1376a1 100644
--- a/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/random/SobolSequenceGeneratorTest.java
+++ b/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/quasirandom/SobolSequenceGeneratorTest.java
@@ -14,15 +14,14 @@
  * See the License for the specific language governing permissions and
  * limitations under the License.
  */
-package org.apache.commons.math4.legacy.random;
-
-import org.junit.Assert;
-
-import java.io.InputStream;
+package org.apache.commons.math4.legacy.quasirandom;
 
+import org.apache.commons.math4.legacy.exception.NotPositiveException;
 import org.apache.commons.math4.legacy.exception.OutOfRangeException;
+import org.junit.Assert;
 import org.junit.Before;
 import org.junit.Test;
+import java.io.InputStream;
 
 public class SobolSequenceGeneratorTest {
     private static final String RESOURCE_NAME = "/assets/org/apache/commons/math4/legacy/random/new-joe-kuo-6.21201";
@@ -52,7 +51,6 @@ public void test3DReference() {
         for (int i = 0; i < referenceValues.length; i++) {
             double[] result = generator.get();
             Assert.assertArrayEquals(referenceValues[i], result, 1e-6);
-            Assert.assertEquals(i + 1, generator.getNextIndex());
         }
     }
 
@@ -92,16 +90,26 @@ public void testConstructor2() throws Exception{
     }
 
     @Test
-    public void testSkip() {
-        double[] result = generator.skipTo(5);
+    public void testJump() {
+        LowDiscrepancySequence copyOfSeq = generator.jump(5);
+        double[] result = copyOfSeq.get();
         Assert.assertArrayEquals(referenceValues[5], result, 1e-6);
-        Assert.assertEquals(6, generator.getNextIndex());
 
         for (int i = 6; i < referenceValues.length; i++) {
-            result = generator.get();
+            result = copyOfSeq.get();
             Assert.assertArrayEquals(referenceValues[i], result, 1e-6);
-            Assert.assertEquals(i + 1, generator.getNextIndex());
         }
     }
 
+    @Test(expected = NotPositiveException.class)
+    public void testJumpNegativeIndex() {
+        LowDiscrepancySequence copyOfSeq = generator.jump(-5);
+    }
+
+    @Test
+    public void testFirstSupplying() {
+        LowDiscrepancySequence sequence = new SobolSequenceGenerator(3);
+        Assert.assertArrayEquals(new double[]{0.0, 0.0, 0.0}, sequence.get(),1e-6);
+    }
+
 }
diff --git a/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/random/StableRandomGeneratorTest.java b/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/random/StableRandomGeneratorTest.java
deleted file mode 100644
index 4f987ebc5d..0000000000
--- a/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/random/StableRandomGeneratorTest.java
+++ /dev/null
@@ -1,130 +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.legacy.random;
-
-import org.apache.commons.math4.legacy.exception.OutOfRangeException;
-import org.apache.commons.rng.UniformRandomProvider;
-import org.apache.commons.rng.simple.RandomSource;
-import org.apache.commons.math4.legacy.stat.StatUtils;
-import org.apache.commons.math4.legacy.stat.descriptive.DescriptiveStatistics;
-import org.junit.Assert;
-import org.junit.Test;
-
-/**
- * Tests for the class {@link StableRandomGenerator}.
- */
-public class StableRandomGeneratorTest {
-    private final UniformRandomProvider rg = RandomSource.create(RandomSource.WELL_19937_C, 100);
-    private static final int sampleSize = 10000;
-
-    /**
-     * Run the double nextDouble() method test Due to leptokurtic property the
-     * acceptance range is widened.
-     *
-     * TODO: verify that tolerance this wide is really OK
-     */
-    @Test
-    public void testNextDouble() {
-        StableRandomGenerator generator = new StableRandomGenerator(rg, 1.3,
-                0.1);
-        double[] sample = new double[2 * sampleSize];
-        for (int i = 0; i < sample.length; ++i) {
-            sample[i] = generator.nextNormalizedDouble();
-        }
-        Assert.assertEquals(0.0, StatUtils.mean(sample), 0.3);
-    }
-
-    /**
-     * If alpha = 2, than it must be Gaussian distribution
-     */
-    @Test
-    public void testGaussianCase() {
-        StableRandomGenerator generator = new StableRandomGenerator(rg, 2d, 0.0);
-
-        double[] sample = new double[sampleSize];
-        for (int i = 0; i < sample.length; ++i) {
-            sample[i] = generator.nextNormalizedDouble();
-        }
-        Assert.assertEquals(0.0, StatUtils.mean(sample), 0.02);
-        Assert.assertEquals(1.0, StatUtils.variance(sample), 0.02);
-    }
-
-    /**
-     * If alpha = 1, than it must be Cauchy distribution
-     */
-    @Test
-    public void testCauchyCase() {
-        StableRandomGenerator generator = new StableRandomGenerator(rg, 1d, 0.0);
-        DescriptiveStatistics summary = new DescriptiveStatistics();
-
-        for (int i = 0; i < sampleSize; ++i) {
-            double sample = generator.nextNormalizedDouble();
-            summary.addValue(sample);
-        }
-
-        // Standard Cauchy distribution should have zero median and mode
-        double median = summary.getPercentile(50);
-        Assert.assertEquals(0.0, median, 0.2);
-    }
-
-    /**
-     * Input parameter range tests
-     */
-    @Test
-    public void testAlphaRangeBelowZero() {
-        try {
-            new StableRandomGenerator(rg,
-                    -1.0, 0.0);
-            Assert.fail("Expected OutOfRangeException");
-        } catch (OutOfRangeException e) {
-            Assert.assertEquals(-1.0, e.getArgument());
-        }
-    }
-
-    @Test
-    public void testAlphaRangeAboveTwo() {
-        try {
-            new StableRandomGenerator(rg,
-                    3.0, 0.0);
-            Assert.fail("Expected OutOfRangeException");
-        } catch (OutOfRangeException e) {
-            Assert.assertEquals(3.0, e.getArgument());
-        }
-    }
-
-    @Test
-    public void testBetaRangeBelowMinusOne() {
-        try {
-            new StableRandomGenerator(rg,
-                    1.0, -2.0);
-            Assert.fail("Expected OutOfRangeException");
-        } catch (OutOfRangeException e) {
-            Assert.assertEquals(-2.0, e.getArgument());
-        }
-    }
-
-    @Test
-    public void testBetaRangeAboveOne() {
-        try {
-            new StableRandomGenerator(rg,
-                    1.0, 2.0);
-            Assert.fail("Expected OutOfRangeException");
-        } catch (OutOfRangeException e) {
-            Assert.assertEquals(2.0, e.getArgument());
-        }
-    }
-}
diff --git a/src/changes/changes.xml b/src/changes/changes.xml
index ddee0b76ff..7cb78f26fa 100644
--- a/src/changes/changes.xml
+++ b/src/changes/changes.xml
@@ -57,7 +57,7 @@ If the output is not quite correct, check for invisible trailing spaces!
 This is a major release that requires source changes in applications.
 
 Lower-level functionality has been split off to new components, often
-  with many bug-fixes, performance enhancements and improved API:
+  with bug-fixes, performance enhancements and improved API:
   Commons RNG,
   Commons Numbers,
   Commons Geometry,
@@ -75,8 +75,7 @@ The minimum version of the Java platform required to compile and use
 Users are encouraged to upgrade to this version, and to the new
  components listed above (they are now Commons Math's dependencies).
  Note that the new artefacts can be used together with v3.6.1 as the
- Java packages are different; ugrading can thus be performed
- incrementally.
+ Java packages are different.
 For codes that are now located in the 'commons-math-legacy' module,
  upgrading will be a matter of changing the top-level package name:
  From 'org.apache.commons.math3' to 'org.apache.commons.math4.legacy'.
@@ -85,6 +84,10 @@ For codes that are now located in the 'commons-math-legacy' module,
  and package 'org.apache.commons.math3.util' (e.g. 'MathArrays') over
  to package 'org.apache.commons.math4.legacy.core'.
 
+Notes:
+ Class 'FastMath' has been renamed 'AccurateMath' (now located in
+ package 'org.apache.commons.math4.legacy.core.jdkmath').
+
 Caveat:
  Although hundreds of issues have been fixed (see 'git log' and the
  reports in the bug-tracking system of all the listed components),
@@ -93,6 +96,12 @@ Caveat:
  nightmare was one of the main reasons for creating more focused
  components.]
 ">
+      <action dev="erans" type="update" issue="MATH-1615">
+        Remove class "StableRandomGenerator" (in package "o.a.c.m.legacy.random").
+      </action>
+      <action dev="erans" type="update" issue="MATH-1614">
+        Refactor "SimplexOptimizer" (in package "o.a.c.m.legacy.optim.nonlinear.scalar.noderiv").
+      </action>
       <action dev="erans" type="fix" issue="MATH-1613">
         Decrease maximum number of iterations in "SimpsonIntegrator" and
         "TrapezoidIntegrator" (due to the available range of the internal