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    }