29:null 2a:["$","title",null,{"children":"Halton Sequence | Jacob Rus | Observable"}] 30:T7f3,{ const svgwidth = Math.min(width, 600); const svgheight = svgwidth; const margin = 20; const svg = d3.select(DOM.svg(svgwidth, svgheight)); const x = d3.scaleLinear() .domain([0, 1]) .range([margin, svgwidth-margin]); const y = d3.scaleLinear() .domain([0, 1]) .range([svgwidth-margin, margin]); const b1 = base1; const b2 = base2; const line_positions1 = [[0,.5], [1, .5]]; for (let bn = b1; width / bn > 5; bn *= b1) { d3.range(bn).filter(d => d % b1) .forEach(d => line_positions1.push([d/bn, 1/bn])); } const vertical_lines = svg.append("g").selectAll("path") .data(line_positions1) .enter().append("path") .attr("d", d => { const xpos = x(d[0]); return `M${xpos},${y(0)}L${xpos},${y(1)}`}) .attr("stroke", "#aaa") .attr("opacity", d => Math.pow(d[1], .6)); const line_positions2 = [[0,.5], [1, .5]]; for (let bn = b2; width / bn > 5; bn *= b2) { d3.range(bn).filter(d => d % b2) .forEach(d => line_positions2.push([d/bn, 1/bn])); } const horizontal_lines = svg.append("g").selectAll("path") .data(line_positions2) .enter().append("path") .attr("d", d => { const ypos = y(d[0]); return `M${x(0)},${ypos}L${x(1)},${ypos}`}) .attr("stroke", "#aaa") .attr("opacity", d => Math.pow(d[1], .6)); const dots = svg.append("g") const update = function () { dots.selectAll("circle").remove() dots.selectAll("circle").data(data.slice(0, viewof npts3.value).reverse()) .enter().append("circle") .attr("cx", d=>x(d[0])) .attr("cy", d=>y(d[1])) .attr("fill", (d,i) => (i == 0) ? "#d03" : coloring(d, npts-i)) .attr("r", svgwidth / 150); // TODO: consider figuring out how to properly use D3 data to reuse dot elements } update(); viewof npts3.addEventListener("input", update); return Generators.disposable( svg.node(), () => {viewof npts3.removeEventListener("input", update);}); }31:T732,Explanation = md` ## Explanation A [Halton sequence][] is a sequence of points within a unit ${$`n`}-cube which have [low discrepancy][] (that is, they appear to be randomish but cover the domain uniformly), while being deterministic, fast to compute, and easy to understand and implement. They generalize the 1-dimensional [Van der Corput sequences][vdc sequence].

[Johannes van der Corput][] (1935). [“Verteilungsfunktionen (Erste Mitteilung).”][vdc paper] *Proceedings of the Koninklijke Akademie van Wetenschappen*, 38: 813–821. (In German.)

[John Halton][] (1960). [“On the efficiency of certain quasi-random sequences of points in evaluating multi-dimensional integrals.”][halton paper] *Numerische Mathematik*, 2(1): 84-90.

These and other low-discrepancy sequences are commonly used in [Monte Carlo simulation][monte carlo], [sampling][supersampling] for computer graphics, and similar applications, where they often result in faster convergence than using completely random sequences. Compare Halton points (left) vs. pseudorandom coordinates (right): [Halton sequence]: https://en.wikipedia.org/wiki/Halton_sequence [low discrepancy]: https://en.wikipedia.org/wiki/Low-discrepancy_sequence [vdc sequence]: https://en.wikipedia.org/wiki/Van_der_Corput_sequence [Johannes van der Corput]: https://en.wikipedia.org/wiki/Johannes_van_der_Corput [vdc paper]: http://www.dwc.knaw.nl/DL/publications/PU00016771.pdf [John Halton]: http://www.cs.unc.edu/~halton/ [halton paper]: https://link.springer.com/article/10.1007%2FBF01386213 [monte carlo]: https://en.wikipedia.org/wiki/Quasi-Monte_Carlo_method [supersampling]: https://en.wikipedia.org/wiki/Supersampling `32:T66d,{ const svgwidth = Math.min(width, 300); const margin = 10; const svg1 = d3.select(DOM.svg(svgwidth, svgwidth)); const svg2 = d3.select(DOM.svg(svgwidth, svgwidth)); const x = d3.scaleLinear() .domain([0, 1]) .range([margin, svgwidth-margin]); const y = d3.scaleLinear() .domain([0, 1]) .range([svgwidth-margin, margin]); const blue = (d, i) => d3.interpolateBlues(Math.pow(i/npts, 5)); const red = (d, i) => d3.interpolateReds(Math.pow(i/npts, 5)); const data_pi = data.slice(0, 2000).map(d => [...d, in_circle(d)]).reverse(); const prngdata_pi = prngdata.slice(0, 2000).map(d => [...d, in_circle(d)]).reverse(); svg2.append("g").selectAll("circle") .data(prngdata_pi) .enter().append("circle") .attr("cx", d=>x(d[0])) .attr("cy", d=>y(d[1])) .attr("fill", (d,i) => d[2] ? blue(d, i+2000) : red(d, i+2000)) .attr("r", svgwidth / 150); svg1.append("g").selectAll("circle") .data(data_pi) .enter().append("circle") .attr("cx", d=>x(d[0])) .attr("cy", d=>y(d[1])) .attr("fill", (d,i) => d[2] ? blue(d, i+2000) : red(d, i+2000)) .attr("r", svgwidth / 150); const pi1 = data_pi.reduce((s, d) => s + d[2], 0) / data_pi.length * 4; const pi2 = prngdata_pi.reduce((s, d) => s + d[2], 0) / data_pi.length * 4; return html`

${svg1.node()}
${$`\pi \approx ${pi1.toFixed(4)}`}

${svg2.node()}
${$`\pi \approx ${pi2.toFixed(4)}`}

`; }33:T5e0,{ const npts2 = 200; const svgwidth = Math.min(width, 600); const svgheight = svgwidth*2; const margin = 20; const svg = d3.select(DOM.svg(svgwidth, svgheight)); const x = d3.scaleLinear() .domain([0, 1]) .range([margin, svgwidth-margin]); const y = d3.scaleLinear() .domain([0, npts2]) .range([margin, svgheight-margin]); const b = base1; const horizontal_stripes = svg.append("g").selectAll("path") .data([[1,b], [b*b, b*b*b], [b*b*b*b, b*b*b*b*b], [b*b*b*b*b*b, b*b*b*b*b*b*b]]) .enter().append("path") .attr("d", (d, i) => { return `M${x(0)},${y(d[0])}L${x(1)},${y(d[0])}L${x(1)},${y(d[1])}L${x(0)},${y(d[1])}Z`}) .attr("fill", "#f3f6ff"); const line_positions = [[0,.5], [1, .5]]; for (let bn = b; width / bn > 5; bn *= b) { d3.range(bn).filter(d => d % b) .forEach(d => line_positions.push([d/bn, 1/bn])); } const vertical_lines = svg.append("g").selectAll("path") .data(line_positions) .enter().append("path") .attr("d", d => { const xpos = x(d[0]); return `M${xpos},${margin}L${xpos},${svgheight-margin}`}) .attr("stroke", "#aaa") .attr("opacity", d => Math.pow(d[1], .6)); const dots = svg.append("g").selectAll("circle") .data(data.slice(0, npts2)) .enter().append("circle") .attr("cx", d=>x(d[0])) .attr("cy", (d,i)=>y(i)) .attr("fill", (d,i) => coloring(d, npts-i)) .attr("r", svgwidth / 150); return svg.node(); }34:[] 2d:{"id":"d94c3e277cd0ff97","slug":"halton","title":"Halton Sequence","thumbnail":"94d366d29cfd4d601bda3ddf834ce651a0a60f87652bc61229d78ef44655f7a9","access_level":"public","create_time":"$D2018-05-25T19:05:18.691Z","update_time":"$D2018-06-04T23:42:37.195Z","view_time":null,"first_public_time":"$D2018-05-25T19:30:16.786Z","fork":false,"stars":15,"starred":false,"owner":{"type":"team","id":"9b8a8e177a49059a","login":"jrus","name":"Jacob Rus","avatar_url":"https://avatars-next.observableusercontent.com/c33dbeed76527c50e64208b0e005e0b6e54a3bef2a14ea6a823f1c7de4867662","tier":"starter_2024"},"version":1131,"publish_version":1131,"role":null,"authored":false,"authors":[{"type":"individual","id":"ddd7235babc79ce0","login":"jrus","name":"Jacob Rus","avatar_url":"https://avatars-next.observableusercontent.com/c33dbeed76527c50e64208b0e005e0b6e54a3bef2a14ea6a823f1c7de4867662","tier":null}],"members":[{"type":"individual","id":"ddd7235babc79ce0","login":"jrus","name":"Jacob Rus","avatar_url":"https://avatars-next.observableusercontent.com/c33dbeed76527c50e64208b0e005e0b6e54a3bef2a14ea6a823f1c7de4867662","tier":null,"role":"editor"}],"comments":[{"id":"4e2812928a6bc670","reply_to_id":null,"content":"This is a really compelling example of how one can be more useful than the other.","create_time":"$D2022-06-23T18:50:33.599Z","update_time":null,"resolved":true,"user":{"type":"individual","id":"dbf69c13a14d2f62","login":"visnup","name":"Visnu Pitiyanuvath","avatar_url":"https://avatars-next.observableusercontent.com/b847c4c49993b1df7599d355f47dc0ad9c2dcf43439ce9e8f3410e6368b3aac9","tier":null}},{"id":"d7f0fe13d7de2dde","reply_to_id":"4e2812928a6bc670","content":"Or perhaps just a demonstration that monte carlo simulation is a really bad way to calculate π, and picking random points can make it still worse. ;-)","create_time":"$D2022-06-23T21:01:59.514Z","update_time":null,"resolved":true,"user":{"type":"individual","id":"ddd7235babc79ce0","login":"jrus","name":"Jacob Rus","avatar_url":"https://avatars-next.observableusercontent.com/c33dbeed76527c50e64208b0e005e0b6e54a3bef2a14ea6a823f1c7de4867662","tier":null}}],"custom_thumbnail":false,"trashed":false,"subscribed":null,"forks":7,"collection":null,"type":"notebook","body":{"title":"Halton Sequence","cells":[{"id":0,"value":"Title = md`# Halton Sequence`","pinned":false,"hidden":false,"mode":"ojs"},{"id":13,"value":"plot()","pinned":false,"hidden":false,"mode":"ojs"},{"id":71,"value":"md`Pick two [coprime](https://en.wikipedia.org/wiki/Coprime_integers) numbers for the bases (or don’t, and admire the poor results):`","pinned":false,"hidden":false,"mode":"ojs"},{"id":27,"value":"viewof base1 = slider({\n title: \"first base\", \n min: 2, \n max: 30, \n step: 1, \n value: 2, \n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":30,"value":"viewof base2 = slider({\n title: \"second base\", \n min: 2, \n max: 30, \n step: 1, \n value: 3, \n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":538,"value":"$30","pinned":false,"hidden":false,"mode":"ojs"},{"id":599,"value":"{\n const max = 200;\n const text = [\"Start Animation\", \"Pause Animation\"];\n let animate = true;\n const button = html``;\n button.addEventListener('click', () => { animate = !animate; button.value = text[+animate]; });\n \n base1, base2; viewof npts3.value = max - 10;\n yield button;\n while (true) {\n if (animate) {\n viewof npts3.value = (viewof npts3.value + 1) % max;\n }\n yield Promises.delay(viewof framerate.value, button);\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":750,"value":"viewof framerate = slider({\n title: \"animation framerate (ms/frame)\", \n min: 40, \n max: 1000, \n step: 20, \n value: 200, \n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":761,"value":"md`* * *`","pinned":false,"hidden":false,"mode":"ojs"},{"id":56,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":834,"value":"{\n const svgwidth = Math.min(width, 300);\n const margin = 10;\n const svg1 = d3.select(DOM.svg(svgwidth, svgwidth));\n const svg2 = d3.select(DOM.svg(svgwidth, svgwidth));\n\n const x = d3.scaleLinear()\n .domain([0, 1])\n .range([margin, svgwidth-margin]);\n const y = d3.scaleLinear()\n .domain([0, 1])\n .range([svgwidth-margin, margin]);\n \n svg2.append(\"g\").selectAll(\"circle\")\n .data(prngdata.slice().reverse())\n .enter().append(\"circle\")\n .attr(\"cx\", d=>x(d[0]))\n .attr(\"cy\", d=>y(d[1]))\n .attr(\"fill\", coloring)\n .attr(\"r\", svgwidth / 150);\n\n svg1.append(\"g\").selectAll(\"circle\")\n .data(data.slice().reverse())\n .enter().append(\"circle\")\n .attr(\"cx\", d=>x(d[0]))\n .attr(\"cy\", d=>y(d[1]))\n .attr(\"fill\", coloring)\n .attr(\"r\", svgwidth / 150);\n\n return html`${svg1.node()}${svg2.node()}`;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":970,"value":"md`\nAs a simple example of a [Monte Carlo integration][] problem, we can try to\ncompute the value of [${$`\\pi`}][pi] by determining what proportion of points in\nthe unit square:\n\n${$$`\\left(x, y\\right) \\in [0, 1] \\times [0, 1]`}\n\nare also in the disk:\n\n${$$`x^2 + y^2 + \\frac{1}{4} < x+y`}\n\nThis proportion should be ${$`\\pi/4`}.\n\n [Monte Carlo integration]: https://en.wikipedia.org/wiki/Monte_Carlo_integration\n [pi]: https://en.wikipedia.org/wiki/Pi#Monte_Carlo_methods\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1115,"value":"in_circle = ([x, y]) => x*x + y*y + 0.25 < x + y","pinned":true,"hidden":false,"mode":"ojs"},{"id":987,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":878,"value":"md`\nAt the beginning of the sequence, the entries follow a straight line\nfrom the corner, so if only a few entries will be used (especially\nwith larger bases) it is common to drop the first few terms and start\nat a higher index, scramble the indices in some way, or sometimes to\nonly use every ${$`n`}th element of the sequence.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1090,"value":"md`And now, without further ado, the code we’ve all been waiting for:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":4,"value":"halton = function halton (index, base) {\n let fraction = 1;\n let result = 0;\n while (index > 0) {\n fraction /= base;\n result += fraction * (index % base);\n index = ~~(index / base); // floor division\n }\n return result;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":274,"value":"md`\nIn prose, what the above function does is to take a natural number ${$`k`},\nexpress it in positional notation in base ${$`b`}, reverse the order of the\ndigits, and then consider that as a fraction (still in base ${$`b`}).\n\nWe can compare the code above to a simple function which expresses a\nnatural number ${$`k`} in base ${$`b`}: \n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":278,"value":"positional = function positional (number, base) {\n const result = [];\n while (number > 0) {\n result.push(number % base);\n number = ~~(number / base); // floor division\n }\n return result.reverse();\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":290,"value":"{\n const pos = (k, b) =>\n md` \\`positional(${k}, ${b})\\` → ${$`[${positional(+k, b).join(',')}]`} `;\n \n const hal = (k, b) => {\n const fracs = positional(+k, b).reverse()\n .map((d, i) => [d, `\\\\frac{${d}}{${Math.pow(b, i+1)}}`])\n .filter(d => d[0]).map(d=>d[1]) // drop zeros\n .join('+');\n return md` \\`halton(${k}, ${b})\\` → ${$`${fracs}`} `;\n }\n \n return md`\nFor example, \n${pos(1234, 10)}\n${hal(1234, 10)}\n${pos(1234, 2)}\n${hal(1234, 2)}\n${pos(1234, 3)}\n${hal(1234, 3)}\n${pos('0x4d2', 16)}\n${hal('0x4d2', 16)}\n`;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":233,"value":"md`\nWhen you want to make a Halton sequence covering e.g. the unit cube, you call\nthe \\`halton\\` function with 3 different coprime bases, e.g. the following\nfunction will produce an array of ${$`n`} low-discrepancy points in the cube:\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":256,"value":"cube_points = n => d3.range(100, 100 + n).map(\n i => [halton(i, 2), halton(i, 3), halton(i, 5)])","pinned":true,"hidden":false,"mode":"ojs"},{"id":348,"value":"md`\nOr we could make a generator\n(cf. [@mbostock/custom-generators](/@mbostock/custom-generators))\nout of it:\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":350,"value":"{\n function* gen() {\n for (let i = 100;; i++)\n yield [halton(i, 2), halton(i, 3), halton(i, 5)];\n }\n \n const g = gen();\n while (true)\n yield Promises.delay(1000, g.next().value);\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":129,"value":"md`\nTo unpack how the \\`halton\\` function works in practice, we can plot\nthe values \\`halton(k, b1)\\` of the first base ${$`b_1`} (horizonal axis)\nvs. index ${$`k`} (vertical axis, from top to bottom). Try dragging the\nslider above to see how the distribution of points changes as the base\nchanges:\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":133,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":119,"value":"md`\n* * *\n\nNow to actually make the plot at the top of the page...\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":15,"value":"npts = 4000","pinned":true,"hidden":false,"mode":"ojs"},{"id":21,"value":"data = [...Array(npts).keys()].map(\n d => [halton(d, base1), halton(d, base2)])","pinned":true,"hidden":false,"mode":"ojs"},{"id":82,"value":"coloring = (d, i) => d3.interpolatePurples(Math.pow(i/npts, 5))","pinned":true,"hidden":false,"mode":"ojs"},{"id":96,"value":"plot = () => {\n const svgwidth = Math.min(width, 600);\n const margin = 20;\n const svg = d3.select(DOM.svg(svgwidth, svgwidth));\n\n const x = d3.scaleLinear()\n .domain([0, 1])\n .range([margin, svgwidth-margin]);\n const y = d3.scaleLinear()\n .domain([0, 1])\n .range([svgwidth-margin, margin]);\n \n const dots = svg.append(\"g\").selectAll(\"circle\")\n .data(data.slice().reverse())\n .enter().append(\"circle\")\n .attr(\"cx\", d=>x(d[0]))\n .attr(\"cy\", d=>y(d[1]))\n .attr(\"fill\", coloring)\n .attr(\"r\", svgwidth / 150);\n \n return svg.node();\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":121,"value":"plot()","pinned":true,"hidden":false,"mode":"ojs"},{"id":434,"value":"md`* * *\n\n### Other Implementation Details and Imports`","pinned":false,"hidden":false,"mode":"ojs"},{"id":891,"value":"md`\nImplement our own so-so pseudorandom number generator here, in case the built-in browser \\`Math.random()\\` happens to be bad (this might be unnecessary in any modern browser, but it doesn’t hurt anything to use our own PRNG):\n\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":810,"value":"randfloat = {\n // Use a Multiply With Carry PRNG, with an XOR-shift successor\n // Both from Numerical Recipes, 3rd Edition [H1, G1]\n const prng = function* prng(seed) {\n if (seed == null) seed = +new Date;\n\n // these two numbers were arbitrarily chosen\n let z, x = 3395989511 ^ seed, y = 1716319410 ^ seed;\n\n while (true) {\n x = 62904 * (x & 0xffff) + (x >>> 16);\n y = 41874 * (y & 0xffff) + (y >>> 16);\n z = (x << 16) + y;\n z ^= z >>> 13; z ^= z << 17; z ^= z >>> 5;\n yield z;\n }\n }\n \n const POW_NEG_26 = Math.pow(2, -26);\n const POW_NEG_27 = Math.pow(2, -27);\n\n return function* randfloat(seed) {\n const p = prng(seed);\n while (true) {\n const low_bits = p.next().value >>> 6;\n const high_bits = p.next().value >>> 5;\n yield (high_bits + low_bits * POW_NEG_26) * POW_NEG_27;\n }\n }\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":801,"value":"prngdata = {\n const r = randfloat();\n return d3.range(npts).map(d => [r.next().value, r.next().value]);\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":607,"value":"viewof npts3 = new View(100)","pinned":false,"hidden":false,"mode":"ojs"},{"id":18,"value":"d3 = require(\"d3@5\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":67,"value":"import {slider} from \"@jashkenas/inputs\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":605,"value":"import {View} from \"@mbostock/synchronized-views\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":522,"value":"import {math_css, $, $$} with {$css as $css} from \"@jrus/misc\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":524,"value":"$$css = html`${math_css}`","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W34","stdlib":"1","titling":"auto"}} 2e:"$2d:body"