lh3
7/19/2013 - 6:00 PM

Given a discrete one-dimention function f(x), fit it with Bernstein polynomial and the find the max. k8 is required.

Given a discrete one-dimention function f(x), fit it with Bernstein polynomial and the find the max. k8 is required.

var getopt = function(args, ostr) {
  var oli; // option letter list index
	if (typeof(getopt.place) == 'undefined')
		getopt.ind = 0, getopt.arg = null, getopt.place = -1;
	if (getopt.place == -1) { // update scanning pointer
		if (getopt.ind >= args.length || args[getopt.ind].charAt(getopt.place = 0) != '-') {
			getopt.place = -1;
			return null;
		}
		if (getopt.place + 1 < args[getopt.ind].length && args[getopt.ind].charAt(++getopt.place) == '-') { // found "--"
			++getopt.ind;
			getopt.place = -1;
			return null;
		}
	}
	var optopt = args[getopt.ind].charAt(getopt.place++); // character checked for validity
	if (optopt == ':' || (oli = ostr.indexOf(optopt)) < 0) {
		if (optopt == '-') return null; //  if the user didn't specify '-' as an option, assume it means null.
		if (getopt.place < 0) ++getopt.ind;
		return '?';
	}
	if (oli+1 >= ostr.length || ostr.charAt(++oli) != ':') { // don't need argument
		getopt.arg = null;
		if (getopt.place < 0 || getopt.place >= args[getopt.ind].length) ++getopt.ind, getopt.place = -1;
	} else { // need an argument
		if (getopt.place >= 0 && getopt.place < args[getopt.ind].length)
			getopt.arg = args[getopt.ind].substr(getopt.place);
		else if (args.length <= ++getopt.ind) { // no arg
			getopt.place = -1;
			if (ostr.length > 0 && ostr.charAt(0) == ':') return ':';
			return '?';
		} else getopt.arg = args[getopt.ind]; // white space
		getopt.place = -1;
		++getopt.ind;
	}
	return optopt;
}

Math.bernstein_poly = function(beta) // Bernstein polynomial with De Casteljau's algorithm
{
	var n = beta.length - 1;
	return function(t) {
		var prev = [], next = [];
		for (var i = 0; i <= n; ++i) prev[i] = beta[i];
		for (var j = 1; j <= n; ++j) {
			for (var i = 0; i <= n - j; ++i)
				next[i] = prev[i] * (1 - t) + prev[i+1] * t;
			var tmp = prev;
			prev = next; next = tmp;
		}
		return prev[0];
	}
}

// ==============> START <================

var c, col_x = 0, col_y = 1, res = 1e-3;

while ((c = getopt(arguments, "x:y:r:")) != null) {
	if (c == 'x') col_x = parseInt(getopt.arg) - 1;
	else if (c == 'y') col_y = parseInt(getopt.arg) - 1;
	else if (c == 'r') res = parseFloat(getopt.arg);
}

if (getopt.ind == arguments.length) {
	print("Usage: k8 max.js [-x 1] [-y 2] [-r 0.001] <table.txt>");
	exit(0);
}

var f = new File(arguments[getopt.ind]);
var s = new Bytes();
var a_x = [], a_y = [];
while (f.readline(s) >= 0) {
	var t = s.toString().split(/[ \t]+/);
	a_x.push(t[col_x]); a_y.push(t[col_y]);
}

var b_x = Math.bernstein_poly(a_x);
var b_y = Math.bernstein_poly(a_y);

var max_y = b_y(0), max_x = b_x(0);
for (var t = res; t <= 1; t += res) {
	var y = b_y(t);
	if (max_y < y) max_y = y, max_x = b_x(t);
}

print(max_x, max_y);