    001    /*
    002     * Licensed to the Apache Software Foundation (ASF) under one or more
    003     * contributor license agreements.  See the NOTICE file distributed with
    004     * this work for additional information regarding copyright ownership.
    005     * The ASF licenses this file to You under the Apache License, Version 2.0
    006     * (the "License"); you may not use this file except in compliance with
    007     * the License.  You may obtain a copy of the License at
    008     *
    009     *      http://www.apache.org/licenses/LICENSE-2.0
    010     *
    011     * Unless required by applicable law or agreed to in writing, software
    012     * distributed under the License is distributed on an "AS IS" BASIS,
    013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    014     * See the License for the specific language governing permissions and
    015     * limitations under the License.
    016     */
    017    package org.apache.commons.math3.analysis.integration.gauss;
    018
    019    import java.util.Map;
    020    import java.util.TreeMap;
    021    import org.apache.commons.math3.util.Pair;
    022    import org.apache.commons.math3.exception.DimensionMismatchException;
    023    import org.apache.commons.math3.exception.NotStrictlyPositiveException;
    024
    025    /**
    026     * Base class for rules that determines the integration nodes and their
    027     * weights.
    028     * Subclasses must implement the {@link #computeRule(int) computeRule} method.
    029     *
    030     * @param <T> Type of the number used to represent the points and weights of
    031     * the quadrature rules.
    032     *
    033     * @since 3.1
    034     * @version $Id: BaseRuleFactory.java 1382197 2012-09-07 22:35:01Z erans $
    035     */
    036    public abstract class BaseRuleFactory<T extends Number> {
    037        /** List of points and weights, indexed by the order of the rule. */
    038        private final Map<Integer, Pair<T[], T[]>> pointsAndWeights
    039            = new TreeMap<Integer, Pair<T[], T[]>>();
    040        /** Cache for double-precision rules. */
    041        private final Map<Integer, Pair<double[], double[]>> pointsAndWeightsDouble
    042            = new TreeMap<Integer, Pair<double[], double[]>>();
    043
    044        /**
    045         * Gets a copy of the quadrature rule with the given number of integration
    046         * points.
    047         *
    048         * @param numberOfPoints Number of integration points.
    049         * @return a copy of the integration rule.
    050         * @throws NotStrictlyPositiveException if {@code numberOfPoints < 1}.
    051         */
    052        public Pair<double[], double[]> getRule(int numberOfPoints)
    053            throws NotStrictlyPositiveException {
    054            // Try to obtain the rule from the cache.
    055            Pair<double[], double[]> cached = pointsAndWeightsDouble.get(numberOfPoints);
    056
    057            if (cached == null) {
    058                // Rule not computed yet.
    059
    060                // Compute the rule.
    061                final Pair<T[], T[]> rule = getRuleInternal(numberOfPoints);
    062                cached = convertToDouble(rule);
    063
    064                // Cache it.
    065                pointsAndWeightsDouble.put(numberOfPoints, cached);
    066            }
    067
    068            // Return a copy.
    069            return new Pair<double[], double[]>(cached.getFirst().clone(),
    070                                                cached.getSecond().clone());
    071        }
    072
    073        /**
    074         * Gets a rule.
    075         * Synchronization ensures that rules will be computed and added to the
    076         * cache at most once.
    077         * The returned rule is a reference into the cache.
    078         *
    079         * @param numberOfPoints Order of the rule to be retrieved.
    080         * @return the points and weights corresponding to the given order.
    081         * @throws NotStrictlyPositiveException if {@code numberOfPoints < 1}.
    082         */
    083        protected synchronized Pair<T[], T[]> getRuleInternal(int numberOfPoints)
    084            throws NotStrictlyPositiveException {
    085            final Pair<T[], T[]> rule = pointsAndWeights.get(numberOfPoints);
    086            if (rule == null) {
    087                addRule(computeRule(numberOfPoints));
    088                // The rule should be available now.
    089                return getRuleInternal(numberOfPoints);
    090            }
    091            return rule;
    092        }
    093
    094        /**
    095         * Stores a rule.
    096         *
    097         * @param rule Rule to be stored.
    098         * @throws DimensionMismatchException if the elements of the pair do not
    099         * have the same length.
    100         */
    101        protected void addRule(Pair<T[], T[]> rule) {
    102            if (rule.getFirst().length != rule.getSecond().length) {
    103                throw new DimensionMismatchException(rule.getFirst().length,
    104                                                     rule.getSecond().length);
    105            }
    106
    107            pointsAndWeights.put(rule.getFirst().length, rule);
    108        }
    109
    110        /**
    111         * Computes the rule for the given order.
    112         *
    113         * @param numberOfPoints Order of the rule to be computed.
    114         * @return the computed rule.
    115         */
    116        protected abstract Pair<T[], T[]> computeRule(int numberOfPoints);
    117
    118        /**
    119         * Converts the from the actual {@code Number} type to {@code double}
    120         *
    121         * @param <T> Type of the number used to represent the points and
    122         * weights of the quadrature rules.
    123         * @param rule Points and weights.
    124         * @return points and weights as {@code double}s.
    125         */
    126        private static <T extends Number> Pair<double[], double[]> convertToDouble(Pair<T[], T[]> rule) {
    127            final T[] pT = rule.getFirst();
    128            final T[] wT = rule.getSecond();
    129
    130            final int len = pT.length;
    131            final double[] pD = new double[len];
    132            final double[] wD = new double[len];
    133
    134            for (int i = 0; i < len; i++) {
    135                pD[i] = pT[i].doubleValue();
    136                wD[i] = wT[i].doubleValue();
    137            }
    138
    139            return new Pair<double[], double[]>(pD, wD);
    140        }
    141    }