}

    /**
     * Returns curvature matrix.
     *
     * @return curvature matrix.
     */
    public Matrix getAlpha() {
        return alpha;
    }

    /**
     * Returns degrees of freedom of computed chi square value.
     * Degrees of freedom is equal to the number of sampled data minus the
     * number of estimated parameters.
     *
     * @return degrees of freedom of computed chi square value.
     */
    public int getChisqDegreesOfFreedom() {
        return ndat - ma;
    }

    /**
     * Gets mean square error produced by estimated parameters respect to
     * provided sample data.
     *
     * @return mean square error.
     */
    public double getMse() {
        return mse;
    }

    /**
     * Gets the probability of finding a smaller chi square value.
     * The smaller the found chi square value is, the better the fit of the estimated
     * parameters to the actual parameter.
     * Thus, the smaller the chance of finding a smaller chi square value, then the
     * better the estimated fit is.
     *
     * @return probability of finding a smaller chi square value (better fit), expressed
     * as a value between 0.0 and 1.0.
     * @throws MaxIterationsExceededException if convergence of incomplete
     *                                        gamma function cannot be reached. This is rarely thrown and happens
     *                                        usually for numerically unstable input values.
     */
    public double getP() throws MaxIterationsExceededException {
        return ChiSqDist.cdf(getChisq(), getChisqDegreesOfFreedom());
    }

    /**
     * Gets a measure of quality of estimated fit as a value between 0.0 and 1.0.
     * The larger the quality value is, the better the fit that has been estimated.
     *
     * @return measure of quality of estimated fit.
     * @throws MaxIterationsExceededException if convergence of incomplete
     *                                        gamma function cannot be reached. This is rarely thrown and happens
     *                                        usually for numerically unstable input values.
     */
    public double getQ() throws MaxIterationsExceededException {
        return 1.0 - getP();
    }

    /**
     * Indicates whether covariance must be adjusted or not.
     * When covariance adjustment is enabled, then covariance is recomputed taking
     * into account input samples, input standard deviations of the samples and
     * jacobians of the model function overestimated parameters using the following
     * expression: Cov = (J'*W*J)^-1 where:
     * Cov is the covariance of estimated parameters
     * J is a matrix containing the Jacobians of the function overestimated parameters
     * for each input parameter x. Each row of J matrix contains an evaluation of
     * the model function Jacobian for i-th input parameter x. Each column of J matrix
     * contains the partial derivative of model function over j-th estimated parameter.
     * W is the inverse of input variances. It's a diagonal matrix containing the
     * reciprocal of the input variances (squared input standard deviations). That is:
     * W = diag(w) where k element of w is wk = 1 / sigmak^2, which corresponds to
     * the k-th standard deviation of input sample k.
     * By default, covariance is adjusted after fitting finishes.
     * <a href="http://people.duke.edu/~hpgavin/ce281/lm.pdf">http://people.duke.edu/~hpgavin/ce281/lm.pdf</a>
     * <a href="https://www8.cs.umu.se/kurser/5DA001/HT07/lectures/lsq-handouts.pdf">https://www8.cs.umu.se/kurser/5DA001/HT07/lectures/lsq-handouts.pdf</a>
     * Numerical Recipes 3rd Ed, page 812
     *
     * @return true if covariance must be adjusted, false otherwise.
     */
    public boolean isCovarianceAdjusted() {
        return adjustCovariance;
    }

    /**
     * Specifies whether covariance must be adjusted or not.
     * When covariance adjustment is enabled, then covariance is recomputed taking
     * into account input samples, input standard deviations of the samples and
     * jacobians of the model function overestimated parameters using the following
     * expression: Cov = (J'*W*J)^-1 where:
     * Cov is the covariance of estimated parameters
     * J is a matrix containing the Jacobians of the function overestimated parameters
     * for each input parameter x. Each row of J matrix contains an evaluation of
     * the model function Jacobian for i-th input parameter x. Each column of J matrix
     * contains the partial derivative of model function over j-th estimated parameter.
     * W is the inverse of input variances. It's a diagonal matrix containing the
     * reciprocal of the input variances (squared input standard deviations). That is:
     * W = diag(w) where k element of w is wk = 1 / sigmak^2, which corresponds to
     * the k-th standard deviation of input sample k.
     * By default, covariance is adjusted after fitting finishes.
     *
     * @param adjustCovariance true if covariance must be adjusted, false otherwise.
     */
    public void setCovarianceAdjusted(final boolean adjustCovariance) {
        this.adjustCovariance = adjustCovariance;
    }

    /**
     * Fits a function to provided data so that parameters associated to that
     * function can be estimated along with their covariance matrix and chi
     * square value.
     * If chi square value is close to 1, the fit is usually good.
     * If it is much larger, then error cannot be properly fitted.
     * If it is close to zero, then the model over-fits the error.
     * Methods {@link #getP()} and {@link #getQ()} can also be used to determine
     * the quality of the fit.
     *
     * @throws FittingException  if fitting fails.
     * @throws NotReadyException if enough input data has not yet been provided.
     */
    @Override
    public void fit() throws FittingException, NotReadyException {
        if (!isReady()) {
            throw new NotReadyException();
        }

        try {
            resultAvailable = false;

            int j;
            int k;
            int l;
            int iter;
            int done = 0;
            double alamda = 0.001;
            double ochisq;
            final var atry = new double[ma];
            final var beta = new double[ma];
            final var da = new double[ma];

            // number of parameters to be fitted
            mfit = 0;
            for (j = 0; j < ma; j++) {
                if (ia[j]) {
                    mfit++;
                }
            }

            final var oneda = new Matrix(mfit, 1);
            final var temp = new Matrix(mfit, mfit);

            // initialization
            mrqcof(a, alpha, beta);
            System.arraycopy(a, 0, atry, 0, ma);

            ochisq = chisq;
            for (iter = 0; iter < itmax; iter++) {

                if (done == ndone) {
                    // last pass. Use zero alamda
                    alamda = 0.0;
                }

                for (j = 0; j < mfit; j++) {
                    // alter linearized fitting matrix, by augmenting diagonal
                    // elements
                    for (k = 0; k < mfit; k++) {
                        covar.setElementAt(j, k, alpha.getElementAt(j, k));
                    }
                    covar.setElementAt(j, j, alpha.getElementAt(j, j) * (1.0 + alamda));
                    for (k = 0; k < mfit; k++) {
                        temp.setElementAt(j, k, covar.getElementAt(j, k));
                    }
                    oneda.setElementAt(j, 0, beta[j]);
                }

                // matrix solution
                GaussJordanElimination.process(temp, oneda);

                for (j = 0; j < mfit; j++) {
                    for (k = 0; k < mfit; k++) {
                        covar.setElementAt(j, k, temp.getElementAt(j, k));
                    }
                    da[j] = oneda.getElementAt(j, 0);
                }

                if (done == ndone) {
                    // Converged. Clean up and return
                    covsrt(covar);
                    covsrt(alpha);

                    if (adjustCovariance) {
                        adjustCovariance();
                    }

                    resultAvailable = true;

                    return;
                }

                // did the trial succeed?
                for (j = 0, l = 0; l < ma; l++) {
                    if (ia[l]) {
                        atry[l] = a[l] + da[j++];
                    }
                }

                mrqcof(atry, covar, da);
                if (Math.abs(chisq - ochisq) < Math.max(tol, tol * chisq)) {
                    done++;
                }

                if (chisq < ochisq) {
                    // success, accept the new solution
                    alamda *= 0.1;
                    ochisq = chisq;
                    for (j = 0; j < mfit; j++) {
                        for (k = 0; k < mfit; k++) {
                            alpha.setElementAt(j, k, covar.getElementAt(j, k));
                        }
                        beta[j] = da[j];
                    }
                    System.arraycopy(atry, 0, a, 0, ma);
                } else {
                    // failure, increase alamda
                    alamda *= 10.0;
                    chisq = ochisq;
                }
            }

            // too many iterations
            throw new FittingException("too many iterations");

        } catch (final AlgebraException | EvaluationException e) {
            throw new FittingException(e);
        }
    }

    /**
     * Prevents parameter at position i of linear combination of basis functions
     * to be modified during function fitting.
     *
     * @param i   position of parameter to be retained.
     * @param val value to be set for parameter at position i.
     */
    public void hold(final int i, final double val) {
        ia[i] = false;
        a[i] = val;
    }

    /**
     * Releases parameter at position i of linear combination of basis functions,
     * so it can be modified again if needed.
     *
     * @param i position of parameter to be released.
     */
    public void free(final int i) {
        ia[i] = true;
    }

    /**
     * Adjusts covariance.
     * Covariance must be adjusted to produce more real results close to the scale
     * of problem, otherwise estimated covariance will just be a measure of
     * goodness similar to chi square value because it will be the inverse of
     * the curvature matrix, which is just a solution of the covariance up to scale.
     * <p>
     * Covariance is adjusted taking into account input samples, input standard
     * deviations of the samples and jacobians of the model function overestimated
     * parameters using the following expression: Cov = (J'*W*J)^-1 where:
     * Cov is the covariance of estimated parameters
     * J is a matrix containing the Jacobians of the function overestimated parameters
     * for each input parameter x. Each row of J matrix contains an evaluation of
     * the model function Jacobian for i-th input parameter x. Each column of J matrix
     * contains the partial derivative of model function over j-th estimated parameter.
     * W is the inverse of input variances. It's a diagonal matrix containing the
     * reciprocal of the input variances (squared input standard deviations). That is:
     * W = diag(w) where k element of w is wk = 1 / sigmak^2, which corresponds to
     * the k-th standard deviation of input sample k.
     *
     * @throws AlgebraException    if there are numerical instabilities.
     * @throws EvaluationException if function evaluation fails.
     */
    private void adjustCovariance() throws AlgebraException, EvaluationException {

        final var nVars = evaluator.getNumberOfVariables();

public void minimize() throws LockedException, NotReadyException, OptimizationException {

        if (isLocked()) {
            throw new LockedException();
        }
        if (!isReady()) {
            throw new NotReadyException();
        }

        locked = true;

        final var n = p.length;

        // set vector of directions
        if (!isDirectionAvailable()) {
            xi = new double[n];
        } else {
            if (xi.length != n) {
                xi = new double[n];
            }
        }

        var validResult = false;
        try {
            double gg;
            double dgg;

            final var g = new double[n];
            final var h = new double[n];

            var fp = listener.evaluate(p);
            gradientListener.evaluateGradient(p, xi);

            for (var j = 0; j < n; j++) {
                g[j] = -xi[j];
                h[j] = g[j];
                xi[j] = h[j];
            }
            for (var its = 0; its < ITMAX; its++) {
                iter = its;
                fret = linmin();
                if (2.0 * Math.abs(fret - fp) <= tolerance * (Math.abs(fret) + Math.abs(fp) + EPS)) {
                    // minimum found
                    validResult = true;

                    if (iterationCompletedListener != null) {
                        iterationCompletedListener.onIterationCompleted(this, its, ITMAX);
                    }
                    break;
                }

                fp = fret;

                gradientListener.evaluateGradient(p, xi);

                var test = 0.0;
                final var den = Math.max(Math.abs(fp), 1.0);
                for (var j = 0; j < n; j++) {
                    final var temp = Math.abs(xi[j]) * Math.max(Math.abs(p[j]), 1.0) / den;

                    if (temp > test) {
                        test = temp;
                    }
                }
                if (test < GTOL) {
                    // minimum found
                    validResult = true;

                    if (iterationCompletedListener != null) {
                        iterationCompletedListener.onIterationCompleted(this, its, ITMAX);
                    }
                    break;
                }

                dgg = gg = 0.0;
                for (var j = 0; j < n; j++) {
                    gg += g[j] * g[j];

                    if (isPolakRibiereEnabled()) {
                        // This statement for Polak-Ribiere
                        dgg += (xi[j] + g[j]) * xi[j];

final int index = alpha.getIndex(j, k++);
                            alphaBuffer[index] += wt * dyda[m];
                        }
                    }
                    beta[j++] += dy * wt;
                }
            }

            // add to mse
            mse += dy * dy / Math.abs(degreesOfFreedom);

            // and find chi square
            chisq += dy * dy * sig2i / degreesOfFreedom;
        }

        // fill in the symmetric side
        for (j = 1; j < mfit; j++) {
            for (k = 0; k < j; k++) {
                alpha.setElementAt(k, j, alpha.getElementAt(j, k));
            }
        }
    }

    /**
     * Expand in storage the covariance matrix covar, to take into account
     * parameters that are being held fixed. (For the latter, return zero
     * covariances).
     *
     * @param covar covariance matrix.
     */
    private void covsrt(final Matrix covar) {
        int i;
        int j;
        int k;
        for (i = mfit; i < ma; i++) {
            for (j = 0; j < i + 1; j++) {
                covar.setElementAt(i, j, 0.0);
                covar.setElementAt(j, i, 0.0);
            }
        }

        k = mfit - 1;
        for (j = ma - 1; j >= 0; j--) {
            if (ia[j]) {
                final var buffer = covar.getBuffer();
                for (i = 0; i < ma; i++) {
                    final var pos1 = covar.getIndex(i, k);
                    final var pos2 = covar.getIndex(i, j);
                    swap(buffer, buffer, pos1, pos2);
                }
                for (i = 0; i < ma; i++) {
                    final var pos1 = covar.getIndex(k, i);
                    final var pos2 = covar.getIndex(j, i);
                    swap(buffer, buffer, pos1, pos2);
                }

                k--;
            }
        }
    }

    /**
     * Swaps values of arrays at provided positions.
     *
     * @param array1 1st array.
     * @param array2 2nd array.
     * @param pos1   1st position.
     * @param pos2   2nd position.
     */
    private static void swap(final double[] array1, final double[] array2, final int pos1, final int pos2) {
        final var value1 = array1[pos1];
        final var value2 = array2[pos2];
        array1[pos1] = value2;
        array2[pos2] = value1;
    }
}

