FieldAbstractRuleFactory

/*
 * Licensed to the Hipparchus project under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The Hipparchus project 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
 *
 *      https://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.hipparchus.analysis.integration.gauss;

import java.util.Arrays;
import java.util.SortedMap;
import java.util.TreeMap;

import org.hipparchus.CalculusFieldElement;
import org.hipparchus.Field;
import org.hipparchus.analysis.CalculusFieldUnivariateFunction;
import org.hipparchus.exception.LocalizedCoreFormats;
import org.hipparchus.exception.MathIllegalArgumentException;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.Incrementor;
import org.hipparchus.util.MathArrays;
import org.hipparchus.util.Pair;

/**
 * Base class for rules that determines the integration nodes and their
 * weights.
 * Subclasses must implement the {@link #computeRule(int) computeRule} method.
 *
 * @param <T> Type of the number used to represent the points and weights of
 * the quadrature rules.
 * @since 2.0
 */
public abstract class FieldAbstractRuleFactory<T extends CalculusFieldElement<T>> implements FieldRuleFactory<T> {

    /** Field to which rule coefficients belong. */
    private final Field<T> field;

    /** List of points and weights, indexed by the order of the rule. */
    private final SortedMap<Integer, Pair<T[], T[]>> pointsAndWeights;

    /** Simple constructor
     * @param field field to which rule coefficients belong
     */
    protected FieldAbstractRuleFactory(final Field<T> field) {
        this.field            = field;
        this.pointsAndWeights = new TreeMap<>();
    }

    /** Get the field to which rule coefficients belong.
     * @return field to which rule coefficients belong
     */
    public Field<T> getField() {
        return field;
    }

    /** {@inheritDoc} */
    @Override
    public Pair<T[], T[]> getRule(int numberOfPoints)
        throws MathIllegalArgumentException {

        if (numberOfPoints <= 0) {
            throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_OF_POINTS,
                                                   numberOfPoints);
        }
        if (numberOfPoints > 1000) {
            throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_LARGE,
                                                   numberOfPoints, 1000);
        }

        Pair<T[], T[]> rule;
        synchronized (pointsAndWeights) {
            // Try to obtain the rule from the cache.
            rule = pointsAndWeights.get(numberOfPoints);

            if (rule == null) {
                // Rule not computed yet.

                // Compute the rule.
                rule = computeRule(numberOfPoints);

                // Cache it.
                pointsAndWeights.put(numberOfPoints, rule);
            }
        }

        // Return a copy.
        return new Pair<>(rule.getFirst().clone(), rule.getSecond().clone());

    }

    /**
     * Computes the rule for the given order.
     *
     * @param numberOfPoints Order of the rule to be computed.
     * @return the computed rule.
     * @throws MathIllegalArgumentException if the elements of the pair do not
     * have the same length.
     */
    protected abstract Pair<T[], T[]> computeRule(int numberOfPoints)
        throws MathIllegalArgumentException;

    /** Computes roots of the associated orthogonal polynomials.
     * <p>
     * The roots are found using the <a href="https://en.wikipedia.org/wiki/Aberth_method">Aberth method</a>.
     * The guess points for initializing search for degree n are fixed for degrees 1 and 2 and are
     * selected from n-1 roots of rule n-1 (the two extreme roots are used, plus the n-1 intermediate
     * points between all roots).
     * </p>
     * @param n number of roots to search for
     * @param ratioEvaluator function evaluating the ratio Pₙ(x)/Pₙ'(x)
     * @return sorted array of roots
     */
    protected T[] findRoots(final int n, final CalculusFieldUnivariateFunction<T> ratioEvaluator) {

        final T[] roots  = MathArrays.buildArray(field, n);

        // set up initial guess
        if (n == 1) {
            // arbitrary guess
            roots[0] = field.getZero();
        } else if (n == 2) {
            // arbitrary guess
            roots[0] = field.getOne().negate();
            roots[1] = field.getOne();
        } else {

            // get roots from previous rule.
            // If it has not been computed yet it will trigger a recursive call
            final T[] previousPoints = getRule(n - 1).getFirst();

            // first guess at previous first root
            roots[0] = previousPoints[0];

            // intermediate guesses between previous roots
            for (int i = 1; i < n - 1; ++i) {
                roots[i] = previousPoints[i - 1].add(previousPoints[i]).multiply(0.5);
            }

            // last guess at previous last root
            roots[n - 1] = previousPoints[n - 2];

        }

        // use Aberth method to find all roots simultaneously
        final T[]         ratio       = MathArrays.buildArray(field, n);
        final Incrementor incrementor = new Incrementor(1000);
        double            tol;
        double            maxOffset;
        do {

            // safety check that triggers an exception if too much iterations are made
            incrementor.increment();

            // find the ratio P(xᵢ)/P'(xᵢ) for all current roots approximations
            for (int i = 0; i < n; ++i) {
                ratio[i] = ratioEvaluator.value(roots[i]);
            }

            // move roots approximations all at once, using Aberth method
            maxOffset = 0;
            for (int i = 0; i < n; ++i) {
                T sum = field.getZero();
                for (int j = 0; j < n; ++j) {
                    if (j != i) {
                        sum = sum.add(roots[i].subtract(roots[j]).reciprocal());
                    }
                }
                final T offset = ratio[i].divide(sum.multiply(ratio[i]).negate().add(1));
                maxOffset = FastMath.max(maxOffset, FastMath.abs(offset).getReal());
                roots[i] = roots[i].subtract(offset);
            }

            // we set tolerance to 1 ulp of the largest root
            tol = 0;
            for (final T r : roots) {
                tol = FastMath.max(tol, FastMath.ulp(r.getReal()));
            }

        } while (maxOffset > tol);

        // sort the roots
        Arrays.sort(roots, (r1, r2) -> Double.compare(r1.getReal(), r2.getReal()));

        return roots;

    }

    /** Enforce symmetry of roots.
     * @param roots roots to process in place
     */
    protected void enforceSymmetry(final T[] roots) {

        final int n = roots.length;

        // enforce symmetry
        for (int i = 0; i < n / 2; ++i) {
            final int idx = n - i - 1;
            final T c = roots[i].subtract(roots[idx]).multiply(0.5);
            roots[i]   = c;
            roots[idx] = c.negate();
        }

        // If n is odd, 0 is a root.
        // Note: as written, the test for oddness will work for negative
        // integers too (although it is not necessary here), preventing
        // a FindBugs warning.
        if (n % 2 != 0) {
            roots[n / 2] = field.getZero();
        }

    }

}