Skip to content

Instantly share code, notes, and snippets.

@lh3
Created July 19, 2013 18:00
Show Gist options
  • Save lh3/6041085 to your computer and use it in GitHub Desktop.
Save lh3/6041085 to your computer and use it in GitHub Desktop.
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);
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment