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 */
017package org.apache.commons.math3.analysis.integration.gauss;
018
019import java.util.Map;
020import java.util.TreeMap;
021import org.apache.commons.math3.util.Pair;
022import org.apache.commons.math3.exception.DimensionMismatchException;
023import org.apache.commons.math3.exception.NotStrictlyPositiveException;
024import org.apache.commons.math3.exception.util.LocalizedFormats;
025
026/**
027 * Base class for rules that determines the integration nodes and their
028 * weights.
029 * Subclasses must implement the {@link #computeRule(int) computeRule} method.
030 *
031 * @param <T> Type of the number used to represent the points and weights of
032 * the quadrature rules.
033 *
034 * @since 3.1
035 */
036public 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     * @throws DimensionMismatchException if the elements of the rule pair do not
052     * have the same length.
053     */
054    public Pair<double[], double[]> getRule(int numberOfPoints)
055        throws NotStrictlyPositiveException, DimensionMismatchException {
056
057        if (numberOfPoints <= 0) {
058            throw new NotStrictlyPositiveException(LocalizedFormats.NUMBER_OF_POINTS,
059                                                   numberOfPoints);
060        }
061
062        // Try to obtain the rule from the cache.
063        Pair<double[], double[]> cached = pointsAndWeightsDouble.get(numberOfPoints);
064
065        if (cached == null) {
066            // Rule not computed yet.
067
068            // Compute the rule.
069            final Pair<T[], T[]> rule = getRuleInternal(numberOfPoints);
070            cached = convertToDouble(rule);
071
072            // Cache it.
073            pointsAndWeightsDouble.put(numberOfPoints, cached);
074        }
075
076        // Return a copy.
077        return new Pair<double[], double[]>(cached.getFirst().clone(),
078                                            cached.getSecond().clone());
079    }
080
081    /**
082     * Gets a rule.
083     * Synchronization ensures that rules will be computed and added to the
084     * cache at most once.
085     * The returned rule is a reference into the cache.
086     *
087     * @param numberOfPoints Order of the rule to be retrieved.
088     * @return the points and weights corresponding to the given order.
089     * @throws DimensionMismatchException if the elements of the rule pair do not
090     * have the same length.
091     */
092    protected synchronized Pair<T[], T[]> getRuleInternal(int numberOfPoints)
093        throws DimensionMismatchException {
094        final Pair<T[], T[]> rule = pointsAndWeights.get(numberOfPoints);
095        if (rule == null) {
096            addRule(computeRule(numberOfPoints));
097            // The rule should be available now.
098            return getRuleInternal(numberOfPoints);
099        }
100        return rule;
101    }
102
103    /**
104     * Stores a rule.
105     *
106     * @param rule Rule to be stored.
107     * @throws DimensionMismatchException if the elements of the pair do not
108     * have the same length.
109     */
110    protected void addRule(Pair<T[], T[]> rule) throws DimensionMismatchException {
111        if (rule.getFirst().length != rule.getSecond().length) {
112            throw new DimensionMismatchException(rule.getFirst().length,
113                                                 rule.getSecond().length);
114        }
115
116        pointsAndWeights.put(rule.getFirst().length, rule);
117    }
118
119    /**
120     * Computes the rule for the given order.
121     *
122     * @param numberOfPoints Order of the rule to be computed.
123     * @return the computed rule.
124     * @throws DimensionMismatchException if the elements of the pair do not
125     * have the same length.
126     */
127    protected abstract Pair<T[], T[]> computeRule(int numberOfPoints)
128        throws DimensionMismatchException;
129
130    /**
131     * Converts the from the actual {@code Number} type to {@code double}
132     *
133     * @param <T> Type of the number used to represent the points and
134     * weights of the quadrature rules.
135     * @param rule Points and weights.
136     * @return points and weights as {@code double}s.
137     */
138    private static <T extends Number> Pair<double[], double[]> convertToDouble(Pair<T[], T[]> rule) {
139        final T[] pT = rule.getFirst();
140        final T[] wT = rule.getSecond();
141
142        final int len = pT.length;
143        final double[] pD = new double[len];
144        final double[] wD = new double[len];
145
146        for (int i = 0; i < len; i++) {
147            pD[i] = pT[i].doubleValue();
148            wD[i] = wT[i].doubleValue();
149        }
150
151        return new Pair<double[], double[]>(pD, wD);
152    }
153}