var length = array.length;
        for (var i = 0; i < length / 2; i++) {
            var temp = array[i];
            var pos = length - 1 - i;
            array[i] = array[pos];
            array[pos] = temp;
        }
    }

    /**
     * Computes number of required iterations to achieve required confidence
     * with current probability of inlier and sample subset size.
     *
     * @param probInlier probability of inlier.
     * @param subsetSize sample subset size.
     * @param confidence required confidence of result.
     * @return number of required iterations.
     */
    private static int computeIterations(final double probInlier, final int subsetSize, final double confidence) {

        // compute number of times the algorithm needs to be executed depending
        // on number of inliers respect total points to achieve with probability
        // confidence that we have all inliers and probability 1 - confidence
        // that we have some outliers
        final var probSubsetAllInliers = Math.pow(probInlier, subsetSize);
        if (Math.abs(probSubsetAllInliers) < Double.MIN_VALUE || Double.isNaN(probSubsetAllInliers)) {
            return Integer.MAX_VALUE;
        } else {
            final var logProbSomeOutliers = Math.log(1.0 - probSubsetAllInliers);
            if (Math.abs(logProbSomeOutliers) < Double.MIN_VALUE || Double.isNaN(logProbSomeOutliers)) {
                return Integer.MAX_VALUE;
            } else {
                return (int) Math.ceil(Math.abs(Math.log(1.0 - confidence) / logProbSomeOutliers));
            }
        }
    }

    /**
     * Non randomness states that i-m (where i is the cardinal of the set of
     * inliers for a wrong model) follows the binomial distribution B(n,beta).
     * For n big enough, B(n,beta) approximates to normal distribution N(mu,
     * sigma^2) by the central limit theorem, with mu = n*beta and sigma =
     * sqrt(n*beta*(1 - beta)).
     * Psi, the probability that In_star out of n_star data points are by chance
     * inliers to an arbitrary incorrect model, is set to 0.05 (5%, as in the
     * original paper), and you must change the Chi2 value if you chose a
     * different value for psi.
     *
     * @param subsetSize sample subset size.
     * @param sampleSize total number of samples.
     * @param beta       beta value.
     * @return i-m.
     */
    private static int imin(final int subsetSize, final int sampleSize, final double beta) {
        final var mu = sampleSize * beta;
        final var sigma = Math.sqrt(sampleSize * beta * (1.0 - beta));

        return (int) Math.ceil(subsetSize + mu + sigma * Math.sqrt(CHI_SQUARED));
    }

thresh = (tol > 0. ? tol * svd.getSingularValues()[0] : -1.0);
            svd.solve(b, thresh, a);
            chisq = 0.0;
            for (i = 0; i < ndat; i++) {
                sum = 0.0;
                for (j = 0; j < ma; j++) {
                    sum += aa.getElementAt(i, j) * a[j];
                }
                chisq += Math.pow(sum - b[i], 2.0);
            }
            for (i = 0; i < ma; i++) {
                for (j = 0; j < i + 1; j++) {
                    sum = 0.0;
                    final var w = svd.getSingularValues();
                    final var tsh = svd.getNegligibleSingularValueThreshold();
                    final var v = svd.getV();
                    for (k = 0; k < ma; k++) {
                        if (w[k] > tsh) {
                            sum += v.getElementAt(i, k) * v.getElementAt(j, k) / Math.pow(w[k], 2.0);
                        }
                    }
                    covar.setElementAt(j, i, sum);
                    covar.setElementAt(i, j, sum);
                }
            }

            resultAvailable = true;

        } catch (final AlgebraException | EvaluationException e) {
            throw new FittingException(e);
        }
    }
}

dgg += (xi[j] + g[j]) * xi[j];
                    } else {
                        // This statement for Fletcher-Reeves
                        dgg += xi[j] * xi[j];
                    }
                }

                if (gg == 0.0) {
                    // minimum found
                    validResult = true;

                    if (iterationCompletedListener != null) {
                        iterationCompletedListener.onIterationCompleted(this, its, ITMAX);
                    }
                    break;
                }

                final var gam = dgg / gg;
                for (var j = 0; j < n; j++) {
                    g[j] = -xi[j];
                    h[j] = g[j] + gam * h[j];
                    xi[j] = h[j];
                }

                if (iterationCompletedListener != null) {
                    iterationCompletedListener.onIterationCompleted(this, its, ITMAX);
                }
            }

            if (!validResult) {
                // too many iterations
                locked = false;
                throw new OptimizationException();
            }
        } catch (final EvaluationException e) {
            throw new OptimizationException(e);
        } finally {
            locked = false;
        }

        // set result
        xmin = p;
        resultAvailable = true;
        fmin = fret;
    }

    /**
     * Returns boolean indicating whether this instance is ready to start the
     * estimation of a local minimum.
     * An instance is ready once a listener, a gradient listener and a start
     * point are provided.
     *
     * @return True if this instance is ready, false otherwise.
     */
    @Override
    public boolean isReady() {
        return isListenerAvailable() && isGradientListenerAvailable() && isStartPointAvailable();
    }

    /**
     * Returns tolerance or accuracy to be expected on estimated local minimum.
     *
     * @return Tolerance or accuracy to be expected on estimated local minimum.
     */
    public double getTolerance() {
        return tolerance;
    }

    /**
     * Sets tolerance or accuracy to be expected on estimated local minimum.
     *
     * @param tolerance Tolerance or accuracy to be expected on estimated local
     *                  minimum.
     * @throws LockedException          Raised if this instance is locked.
     * @throws IllegalArgumentException Raised if provided tolerance is
     *                                  negative.
     */
    public void setTolerance(final double tolerance) throws LockedException {
        if (isLocked()) {
            throw new LockedException();
        }
        internalSetTolerance(tolerance);
    }

    /**
     * Returns gradient listener in charge of obtaining gradient values for the
     * function to be evaluated.
     *
     * @return Gradient listener.
     * @throws NotAvailableException Raised if gradient listener has not yet
     *                               been provided.
     */
    public GradientFunctionEvaluatorListener getGradientListener() throws NotAvailableException {

}

    /**
     * Returns maximum allowed outliers proportion in the input data. This is
     * used to compute number of iterations to be done (nIters). It typically
     * can be as high as 0.95. Higher values, up to 1 are possible but not
     * recommended.
     * In this implementation, PROSAC won't stop before having reached the
     * corresponding inliers rate on the complete data set.
     *
     * @return maximum allowed outliers proportion in the input data.
     */
    public double getMaxOutliersProportion() {
        return maxOutliersProportion;
    }

    /**
     * Sets maximum allowed outliers proportion in the input data. This is used
     * to compute number of iterations to be done (nIters). It typically can be
     * as high as 0.95. Higher values, up to 1 are possible but not recommended.
     * In this implementation, PROSAC won't stop before having reached the
     * corresponding inliers rate on the complete data set.
     *
     * @param maxOutliersProportion maximum allowed outliers proportion in the
     *                              input data.
     * @throws IllegalArgumentException if provided value is less than 0.0 or
     *                                  greater than 1.0.
     * @throws LockedException          if this estimator is locked because an estimation
     *                                  is being computed.
     */
    public void setMaxOutliersProportion(final double maxOutliersProportion) throws LockedException {
        if (isLocked()) {
            throw new LockedException();
        }
        if (maxOutliersProportion < MIN_MAX_OUTLIERS_PROPORTION
                || maxOutliersProportion > MAX_MAX_OUTLIERS_PROPORTION) {
            throw new IllegalArgumentException();
        }

        this.maxOutliersProportion = maxOutliersProportion;
    }

    /**
     * Return eta0, which is the maximum probability that a solution with more
     * than inliersNStar inliers in U_nStar exists and was not found after k
     * samples (typically set to 5%).
     *
     * @return eta0 value.
     */
    public double getEta0() {
        return eta0;
    }

    /**
     * Sets eta0, which is the maximum probability that a solution with more
     * than inliersNStar inliers in U_nStar exists and was not found after k
     * samples (typically set to 5%).
     *
     * @param eta0 eta0 value to be set.
     * @throws IllegalArgumentException if provided value is less than 0.0 or
     *                                  greater than 1.0.
     * @throws LockedException          if this estimator is locked because an estimation
     *                                  is being computed.
     */
    public void setEta0(final double eta0) throws LockedException {
        if (isLocked()) {
            throw new LockedException();
        }
        if (eta0 < MIN_ETA0 || eta0 > MAX_ETA0) {
            throw new IllegalArgumentException();
        }

        this.eta0 = eta0;
    }

    /**
     * Returns beta, which is the probability that a match is declared inlier by
     * mistake, i.e. the ratio of the "inlier" surface by the total surface. The
     * inlier surface is a disc with radius 1.96s for homography/displacement
     * computation, or a band with width 1.96s*2 for epipolar geometry (s is
     * the detection noise), and the total surface is the surface of the image
     * YOU MUST ADJUST THIS VALUE, DEPENDING ON YOUR PROBLEM!
     *
     * @return beta value.
     */
    public double getBeta() {
        return beta;
    }

    /**
     * Sets beta, which is the probability that a match is declared inlier by
     * mistake, i.e. the ratio of the "inlier" surface by the total surface. The
     * inlier surface is a disc with radius 1.96s for homography/displacement
     * computation, or a band with width 1.96s*2 for epipolar geometry (s is
     * the detection noise), and the total surface is the surface of the image
     * YOU MUST ADJUST THIS VALUE, DEPENDING ON YOUR PROBLEM!
     *
     * @param beta beta value to be set.
     * @throws IllegalArgumentException if provided value is less than 0.0 or
     *                                  greater than 1.0.
     * @throws LockedException          if this estimator is locked because an estimation
     *                                  is being computed.
     */
    public void setBeta(final double beta) throws LockedException {
        if (isLocked()) {
            throw new LockedException();
        }
        if (beta < MIN_BETA || beta > MAX_BETA) {
            throw new IllegalArgumentException();
        }

        this.beta = beta;
    }

    /**
     * Returns number of iterations to be done to obtain required confidence.
     * This does not need to be equal to the actual number of iterations the
     * algorithm finally required to obtain a solution.
     *
     * @return number of iterations to be done to obtain required confidence.
     */
    public int getNIters() {
        return iters;

var sampleSizeBest = totalSamples;
                        var inliersSampleSizeBest = inliersCurrent;
                        // test value for the termination length
                        int sampleSizeTest;
                        // number of inliers for that test value
                        int inliersSampleSizeTest;
                        var epsilonSampleSizeBest = (double) inliersSampleSizeBest / (double) sampleSizeBest;

                        for (sampleSizeTest = totalSamples, inliersSampleSizeTest = inliersCurrent;
                             sampleSizeTest > subsetSize; sampleSizeTest--) {
                            // Loop invariants:
                            // - inliersSampleSizeTest is the number of inliers
                            //   for the sampleSizeTest first correspondences
                            // - sampleSizeBest is the value between
                            //   sampleSizeTest+1 and totalSamples that maximizes
                            //   the ratio inliersSampleSizeBest/sampleSizeBest

                            // - Non-randomness: In >= imin(n*)
                            // - Maximality: the number of samples that were drawn
                            //   so far must be enough so that the probability of
                            //   having missed a set of inliers is below eta=0.01.
                            //   This is the classical RANSAC termination criterion,
                            //   except that it takes into account only the
                            //   sampleSize first samples (not the total number
                            //   of samples)
                            //   kNStar = log(eta0) / log(1 - (inliersNStar/
                            //   sampleSizeStar)^subsetSize
                            //   We have to minimize kNStar, e.g. maximize
                            //   inliersNStar/sampleSizeStar, a straightforward
                            //   implementation would use the following test:
                            //   if(inliersSampleSizeTest > epsilonSampleSizeBest *
                            //   sampleSizeTest){ ... blah blah blah
                            //   However, since In is binomial, and in the case of
                            //   evenly distributed inliers, a better test would be
                            //   to reduce sampleSizeStar only if there's a
                            //   significant improvement in epsilon. Thus, we use a
                            //   Chi-squared test (P=0.10), together with the normal
                            //   approximation to the binomial (mu =
                            //   epsilonSampleSizeStart * sampleSizeTest, sigma =
                            //   sqrt(sampleSizeTest * epsilonSampleSizeStar * (1 -
                            //   epsilonSampleSizeStar))).
                            //   There is a significant difference between the two
                            //   tests (e.g. with the computeInliers function
                            //   provided)
                            //   We do the cheap test first, and the expensive test
                            //   only if the cheap one passes
                            if ((inliersSampleSizeTest * sampleSizeBest > inliersSampleSizeBest * sampleSizeTest)
                                    && (inliersSampleSizeTest > epsilonSampleSizeBest * sampleSizeTest
                                    + Math.sqrt(sampleSizeTest * epsilonSampleSizeBest
                                    * (1.0 - epsilonSampleSizeBest) * CHI_SQUARED))) {

                                if (inliersSampleSizeTest < imin(subsetSize, sampleSizeTest, beta)) {
                                    // equation not satisfied, no need to test for
                                    // smaller sampleSizeTest values anyway

                                    // jump out of the for sampleSizeTest loop
                                    break;
                                }
                                sampleSizeBest = sampleSizeTest;
                                inliersSampleSizeBest = inliersSampleSizeTest;
                                epsilonSampleSizeBest = (double) inliersSampleSizeBest / (double) sampleSizeBest;
                            }

                            // prepare for next loop iteration
                            inliersSampleSizeTest -= inliers.get(sortedIndices[sampleSizeTest - 1]) ? 1 : 0;
                        }

                        // is the best one we found even better than sampleSizeStar?
                        if (inliersSampleSizeBest * sampleSizeStar > inliersNStar * sampleSizeBest) {

                            // update all values
                            sampleSizeStar = sampleSizeBest;
                            inliersNStar = inliersSampleSizeBest;
                            kNStar = computeIterations((double) inliersNStar
                                            / (double) sampleSizeStar, subsetSize, 1.0 - eta0);
                        }
                    }
                }

