/*
* 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);
}
}
}