Skip to content

Instantly share code, notes, and snippets.

@benallard

benallard/.block

Last active February 1, 2022 12:37
Show Gist options
  • Save benallard/7d66ba5fcdfe0df2969b8a0baf658fa9 to your computer and use it in GitHub Desktop.
Save benallard/7d66ba5fcdfe0df2969b8a0baf658fa9 to your computer and use it in GitHub Desktop.
Download ZIP
Darts
Raw
.block
license: gpl-3.0
height: 800
scrolling: no
border: no
Raw
index.html
<!DOCTYPE html>
<html>
<head>
<meta charset="utf-8" />
<meta viewport="width=460, initial-scale=1.0" />
<script src="https://d3js.org/d3.v7.min.js"></script>
<script src="stats.js"></script>
</head>
<body>
<h1>Darts</h1>
<p>Throw a few darts, aim at the center of the board (Double-Bull), click on the board where the darts lands, and get a personalized heatmap where to aim to maximize your expected score.</p>
<div id="dart_board"></div>
<label for="score" hidden>Throw a dart</label>
<input id="score" placeholder="score" hidden>
<label for="resolution">Resolution: </label>
<input type="range" min="1" max="100" value="50" class="slider" id="resolution">
<button id="reset">Reset</button>
<script>
// set the dimensions and margins of the graph
var margin = {
top: 10,
right: 10,
bottom: 10,
left: 10
},
width = 460 - margin.left - margin.right,
height = 460 - margin.top - margin.bottom;
// append the svg object to the body of the page
var svg = d3.select("#dart_board")
.append("svg")
.attr("width", width + margin.left + margin.right)
.attr("height", height + margin.top + margin.bottom)
.append("g")
.attr("transform",
"translate(" + margin.left + "," + margin.top + ")");
var heatmapLayer = svg.append("g").attr("id", "heatmapLayer");
var boardLayer = svg.append("g").attr("id", "boardLayer");
var dartsLayer = svg.append("g").attr("id", "dartsLayer");
var crossLayer = svg.append("g").attr("id", "crossLayer");
// Add X axis
var x = d3.scaleLinear()
.domain([-R, R])
.range([margin.left, width - margin.right]);
// Add Y axis
var y = d3.scaleLinear()
.domain([-R, R])
.range([height - margin.bottom, margin.top]);
[R1, R2, R3, R4, R5, R].forEach(r =>
boardLayer.append("circle")
.attr("cx", x(0))
.attr("cy", y(0))
.attr("r", x(r) - x(0))
.style("fill", "none")
.style("stroke", "black"));
for (var i = 0; i < 10; i++) {
var th = Math.PI * 11 / 20 - i * Math.PI / 10;
boardLayer.append("line")
.attr("x1", x(R * Math.cos(th)))
.attr("y1", y(R * Math.sin(th)))
.attr("x2", x(R2 * Math.cos(th)))
.attr("y2", y(R2 * Math.sin(th)))
.style("stroke", "black");
boardLayer.append("line")
.attr("x1", x(R * Math.cos(th + Math.PI)))
.attr("y1", y(R * Math.sin(th + Math.PI)))
.attr("x2", x(R2 * Math.cos(th + Math.PI)))
.attr("y2", y(R2 * Math.sin(th + Math.PI)))
.style("stroke", "black");
}
for (var i = 0; i < d.length; i++) {
var th = Math.PI / 2 - i * Math.PI / 10;
boardLayer.append("text")
.attr("x", x((R + 10) * Math.cos(th)))
.attr("y", y((R + 10) * Math.sin(th)))
.attr("text-anchor", "middle")
.attr("alignment-baseline", "middle")
.text(d[i].toString())
.style("font-size", "30px")
.style("font-family", "sans-serif")
.style("font-weight", "bold")
.style("fill", "black");
}
var maxCross = crossLayer
.append('g')
.append('circle')
.style("fill", "none")
.attr("stroke", "blue")
.attr('r', 4)
.style("opacity", 0)
// Create the text that travels along the curve of chart
var maxCrossText = crossLayer
.append('g')
.append('text')
.style("opacity", 0)
.attr("text-anchor", "left")
.attr("alignment-baseline", "middle")
function add_dart(pt) {
darts.push(pt)
heatmap(darts);
}
d3.select("#score").on("change", function() {
var pt = randomPt(+this.value);
add_dart(pt);
this.value = "";
});
d3.select("#dart_board").on("click", function(event) {
var pointer = d3.pointer(event);
pt = [x.invert(pointer[0] - margin.left), y.invert(pointer[1] - margin.top)];
console.log(pt, score(pt[0], pt[1]), getRegionXy(pt[0], pt[1]));
add_dart([x.invert(pointer[0] - margin.left), y.invert(pointer[1] - margin.top)]);
})
var nn = Math.round(2 * R * 50 / 100);
d3.select("#resolution").on("change", function() {
nn = Math.round(2 * R * this.value / 100);
heatmap(darts);
});
d3.select("#reset").on("click", function() {
darts = [];
heatmap(darts);
});
var darts = [];
numImpSamples = 1;
const easy = false;
function heatmap(darts) {
// First display the darts.
dartsLayer.selectAll("circle")
.data(darts)
.join("circle")
.attr("cx", d => x(d[0]))
.attr("cy", d => y(d[1]))
.attr("r", x(1) - x(0))
.style("stroke", "darkgreen")
.style("stroke-width", "2px");
var variance = simpleEM(darts.map(d => score(d[0], d[1])));
console.log("variance:", Math.sqrt(variance));
var res = generalEM(darts, 100, 100, 0, 1, false);
console.log(Math.sqrt(res[0]));
console.log(Math.sqrt(res[1]));
console.log(res[2] / Math.sqrt(res[0] * res[1]));
var map = [];
if (easy) {
map = computeExpScores(variance, nn);
} else {
map = computeExpScores(res[0], res[1], res[2], nn);
}
const args = getMaxAndArgMax(map);
const color = d3.scaleLinear()
.domain([0, args[0] / 2, args[0]])
.range(["red", "yellow", "white"]);
const xMap = d3.scaleBand()
.domain([...Array(nn).keys()])
.range([margin.left, width - margin.right]);
const yMap = d3.scaleBand()
.domain([...Array(nn).keys()])
.range([height - margin.top, margin.bottom]);
heatmapLayer.selectAll(".row")
.data(map)
.join("g")
.attr("class", "row")
.selectAll(".cell")
.data(function(d, i) {
return d.map(function(a) {
return {
value: a,
row: i
};
})
})
.join("rect")
.attr("class", "cell")
.attr("x", (d, i) => xMap(d.row))
.attr("y", (d, i) => yMap(i))
.attr("width", xMap.bandwidth())
.attr("height", yMap.bandwidth())
.style("fill", d => color(d.value));
maxCross.attr("cx", xMap(args[1]))
.attr("cy", yMap(args[2]))
.style("opacity", 1);
maxCrossText.attr("x", xMap(args[1]))
.attr("y", yMap(args[2]))
.text(args[0].toString())
.style("opacity", 1);
}
</script>
<p>Reference: </p><a href="http://www.stat.cmu.edu/~ryantibs/darts/">A Statistician Plays Darts</a>
</body>
</html>
Raw
stats.js
// Rewrite of Stats.java
// hope I'll get it.
// board dimensions in mm (measured from center of board)
const R1 = 6.35; // distance to double bullseye
const R2 = 15.9; // distance to single bullseye
const R3 = 99; // distance to inside triple ring
const R4 = 107; // distance to outside triple ring
const R5 = 162; // distance to inside double ring
const R = 170; // distance to outside double ring
// ordering of the numbers
const d = [20, 1, 18, 4, 13, 6, 10, 15, 2, 17, 3, 19, 7, 16, 8, 11, 14, 9, 12, 5];
// index ordering of the numbers (unit-based)
const ii = [2, 9, 11, 4, 20, 6, 13, 15, 18, 7, 16, 19, 5, 17, 8, 14, 10, 3, 12, 1];
// Parameters
// scores: an array of scores
// sInit: initial guess for the variance
// numIter: the number of iterations to run EM algorithm
//
// Returns an estimate of the variance
function simpleEM(scores, sInit = 100, numIter = 100) {
var s = sInit;
for (var i = 0; i < numIter; i++) {
s = simpleStep(scores, s);
}
return s;
}
// Parameters
// scores: an array of scores
// s: current estimate of variance
function simpleStep(scores, s) {
var a = [];
var b = [];
// calculate the integral (pg.6 in paper) over each region
// note: the regions are defined as follows
// 0 - double bullseye
// 1 - single bullseye
// 2 - single
// 3 - double ring
// 4 - triple ring
// 5 - off the board
a[0] = integNumerator(s, 0, R1);
a[1] = integNumerator(s, R1, R2);
a[2] = integNumerator(s, R2, R3) / 20 + integNumerator(s, R4, R5) / 20;
a[3] = integNumerator(s, R5, R) / 20;
a[4] = integNumerator(s, R3, R4) / 20;
a[5] = integNumerator(s, R, -1);
// calculate the integral over each region
b[0] = integDenominator(s, 0, R1);
b[1] = integDenominator(s, R1, R2);
b[2] = integDenominator(s, R2, R3) / 20 + integDenominator(s, R4, R5) / 20;
b[3] = integDenominator(s, R5, R) / 20;
b[4] = integDenominator(s, R3, R4) / 20;
b[5] = integDenominator(s, R, -1);
var n = scores.length;
var e = 0;
for (var i = 0; i < n; i++) {
e += computeExp(scores[i], a, b);
}
return e / (2 * n);
}
// Parameters
// x: the score we observe
// a: the numerators we computed for each region
// b: the denominators for each region
function computeExp(x, a, b) {
// return the appropriate expectation, based on how the score can be achieved
if ([1, 5, 7, 11, 13, 17, 19].includes(x)) {
// these numbers can only be achieved by hitting a single
return a[2] / b[2];
} else if ([2, 4, 8, 10, 14, 16, 20].includes(x)) {
// single or double
return (a[2] + a[3]) / (b[2] + b[3]);
} else if ([3, 9, 15].includes(x)) {
// single or triple
return (a[2] + a[4]) / (b[2] + b[4]);
} else if ([6, 12, 18].includes(x)) {
// single, double or triple
return (a[2] + a[3] + a[4]) / (b[2] + b[3] + b[4]);
} else if ([24, 30, 36].includes(x)) {
// double or triple
return (a[3] + a[4]) / (b[3] + b[4]);
} else if ([22, 26, 28, 32, 34, 38, 40].includes(x)) {
// double only
return a[3] / b[3];
} else if ([21, 27, 33, 39, 42, 45, 48, 51, 54, 57, 60].includes(x)) {
// triple only
return a[4] / b[4];
} else if (x == 25) {
// single bullseye
return a[1] / b[1];
} else if (x == 50) {
// double bullseye
return a[0] / b[0];
} else {
// off the board
return a[5] / b[5];
}
}
// Parameters
// s: variance
// r1: lower bound of integration
// r2: upper bound of integration
function integNumerator(s, r1, r2) {
if (r2 == -1) { // r2 is assumed to be infinity
return (r1 * r1 + 2 * s) * Math.exp(-r1 * r1 / (2 * s));
} else {
return (r1 * r1 + 2 * s) * Math.exp(-r1 * r1 / (2 * s)) - (r2 * r2 + 2 * s) * Math.exp(-r2 * r2 / (2 * s));
}
}
// Parameters
// s: variance
// r1: lower bound of integration
// r2: upper bound of integration
function integDenominator(s, r1, r2) {
if (r2 == -1) { // r2 is assumed to be infinity
return Math.exp(-r1 * r1 / (2 * s));
} else {
return Math.exp(-r1 * r1 / (2 * s)) - Math.exp(-r2 * r2 / (2 * s));
}
}
/*
* EM for the general model
*/
// number of importance samples
var numImpSamples = 1000;
// random number generator for importance samplce
const rand = Math.random;
// board areas used for importance sampling
const ar1 = Math.PI * R1 * R1; // double bullseye
const ar2 = Math.PI * (R2 * R2 - R1 * R1); // single bullseye
const ar3 = Math.PI * (R3 * R3 - R2 * R2) / 20; // lower single
const ar4 = Math.PI * (R5 * R5 - R4 * R4) / 20; // upper single
const ar5 = Math.PI * (R * R - R5 * R5) / 20; // double ring
const ar6 = Math.PI * (R4 * R4 - R3 * R3) / 20; // triple ring
// Parameters
// scores: scores array
// s1Init, s2Init, s3Init: initial values for s1 (x variance),
// s2 (y variance), and s3 (correlation)
// numIter: number of steps to run the EM algorithm
//
// Returns an array of the estimates: {s1, s2, s3}
function generalEM(scores, s1Init = 100, s2Init = 100, s3Init = 0, numIter = 150, random = true) {
var s1 = s1Init;
var s2 = s2Init;
var s3 = s3Init;
for (var i = 0; i < numIter; i++) {
var A = generalStep(scores, s1, s2, s3, random);
s1 = A[0];
s2 = A[1];
s3 = A[2];
}
return [s1, s2, s3];
}
// Parameters
// scores: scores array
// s1, s2, s3: current estimate of covariance matrix
//
// Returns an array of updated estimates {s1, s2, s3}
function generalStep(scores, s1, s2, s3, random = true) {
var n = scores.length;
var A = [0, 0, 0]
for (var i = 0; i < n; i++) {
var B = simulateExp(scores[i], s1, s2, s3, random);
A[0] += B[0];
A[1] += B[1];
A[2] += B[2];
}
return A.map((x) => x / n);
}
// Parameters
// x: the score we observe
// s1, s2, s3: current estimate of covariance
function simulateExp(x, s1, s2, s3, random = true) {
var det = s1 * s2 - s3 * s3;
var A = [0, 0, 0];
var W = 0;
for (var i = 0; i < numImpSamples; i++) {
var z;
if (random) {
z = randomPt(x);
} else {
z = x;
}
var w = Math.exp(-(s2 * z[0] * z[0] -
2 * s3 * z[0] * z[1] +
s1 * z[1] * z[1]) /
(2 * det));
A[0] += w * z[0] * z[0];
A[1] += w * z[1] * z[1];
A[2] += w * z[0] * z[1];
W += w;
}
return A.map((x) => x / W);
}
// Parameters
// x: the score we observe
function randomPt(x) {
var u = rand();
if ([1, 5, 7, 11, 13, 17, 19].includes(x)) {
// single only
if (u < ar3 / (ar3 + ar4)) {
return randomSlicePt(x, R2, R3);
} else {
return randomSlicePt(x, R4, R5);
}
} else if ([2, 4, 8, 10, 14, 16, 20].includes(x)) {
// single or double
if (u < ar5 / (ar3 + ar4 + ar5)) {
return randomSlicePt(x / 2, R5, R);
} else if (u < (ar5 + ar3) / (ar3 + ar4 + ar5)) {
return randomSlicePt(x, R2, R3);
} else {
return randomSlicePt(x, R4, R5);
}
} else if ([3, 9, 15].includes(x)) {
// single or triple
if (u < ar6 / (ar3 + ar4 + ar6)) {
return randomSlicePt(x / 3, R3, R4);
} else if (u < (ar6 + ar3) / (ar3 + ar4 + ar6)) {
return randomSlicePt(x, R2, R3);
} else {
return randomSlicePt(x, R4, R5);
}
} else if ([6, 12, 18].includes(x)) {
// single, double or triple
if (u <= ar6 / (ar6 + ar5 + ar3 + ar4)) {
return randomSlicePt(x / 3, R3, R4);
} else if (u <= (ar6 + ar5) /
(ar6 + ar5 + ar3 + ar4)) {
return randomSlicePt(x / 2, R5, R);
} else if (u <= (ar6 + ar5 + ar3) /
(ar6 + ar5 + ar3 + ar4)) {
return randomSlicePt(x, R2, R3);
} else {
return randomSlicePt(x, R4, R5);
}
} else if ([24, 30, 36].includes(x)) {
// double or triple
if (u <= ar6 / (ar6 + ar5)) {
return randomSlicePt(x / 3, R3, R4);
} else {
return randomSlicePt(x / 2, R5, R);
}
} else if ([22, 26, 28, 32, 34, 38, 40].includes(x)) {
// double only
return randomSlicePt(x / 2, R5, R);
} else if ([21, 27, 33, 39, 42, 45, 48, 51, 54, 57, 60].includes(x)) {
// triple only
return randomSlicePt(x / 3, R3, R4);
} else if (x == 25) {
// single bullseye
return randomCirclePt(R1, R2);
} else if (x == 50) {
// double bullseye
return randomCirclePt(0, R1);
} else {
// outside the board
// a bit of a cheat, the second number should be infinity
return randomCirclePt(R, 10 * R);
}
}
// Parameters
// x: the score we observe
// r1, r2: lower and upper limits
function randomSlicePt(x, r1, r2) {
var k = ii[x - 1];
var u = rand();
var th = -2 * Math.PI / 40 + (k - 1) * 2 * Math.PI / 20 + 2 * Math.PI / 20 * u;
th = Math.PI / 2 - th;
var r = randomR(r1, r2);
return [r * Math.cos(th), r * Math.sin(th)];
}
// Parameters:
// r1, r2: lower and upper limits
function randomCirclePt(r1, r2) {
var u = rand();
var th = 2 * Math.PI * u;
var r = randomR(r1, r2);
return [r * Math.cos(th), r * Math.sin(th)];
}
// Parameters
// r1, r2: lower and upper limits
function randomR(r1, r2) {
var u = rand();
return Math.sqrt(r1 * r1 + (r2 * r2 - r1 * r1) * u);
}
/*
* Funtions to compute grid of expected scores
*/
// powers of 2
const pow2 = [1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024];
// Parameters
// s: variance
// n: size of grid (assuming square, and a power of 2)
//
// Returns an n x n grid of expected scores
function computeExpScores1(s, n) {
// round up to the nearest power of 2
var m = 1;
for (i = 0; i < pow2.length; i++) {
m = pow2[i];
if (n <= m) {
break;
}
}
// note: due to the design of all the FFT code (taken from Numerical
// Recipes in C), all of our arrays must be unit-based!
// number of millimeters per pixel
var c = 2 * R / n;
var S = Array.from(Array(2 * m + 1), () => new Array(2 * m + 1).fill(0));
for (var i = 1; i <= 2 * m; i++) {
for (var j = 1; j <= 2 * m; j++) {
S[i][j] = score((i - m) * c, (j - m) * c);
}
}
// helpful constants
var n1 = Math.floor(n / 2)
var n2 = Math.ceil(n / 2)
// if the variance is 0, just return the scores array!
if (s == 0) {
var E = Array.from(Array(n), () => new Array(n).fill(0));
for (var i = m - n1 + 1; i < m + n2; i++) {
for (var j = m - n1 + 1; j < m + n2; j++) {
E[i - m + n1 - 1][j - m + n1 - 1] = S[i][j];
}
}
return E;
}
// build the Gaussian density arrays
var g1 = new Array(2 * m + 1).fill(0);
var g2 = new Array(2 * m + 1).fill(0);
for (var i = 1; i <= 2 * m; i++) {
g1[i] = Math.exp(-(i - m) * (i - m) * c * c / (2 * s)) / Math.sqrt(2 * Math.PI * s) * c;
g2[i] = g1[i];
}
// compute the FT of g1 and g2
var g1f = new Array(4 * m + 1).fill(0);
var g2f = new Array(4 * m + 1).fill(0);
twofft(g1, g2, g1f, g2f, 2 * m);
// compute the FT of each row of S
var A = Array.from(Array(2 * m + 1), () => new Array(4 * m + 1).fill(0));
for (var i = 1; i <= 2 * m - 1; i += 2) {
twofftrow(S, i, i + 1, A, 2 * m);
}
// multiply every row of A by g2f
var re, im;
for (var i = 1; i <= 2 * m; i += 1) {
for (var j = 1; j <= 4 * m - 1; j += 2) {
re = A[i][j] * g2f[j] - A[i][j + 1] * g2f[j + 1];
im = A[i][j] * g2f[j + 1] + A[i][j + 1] * g2f[j];
A[i][j] = re;
A[i][j + 1] = im;
}
}
// compute the inverse FT of every row of A
for (var i = 1; i <= 2 * m; i++) {
four1row(A, i, 2 * m, -1);
}
// take the real part of each row, and switch around some columns
var AA = Array.from(Array(2 * m + 1), () => new Array(2 * m + 1).fill(0));
for (var i = 1; i <= 2 * m; i++) {
for (var j = 1; j <= m; j++) {
AA[i][j + m] = A[i][2 * j - 1] / (2 * m);
}
for (var j = m + 1; j <= 2 * m; j++) {
AA[i][j - m] = A[i][2 * j - 1] / (2 * m);
}
}
// compute the FT of every column of AA
var AAA = Array.from(Array(4 * m + 1), () => new Array(2 * m + 1).fill(0));
for (var j = 1; j <= 2 * m - 1; j += 2) {
twofftcol(AA, j, j + 1, AAA, 2 * m);
}
// multiply every column of AAA by g1f
for (var j = 1; j <= 2 * m; j++) {
for (var i = 1; i <= 4 * m - 1; i += 2) {
re = AAA[i][j] * g1f[i] - AAA[i + 1][j] * g1f[i + 1];
im = AAA[i][j] * g1f[i + 1] + AAA[i + 1][j] * g1f[i];
AAA[i][j] = re;
AAA[i + 1][j] = im;
}
}
// compute the inverse FT of every column of AAA
for (var j = 1; j <= 2 * m; j++) {
four1col(AAA, j, 2 * m, -1);
}
// take the real part of each column, switch around some rows,
// and strip away some part of the array on each side
var E = Array.from(Array(n), () => new Array(n).fill(0));
for (var i = 1; i <= n1; i++) {
for (var j = m - n1 + 1; j <= m + n2; j++) {
E[i + n2 - 1][j - m + n1 - 1] = AAA[2 * i - 1][j] / (2 * m);
}
}
for (var i = 2 * m - n2 + 1; i <= 2 * m; i++) {
for (var j = m - n1 + 1; j <= m + n2; j++) {
E[i - 2 * m + n2 - 1][j - m + n1 - 1] = AAA[2 * i - 1][j] / (2 * m);
}
}
// help the gc out
S = null;
g1 = null;
g2 = null;
g1f = null;
g2f = null;
A = null;
AA = null;
AAA = null;
return E;
}
// Parameters
// s1, s2, s3: x variance, y variance, covariance
// n: size of grid (assuming square, and a power of 2)
//
// Returns an n x n grid of expected scores
//
// Note: as opposed to the simple version of this function:
// computeExpScores(double s, int n)
// this function is not really optimized for speed. The rationale
// here is that it will take long enough to run the general EM
// algorithm anyway.
function computeExpScores(s1, s2, s3, n) {
if (s3 === undefined) {
// poor man function overloading
return computeExpScores1(s1, s2);
}
// round up to the nearest power of 2
var m = 1;
for (var i = 0; i < pow2.length; i++) {
m = pow2[i];
if (n <= pow2[i]) {
break;
}
}
// note: due to the design of all the FFT code (taken from Numerical
// Recipes in C), all of our arrays must be unit-based!
// another note: all of our matrices here are complex. the layout is that
// there are twice as many columns as a corresponding real array representing
// the same dimension, so that m[i][j] and m[i][j+1] represent corresponding
// real and complex parts of the same number
// number of millimeters per pixel
var c = 2 * R / n;
// build the score array
var S = Array.from(Array(2 * m + 1), () => new Array(4 * m + 1).fill(0));
for (var i = 1; i <= 2 * m; i++) {
for (var j = 1; j <= 2 * m; j++) {
S[i][2 * j - 1] = score((i - m) * c, (j - m) * c);
S[i][2 * j] = 0;
}
}
// helpful constants
var n1 = Math.floor(n / 2);
var n2 = Math.ceil(n / 2);
// build the Gaussian density array
var G = Array.from(Array(2 * m + 1), () => new Array(4 * m + 1).fill(0));
var det = s1 * s2 - s3 * s3;
for (var i = 1; i <= 2 * m; i++) {
for (var j = 1; j <= 2 * m; j++) {
G[i][2 * j - 1] = Math.exp(-(s2 * (i - m) * (i - m) -
2 * s3 * (i - m) * (j - m) +
s1 * (j - m) * (j - m)) * c * c / (2 * det)) /
(2 * Math.PI * Math.sqrt(det)) * c * c;
G[i][2 * j] = 0;
}
}
// compute the FT of each matrix (in-place)
four2(S, 2 * m, 1);
four2(G, 2 * m, 1);
// multiply S and G together (store the result in S)
var re, im;
for (var i = 1; i <= 2 * m; i += 1) {
for (var j = 1; j <= 4 * m - 1; j += 2) {
re = S[i][j] * G[i][j] - S[i][j + 1] * G[i][j + 1];
im = S[i][j] * G[i][j + 1] + S[i][j + 1] * G[i][j];
S[i][j] = re;
S[i][j + 1] = im;
}
}
// compute the inverse FT of S (in-place)
four2(S, 2 * m, -1);
// switch around some rows and some columns, take the real part
var EE = Array.from(Array(2 * m + 1), () => new Array(2 * m + 1).fill(0));
for (var i = 1; i <= m; i++) {
for (var j = 1; j <= m; j++) {
EE[i + m][j + m] = S[i][2 * j - 1] / (4 * m * m);
}
for (var j = m + 1; j <= 2 * m; j++) {
EE[i + m][j - m] = S[i][2 * j - 1] / (4 * m * m);
}
}
for (var i = m + 1; i <= 2 * m; i++) {
for (var j = 1; j <= m; j++) {
EE[i - m][j + m] = S[i][2 * j - 1] / (4 * m * m);
}
for (var j = m + 1; j <= 2 * m; j++) {
EE[i - m][j - m] = S[i][2 * j - 1] / (4 * m * m);
}
}
// trim away on each side
var E = Array.from(Array(n), () => new Array(n).fill(0));
for (var i = m - n1; i <= m + n2 - 1; i++) {
for (var j = m - n1; j <= m + n2 - 1; j++) {
E[i - m + n1][j - m + n1] = EE[i][j];
}
}
return E;
}
function score(x, y) {
// compute the radius
var r = Math.sqrt(x * x + y * y);
// check if it's off the board (do this for speed)
if (r > R) {
return 0;
}
// check if it's a double bullseye
if (r <= R1) {
return 50;
}
// check for a single bullseye
if (r <= R2) {
return 25;
}
// get the angle
var th = Math.atan2(y, x);
var phi = (Math.PI / 2 + Math.PI / 20 - th) % (2 * Math.PI);
if (phi < 0) {
phi += 2 * Math.PI;
}
// now get the number
var i = Math.floor(phi / (2 * Math.PI / 20));
if (i < 0) {
i = 0;
}
if (i >= 19) {
i = 19;
}
var n = d[i]
// check for a single
if (r <= R3) {
return n;
}
// check for a triple
if (r <= R4) {
return 3 * n;
}
// check for a single
if (r <= R5) {
return n;
}
// now it must be a double
return 2 * n;
}
function getRegionXy(x, y) {
// compute the radius
var r = Math.sqrt(x * x + y * y);
// check if it's off the board (do this for speed)
if (r > R) {
return "Off";
}
// check if it's a double bullseye
if (r <= R1) {
return "DB";
}
// check for a single bullseye
if (r <= R2) {
return "SB";
}
// get the angle
var th = Math.atan2(y, x);
var phi = (Math.PI / 2 + Math.PI / 20 - th) % (2 * Math.PI);
if (phi < 0) {
phi += 2 * Math.PI;
}
// now get the number
var i = Math.floor(phi / (2 * Math.PI / 20));
if (i < 0) {
i = 0;
}
if (i >= 19) {
i = 19;
}
var n = d[i]
// check for a single
if (r <= R3) {
return "S" + n;
}
// check for a triple
if (r <= R4) {
return "T" + n;
}
// check for a single
if (r <= R5) {
return "S" + n;
}
// now it must be a double
return "D" + n;
}
function getRegion(i, j, n1, n2) {
if (n1 === undefined) {
return getRegionXy(i, j);
}
return getRegion((i + 1 - n1 / 2) * 2 * R / n1, (j + 2 - n2 / 2) * 2 * R / n2);
}
function getMaxAndArgMax(data) {
var max = data[0][0];
var ii = 0,
jj = 0;
for (var i = 0; i < data.length; i++) {
for (var j = 0; j < data[i].length; j++) {
if (data[i][j] > max) {
max = data[i][j];
ii = i;
jj = j;
}
}
}
return [max, ii, jj];
}
/*
* FFT Funtions
* Adapted from Numerical Recipes in C
*/
// Replaces data[1..2*nn] by its discrete Fourier transform, if isign is input
// as 1;
// or replaces data[1..2*nn] by nn times its inverse discrete Fourier transform,
// if
// isign is input as -1.
// data is a complex array of length nn or, equivalently, a real array of length
// 2*nn.
// nn MUST be an integer power of 2 (this is not checked for!).
function four1(data, nn, isign) {
var n, mmax, m, j, istep, i;
var wtemp, wr, wpr, wpi, wi, theta;
var tempr, tempi;
n = nn << 1;
j = 1;
for (i = 1; i < n; i += 2) {
if (j > i) {
// SWAP(data[j], data[i]);
var temp = data[j];
data[j] = data[i];
data[i] = temp;
// SWAP(data[j+1], data[i+1]);
temp = data[j + 1];
data[j + 1] = data[i + 1];
data[i + 1] = temp;
}
m = n >>> 1;
while (m >= 2 && j > m) {
j -= m;
m = m >>> 1;
}
j += m;
}
mmax = 2;
while (n > mmax) {
istep = mmax << 1;
theta = isign * (6.28318530717959 / mmax);
wtemp = Math.sin(0.5 * theta);
wpr = -2.0 * wtemp * wtemp;
wpi = Math.sin(theta);
wr = 1.0;
wi = 0.0;
for (var m = 1; m < mmax; m += 2) {
for (var i = m; i <= n; i += istep) {
j = i + mmax;
tempr = wr * data[j] - wi * data[j + 1];
tempi = wr * data[j + 1] + wi * data[j];
data[j] = data[i] - tempr;
data[j + 1] = data[i + 1] - tempi;
data[i] += tempr;
data[i + 1] += tempi;
}
wr = (wtemp = wr) * wpr - wi * wpi + wr;
wi = wi * wpr + wtemp * wpi + wi;
}
mmax = istep;
}
}
// four1 applied to the kth row of data (a matrix)
function four1row(data, k, nn, isign) {
var n, mmax, m, j, istep, i;
var wtemp, wr, wpr, wpi, wi, theta;
var tempr, tempi;
n = nn << 1;
j = 1;
for (i = 1; i < n; i += 2) {
if (j > i) {
// SWAP(data[k][j], data[k][i]);
var temp = data[k][j];
data[k][j] = data[k][i];
data[k][i] = temp;
// SWAP(data[k][j+1], data[k][i+1]);
temp = data[k][j + 1];
data[k][j + 1] = data[k][i + 1];
data[k][i + 1] = temp;
}
m = n >>> 1;
while (m >= 2 && j > m) {
j -= m;
m = m >>> 1;
}
j += m;
}
mmax = 2;
while (n > mmax) {
istep = mmax << 1;
theta = isign * (6.28318530717959 / mmax);
wtemp = Math.sin(0.5 * theta);
wpr = -2.0 * wtemp * wtemp;
wpi = Math.sin(theta);
wr = 1.0;
wi = 0.0;
for (var m = 1; m < mmax; m += 2) {
for (var i = m; i <= n; i += istep) {
j = i + mmax;
tempr = wr * data[k][j] - wi * data[k][j + 1];
tempi = wr * data[k][j + 1] + wi * data[k][j];
data[k][j] = data[k][i] - tempr;
data[k][j + 1] = data[k][i + 1] - tempi;
data[k][i] += tempr;
data[k][i + 1] += tempi;
}
wr = (wtemp = wr) * wpr - wi * wpi + wr;
wi = wi * wpr + wtemp * wpi + wi;
}
mmax = istep;
}
}
// four1 applied to the kth column of data (a matrix)
function four1col(data, k, nn, isign) {
var n, mmax, m, j, istep, i;
var wtemp, wr, wpr, wpi, wi, theta;
var tempr, tempi;
n = nn << 1;
j = 1;
for (var i = 1; i < n; i += 2) {
if (j > i) {
// SWAP(data[j][k], data[i][k]);
var temp = data[j][k];
data[j][k] = data[i][k];
data[i][k] = temp;
// SWAP(data[j+1][k], data[i+1][k]);
temp = data[j + 1][k];
data[j + 1][k] = data[i + 1][k];
data[i + 1][k] = temp;
}
m = n >>> 1;
while (m >= 2 && j > m) {
j -= m;
m = m >>> 1;
}
j += m;
}
mmax = 2;
while (n > mmax) {
istep = mmax << 1;
theta = isign * (6.28318530717959 / mmax);
wtemp = Math.sin(0.5 * theta);
wpr = -2.0 * wtemp * wtemp;
wpi = Math.sin(theta);
wr = 1.0;
wi = 0.0;
for (var m = 1; m < mmax; m += 2) {
for (var i = m; i <= n; i += istep) {
j = i + mmax;
tempr = wr * data[j][k] - wi * data[j + 1][k];
tempi = wr * data[j + 1][k] + wi * data[j][k];
data[j][k] = data[i][k] - tempr;
data[j + 1][k] = data[i + 1][k] - tempi;
data[i][k] += tempr;
data[i + 1][k] += tempi;
}
wr = (wtemp = wr) * wpr - wi * wpi + wr;
wi = wi * wpr + wtemp * wpi + wi;
}
mmax = istep;
}
}
// Given two real input arrays data1[1..n] and data2[1..n], this routine calls
// four1
// and returns two complex output arrays, fft1[1..2n] and fft2[1..2n], each of
// complex length n (i.e., real length 2*n), which contain the discrete Fourier
// transforms of the respective data arrays. n MUST be an integer power of 2.
function twofft(data1, data2, fft1, fft2, n) {
var nn3, nn2;
var rep, rem, aip, aim;
nn3 = 1 + (nn2 = 2 + n + n);
for (var j = 1, jj = 2; j <= n; j++, jj += 2) {
fft1[jj - 1] = data1[j];
fft1[jj] = data2[j];
}
four1(fft1, n, 1);
fft2[1] = fft1[2];
fft1[2] = fft2[2] = 0.0;
for (var j = 3; j <= n + 1; j += 2) {
rep = 0.5 * (fft1[j] + fft1[nn2 - j]);
rem = 0.5 * (fft1[j] - fft1[nn2 - j]);
aip = 0.5 * (fft1[j + 1] + fft1[nn3 - j]);
aim = 0.5 * (fft1[j + 1] - fft1[nn3 - j]);
fft1[j] = rep;
fft1[j + 1] = aim;
fft1[nn2 - j] = rep;
fft1[nn3 - j] = -aim;
fft2[j] = aip;
fft2[j + 1] = -rem;
fft2[nn2 - j] = aip;
fft2[nn3 - j] = rem;
}
}
// twofft applied to the k1st and k2nd rows of data, results are stored
// in the corresponding rows of fft
function twofftrow(data, k1, k2, fft, n) {
var nn3, nn2, jj, j;
var rep, rem, aip, aim;
nn3 = 1 + (nn2 = 2 + n + n);
for (var j = 1, jj = 2; j <= n; j++, jj += 2) {
fft[k1][jj - 1] = data[k1][j];
fft[k1][jj] = data[k2][j];
}
four1row(fft, k1, n, 1);
fft[k2][1] = fft[k1][2];
fft[k1][2] = fft[k2][2] = 0.0;
for (var j = 3; j <= n + 1; j += 2) {
rep = 0.5 * (fft[k1][j] + fft[k1][nn2 - j]);
rem = 0.5 * (fft[k1][j] - fft[k1][nn2 - j]);
aip = 0.5 * (fft[k1][j + 1] + fft[k1][nn3 - j]);
aim = 0.5 * (fft[k1][j + 1] - fft[k1][nn3 - j]);
fft[k1][j] = rep;
fft[k1][j + 1] = aim;
fft[k1][nn2 - j] = rep;
fft[k1][nn3 - j] = -aim;
fft[k2][j] = aip;
fft[k2][j + 1] = -rem;
fft[k2][nn2 - j] = aip;
fft[k2][nn3 - j] = rem;
}
}
// twofft applied to the k1st and k2nd columns of data, results are stored
// in the corresponding columns of fft
function twofftcol(data, k1, k2, fft, n) {
var nn3, nn2, jj, j;
var rep, rem, aip, aim;
nn3 = 1 + (nn2 = 2 + n + n);
for (var j = 1, jj = 2; j <= n; j++, jj += 2) {
fft[jj - 1][k1] = data[j][k1];
fft[jj][k1] = data[j][k2];
}
four1col(fft, k1, n, 1);
fft[1][k2] = fft[2][k1];
fft[2][k1] = fft[2][k2] = 0.0;
for (var j = 3; j <= n + 1; j += 2) {
rep = 0.5 * (fft[j][k1] + fft[nn2 - j][k1]);
rem = 0.5 * (fft[j][k1] - fft[nn2 - j][k1]);
aip = 0.5 * (fft[j + 1][k1] + fft[nn3 - j][k1]);
aim = 0.5 * (fft[j + 1][k1] - fft[nn3 - j][k1]);
fft[j][k1] = rep;
fft[j + 1][k1] = aim;
fft[nn2 - j][k1] = rep;
fft[nn3 - j][k1] = -aim;
fft[j][k2] = aip;
fft[j + 1][k2] = -rem;
fft[nn2 - j][k2] = aip;
fft[nn3 - j][k2] = rem;
}
}
// Replaces data[1..n][1..2*nn] by its discrete 2D Fourier transform, if isign
// is 1; or replaces data by nn times its inverse discrete Fourier transform, if
// isign is -1. note: data is a complex matrix.
// n MUST be an integer power of 2 (this is not checked for!).
function four2(data, nn, isign) {
// compute the 1D FT of every row
for (var i = 1; i <= nn; i++) {
four1row(data, i, nn, isign);
}
// compute the 1D FT of every column
var a = new Array(2 * nn + 1).fill(0);
for (var j = 1; j <= 2 * nn - 1; j += 2) {
//copy into the array a
for (var i = 1; i <= nn; i++) {
a[2 * i - 1] = data[i][j];
a[2 * i] = data[i][j + 1];
}
// compute the 1D FT of a
four1(a, nn, isign);
// copy back into data
for (var i = 1; i <= nn; i++) {
data[i][j] = a[2 * i - 1];
data[i][j + 1] = a[2 * i];
}
}
}
function main() {
const benScores = [
1, 1,
2, 2, 2, 2, 2,
3, 3, 3, 3,
4, 4,
5, 5, 5,
6, 6, 6,
7, 7, 7, 7, 7, 7,
8, 8, 8, 8, 8, 8,
9, 9, 9, 9, 9,
10, 10, 10, 10, 10,
11, 11, 11, 11, 11, 11,
12,
13, 13, 13, 13, 13,
14, 14, 14, 14,
15, 15, 15,
16, 16, 16, 16, 16, 16, 16,
17, 17, 17, 17, 17, 17, 17, 17,
18, 18, 18, 18, 18, 18,
19, 19, 19, 19, 19, 19,
20, 20, 20,
25, 25, 25, 25, 25, 25, 25, 25, 25, 25,
50, 50,
];
console.log(Math.sqrt(simpleEM(benScores)));
var p = generalEM(benScores);
console.log(Math.sqrt(p[0]));
console.log(Math.sqrt(p[1]));
console.log(p[2] / Math.sqrt(p[0] * p[1]));
}
//main();
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment