/* * This file is part of the X10 project (http://x10-lang.org). * * This file is licensed to You under the Eclipse Public License (EPL); * You may not use this file except in compliance with the License. * You may obtain a copy of the License at * http://www.opensource.org/licenses/eclipse-1.0.php * * (C) Copyright IBM Corporation 2006-2010. */ import x10.util.HashMap; import utils.CmdLnArgs; /** * See section 1.5 of the companion for a discussion of Simpson's rule and this code. */ public class Simpson { /** * If the subintervals are too small we run into precision problems because we only * have so many sig figs to work with. Here we choose a "resaonable" default minimum. * Real numerical analysts would tune this together with the relative precision * expected for the answer. */ public static DEFAULT_MIN_H = 1.0e-240; /** * We want a reasonable default relative precision so that we know when to terminate. * Stopping when an iteration yields a relative change of at most 1 part in 10 billion * seemed reasonable. */ public static DEFAULT_PRECISION = 1.0E-16; /** * We also supply a reasonable stopping point in terms of at most how many subintervals * to use, independent of the length b-a of the interval over which we are integrating. */ public static DEFAULT_LAST_N = 512*1024; /** * It makes no sense to add extra Places if each Place does not have a reasonable * amount of work to do. This number is the minimum number of subintervals to * require per Place */ public static DEFAULT_MIN_N = 250.0; // fields for the cleaned up command line parameters /** the relative or absolute improvement we require per iteration on h */ public val minDelta: Double; /** the left endpoint of the interval integrated over */ public val x0: Double; /** the right endpoint of the interval integrated over */ public val xn: Double; /** the minimum subinterval length to allow */ public val minH: Double; /** how many subintervals to use in the first iteration */ public val firstN: Double; /** upper bound on the number of subintervals to use per iteration */ public val lastN: Double; /** after each iteration, multiply the number of subintervals by this factor */ public val nfactor: Double; /** do not use more places than will guarantee each uses at least this many subintervals */ public val minN: Double; /** use a relative or absolute improvement check? relative by default */ public val isRelative: Boolean; /** if we are using a relative check, signal an error if 0 is the integral */ public val beStrict: Boolean; public val out: Out; // fields for the output public static class Out { /** the approximation to the integral via Simpson's Rule */ public var answer: Double = 0.0; /** * the absolute error at the last iteration; when combining the contributions * of multiple places, we form the sum. */ public var delta: Double = 0.0; /** * the subinterval length at the last iteration: when combining the * contributions of multiple places, we use the max over all places, * because we want to drive the max down to the min for every place. */ public var h: Double = 0.0; /** * the number of subintervals at the last iteration, or when combining * the contributions of multiple places, the sum over all places of the * number per place. */ public var n: Int = 0; /** returns a JSON literal */ public def toString() { return "{\r\n\t answer: " +answer +",\r\n\t delta: " +delta + ",\r\n\t h: "+h+",\r\n\t points_sampled: "+n+"\r\n\t}"; } /** * combines the cumulutive results from some set of subintervals with another * subinterval's result */ public def add(other: Out): Out { //Console.OUT.println("Add "+answer+" to "+other.answer); answer +=other.answer; delta += other.delta; h = Math.max(h, other.h); n += other.n; return this; } } /** * constructor that uses the command line to set the parameters * See the comments on main() for more. * @param args the command line arguments as a hash */ public def this(args: CmdLnArgs) { minDelta = args.getDouble("mindelta", DEFAULT_PRECISION); x0 = args.getDouble("x0", 1.0e260); xn = args.getDouble("xn", 1.0e-250); minH = args.getDouble("minh", DEFAULT_MIN_H); firstN = args.getInt("firstn", 2) as Double; lastN = args.getInt("lastn", DEFAULT_LAST_N) as Double; nfactor = args.getDouble("nfactor", 2.0); minN = args.getDouble("minn", DEFAULT_MIN_N); isRelative = args.getBoolean("relative", true); beStrict = args.getBoolean("strict", false); out = new Out(); } /** * create local parameters from a global one: * resets x0 and xn to suitable values for the k-th subinterval out of the total * number, which is placeCount. * @param oldStuff supplies the defaults missing from newStuff * @param newStuff a hash using the same keys as the command line to supply * new values for some of the fields in oldStuff. */ public def this(oldStuff: Simpson, k: Int, placeCount: Int, newTolerance: Double) { val newStuff = new HashMap[String, Double](); val subintervalLength: Double = (oldStuff.xn - oldStuff.x0)/(placeCount as Double); val kAsDouble = k as Double; val offset = k*subintervalLength; var newFirstN: Int = Math.ceil(oldStuff.firstN/placeCount) as Int; if(newFirstN % 2 == 1) newFirstN += 1; else if(newFirstN == 0) newFirstN = 2; var newLastN: Int = Math.ceil(oldStuff.lastN/placeCount) as Int; if(newLastN % 2 == 1) newLastN += 1; else if(newLastN == 0) newLastN = 2; minDelta = oldStuff.minDelta; x0 = oldStuff.x0+offset; xn = oldStuff.x0+offset+subintervalLength; minH = oldStuff.minH; firstN = newFirstN; lastN = newLastN; nfactor = oldStuff.nfactor; minN = oldStuff.minN; isRelative = (1.0 == newStuff.getOrElse("relative", oldStuff.isRelative?1.0:0.0)); beStrict = (1.0 == newStuff.getOrElse("strict", oldStuff.beStrict?1.0:0.0)); out = new Out(); } /** * computes the actual minimal subinterval length, the "target h", that is * allowed by these input parameters * @return the length computed */ public def getTargetH(): Double { val naiveValue = Math.max(minH, (xn - x0)/lastN); var lastH: Double = Double.MAX_VALUE; for (var h: Double = (xn - x0)/firstN; h >= naiveValue; h /= nfactor) lastH = h; return lastH; } /** * checks whether the input parameters make sense. */ public def isLegal() { return (firstN as Int) % 2 == 0 && (lastN as Int) % 2 == 0 && xn > x0 && xn-x0 > minH; } /** returns a JSON literal for the entire instance */ public def toString() { return "{\r\n\tinput: "+inputToString()+",\r\n\tresult: "+out.toString()+"\r\n}"; } /** returns the paramters as a JSON literal */ public def inputToString() { return "{\r\n\t minDelta: "+ minDelta+",\r\n\t " + "x0: "+x0+",\r\n\t " + "xn: "+xn+",\r\n\t " + "minH: "+minH+",\r\n\t " + "firstN: "+firstN+",\r\n\t " + "lastN: "+lastN+",\r\n\t " + "nfactor:"+nfactor+",\r\n\t " + "minN: "+minN+",\r\n\t " + "isRelative: "+isRelative+",\r\n\t " + "beStrict: "+beStrict+"\r\n\t}"; } /** * This is an iterative reapplication of Simpson's rule, with smaller and smaller * subintervals being used until the answer we get is "good enough", or until the * subinterval lengths reach an agreed minimum. The parameters specify both the * number of subintervals to start with, "firstN", and the maximum number, "lastN". * If firstN == lastN, exactly one iteration will be made. Otherwise, if on the * current iteration we don't get a good enough approximation, we multiply the * number of subintervals by the factor "nfactor", which is 2 by default and try * again... */ /* * PROGRAMMING NOTE * * The outer loop here is the loop over the subinterval lengths, and the inner loop * is the sum specified in Simpson's rule. There are two ways one might want to * measure one's progress: if oldA is the last answer we computed, and newA is the * one computed now, we can look for either an absolute change |newA - oldA|, * or we can bound the relative change|(newA - oldA)/newA|. We assume that newA * is always a better approximation to the real answer than oldA, so it makes sense * to compare the delta to newA, rather than to oldA or the averge of the two. * The input parameter isRelative is a Boolean that allows you to control * which of the two tests, relative or absolute, to use. * @param f the function we are integrating. * @param parms the parameters that describe the whole interval, the sizes of * subinterval we should use, the amount of real or absolute progress we should * make at each iteration, and how that check is to be made. * @return the answer, together with difference between the last answer and this one, * if we did more than 1 iteration, the size of last subinterval we used and the * number of subintervals. */ public def serial(f: (d:Double)=>Double) { val intervalLength = xn - x0; val boundaryTerms = f(x0) - f(xn); val firstH = intervalLength / firstN; var nextN:Double = firstN; val lastH = Math.max(minH, intervalLength / lastN); out.answer = out.delta = 0.0; out.h = getTargetH(); // Console.OUT.println(""+here.id+": "+this); for (var hToTry: Double = firstH; hToTry >= lastH; hToTry = intervalLength/nextN) { var odds: Double = 0.0; var evens: Double = 0.0; for(var xOdd: Double = x0 + hToTry; xOdd < xn; xOdd += 2.0*hToTry) { odds += f(xOdd); evens += f(xOdd + hToTry); } val bestGuess = ((boundaryTerms + 4.0*odds + 2.0*evens)*hToTry)/3.0; out.delta = out.answer - bestGuess; if (out.delta < 0.0) out.delta = -out.delta; out.answer = bestGuess; if (hToTry != firstH ) { // enough progress to try again? if(isRelative) { if(out.answer != 0.0) { if(out.delta < minDelta*out.answer) break; } else if (beStrict) { throw new IllegalArgumentException("Cannot use relative errors on a 0 solution."); } } else if (out.delta <= minDelta) break; } var nextNAsInt: Int = Math.floor(nextN*nfactor) as Int; if (nextNAsInt %2 == 1) nextNAsInt++; nextN = nextNAsInt as Double; out.h = hToTry; } out.n = Math.floor(intervalLength/out.h) as Int; out.delta = out.h==firstH?0:out.delta; return this; } /** * About the only thing worth breaking up is the interval [x0, xn] over which we are * integrating. If there are vast numbers of processers, assigning each processor an * equal share of the interval may not be smart: the processor may have too little to * work with. So step 1 is to decide how many processors to use. Step 2 is to let * each processor do its thing. There is no magic here: each processor just executes * serial(). The last step is to collect the results. * @param f the function we are integrating. * @return this instance of Simpson with the results in the appropriate instance * fields. */ /* * PROGRAMMING NOTE: When is the number of places too big? Answer: when the number * of subintervals each place is responsible for is too small. So we begin by * making sure that for our minimal subinterval length, there are enough * subintervals for each place. */ public def parallel(val f: (d:Double)=>Double) { out.answer = out.delta = 0.0; out.h = getTargetH(); // smallest "h" value to shoot for. var placeCount: Double = Place.MAX_PLACES as Double; for(; placeCount > 1.0; placeCount -= 1.0) { var intervalAtEachPlace: Double = (xn-x0) / placeCount; val targetSubintervalCount = Math.floor(intervalAtEachPlace/out.h); if (targetSubintervalCount >= minN) break; } val placeCountValue = placeCount as Int; Console.OUT.println("There are "+ placeCountValue + " places."); val dist = Dist.makeUnique() | (0..(placeCountValue-1)); val localMinDelta: Double = minDelta/placeCount; val outcomes = new Array[Out](placeCountValue, new Out()); finish for (var k: Int = 0; k < placeCount; ++k) { val kval = k; async outcomes(kval) = at(Place.place(kval)) new Simpson(this, kval, placeCountValue, localMinDelta).serial(f).out; } // you could use outcomes.reduce((a: Out, b: Out)=>a.add(b), this.out) here for(var k: Int = 0; k < placeCountValue; ++k) this.out.add(outcomes(k)); return this; } // given the coefficients num for the numerator and den for the denominator, // create the rational function that is their quotient. The return value // is the function (x: Double) => num(x)/den(x). private static def makeFunction(num: Array[Double](1), den: Array[Double](1)) { var f : (x: Double) => Double; if (den.size == 1) { if (num.size == 1) { val answer = num(0)/den(0); f = (x: Double): Double => answer; } else if (den(0) != 1.0) f = (x: Double): Double => poly(num,x)/den(0); else f = (x : Double): Double => poly(num,x); } else if (num.size == 1 ) { f = (x: Double): Double => num(0)/poly(den,x); } else f = (x: Double): Double => poly(num,x)/poly(den,x); return f; } // evaluates the polynomial with coefficients p at x. p(j) is the coefficient // of the term with exponent p.size - j. The reason is to make the input from // the command line translate sensibly: 1.0 2.0 3.0 => 1.0*x*x + 2.0*x + 3.0 // The point is that polynomials are usually written with the higher degree terms // to the left. private static def poly(p: Array[Double](1), x: Double) { var answer: Double = 0.0; for(v in p.values()) { answer = answer*x + v; } return answer; } // transforms an array of Double literal Strings to an array of Doubles private static def convertToDoubles(s: Array[String](1)) { return new Array[Double](s.size, (n: Int)=>Double.parse(s(n))); } /** * Carries out Simpson's rule for a polynomial or rational function. The command * line arguments specifying the function look like * -num cn ... c0 -den c'm ... c'0 * in which ck and c'k are the coefficients of x to the k-th in the numerator and * denominator, respectively. Both are optional: by default both numerator and * denominator are 1.0, so you need'nt supply a denominator if all you want is a * polynomial. * * The remaining command line arguments are: * x0 xxx the left endpoint of the interval * xn xxx the right endpoint of the interval * minh xxx minimum subinterval length to use * firstn nnn number of subintervals to use initially * lastn nnn maximum number of subintervals to use * nfactor xxx we multiply n by this factor after each iteration * mindelta xxx the absolute or relative error in the result that's good enough * minN xxx minimum number of subintervals for each Place * relative bbb a boolean: if true, mindelta is a relative error * strict bbb a boolean: if true, and if we are looking at relative errors, * report failure if the approximate answer is 0.0. * In each case, xxx is a Double and nnn is an even Int: note carefully: we want * even integers. The default values are * minh = 1.0E10*Double.MIN_NORMAL * mindelta = 1.0E-10 * firstN = 2 and lastN = 512*1024 * nfactor = 2.0 * relative = strict = true * There are no defaults for x0 and x1. Try * x10 Simpson -num 4.0 -den 1.0 0.0 1.0 -x0 0 -xn 1 * This is the integral from 0 to 1 of 4/(x*x + 1) == 4*arctan(1)-4*arctan(0) == * 4*pi/4 - 0 = pi. We do love pi. */ public static def main(args: Array[String](1)): void { try { val argsHash = new CmdLnArgs(args); // converts args to a hash val numArg = argsHash.get("num", "1.0"); // an array of strings val num = convertToDoubles(numArg); // the array of doubles we want val denArg = argsHash.get("den", "1.0"); // ditto val den = convertToDoubles(denArg); val serialData = new Simpson(argsHash); Console.OUT.println("Numerator "+num+"\r\nDenominator "+den); if (!serialData.isLegal()) { Console.OUT.println("The parameters are not consistent."); } else { val f = makeFunction(num, den); val sOut = serialData.serial(f); val pOut = new Simpson(argsHash).parallel(f); Console.OUT.println("Results:\r\n parallel\r\n "+pOut+",\r\n serial\r\n "+sOut); } } catch(e: Exception) { Console.ERR.println(e); } } }