for (int j = 1; j <= JMAX; j++) {
                q.next(s[j - 1]);
                // update sInterp
                for (var i = 0; i < elems; i++) {
                    sInterp[i][j - 1] = s[j - 1].getElementAtIndex(i);
                }
                if (j >= K) {
                    var finished = true;
                    for (var i = 0; i < elems; i++) {
                        final var ss = interpolators[i].rawinterp(j - K, 0.0);
                        if (Double.isNaN(ss)) {
                            throw new IntegrationException("NaN was found");
                        }
                        result.setElementAtIndex(i, ss);
                        if (Math.abs(interpolators[i].getDy()) > eps * Math.abs(ss)) {
                            finished = false;
                        }
                    }

                    if (finished) {
                        return;
                    }
                }
                h[j] = h[j - 1] / 9.0;

public BrentSingleRootEstimator(
            final SingleDimensionFunctionEvaluatorListener listener, final double minEvalPoint,
            final double maxEvalPoint, final double tolerance) throws InvalidBracketRangeException {
        super(listener, minEvalPoint, maxEvalPoint);
        internalSetTolerance(tolerance);
    }

    /**
     * Returns tolerance value.
     * Tolerance is the accuracy to be achieved when estimating a root.
     * If a root is found by this class, it is ensured to have an accuracy below
     * the tolerance value.
     *
     * @return Tolerance value.
     */
    public double getTolerance() {
        return tolerance;
    }

    /**
     * Internal method to set tolerance value.
     * Tolerance is the accuracy to be achieved when estimating a root.
     * If a root is found by this class, it is ensured to have an accuracy below
     * provided tolerance value.
     * This method does not check whether this instance is locked or not.
     *
     * @param tolerance Tolerance value.
     * @throws IllegalArgumentException Raised if provided tolerance value is
     *                                  negative.
     */
    private void internalSetTolerance(final double tolerance) {
        if (tolerance < MIN_TOLERANCE) {
            throw new IllegalArgumentException();
        }
        this.tolerance = tolerance;
    }

    /**
     * Sets tolerance value.
     * Tolerance is the accuracy to be achieved when estimating a root.
     * If a root is found by this class, it is ensured to have an accuracy below
     * provided tolerance value.
     *
     * @param tolerance Tolerance value.
     * @throws LockedException          Raised if this instance is locked.
     * @throws IllegalArgumentException Raised if provided tolerance value is
     *                                  negative.
     */
    public void setTolerance(final double tolerance) throws LockedException {
        if (isLocked()) {
            throw new LockedException();
        }
        internalSetTolerance(tolerance);
    }

    /**
     * Estimates a local root for a given single dimension function being
     * evaluated by provided listener.
     *
     * @throws LockedException         Exception raised if this instance is already
     *                                 locked.
     * @throws NotReadyException       Exception raised if either a listener has not
     *                                 yet been provided or a bracket has not been provided or computed.
     * @throws RootEstimationException Raised if the root estimation failed for
     *                                 some other reason (usually inability to evaluate the function,
     *                                 numerical instability or convergence problems, or no roots are found).
     */
    @Override
    @SuppressWarnings("Duplicates")
    public void estimate() throws LockedException, NotReadyException, RootEstimationException {
        if (isLocked()) {
            throw new LockedException();
        }
        if (!isReady()) {
            throw new NotReadyException();
        }

        locked = true;

        final var x1 = minEvalPoint;
        final var x2 = maxEvalPoint;
        final var tol = tolerance;

@Override
                    public int getTotalSamples() {
                        return evaluations.size();
                    }

                    @Override
                    public int getSubsetSize() {
                        return polynomialEstimator.getMinNumberOfEvaluations();
                    }

                    @Override
                    public void estimatePreliminarSolutions(
                            final int[] samplesIndices, final List<Polynomial> solutions) {
                        subsetEvaluations.clear();
                        for (final var samplesIndex : samplesIndices) {
                            subsetEvaluations.add(evaluations.get(samplesIndex));
                        }

                        try {
                            polynomialEstimator.setLMSESolutionAllowed(false);
                            polynomialEstimator.setEvaluations(subsetEvaluations);

                            final var polynomial = polynomialEstimator.estimate();
                            solutions.add(polynomial);
                        } catch (final Exception e) {
                            // if anything fails, no solution is added
                        }
                    }

                    @Override
                    public double computeResidual(final Polynomial currentEstimation, final int i) {
                        final var eval = evaluations.get(i);
                        return getDistance(eval, currentEstimation);
                    }

                    @Override
                    public boolean isReady() {
                        return LMedSPolynomialRobustEstimator.this.isReady();

final Matrix x, final double[] y, final double sig) throws FittingException {
        this(x, y, sig);
        setFunctionEvaluator(evaluator);
    }

    /**
     * Returns convergence parameter.
     *
     * @return convergence parameter.
     */
    public int getNdone() {
        return ndone;
    }

    /**
     * Sets convergence parameter.
     *
     * @param ndone convergence parameter.
     * @throws IllegalArgumentException if provided value is less than 1.
     */
    public void setNdone(final int ndone) {
        if (ndone < 1) {
            throw new IllegalArgumentException();
        }
        this.ndone = ndone;
    }

    /**
     * Returns maximum number of iterations.
     *
     * @return maximum number of iterations.
     */
    public int getItmax() {
        return itmax;
    }

    /**
     * Sets maximum number of iterations.
     *
     * @param itmax maximum number of iterations.
     * @throws IllegalArgumentException if provided value is zero or negative.
     */
    public void setItmax(final int itmax) {
        if (itmax <= 0) {
            throw new IllegalArgumentException();
        }
        this.itmax = itmax;
    }

    /**
     * Returns tolerance to reach convergence.
     *
     * @return tolerance to reach convergence.
     */
    public double getTol() {
        return tol;
    }

    /**
     * Sets tolerance to reach convergence.
     *
     * @param tol tolerance to reach convergence.
     * @throws IllegalArgumentException if provided value is zero or negative.
     */
    public void setTol(final double tol) {
        if (tol <= 0.0) {
            throw new IllegalArgumentException();
        }
        this.tol = tol;
    }

    /**
     * Returns function evaluator to evaluate function at a given point and
     * obtain function derivatives respect to each provided parameter.
     *
     * @return function evaluator.
     */
    public LevenbergMarquardtMultiDimensionFunctionEvaluator getFunctionEvaluator() {

&& x.getColumns() == evaluator.getNumberOfDimensions();
    }

    /**
     * Returns curvature matrix.
     *
     * @return curvature matrix.
     */
    public Matrix getAlpha() {
        return alpha;
    }

    /**
     * Returns degrees of freedom of computed chi square value.
     * Degrees of freedom is equal to the number of sampled data minus the
     * number of estimated parameters.
     *
     * @return degrees of freedom of computed chi square value.
     */
    public int getChisqDegreesOfFreedom() {
        return ndat - ma;
    }

    /**
     * Gets mean square error produced by estimated parameters respect to
     * provided sample data.
     *
     * @return mean square error.
     */
    public double getMse() {
        return mse;
    }

    /**
     * Gets the probability of finding a smaller chi square value.
     * The smaller the found chi square value is, the better the fit of the estimated
     * parameters to the actual parameter.
     * Thus, the smaller the chance of finding a smaller chi square value, then the
     * better the estimated fit is.
     *
     * @return probability of finding a smaller chi square value (better fit), expressed
     * as a value between 0.0 and 1.0.
     * @throws MaxIterationsExceededException if convergence of incomplete
     *                                        gamma function cannot be reached. This is rarely thrown and happens
     *                                        usually for numerically unstable input values.
     */
    public double getP() throws MaxIterationsExceededException {
        return ChiSqDist.cdf(getChisq(), getChisqDegreesOfFreedom());
    }

    /**
     * Gets a measure of quality of estimated fit as a value between 0.0 and 1.0.
     * The larger the quality value is, the better the fit that has been estimated.
     *
     * @return measure of quality of estimated fit.
     * @throws MaxIterationsExceededException if convergence of incomplete
     *                                        gamma function cannot be reached. This is rarely thrown and happens
     *                                        usually for numerically unstable input values.
     */
    public double getQ() throws MaxIterationsExceededException {
        return 1.0 - getP();
    }

    /**
     * Indicates whether covariance must be adjusted or not.
     * When covariance adjustment is enabled, then covariance is recomputed taking
     * into account input samples, input standard deviations of the samples and
     * jacobians of the model function overestimated parameters using the following
     * expression: Cov = (J'*W*J)^-1 where:
     * Cov is the covariance of estimated parameters
     * J is a matrix containing the Jacobians of the function overestimated parameters
     * for each input parameter x. Each row of J matrix contains an evaluation of
     * the model function Jacobian for i-th input parameter x. Each column of J matrix
     * contains the partial derivative of model function over j-th estimated parameter.
     * W is the inverse of input variances. It's a diagonal matrix containing the
     * reciprocal of the input variances (squared input standard deviations). That is:
     * W = diag(w) where k element of w is wk = 1 / sigmak^2, which corresponds to
     * the k-th standard deviation of input sample k.
     * By default, covariance is adjusted after fitting finishes.
     * <a href="http://people.duke.edu/~hpgavin/ce281/lm.pdf">http://people.duke.edu/~hpgavin/ce281/lm.pdf</a>
     * <a href="https://www8.cs.umu.se/kurser/5DA001/HT07/lectures/lsq-handouts.pdf">https://www8.cs.umu.se/kurser/5DA001/HT07/lectures/lsq-handouts.pdf</a>
     * Numerical Recipes 3rd Ed, page 812
     *
     * @return true if covariance must be adjusted, false otherwise.
     */
    public boolean isCovarianceAdjusted() {
        return adjustCovariance;
    }

    /**
     * Specifies whether covariance must be adjusted or not.
     * When covariance adjustment is enabled, then covariance is recomputed taking
     * into account input samples, input standard deviations of the samples and
     * jacobians of the model function overestimated parameters using the following
     * expression: Cov = (J'*W*J)^-1 where:
     * Cov is the covariance of estimated parameters
     * J is a matrix containing the Jacobians of the function overestimated parameters
     * for each input parameter x. Each row of J matrix contains an evaluation of
     * the model function Jacobian for i-th input parameter x. Each column of J matrix
     * contains the partial derivative of model function over j-th estimated parameter.
     * W is the inverse of input variances. It's a diagonal matrix containing the
     * reciprocal of the input variances (squared input standard deviations). That is:
     * W = diag(w) where k element of w is wk = 1 / sigmak^2, which corresponds to
     * the k-th standard deviation of input sample k.
     * By default, covariance is adjusted after fitting finishes.
     *
     * @param adjustCovariance true if covariance must be adjusted, false otherwise.
     */
    public void setCovarianceAdjusted(final boolean adjustCovariance) {
        this.adjustCovariance = adjustCovariance;
    }

    /**
     * Fits a function to provided data so that parameters associated to that
     * function can be estimated along with their covariance matrix and chi
     * square value.
     * If chi square value is close to 1, the fit is usually good.
     * If it is much larger, then error cannot be properly fitted.
     * If it is close to zero, then the model over-fits the error.
     * Methods {@link #getP()} and {@link #getQ()} can also be used to determine
     * the quality of the fit.
     *
     * @throws FittingException  if fitting fails.
     * @throws NotReadyException if enough input data has not yet been provided.
     */
    @Override
    public void fit() throws FittingException, NotReadyException {
        if (!isReady()) {
            throw new NotReadyException();
        }

        try {
            resultAvailable = false;

            int j;
            int k;
            int l;
            int iter;
            int done = 0;

}

    /**
     * Returns curvature matrix.
     *
     * @return curvature matrix.
     */
    public Matrix getAlpha() {
        return alpha;
    }

    /**
     * Returns degrees of freedom of computed chi square value.
     * Degrees of freedom is equal to the number of sampled data minus the
     * number of estimated parameters.
     *
     * @return degrees of freedom of computed chi square value.
     */
    public int getChisqDegreesOfFreedom() {
        return ndat - ma;
    }

    /**
     * Gets mean square error produced by estimated parameters respect to
     * provided sample data.
     *
     * @return mean square error.
     */
    public double getMse() {
        return mse;
    }

    /**
     * Gets the probability of finding a smaller chi square value.
     * The smaller the found chi square value is, the better the fit of the estimated
     * parameters to the actual parameter.
     * Thus, the smaller the chance of finding a smaller chi square value, then the
     * better the estimated fit is.
     *
     * @return probability of finding a smaller chi square value (better fit), expressed
     * as a value between 0.0 and 1.0.
     * @throws MaxIterationsExceededException if convergence of incomplete
     *                                        gamma function cannot be reached. This is rarely thrown and happens
     *                                        usually for numerically unstable input values.
     */
    public double getP() throws MaxIterationsExceededException {
        return ChiSqDist.cdf(getChisq(), getChisqDegreesOfFreedom());
    }

    /**
     * Gets a measure of quality of estimated fit as a value between 0.0 and 1.0.
     * The larger the quality value is, the better the fit that has been estimated.
     *
     * @return measure of quality of estimated fit.
     * @throws MaxIterationsExceededException if convergence of incomplete
     *                                        gamma function cannot be reached. This is rarely thrown and happens
     *                                        usually for numerically unstable input values.
     */
    public double getQ() throws MaxIterationsExceededException {
        return 1.0 - getP();
    }

    /**
     * Indicates whether covariance must be adjusted or not.
     * When covariance adjustment is enabled, then covariance is recomputed taking
     * into account input samples, input standard deviations of the samples and
     * jacobians of the model function overestimated parameters using the following
     * expression: Cov = (J'*W*J)^-1 where:
     * Cov is the covariance of estimated parameters
     * J is a matrix containing the Jacobians of the function overestimated parameters
     * for each input parameter x. Each row of J matrix contains an evaluation of
     * the model function Jacobian for i-th input parameter x. Each column of J matrix
     * contains the partial derivative of model function over j-th estimated parameter.
     * W is the inverse of input variances. It's a diagonal matrix containing the
     * reciprocal of the input variances (squared input standard deviations). That is:
     * W = diag(w) where k element of w is wk = 1 / sigmak^2, which corresponds to
     * the k-th standard deviation of input sample k.
     * By default, covariance is adjusted after fitting finishes.
     * <a href="http://people.duke.edu/~hpgavin/ce281/lm.pdf">http://people.duke.edu/~hpgavin/ce281/lm.pdf</a>
     * <a href="https://www8.cs.umu.se/kurser/5DA001/HT07/lectures/lsq-handouts.pdf">https://www8.cs.umu.se/kurser/5DA001/HT07/lectures/lsq-handouts.pdf</a>
     * Numerical Recipes 3rd Ed, page 812
     *
     * @return true if covariance must be adjusted, false otherwise.
     */
    public boolean isCovarianceAdjusted() {
        return adjustCovariance;
    }

    /**
     * Specifies whether covariance must be adjusted or not.
     * When covariance adjustment is enabled, then covariance is recomputed taking
     * into account input samples, input standard deviations of the samples and
     * jacobians of the model function overestimated parameters using the following
     * expression: Cov = (J'*W*J)^-1 where:
     * Cov is the covariance of estimated parameters
     * J is a matrix containing the Jacobians of the function overestimated parameters
     * for each input parameter x. Each row of J matrix contains an evaluation of
     * the model function Jacobian for i-th input parameter x. Each column of J matrix
     * contains the partial derivative of model function over j-th estimated parameter.
     * W is the inverse of input variances. It's a diagonal matrix containing the
     * reciprocal of the input variances (squared input standard deviations). That is:
     * W = diag(w) where k element of w is wk = 1 / sigmak^2, which corresponds to
     * the k-th standard deviation of input sample k.
     * By default, covariance is adjusted after fitting finishes.
     *
     * @param adjustCovariance true if covariance must be adjusted, false otherwise.
     */
    public void setCovarianceAdjusted(final boolean adjustCovariance) {
        this.adjustCovariance = adjustCovariance;
    }

    /**
     * Fits a function to provided data so that parameters associated to that
     * function can be estimated along with their covariance matrix and chi
     * square value.
     * If chi square value is close to 1, the fit is usually good.
     * If it is much larger, then error cannot be properly fitted.
     * If it is close to zero, then the model over-fits the error.
     * Methods {@link #getP()} and {@link #getQ()} can also be used to determine
     * the quality of the fit.
     *
     * @throws FittingException  if fitting fails.
     * @throws NotReadyException if enough input data has not yet been provided.
     */
    @Override
    public void fit() throws FittingException, NotReadyException {
        if (!isReady()) {
            throw new NotReadyException();
        }

        try {
            resultAvailable = false;

            int j;
            int k;
            int l;
            int iter;
            int done = 0;

}

                    @Override
                    public int getTotalSamples() {
                        return evaluations.size();
                    }

                    @Override
                    public int getSubsetSize() {
                        return polynomialEstimator.getMinNumberOfEvaluations();
                    }

                    @SuppressWarnings("DuplicatedCode")
                    @Override
                    public void estimatePreliminarSolutions(
                            final int[] samplesIndices, final List<Polynomial> solutions) {
                        subsetEvaluations.clear();
                        for (var samplesIndex : samplesIndices) {
                            subsetEvaluations.add(evaluations.get(samplesIndex));
                        }

                        try {
                            polynomialEstimator.setLMSESolutionAllowed(false);
                            polynomialEstimator.setEvaluations(subsetEvaluations);

                            final var polynomial = polynomialEstimator.estimate();
                            solutions.add(polynomial);
                        } catch (final Exception e) {
                            // if anything fails, no solution is added
                        }
                    }

                    @Override
                    public double computeResidual(final Polynomial currentEstimation, int i) {

final var inliers = inliersData.getInliers();
        var bestMedianResidual = inliersData.getBestMedianResidual();
        var medianResidualImproved = false;

        final var totalSamples = residuals.length;

        for (var i = 0; i < totalSamples; i++) {
            residuals[i] = Math.abs(listener.computeResidual(iterResult, i));
        }
        System.arraycopy(residuals, 0, residualsTemp, 0, residuals.length);
        final var medianResidual = sorter.median(residualsTemp);
        if (medianResidual < bestMedianResidual) {
            bestMedianResidual = medianResidual;
            medianResidualImproved = true;
        }

        final var standardDeviation = STD_CONSTANT * (1.0 + 5.0 / (totalSamples - subsetSize))
                * Math.sqrt(medianResidual);
        final var normEstimatedThreshold = inlierFactor * medianResidual;

        // determine which points are inliers
        var numInliers = 0;

nIters = maxIterations;
        bestResult = null;
        bestInliersData = null;
        computeAndKeepInliers = DEFAULT_COMPUTE_AND_KEEP_INLIERS;
        computeAndKeepResiduals = DEFAULT_COMPUTE_AND_KEEP_RESIDUALS;
    }

    /**
     * Returns amount of confidence expressed as a value between 0 and 1.0
     * (which is equivalent to 100%). The amount of confidence indicates the
     * probability that the estimated result is correct. Usually this value will
     * be close to 1.0, but not exactly 1.0.
     *
     * @return amount of confidence as a value between 0.0 and 1.0.
     */
    public double getConfidence() {
        return confidence;
    }

    /**
     * Sets amount of confidence expressed as a value between 0 and 1.0 (which
     * is equivalent to 100%). The amount of confidence indicates the
     * probability that the estimated result is correct. Usually this value will
     * be close to 1.0, but not exactly 1.0.
     *
     * @param confidence confidence to be set as a value between 0.0 and 1.0.
     * @throws IllegalArgumentException if provided value is not between 0.0 and
     *                                  1.0.
     * @throws LockedException          if this estimator is locked because an estimation
     *                                  is being computed.
     */
    public void setConfidence(final double confidence) throws LockedException {
        if (isLocked()) {
            throw new LockedException();
        }
        if (confidence < MIN_CONFIDENCE || confidence > MAX_CONFIDENCE) {
            throw new IllegalArgumentException();
        }
        this.confidence = confidence;
    }

    /**
     * Maximum allowed number of iterations. When the maximum number of
     * iterations is exceeded, result will not be available, however an
     * approximate result will be available for retrieval.
     *
     * @return maximum allowed number of iterations.
     */
    public int getMaxIterations() {
        return maxIterations;
    }

    /**
     * Sets maximum allowed number of iterations. When the maximum number of
     * iterations is exceeded, result will not be available, however an
     * approximate result will be available for retrieval.
     *
     * @param maxIterations maximum allowed number of iterations to be set.
     * @throws IllegalArgumentException if provided value is less than 1.
     * @throws LockedException          if this estimator is locked because an estimation
     *                                  is being computed.
     */
    public void setMaxIterations(final int maxIterations) throws LockedException {
        if (isLocked()) {
            throw new LockedException();
        }
        if (maxIterations < MIN_ITERATIONS) {
            throw new IllegalArgumentException();
        }
        this.maxIterations = maxIterations;
    }

    /**
     * Returns maximum allowed outliers proportion in the input data. This is
     * used to compute number of iterations to be done (nIters). It typically
     * can be as high as 0.95. Higher values, up to 1 are possible but not
     * recommended.
     * In this implementation, PROSAC won't stop before having reached the
     * corresponding inliers rate on the complete data set.
     *
     * @return maximum allowed outliers proportion in the input data.
     */
    public double getMaxOutliersProportion() {

super(listener, minEvalPoint, maxEvalPoint);
        internalSetTolerance(tolerance);
    }

    /**
     * Returns tolerance value.
     * Tolerance is the accuracy to be achieved when estimating a root.
     * If a root is found by this class, it is ensured to have an accuracy below
     * the tolerance value.
     *
     * @return Tolerance value.
     */
    public double getTolerance() {
        return tolerance;
    }

    /**
     * Internal method to set tolerance value.
     * Tolerance is the accuracy to be achieved when estimating a root.
     * If a root is found by this class, it is ensured to have an accuracy below
     * provided tolerance value.
     * This method does not check whether this instance is locked or not.
     *
     * @param tolerance Tolerance value.
     * @throws IllegalArgumentException Raised if provided tolerance value is
     *                                  negative.
     */
    private void internalSetTolerance(final double tolerance) {
        if (tolerance < MIN_TOLERANCE) {
            throw new IllegalArgumentException();
        }
        this.tolerance = tolerance;
    }

    /**
     * Sets tolerance value.
     * Tolerance is the accuracy to be achieved when estimating a root.
     * If a root is found by this class, it is ensured to have an accuracy below
     * provided tolerance value.
     *
     * @param tolerance Tolerance value.
     * @throws LockedException          Raised if this instance is locked.
     * @throws IllegalArgumentException Raised if provided tolerance value is
     *                                  negative.
     */
    public void setTolerance(final double tolerance) throws LockedException {
        if (isLocked()) {
            throw new LockedException();
        }
        internalSetTolerance(tolerance);
    }

    /**
     * Estimates a local root for a given single dimension function being
     * evaluated by provided listener.
     *
     * @throws LockedException         Exception raised if this instance is already
     *                                 locked.
     * @throws NotReadyException       Exception raised if either a listener has not
     *                                 yet been provided or a bracket has not been provided or computed.
     * @throws RootEstimationException Raised if the root estimation failed for
     *                                 some other reason (usually inability to evaluate the function,
     *                                 numerical instability or convergence problems, or no roots are found).
     */
    @Override
    @SuppressWarnings("Duplicates")
    public void estimate() throws LockedException, NotReadyException, RootEstimationException {
        if (isLocked()) {
            throw new LockedException();
        }
        if (!isReady()) {
            throw new NotReadyException();
        }

        locked = true;

super(listener, minEvalPoint, maxEvalPoint);
        internalSetTolerance(tolerance);
    }

    /**
     * Returns tolerance value.
     * Tolerance is the accuracy to be achieved when estimating a root.
     * If a root is found by this class, it is ensured to have an accuracy below
     * the tolerance value.
     *
     * @return Tolerance value.
     */
    public double getTolerance() {
        return tolerance;
    }

    /**
     * Internal method to set tolerance value.
     * Tolerance is the accuracy to be achieved when estimating a root.
     * If a root is found by this class, it is ensured to have an accuracy below
     * provided tolerance value.
     * This method does not check whether this instance is locked or not.
     *
     * @param tolerance Tolerance value.
     * @throws IllegalArgumentException Raised if provided tolerance value is
     *                                  negative.
     */
    private void internalSetTolerance(final double tolerance) {
        if (tolerance < MIN_TOLERANCE) {
            throw new IllegalArgumentException();
        }
        this.tolerance = tolerance;
    }

    /**
     * Sets tolerance value.
     * Tolerance is the accuracy to be achieved when estimating a root.
     * If a root is found by this class, it is ensured to have an accuracy below
     * provided tolerance value.
     *
     * @param tolerance Tolerance value.
     * @throws LockedException          Raised if this instance is locked.
     * @throws IllegalArgumentException Raised if provided tolerance value is
     *                                  negative.
     */
    public void setTolerance(final double tolerance) throws LockedException {
        if (isLocked()) {
            throw new LockedException();
        }
        internalSetTolerance(tolerance);
    }

    /**
     * Estimates a local root for a given single dimension function being
     * evaluated by provided listener.
     *
     * @throws LockedException         Exception raised if this instance is already
     *                                 locked.
     * @throws NotReadyException       Exception raised if either a listener has not
     *                                 yet been provided or a bracket has not been provided or computed.
     * @throws RootEstimationException Raised if the root estimation failed for
     *                                 some other reason (usually inability to evaluate the function,
     *                                 numerical instability or convergence problems, or no roots are found).
     */
    @Override
    @SuppressWarnings("Duplicates")
    public void estimate() throws LockedException, NotReadyException, RootEstimationException {
        if (isLocked()) {
            throw new LockedException();
        }
        if (!isReady()) {
            throw new NotReadyException();
        }

        locked = true;

public MSACPolynomialRobustEstimator(
            final int degree, final List<PolynomialEvaluation> evaluations,
            final PolynomialRobustEstimatorListener listener) {
        super(degree, evaluations, listener);
        threshold = DEFAULT_THRESHOLD;
    }

    /**
     * Returns threshold to determine whether polynomials are inliers or not
     * when testing possible estimation solutions.
     *
     * @return threshold to determine whether polynomials are inliers or not
     * when testing possible estimation solutions.
     */
    public double getThreshold() {
        return threshold;
    }

    /**
     * Sets threshold to determine whether polynomials are inliers or not when
     * testing possible estimation solutions.
     *
     * @param threshold threshold to determine whether polynomials are inliers
     *                  or not when testing possible estimation solutions.
     * @throws IllegalArgumentException if provided value is equal or less than
     *                                  zero.
     * @throws LockedException          if robust estimator is locked.
     */
    public void setThreshold(final double threshold) throws LockedException {
        if (isLocked()) {
            throw new LockedException();
        }
        if (threshold <= MIN_THRESHOLD) {
            throw new IllegalArgumentException();
        }
        this.threshold = threshold;
    }


    /**
     * Estimates polynomial.
     *
     * @return estimated polynomial.
     * @throws LockedException          if robust estimator is locked because an
     *                                  estimation is already in progress.
     * @throws NotReadyException        if provided input data is not enough to start
     *                                  the estimation.
     * @throws RobustEstimatorException if estimation fails for any other reason
     *                                  (i.e. numerical instability, no solution available, etc).
     */
    @Override
    public Polynomial estimate() throws LockedException, NotReadyException, RobustEstimatorException {
        if (isLocked()) {
            throw new LockedException();
        }
        if (!isReady()) {
            throw new NotReadyException();
        }

        final MSACRobustEstimator<Polynomial> innerEstimator =

public FalsePositionSingleRootEstimator(
            final SingleDimensionFunctionEvaluatorListener listener, final double minEvalPoint,
            final double maxEvalPoint, final double tolerance) throws InvalidBracketRangeException {
        super(listener, minEvalPoint, maxEvalPoint);
        internalSetTolerance(tolerance);
    }

    /**
     * Returns tolerance to find a root. Whenever the variation of the estimated
     * root is smaller than returned tolerance, then the algorithm is assumed to
     * be converged, and the estimated root is ensured to have an accuracy that
     * equals the returned tolerance.
     *
     * @return Tolerance to find a root.
     */
    public double getTolerance() {
        return tolerance;
    }

    /**
     * Sets tolerance to find a root. Whenever the variation of the estimated
     * root is smaller than provided tolerance, then the algorithm is assumed to
     * be converged, and the estimated root is ensured to have an accuracy that
     * equals provided tolerance.
     *
     * @param tolerance Tolerance to find a root.
     * @throws LockedException          Raised if this instance is locked while doing
     *                                  some computations.
     * @throws IllegalArgumentException Raised if provided tolerance is negative.
     */
    public void setTolerance(final double tolerance) throws LockedException {
        if (isLocked()) {
            throw new LockedException();
        }
        internalSetTolerance(tolerance);
    }

    /**
     * Estimates a single root of the provided single dimension function
     * contained within a given bracket of values.
     *
     * @throws LockedException         Exception raised if this instance is already
     *                                 locked.
     * @throws NotReadyException       Exception raised if not enough parameters have
     *                                 been provided in order to start the estimation.
     * @throws RootEstimationException Raised if the root estimation failed for
     *                                 some other reason (usually inability to evaluate the function,
     *                                 numerical instability or convergence problems, or no roots are found).
     */
    @Override
    public void estimate() throws LockedException, NotReadyException, RootEstimationException {
        if (isLocked()) {
            throw new LockedException();
        }
        if (!isReady()) {
            throw new NotReadyException();
        }

        locked = true;
        rootAvailable = false;

        try {
            final var v1 = new double[1];

final var probSubsetAllInliers = Math.pow((double) bestNumInliers / (double) totalSamples,
                                subsetSize);

                        if (Math.abs(probSubsetAllInliers) < Double.MIN_VALUE || Double.isNaN(probSubsetAllInliers)) {
                            newNIters = Integer.MAX_VALUE;
                        } else {
                            final var logProbSomeOutliers = Math.log(1.0 - probSubsetAllInliers);
                            if (Math.abs(logProbSomeOutliers) < Double.MIN_VALUE || Double.isNaN(logProbSomeOutliers)) {
                                newNIters = Integer.MAX_VALUE;
                            } else {
                                newNIters = (int) Math.ceil(Math.abs(Math.log(1.0 - confidence) / logProbSomeOutliers));
                            }
                        }
                        if (newNIters < iters) {

}

    /**
     * Returns quality scores corresponding to each provided point.
     * The larger the score value the better the quality of the sampled point
     *
     * @return quality scores corresponding to each point
     */
    @Override
    public double[] getQualityScores() {
        return qualityScores;
    }

    /**
     * Sets quality scores corresponding to each provided point.
     * The larger the score value the better the quality of the sampled point.
     *
     * @param qualityScores quality scores corresponding to each point.
     * @throws LockedException          if robust estimator is locked because an
     *                                  estimation is already in progress.
     * @throws IllegalArgumentException if provided quality scores length is
     *                                  smaller than required minimum size.
     */
    @Override
    public void setQualityScores(final double[] qualityScores) throws LockedException {
        if (isLocked()) {
            throw new LockedException();
        }
        internalSetQualityScores(qualityScores);
    }

    /**
     * Indicates if estimator is ready to start the polynomial estimation.
     * This is true when input data (i.e. polynomial evaluations and quality
     * scores) are provided and enough data is available.
     *
     * @return true if estimator is ready, false otherwise
     */
    @Override
    public boolean isReady() {
        return super.isReady() && qualityScores != null && qualityScores.length == evaluations.size();
    }

    /**
     * Estimates polynomial.
     *
     * @return estimated polynomial.
     * @throws LockedException          if robust estimator is locked because an
     *                                  estimation is already in progress.
     * @throws NotReadyException        if provided input data is not enough to start
     *                                  the estimation.
     * @throws RobustEstimatorException if estimation fails for any other reason
     *                                  (i.e. numerical instability, no solution available, etc).
     */
    @Override
    public Polynomial estimate() throws LockedException, NotReadyException, RobustEstimatorException {
        if (isLocked()) {
            throw new LockedException();
        }
        if (!isReady()) {
            throw new NotReadyException();
        }

        final PROMedSRobustEstimator<Polynomial> innerEstimator = new PROMedSRobustEstimator<>(

private void internalSetFunctionEvaluator(final LevenbergMarquardtMultiDimensionFunctionEvaluator evaluator)
            throws FittingException {

        try {
            this.evaluator = evaluator;

            if (evaluator != null) {
                a = evaluator.createInitialParametersArray();
                ma = a.length;
                covar = new Matrix(ma, ma);
                alpha = new Matrix(ma, ma);
                ia = new boolean[ma];
                Arrays.fill(ia, true);
            }
        } catch (final AlgebraException e) {
            throw new FittingException(e);
        }
    }

    /**
     * Used by fit to evaluate the linearized fitting matrix alpha, and vector
     * beta to calculate chi square.
     *
     * @param a     estimated parameters so far.
     * @param alpha curvature (i.e. fitting) matrix.
     * @param beta  array where derivative increments for each parameter are
     *              stored.
     * @throws AlgebraException    if there are numerical instabilities.
     * @throws EvaluationException if function evaluation fails.
     */
    private void mrqcof(final double[] a, final Matrix alpha, final double[] beta)
            throws AlgebraException, EvaluationException {

        int i;
        int j;
        int k;
        int l;
        int m;
        double ymod;

}

    /**
     * Returns tolerance value. Estimated result will be found with an accuracy
     * below or equal to provided tolerance value.
     *
     * @return Tolerance value.
     */
    public double getTolerance() {
        return tolerance;
    }

    /**
     * Sets tolerance value. Estimated result will be found with an accuracy
     * below or equal to provided tolerance value.
     *
     * @param tolerance Tolerance value.
     * @throws LockedException          Raised if this instance is locked.
     * @throws IllegalArgumentException Raised if tolerance is negative.
     */
    public void setTolerance(final double tolerance) throws LockedException {
        if (isLocked()) {
            throw new LockedException();
        }
        internalSetTolerance(tolerance);
    }

    /**
     * This function estimates a function minimum within provided or computed
     * bracket of values.
     * Given a function f that computes a function and also its derivative
     * function df, and given a bracketing triplet of abscissas "ax", "bx", "cx" (such
     * that bx is between ax and cx, and f(bx) is less than both f(ax) and
     * f(cx), this routine isolates the minimum to a fractional precision of
     * about tolerance using a modification of Brent's method that uses
     * derivatives. The abscissa of the minimum is returned as "xmin" and the
     * minimum function value is returned as "fmin".
     *
     * @throws LockedException       Raised if this instance is locked, because
     *                               estimation is being computed.
     * @throws NotReadyException     Raised if this instance is not ready because
     *                               either a listener or a bracket has not yet been provided or computed.
     * @throws OptimizationException Raised if the algorithm failed because of
     *                               lack of convergence or because function couldn't be evaluated.
     */
    @SuppressWarnings("DuplicatedCode")
    @Override
    public void minimize() throws LockedException, NotReadyException, OptimizationException {
        if (isLocked()) {
            throw new LockedException();
        }
        if (!isReady()) {
            throw new NotReadyException();
        }

        locked = true;

        final var v1 = new double[1];
        final var v2 = new double[2];
        final var v3 = new double[3];

        try {

}

    /**
     * Returns array of second degree polynomial parameters.
     * A second degree polynomial is defined by p(x) = a * x^2 + b * x + c, and
     * the array is returned as [c, b, a].
     * Note: This class only supports real polynomial parameters
     *
     * @return Array of first degree polynomial parameters
     * @throws NotAvailableException Raised if polynomial parameter have not yet
     *                               been provided
     */
    public double[] getRealPolynomialParameters() throws NotAvailableException {
        if (!arePolynomialParametersAvailable()) {
            throw new NotAvailableException();
        }
        return realPolyParams;
    }

    /**
     * Returns boolean indicating whether REAL polynomial parameters have been
     * provided and is available for retrieval.
     * Note: This class only supports real polynomial parameters
     *
     * @return True if available, false otherwise
     */
    @Override
    public boolean arePolynomialParametersAvailable() {
        return realPolyParams != null;
    }

    /**
     * This method will always raise a NotAvailableException because this class
     * only supports REAL polynomial parameters
     *
     * @return always throws NotAvailableException
     * @throws NotAvailableException always thrown
     */
    @Override
    public Complex[] getPolynomialParameters() throws NotAvailableException {
        throw new NotAvailableException();
    }

    /**
     * Estimates the roots of provided polynomial.
     *
     * @throws LockedException         Raised if this instance is locked estimating
     *                                 roots.
     * @throws NotReadyException       Raised if this instance is not ready because
     *                                 polynomial parameters have not been provided
     * @throws RootEstimationException Raised if roots cannot be estimated for
     *                                 some reason
     */
    @Override
    public void estimate() throws LockedException, NotReadyException, RootEstimationException {

        if (isLocked()) {
            throw new LockedException();
        }
        if (!isReady()) {
            throw new NotReadyException();
        }

        locked = true;

        roots = new Complex[VALID_POLY_PARAMS_LENGTH - 1];

        final var c = realPolyParams[0];

}

    /**
     * Returns amount of confidence expressed as a value between 0.0 and 1.0
     * (which is equivalent to 100%). The amount of confidence indicates the
     * probability that the estimated result is correct. Usually this value will
     * be close to 1.0, but not exactly 1.0.
     *
     * @return amount of confidence as a value between 0.0 and 1.0.
     */
    public double getConfidence() {
        return confidence;
    }

    /**
     * Sets amount of confidence expressed as a value between 0.0 and 1.0 (which
     * is equivalent to 100%). The amount of confidence indicates the
     * probability that the estimated result is correct. Usually this value will
     * be close to 1.0 but not exactly 1.0.
     *
     * @param confidence confidence to be set as a value between 0.0 and 1.0.
     * @throws IllegalArgumentException if provided value is not between 0.0 and
     *                                  1.0.
     * @throws LockedException          if this estimator is locked because an estimator
     *                                  is being computed.
     */
    public void setConfidence(final double confidence) throws LockedException {
        if (isLocked()) {
            throw new LockedException();
        }
        if (confidence < MIN_CONFIDENCE || confidence > MAX_CONFIDENCE) {
            throw new IllegalArgumentException();
        }
        this.confidence = confidence;
    }

    /**
     * Returns maximum allowed number of iterations. If maximum allowed number
     * of iterations is achieved without converging to a result when calling
     * estimate(), a RobustEstimatorException will be raised.
     *
     * @return maximum allowed number of iterations.
     */
    public int getMaxIterations() {
        return maxIterations;
    }

    /**
     * Sets maximum allowed number of iterations. When the maximum number of
     * iterations is exceeded, result will not be available, however an
     * approximate result will be available for retrieval.
     *
     * @param maxIterations maximum allowed number of iterations to be set.
     * @throws IllegalArgumentException if provided value is less than 1.
     * @throws LockedException          if this estimator is locked because an estimation
     *                                  is being computed.
     */
    public void setMaxIterations(final int maxIterations) throws LockedException {
        if (isLocked()) {
            throw new LockedException();
        }
        if (maxIterations < MIN_ITERATIONS) {
            throw new IllegalArgumentException();
        }
        this.maxIterations = maxIterations;
    }

    /**
     * Indicates whether geometric distance will be used to find outliers or
     * algebraic distance will be used instead.
     *
     * @return true if geometric distance is used, false otherwise.
     */
    public boolean isGeometricDistanceUsed() {

final var probSubsetAllInliers = Math.pow(probInlier, subsetSize);

                        if (Math.abs(probSubsetAllInliers) < Double.MIN_VALUE || Double.isNaN(probSubsetAllInliers)) {
                            newNIters = Integer.MAX_VALUE;
                        } else {
                            final var logProbSomeOutliers = Math.log(1.0 - probSubsetAllInliers);
                            if (Math.abs(logProbSomeOutliers) < Double.MIN_VALUE || Double.isNaN(logProbSomeOutliers)) {
                                newNIters = Integer.MAX_VALUE;
                            } else {
                                newNIters = (int) Math.ceil(Math.abs(Math.log(1.0 - confidence) / logProbSomeOutliers));
                            }
                        }

                        if (newNIters < iters) {
                            iters = newNIters;
                        }

}

    /**
     * Returns amount of confidence expressed as a value between 0.0 and 1.0
     * (which is equivalent to 100%). The amount of confidence indicates the
     * probability that the estimated result is correct. Usually this value will
     * be close to 1.0, but not exactly 1.0.
     *
     * @return amount of confidence as a value between 0.0 and 1.0.
     */
    public double getConfidence() {
        return confidence;
    }

    /**
     * Sets amount of confidence expressed as a value between 0.0 and 1.0 (which
     * is equivalent to 100%). The amount of confidence indicates the
     * probability that the estimated result is correct. Usually this value will
     * be close to 1.0 but not exactly 1.0.
     *
     * @param confidence confidence to be set as a value between 0.0 and 1.0.
     * @throws IllegalArgumentException if provided value is not between 0.0 and
     *                                  1.0.
     * @throws LockedException          if this estimator is locked because an estimator
     *                                  is being computed.
     */
    public void setConfidence(final double confidence) throws LockedException {
        if (isLocked()) {
            throw new LockedException();
        }
        if (confidence < MIN_CONFIDENCE || confidence > MAX_CONFIDENCE) {
            throw new IllegalArgumentException();
        }
        this.confidence = confidence;
    }

    /**
     * Returns maximum allowed number of iterations. If maximum allowed number
     * of iterations is achieved without converging to a result when calling
     * estimate(), a RobustEstimatorException will be raised.
     *
     * @return maximum allowed number of iterations.
     */
    public int getMaxIterations() {
        return maxIterations;
    }

    /**
     * Sets maximum allowed number of iterations. When the maximum number of
     * iterations is exceeded, result will not be available, however an
     * approximate result will be available for retrieval.
     *
     * @param maxIterations maximum allowed number of iterations to be set.
     * @throws IllegalArgumentException if provided value is less than 1.
     * @throws LockedException          if this estimator is locked because an estimation
     *                                  is being computed.
     */
    public void setMaxIterations(final int maxIterations) throws LockedException {
        if (isLocked()) {
            throw new LockedException();
        }
        if (maxIterations < MIN_ITERATIONS) {
            throw new IllegalArgumentException();
        }
        this.maxIterations = maxIterations;
    }

    /**
     * Indicates whether geometric distance will be used to find outliers or
     * algebraic distance will be used instead.
     *
     * @return true if geometric distance is used, false otherwise.
     */
    public boolean isGeometricDistanceUsed() {

public PROSACInliersData getBestInliersData() {
        return bestInliersData;
    }

    /**
     * Indicates whether inliers must be computed and kept.
     *
     * @return true if inliers must be computed and kept, false if inliers
     * only need to be computed but not kept.
     */
    public boolean isComputeAndKeepInliersEnabled() {
        return computeAndKeepInliers;
    }

    /**
     * Specifies whether inliers must be computed and kept.
     *
     * @param computeAndKeepInliers true if inliers must be computed and kept,
     *                              false if inliers only need to be computed but not kept.
     * @throws LockedException if estimator is locked.
     */
    public void setComputeAndKeepInliersEnabled(final boolean computeAndKeepInliers) throws LockedException {
        if (isLocked()) {
            throw new LockedException();
        }
        this.computeAndKeepInliers = computeAndKeepInliers;
    }

    /**
     * Indicates whether residuals must be computed and kept.
     *
     * @return true if residuals must be computed and kept, false if residuals
     * only need to be computed but not kept.
     */
    public boolean isComputeAndKeepResidualsEnabled() {
        return computeAndKeepResiduals;
    }

    /**
     * Specifies whether residuals must be computed and kept.
     *
     * @param computeAndKeepResiduals true if residuals must be computed and
     *                                kept, false if residuals only need to be computed but not kept.
     * @throws LockedException if estimator is locked.
     */
    public void setComputeAndKeepResidualsEnabled(final boolean computeAndKeepResiduals) throws LockedException {
        if (isLocked()) {
            throw new LockedException();
        }
        this.computeAndKeepResiduals = computeAndKeepResiduals;
    }

    /**
     * Indicates if estimator is ready to start the estimation process.
     *
     * @return true if ready, false otherwise.
     */
    @Override
    public boolean isReady() {
        if (!super.isReady()) {
            return false;
        }
        return (listener instanceof PROSACRobustEstimatorListener);

final var probSubsetAllInliers = Math.pow(probInlier, subsetSize);

                        if (Math.abs(probSubsetAllInliers) < Double.MIN_VALUE || Double.isNaN(probSubsetAllInliers)) {
                            newNIters = Integer.MAX_VALUE;
                        } else {
                            final var logProbSomeOutliers = Math.log(1.0 - probSubsetAllInliers);
                            if (Math.abs(logProbSomeOutliers) < Double.MIN_VALUE || Double.isNaN(logProbSomeOutliers)) {
                                newNIters = Integer.MAX_VALUE;
                            } else {
                                newNIters = (int) Math.ceil(Math.abs(Math.log(1.0 - confidence) / logProbSomeOutliers));
                            }
                        }

                        if (newNIters < iters) {

private void internalSetFunctionEvaluator(final LevenbergMarquardtMultiDimensionFunctionEvaluator evaluator)
            throws FittingException {

        try {
            this.evaluator = evaluator;

            if (evaluator != null) {
                a = evaluator.createInitialParametersArray();
                ma = a.length;
                covar = new Matrix(ma, ma);
                alpha = new Matrix(ma, ma);
                ia = new boolean[ma];
                Arrays.fill(ia, true);
            }
        } catch (final AlgebraException e) {
            throw new FittingException(e);
        }
    }

    /**
     * Used by fit to evaluate the linearized fitting matrix alpha, and vector
     * beta to calculate chi square.
     *
     * @param a     estimated parameters so far.
     * @param alpha curvature (i.e. fitting) matrix.
     * @param beta  array where derivative increments for each parameter are
     *              stored.
     * @throws AlgebraException    if there are numerical instabilities.
     * @throws EvaluationException if function evaluation fails.
     */
    private void mrqcof(final double[] a, final Matrix alpha, final double[] beta)
            throws AlgebraException, EvaluationException {