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.math4.legacy.analysis.integration.gauss; 018 019import java.util.Map; 020import java.util.TreeMap; 021 022import org.apache.commons.math4.legacy.exception.DimensionMismatchException; 023import org.apache.commons.math4.legacy.exception.NotStrictlyPositiveException; 024import org.apache.commons.math4.legacy.exception.util.LocalizedFormats; 025import org.apache.commons.math4.legacy.core.Pair; 026 027/** 028 * Base class for rules that determines the integration nodes and their 029 * weights. 030 * Subclasses must implement the {@link #computeRule(int) computeRule} method. 031 * 032 * @param <T> Type of the number used to represent the points and weights of 033 * the quadrature rules. 034 * 035 * @since 3.1 036 */ 037public abstract class BaseRuleFactory<T extends Number> { 038 /** List of points and weights, indexed by the order of the rule. */ 039 private final Map<Integer, Pair<T[], T[]>> pointsAndWeights 040 = new TreeMap<>(); 041 /** Cache for double-precision rules. */ 042 private final Map<Integer, Pair<double[], double[]>> pointsAndWeightsDouble 043 = new TreeMap<>(); 044 045 /** 046 * Gets a copy of the quadrature rule with the given number of integration 047 * points. 048 * 049 * @param numberOfPoints Number of integration points. 050 * @return a copy of the integration rule. 051 * @throws NotStrictlyPositiveException if {@code numberOfPoints < 1}. 052 * @throws DimensionMismatchException if the elements of the rule pair do not 053 * have the same length. 054 */ 055 public Pair<double[], double[]> getRule(int numberOfPoints) { 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<>(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 final Pair<T[], T[]> rule = pointsAndWeights.get(numberOfPoints); 094 if (rule == null) { 095 addRule(computeRule(numberOfPoints)); 096 // The rule should be available now. 097 return getRuleInternal(numberOfPoints); 098 } 099 return rule; 100 } 101 102 /** 103 * Stores a rule. 104 * 105 * @param rule Rule to be stored. 106 * @throws DimensionMismatchException if the elements of the pair do not 107 * have the same length. 108 */ 109 protected void addRule(Pair<T[], T[]> rule) { 110 if (rule.getFirst().length != rule.getSecond().length) { 111 throw new DimensionMismatchException(rule.getFirst().length, 112 rule.getSecond().length); 113 } 114 115 pointsAndWeights.put(rule.getFirst().length, rule); 116 } 117 118 /** 119 * Computes the rule for the given order. 120 * 121 * @param numberOfPoints Order of the rule to be computed. 122 * @return the computed rule. 123 * @throws DimensionMismatchException if the elements of the pair do not 124 * have the same length. 125 */ 126 protected abstract Pair<T[], T[]> computeRule(int numberOfPoints); 127 128 /** 129 * Converts the from the actual {@code Number} type to {@code double}. 130 * 131 * @param <T> Type of the number used to represent the points and 132 * weights of the quadrature rules. 133 * @param rule Points and weights. 134 * @return points and weights as {@code double}s. 135 */ 136 private static <T extends Number> Pair<double[], double[]> convertToDouble(Pair<T[], T[]> rule) { 137 final T[] pT = rule.getFirst(); 138 final T[] wT = rule.getSecond(); 139 140 final int len = pT.length; 141 final double[] pD = new double[len]; 142 final double[] wD = new double[len]; 143 144 for (int i = 0; i < len; i++) { 145 pD[i] = pT[i].doubleValue(); 146 wD[i] = wT[i].doubleValue(); 147 } 148 149 return new Pair<>(pD, wD); 150 } 151}