Skip to content

Instantly share code, notes, and snippets.

@Kcnarf
Last active October 1, 2018 08:12
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save Kcnarf/096cb061195863d2021e291c56daf16e to your computer and use it in GitHub Desktop.
Save Kcnarf/096cb061195863d2021e291c56daf16e to your computer and use it in GitHub Desktop.
weighted KDE & Voronoï maps
license: gpl-3.0

This block experiments two things:

  1. computation of Weighted KDE, an extension of standard KDE; in standard KDE, each datum counts for the same amount (i.e. 1), whereas in weighted KDE, each data counts for a specific amount (i.e. its weight); for example, one can make a standard KDE of the number of sales per day, or a weighted KDE of the total sales' profit per day
  2. fill the weighted KDE curve/area with cells encoding data's weights (light weight -> small cell, heavy weight -> large cell); the objective is to give a sens of the underlying distribution that produces the weighted KDE; I use the d3-voronoï-map plugin to do so, but:
  3. tweek it so that each site's x-coord remains unchanged durring the voronoï map computation (see file d"-voronoi-map-fixed-x.js)
  4. define a specific initial positioning function (exactly encodes x-coord, and computes a random y-coord)

Usage : use the controller to hide/show objects, and hover a cell or a bin for details.

Indeed, the underlying dataset does not suit the experimentation. I have to find another one.

==original README==

Kernel density estimation is a method of estimating the probability distribution of a random variable based on a random sample. In contrast to a histogram, kernel density estimation produces a smooth estimate. The smoothness can be tuned via the kernel’s bandwidth parameter. With the correct choice of bandwidth, important features of the distribution can be seen, while an incorrect choice results in undersmoothing or oversmoothing and obscured features.

This example shows a histogram and a kernel density estimation for times between eruptions of Old Faithful Geyser in Yellowstone National Park, taken from R’s faithful dataset. The data follow a bimodal distribution; short eruptions are followed by a wait time averaging about 55 minutes, and long eruptions by a wait time averaging about 80 minutes. In recent years, wait times have been increasing, possibly due to the effects of earthquakes on the geyser’s geohydrology.

This example is based on a Protovis version by John Firebaugh. See also a two-dimensional density estimation of this dataset using d3-contour.

// Code from github.com/Kcnarf/d3-voronoi-map/v1.2.0/build/d3-voronoi-map.js
// Adapted to fix x-coord of sites during the entire process
// look for 'adaptPosition' function and '//(fixed-x)' prefixed loc.
(function (global, factory) {
typeof exports === 'object' && typeof module !== 'undefined' ? factory(exports, require('d3-polygon'), require('d3-weighted-voronoi')) :
typeof define === 'function' && define.amd ? define(['exports', 'd3-polygon', 'd3-weighted-voronoi'], factory) :
(factory((global.d3 = global.d3 || {}),global.d3,global.d3));
}(this, function (exports,d3Polygon,d3WeightedVoronoi) { 'use strict';
function FlickeringMitigation () {
/////// Inputs ///////
this.growthChangesLength = DEFAULT_LENGTH;
this.totalAvailableArea = NaN;
//begin: internals
this.lastAreaError = NaN;
this.lastGrowth = NaN;
this.growthChanges = [];
this.growthChangeWeights = generateGrowthChangeWeights(this.growthChangesLength); //used to make recent changes weighter than older changes
this.growthChangeWeightsSum = computeGrowthChangeWeightsSum(this.growthChangeWeights);
//end: internals
}
var DEFAULT_LENGTH = 10;
function direction(h0, h1) {
return (h0 >= h1)? 1 : -1;
}
function generateGrowthChangeWeights(length) {
var initialWeight = 3; // a magic number
var weightDecrement = 1; // a magic number
var minWeight = 1;
var weightedCount = initialWeight;
var growthChangeWeights = [];
for (var i=0; i<length; i++) {
growthChangeWeights.push(weightedCount);
weightedCount -= weightDecrement;
if (weightedCount<minWeight) { weightedCount = minWeight; }
}
return growthChangeWeights;
}
function computeGrowthChangeWeightsSum (growthChangeWeights) {
var growthChangeWeightsSum = 0;
for (var i=0; i<growthChangeWeights.length; i++) {
growthChangeWeightsSum += growthChangeWeights[i];
}
return growthChangeWeightsSum;
}
///////////////////////
///////// API /////////
///////////////////////
FlickeringMitigation.prototype.reset = function () {
this.lastAreaError = NaN;
this.lastGrowth = NaN;
this.growthChanges = [];
this.growthChangesLength = DEFAULT_LENGTH;
this.growthChangeWeights = generateGrowthChangeWeights(this.growthChangesLength);
this.growthChangeWeightsSum = computeGrowthChangeWeightsSum(this.growthChangeWeights);
this.totalAvailableArea = NaN;
return this;
};
FlickeringMitigation.prototype.clear = function () {
this.lastAreaError = NaN;
this.lastGrowth = NaN;
this.growthChanges = [];
return this;
};
FlickeringMitigation.prototype.length = function (_) {
if (!arguments.length) { return this.growthChangesLength; }
if (parseInt(_)>0) {
this.growthChangesLength = Math.floor(parseInt(_));
this.growthChangeWeights = generateGrowthChangeWeights(this.growthChangesLength);
this.growthChangeWeightsSum = computeGrowthChangeWeightsSum(this.growthChangeWeights);
} else {
console.warn("FlickeringMitigation.length() accepts only positive integers; unable to handle "+_);
}
return this;
};
FlickeringMitigation.prototype.totalArea = function (_) {
if (!arguments.length) { return this.totalAvailableArea; }
if (parseFloat(_)>0) {
this.totalAvailableArea = parseFloat(_);
} else {
console.warn("FlickeringMitigation.totalArea() accepts only positive numbers; unable to handle "+_);
}
return this;
};
FlickeringMitigation.prototype.add = function (areaError) {
var secondToLastAreaError, secondToLastGrowth;
secondToLastAreaError = this.lastAreaError;
this.lastAreaError = areaError;
if (!isNaN(secondToLastAreaError)) {
secondToLastGrowth = this.lastGrowth;
this.lastGrowth = direction(this.lastAreaError, secondToLastAreaError);
}
if (!isNaN(secondToLastGrowth)) {
this.growthChanges.unshift(this.lastGrowth!=secondToLastGrowth);
}
if (this.growthChanges.length>this.growthChangesLength) {
this.growthChanges.pop();
}
return this;
};
FlickeringMitigation.prototype.ratio = function () {
var weightedChangeCount = 0;
var ratio;
if (this.growthChanges.length < this.growthChangesLength) { return 0; }
if (this.lastAreaError > this.totalAvailableArea/10) { return 0; }
for(var i=0; i<this.growthChangesLength; i++) {
if (this.growthChanges[i]) {
weightedChangeCount += this.growthChangeWeights[i];
}
}
ratio = weightedChangeCount/this.growthChangeWeightsSum;
if (ratio>0) {
console.log("flickering mitigation ratio: "+Math.floor(ratio*1000)/1000);
}
return ratio;
};
function randomInitialPosition () {
//begin: internals
var clippingPolygon,
extent,
minX, maxX,
minY, maxY,
dx, dy;
//end: internals
///////////////////////
///////// API /////////
///////////////////////
function _random(d, i, arr, voronoiMap) {
var shouldUpdateInternals = false;
var x, y;
if (clippingPolygon !== voronoiMap.clip()) {
clippingPolygon = voronoiMap.clip();
extent = voronoiMap.extent();
shouldUpdateInternals = true;
}
if (shouldUpdateInternals) {
updateInternals();
}
x = minX + dx * voronoiMap.prng()();
y = minY + dy * voronoiMap.prng()();
while (!d3Polygon.polygonContains(clippingPolygon, [x, y])) {
x = minX + dx * voronoiMap.prng()();
y = minY + dy * voronoiMap.prng()();
}
return [x, y];
};
///////////////////////
/////// Private ///////
///////////////////////
function updateInternals() {
minX = extent[0][0];
maxX = extent[1][0];
minY = extent[0][1];
maxY = extent[1][1];
dx = maxX - minX;
dy = maxY - minY;
};
return _random;
};
function pie () {
//begin: internals
var startAngle = 0;
var clippingPolygon,
dataArray,
dataArrayLength,
clippingPolygonCentroid,
halfIncircleRadius,
angleBetweenData;
//end: internals
///////////////////////
///////// API /////////
///////////////////////
function _pie(d, i, arr, voronoiMap) {
var shouldUpdateInternals = false;
if (clippingPolygon !== voronoiMap.clip()) {
clippingPolygon = voronoiMap.clip();
shouldUpdateInternals |= true;
}
if (dataArray !== arr) {
dataArray = arr;
shouldUpdateInternals |= true;
}
if (shouldUpdateInternals) {
updateInternals();
}
// add some randomness to prevent colinear/cocircular points
// substract -0.5 so that the average jitter is still zero
return [
clippingPolygonCentroid[0] + Math.cos(startAngle + i * angleBetweenData) * halfIncircleRadius + (voronoiMap.prng()() - 0.5) * 1E-3,
clippingPolygonCentroid[1] + Math.sin(startAngle + i * angleBetweenData) * halfIncircleRadius + (voronoiMap.prng()() - 0.5) * 1E-3
];
};
_pie.startAngle = function (_) {
if (!arguments.length) {
return startAngle;
}
startAngle = _;
return _pie;
};
///////////////////////
/////// Private ///////
///////////////////////
function updateInternals() {
clippingPolygonCentroid = d3Polygon.polygonCentroid(clippingPolygon);
halfIncircleRadius = computeMinDistFromEdges(clippingPolygonCentroid, clippingPolygon) / 2;
dataArrayLength = dataArray.length;
angleBetweenData = 2 * Math.PI / dataArrayLength;
};
function computeMinDistFromEdges(vertex, clippingPolygon) {
var minDistFromEdges = Infinity,
edgeIndex = 0,
edgeVertex0 = clippingPolygon[clippingPolygon.length - 1],
edgeVertex1 = clippingPolygon[edgeIndex];
var distFromCurrentEdge;
while (edgeIndex < clippingPolygon.length) {
distFromCurrentEdge = vDistance(vertex, edgeVertex0, edgeVertex1);
if (distFromCurrentEdge < minDistFromEdges) {
minDistFromEdges = distFromCurrentEdge;
}
edgeIndex++;
edgeVertex0 = edgeVertex1;
edgeVertex1 = clippingPolygon[edgeIndex];
}
return minDistFromEdges;
}
//from https://stackoverflow.com/questions/849211/shortest-distance-between-a-point-and-a-line-segment
function vDistance(vertex, edgeVertex0, edgeVertex1) {
var x = vertex[0],
y = vertex[1],
x1 = edgeVertex0[0],
y1 = edgeVertex0[1],
x2 = edgeVertex1[0],
y2 = edgeVertex1[1];
var A = x - x1,
B = y - y1,
C = x2 - x1,
D = y2 - y1;
var dot = A * C + B * D;
var len_sq = C * C + D * D;
var param = -1;
if (len_sq != 0) //in case of 0 length line
param = dot / len_sq;
var xx, yy;
if (param < 0) { // this should not arise as clippingpolygon is convex
xx = x1;
yy = y1;
} else if (param > 1) { // this should not arise as clippingpolygon is convex
xx = x2;
yy = y2;
} else {
xx = x1 + param * C;
yy = y1 + param * D;
}
var dx = x - xx;
var dy = y - yy;
return Math.sqrt(dx * dx + dy * dy);
}
return _pie;
}
function halfAverageAreaInitialWeight () {
//begin: internals
var clippingPolygon,
dataArray,
siteCount,
totalArea,
halfAverageArea;
//end: internals
///////////////////////
///////// API /////////
///////////////////////
function _halfAverageArea(d, i, arr, voronoiMap) {
var shouldUpdateInternals = false;
if (clippingPolygon !== voronoiMap.clip()) {
clippingPolygon = voronoiMap.clip();
shouldUpdateInternals |= true;
}
if (dataArray !== arr) {
dataArray = arr;
shouldUpdateInternals |= true;
}
if (shouldUpdateInternals) {
updateInternals();
}
return halfAverageArea;
};
///////////////////////
/////// Private ///////
///////////////////////
function updateInternals() {
siteCount = dataArray.length;
totalArea = d3Polygon.polygonArea(clippingPolygon);
halfAverageArea = totalArea / siteCount / 2; // half of the average area of the the clipping polygon
}
return _halfAverageArea;
};
function voronoiMap() {
//begin: constants
var DEFAULT_CONVERGENCE_RATIO = 0.01;
var DEFAULT_MAX_ITERATION_COUNT = 50;
var DEFAULT_MIN_WEIGHT_RATIO = 0.01;
var DEFAULT_PRNG = Math.random;
var DEFAULT_INITIAL_POSITION = randomInitialPosition();
var DEFAULT_INITIAL_WEIGHT = halfAverageAreaInitialWeight();
var RANDOM_INITIAL_POSITION = randomInitialPosition();
var epsilon = 1;
//end: constants
/////// Inputs ///////
var weight = function (d) {
return d.weight;
}; // accessor to the weight
var convergenceRatio = DEFAULT_CONVERGENCE_RATIO; // targeted allowed error ratio; default 0.01 stops computation when cell areas error <= 1% clipping polygon's area
var maxIterationCount = DEFAULT_MAX_ITERATION_COUNT; // maximum allowed iteration; stops computation even if convergence is not reached; use a large amount for a sole converge-based computation stop
var minWeightRatio = DEFAULT_MIN_WEIGHT_RATIO; // used to compute the minimum allowed weight; default 0.01 means 1% of max weight; handle near-zero weights, and leaves enought space for cell hovering
var prng = DEFAULT_PRNG; // pseudorandom number generator
var initialPosition = DEFAULT_INITIAL_POSITION; // accessor to the initial position; defaults to a random position inside the clipping polygon
var initialWeight = DEFAULT_INITIAL_WEIGHT; // accessor to the initial weight; defaults to the average area of the clipping polygon
var tick = function (polygons, i) {
return true;
}; // hook called at each iteration's end (i = iteration count)
//begin: internals
var weightedVoronoi = d3WeightedVoronoi.weightedVoronoi();
var siteCount,
totalArea,
areaErrorTreshold,
flickeringMitigation = new FlickeringMitigation();
//end: internals
//begin: algorithm conf.
var handleOverweightedVariant = 1; // this option still exists 'cause for further experiments
var handleOverweighted;
//end: algorithm conf.
//begin: utils
function sqr(d) {
return Math.pow(d, 2);
}
function squaredDistance(s0, s1) {
return sqr(s1.x - s0.x) + sqr(s1.y - s0.y);
}
//end: utils
///////////////////////
///////// API /////////
///////////////////////
function _voronoiMap(data) {
//begin: handle algorithm's variants
setHandleOverweighted();
//end: handle algorithm's variants
siteCount = data.length;
(totalArea = Math.abs(d3Polygon.polygonArea(weightedVoronoi.clip()))), (areaErrorTreshold = convergenceRatio * totalArea);
flickeringMitigation.clear().totalArea(totalArea);
var iterationCount = 0,
polygons = initialize(data),
converged = false;
var areaError;
tick(polygons, iterationCount);
while (!(converged || iterationCount >= maxIterationCount)) {
polygons = adapt(polygons, flickeringMitigation.ratio());
iterationCount++;
areaError = computeAreaError(polygons);
flickeringMitigation.add(areaError);
converged = areaError < areaErrorTreshold;
// console.log("error %: "+Math.round(areaError*100*1000/totalArea)/1000);
tick(polygons, iterationCount);
}
return {
polygons: polygons,
iterationCount: iterationCount,
convergenceRatio: areaError / totalArea
};
}
_voronoiMap.weight = function (_) {
if (!arguments.length) {
return weight;
}
weight = _;
return _voronoiMap;
};
_voronoiMap.convergenceRatio = function (_) {
if (!arguments.length) {
return convergenceRatio;
}
convergenceRatio = _;
return _voronoiMap;
};
_voronoiMap.maxIterationCount = function (_) {
if (!arguments.length) {
return maxIterationCount;
}
maxIterationCount = _;
return _voronoiMap;
};
_voronoiMap.minWeightRatio = function (_) {
if (!arguments.length) {
return minWeightRatio;
}
minWeightRatio = _;
return _voronoiMap;
};
_voronoiMap.tick = function (_) {
if (!arguments.length) {
return tick;
}
tick = _;
return _voronoiMap;
};
_voronoiMap.clip = function (_) {
if (!arguments.length) {
return weightedVoronoi.clip();
}
weightedVoronoi.clip(_);
return _voronoiMap;
};
_voronoiMap.extent = function (_) {
if (!arguments.length) {
return weightedVoronoi.extent();
}
weightedVoronoi.extent(_);
return _voronoiMap;
};
_voronoiMap.size = function (_) {
if (!arguments.length) {
return weightedVoronoi.size();
}
weightedVoronoi.size(_);
return _voronoiMap;
};
_voronoiMap.prng = function (_) {
if (!arguments.length) {
return prng;
}
prng = _;
return _voronoiMap;
};
_voronoiMap.initialPosition = function (_) {
if (!arguments.length) {
return initialPosition;
}
initialPosition = _;
return _voronoiMap;
};
_voronoiMap.initialWeight = function (_) {
if (!arguments.length) {
return initialWeight;
}
initialWeight = _;
return _voronoiMap;
};
///////////////////////
/////// Private ///////
///////////////////////
function adapt(polygons, flickeringMitigationRatio) {
var adaptedMapPoints;
adaptPositions(polygons, flickeringMitigationRatio);
adaptedMapPoints = polygons.map(function (p) {
return p.site.originalObject;
});
polygons = weightedVoronoi(adaptedMapPoints);
if (polygons.length < siteCount) {
console.log('at least 1 site has no area, which is not supposed to arise');
debugger;
}
adaptWeights(polygons, flickeringMitigationRatio);
adaptedMapPoints = polygons.map(function (p) {
return p.site.originalObject;
});
polygons = weightedVoronoi(adaptedMapPoints);
if (polygons.length < siteCount) {
console.log('at least 1 site has no area, which is not supposed to arise');
debugger;
}
return polygons;
}
function adaptPositions(polygons, flickeringMitigationRatio) {
var newMapPoints = [],
flickeringInfluence = 0.5;
var flickeringMitigation, d, polygon, mapPoint, centroid, dx, dy;
flickeringMitigation = flickeringInfluence * flickeringMitigationRatio;
d = 1 - flickeringMitigation; // in [0.5, 1]
for (var i = 0; i < siteCount; i++) {
polygon = polygons[i];
mapPoint = polygon.site.originalObject;
centroid = d3Polygon.polygonCentroid(polygon);
//(fixed-x)dx = centroid[0] - mapPoint.x;
dy = centroid[1] - mapPoint.y;
//begin: handle excessive change;
//(fixed-x)dx *= d;
dy *= d;
//end: handle excessive change;
//(fixed-x)mapPoint.x += dx;
mapPoint.y += dy;
newMapPoints.push(mapPoint);
}
handleOverweighted(newMapPoints);
}
function adaptWeights(polygons, flickeringMitigationRatio) {
var newMapPoints = [],
flickeringInfluence = 0.1;
var flickeringMitigation, polygon, mapPoint, currentArea, adaptRatio, adaptedWeight;
flickeringMitigation = flickeringInfluence * flickeringMitigationRatio;
for (var i = 0; i < siteCount; i++) {
polygon = polygons[i];
mapPoint = polygon.site.originalObject;
currentArea = d3Polygon.polygonArea(polygon);
adaptRatio = mapPoint.targetedArea / currentArea;
//begin: handle excessive change;
adaptRatio = Math.max(adaptRatio, 1 - flickeringInfluence + flickeringMitigation); // in [(1-flickeringInfluence), 1]
adaptRatio = Math.min(adaptRatio, 1 + flickeringInfluence - flickeringMitigation); // in [1, (1+flickeringInfluence)]
//end: handle excessive change;
adaptedWeight = mapPoint.weight * adaptRatio;
adaptedWeight = Math.max(adaptedWeight, epsilon);
mapPoint.weight = adaptedWeight;
newMapPoints.push(mapPoint);
}
handleOverweighted(newMapPoints);
}
// heuristics: lower heavy weights
function handleOverweighted0(mapPoints) {
var fixCount = 0;
var fixApplied, tpi, tpj, weightest, lightest, sqrD, adaptedWeight;
do {
fixApplied = false;
for (var i = 0; i < siteCount; i++) {
tpi = mapPoints[i];
for (var j = i + 1; j < siteCount; j++) {
tpj = mapPoints[j];
if (tpi.weight > tpj.weight) {
weightest = tpi;
lightest = tpj;
} else {
weightest = tpj;
lightest = tpi;
}
sqrD = squaredDistance(tpi, tpj);
if (sqrD < weightest.weight - lightest.weight) {
// adaptedWeight = sqrD - epsilon; // as in ArlindNocaj/Voronoi-Treemap-Library
// adaptedWeight = sqrD + lightest.weight - epsilon; // works, but below heuristics performs better (less flickering)
adaptedWeight = sqrD + lightest.weight / 2;
adaptedWeight = Math.max(adaptedWeight, epsilon);
weightest.weight = adaptedWeight;
fixApplied = true;
fixCount++;
break;
}
}
if (fixApplied) {
break;
}
}
} while (fixApplied);
/*
if (fixCount>0) {
console.log("# fix: "+fixCount);
}
*/
}
// heuristics: increase light weights
function handleOverweighted1(mapPoints) {
var fixCount = 0;
var fixApplied, tpi, tpj, weightest, lightest, sqrD, overweight;
do {
fixApplied = false;
for (var i = 0; i < siteCount; i++) {
tpi = mapPoints[i];
for (var j = i + 1; j < siteCount; j++) {
tpj = mapPoints[j];
if (tpi.weight > tpj.weight) {
weightest = tpi;
lightest = tpj;
} else {
weightest = tpj;
lightest = tpi;
}
sqrD = squaredDistance(tpi, tpj);
if (sqrD < weightest.weight - lightest.weight) {
overweight = weightest.weight - lightest.weight - sqrD;
lightest.weight += overweight + epsilon;
fixApplied = true;
fixCount++;
break;
}
}
if (fixApplied) {
break;
}
}
} while (fixApplied);
/*
if (fixCount>0) {
console.log("# fix: "+fixCount);
}
*/
}
function computeAreaError(polygons) {
//convergence based on summation of all sites current areas
var areaErrorSum = 0;
var polygon, mapPoint, currentArea;
for (var i = 0; i < siteCount; i++) {
polygon = polygons[i];
mapPoint = polygon.site.originalObject;
currentArea = d3Polygon.polygonArea(polygon);
areaErrorSum += Math.abs(mapPoint.targetedArea - currentArea);
}
return areaErrorSum;
}
function setHandleOverweighted() {
switch (handleOverweightedVariant) {
case 0:
handleOverweighted = handleOverweighted0;
break;
case 1:
handleOverweighted = handleOverweighted1;
break;
default:
console.log("Variant of 'handleOverweighted' is unknown");
}
}
function initialize(data) {
var maxWeight = data.reduce(function (max, d) {
return Math.max(max, weight(d));
}, -Infinity),
minAllowedWeight = maxWeight * minWeightRatio;
var weights, mapPoints;
//begin: extract weights
weights = data.map(function (d, i, arr) {
return {
index: i,
weight: Math.max(weight(d), minAllowedWeight),
initialPosition: initialPosition(d, i, arr, _voronoiMap),
initialWeight: initialWeight(d, i, arr, _voronoiMap),
originalData: d
};
});
//end: extract weights
// create map-related points
// (with targetedArea, initial position and initialWeight)
mapPoints = createMapPoints(weights);
return weightedVoronoi(mapPoints);
}
function createMapPoints(basePoints) {
var totalWeight = basePoints.reduce(function (acc, bp) {
return (acc += bp.weight);
}, 0);
var initialPosition;
return basePoints.map(function (bp, i, bps) {
initialPosition = bp.initialPosition;
if (!d3Polygon.polygonContains(weightedVoronoi.clip(), initialPosition)) {
initialPosition = RANDOM_INITIAL_POSITION(bp, i, bps, _voronoiMap);
}
return {
index: bp.index,
targetedArea: totalArea * bp.weight / totalWeight,
data: bp,
x: initialPosition[0],
y: initialPosition[1],
weight: bp.initialWeight // ArlindNocaj/Voronoi-Treemap-Library uses an epsilonesque initial weight; using heavier initial weights allows faster weight adjustements, hence faster stabilization
};
});
}
return _voronoiMap;
}
exports.voronoiMap = voronoiMap;
exports.voronoiMapInitialPositionRandom = randomInitialPosition;
exports.voronoiMapInitialPositionPie = pie;
Object.defineProperty(exports, '__esModule', { value: true });
}));
[79,54,74,62,85,55,88,85,51,85,54,84,78,47,83,52,62,84,52,79,51,47,78,69,74,83,55,76,78,79,73,77,66,80,74,52,48,80,59,90,80,58,84,58,73,83,64,53,82,59,75,90,54,80,54,83,71,64,77,81,59,84,48,82,60,92,78,78,65,73,82,56,79,71,62,76,60,78,76,83,75,82,70,65,73,88,76,80,48,86,60,90,50,78,63,72,84,75,51,82,62,88,49,83,81,47,84,52,86,81,75,59,89,79,59,81,50,85,59,87,53,69,77,56,88,81,45,82,55,90,45,83,56,89,46,82,51,86,53,79,81,60,82,77,76,59,80,49,96,53,77,77,65,81,71,70,81,93,53,89,45,86,58,78,66,76,63,88,52,93,49,57,77,68,81,81,73,50,85,74,55,77,83,83,51,78,84,46,83,55,81,57,76,84,77,81,87,77,51,78,60,82,91,53,78,46,77,84,49,83,71,80,49,75,64,76,53,94,55,76,50,82,54,75,78,79,78,78,70,79,70,54,86,50,90,54,54,77,79,64,75,47,86,63,85,82,57,82,67,74,54,83,73,73,88,80,71,83,56,79,78,84,58,83,43,60,75,81,46,90,46,74]
<!DOCTYPE html>
<html>
<head>
<meta charset="utf-8">
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<title>weighted KDE + Voronoï maps</title>
<meta name="description" content="Weighted KDE filled with Voronoï maps in D3.js">
<style>
.hide {
display: none;
}
.standard {
fill: #bbb;
stroke: #bbb;
}
.weighted {
fill: #ddd;
stroke: #ddd;
}
.bin {
stroke: white;
}
.kde {
stroke-width: 1.5;
stroke-linejoin: round;
fill: none;
stroke: black;
}
text.standard {
stroke: none;
}
text.weighted {
stroke: none;
}
.cell {
fill: transparent;
stroke-width: 1;
stroke: grey;
}
.axis--y .domain {
display: none;
}
.dg .property-name {
width: 80% !important;
}
.dg .c {
width: 20% !important;
}
</style>
</head>
<body>
<svg width="960" height="500"></svg>
<script src="https://d3js.org/d3.v5.min.js"></script>
<script src="https://rawcdn.githack.com/Kcnarf/d3-weighted-voronoi/v1.0.0/build/d3-weighted-voronoi.js"></script>
<script src="d3-voronoi-map-fixed-x.js"></script>
<script src="//cdnjs.cloudflare.com/ajax/libs/seedrandom/2.4.3/seedrandom.min.js"></script>
<script src="//cdnjs.cloudflare.com/ajax/libs/dat-gui/0.6.5/dat.gui.min.js"></script>
<script>
var settings = {
standard: {
showBins: true,
showKDE: true
},
weighted: {
showWeightedBins: true,
showWeightedKDE: true,
showVoronoiMaps: true
}
};
var binCount = 40;
var svg = d3.select("svg"),
width = +svg.attr("width"),
height = +svg.attr("height"),
margin = {top: 20, right: 30, bottom: 30, left: 40};
var x = d3.scaleLinear()
.domain([30, 110])
.range([margin.left, width - margin.right]);
var y = d3.scaleLinear()
.range([height - margin.bottom, margin.top]);
var valueAccessor = d => d.value,
weightAccessor = d => d.weight;
var bins, // will store histograms data
density, // will store KDE data
weightedDensity, // will store weighted KDE data
maps; // will store each Voronoi map's data of each bin
d3.json("faithful.json").then(function(data) {
// I know the data is not suitable for such an experimentation :-(
// Data dealing with sales (with the possibility to count the number the sales, or to count the total sales' prices), would be much more better
faithful = data.map((d)=>{
return {
value: d,
weight: Math.round(d/10) // weird weight feature :-(
}
});
//begin: compute histograms data
bins = d3.histogram()
.domain(x.domain())
.value(valueAccessor)
.thresholds(binCount)
(faithful);
bins.forEach(b=>{
b.totalWeight = d3.sum(b, weightAccessor);
})
//end: compute histograms data
y.domain([0, d3.max(bins, b=>b.totalWeight)*1.1]);
//compute KDE's data, using a weight of 1 for each datum
density = kernelDensityEstimator(x.ticks(binCount), scaledEpanechnikov(7, valueAccessor), d=>1)(faithful);
//compute weighted KDE's data, using a particular weight for each datum
weightedDensity = kernelDensityEstimator(x.ticks(binCount), scaledEpanechnikov(7, valueAccessor), weightAccessor)(faithful);
//begin: weighted voronoï map for each right trapeze bin
var seededprng = new Math.seedrandom('my seed'),
voronoiMap = d3.voronoiMap().weight(weightAccessor).prng(seededprng);
maps = [];
weightedDensity.forEach(function (d,i) {
if (i > weightedDensity.length-2) return;
//begin: retrieve weighted KDE of bin's extent for right trapeze bins
var e = weightedDensity[i+1],
x0 = x(d[0]),
y0 = y(d[1]),
x1 = x(e[0]),
y1 = y(e[1]),
baseY = y(0),
data = bins[i];
//begin: retrieve weighted KDE of bin's extent for right trapeze bins
var map;
if (data.length <= 0) return; // no real data
//begin: voronoï map computation of a specific bin, considering weighted KDE
map = voronoiMap
.clip([[x0,baseY], [x0,y0], [x1,y1], [x1,baseY]])
.weight(weightAccessor)
.initialPosition(randomYInitialPosition())
(data).polygons;
//begin: voronoï map computation of a specific bin, considering weighted KDE
maps.push(map);
})
//end: weighted voronoï map for each right trapeze bin
drawBins();
drawVoronoiMaps();
drawKDEs();
drawAxes();
drawDatGui();
});
/******************************/
/* Weighted KDE */
/******************************/
function kernelDensityEstimator(samples, kernel, weightAccessor) {
return function(data) {
return samples.map(function(sample) {
return [sample, d3.sum(data, function(d) {
//return a value depending on kernel function AND weight
return kernel(sample, d) * weightAccessor(d);
})];
});
};
}
function scaledEpanechnikov(bw, valueAccessor) {
// 'bd' states for 'bandwidth'
return function(sample, d) {
var diff = sample - valueAccessor(d);
return epanechnikov(diff/bw)/bw;
};
}
function epanechnikov(diff) {
return Math.abs(diff) <= 1 ? 0.75 * (1 - diff * diff) : 0;
};
/******************************************************/
/* Voronoï map initial positioning with fixed x-coord */
/******************************************************/
function randomYInitialPosition () {
var clippingPolygon,
extent,
minY, maxY,
dy;
function _randomY(d, i, arr, voronoiMap) {
var shouldUpdateInternals = false;
var y;
if (clippingPolygon !== voronoiMap.clip()) {
clippingPolygon = voronoiMap.clip();
extent = voronoiMap.extent();
shouldUpdateInternals = true;
}
if (shouldUpdateInternals) {
updateInternals();
}
y = minY + dy * voronoiMap.prng()();
while (!d3.polygonContains(clippingPolygon, [x(valueAccessor(d)), y])) {
y = minY + dy * voronoiMap.prng()();
}
return [x(valueAccessor(d)), y];
};
///////////////////////
/////// Private ///////
///////////////////////
function updateInternals() {
minY = extent[0][1];
maxY = extent[1][1];
dy = maxY - minY;
};
return _randomY;
};
/******************************/
/* Drawings */
/******************************/
function drawBins() {
//begin: weighted bin
svg.append("g")
.attr("id", "weighted-histogram")
.classed("hide", !settings.weighted.showWeightedBins)
.selectAll("rect")
.data(bins)
.enter().append("rect")
.classed("weighted bin", true)
.attr("x", function(d) { return x(d.x0); })
.attr("y", function(d) { return y(d.totalWeight); })
.attr("width", function(d) { return x(d.x1) - x(d.x0); })
.attr("height", function(d) { return y(0) - y(d.totalWeight); })
.append("title")
.text(d => { return "weighted count: "+d.totalWeight+"\ncount: "+d.length; });;
//end: weighted bin
//begin: standard bin
svg.append("g")
.attr("id", "standard-histogram")
.classed("hide", !settings.standard.showBins)
.selectAll("rect")
.data(bins)
.enter().append("rect")
.classed("standard bin", true)
.attr("x", function(d) { return x(d.x0); })
.attr("y", function(d) { return y(d.length); })
.attr("width", function(d) { return x(d.x1) - x(d.x0); })
.attr("height", function(d) { return y(0) - y(d.length); })
.append("title")
.text(d => { return "count: " + d.length; });
//end: standard bin
}
function drawKDEs () {
//begin: weighted KDE
svg.append("path")
.attr("id", "weighted-kde")
.classed("weighted kde", true)
.classed("hide", !settings.weighted.showWeightedKDE)
.datum(weightedDensity)
.attr("d", d3.line()
.curve(d3.curveBasis)
.x(function(d) { return x(d[0]); })
.y(function(d) { return y(d[1]); }));
//begin: weighted KDE
//begin: standard KDE
svg.append("path")
.attr("id", "standard-kde")
.classed("standard kde", true)
.classed("hide", !settings.standard.showKDE)
.datum(density)
.attr("d", d3.line()
.curve(d3.curveLinear)
.x(function(d) { return x(d[0]); })
.y(function(d) { return y(d[1]); }));
//end: standard KDE
}
function drawVoronoiMaps () {
//begin: draw each cells of each voronoï map
var mapContainer = svg.append("g").attr("id", "voronoi-map-container").classed("hide", !settings.weighted.showVoronoiMaps);
maps.forEach(function(d){
mapContainer.append("g").classed("map", true)
.selectAll("path")
.data(d)
.enter()
.append("path")
.classed("cell", true)
.attr("d", d => { return d3.line().curve(d3.curveLinear)(d) + "z"; })
.append("title")
.text(d => { return "weight: " + weightAccessor(d.site.originalObject.data.originalData); })
})
//end: draw each cells of each voronoï map
}
function drawAxes () {
svg.append("g")
.attr("class", "axis axis--x")
.attr("transform", "translate(0," + (height - margin.bottom) + ")")
.call(d3.axisBottom(x))
.append("text")
.attr("x", width - margin.right)
.attr("y", -6)
.attr("fill", "#000")
.attr("text-anchor", "end")
.attr("font-weight", "bold")
.text("Time between eruptions (min.)");
var yAxis = svg.append("g")
.attr("class", "axis axis--y")
.attr("transform", "translate(" + margin.left + ",0)")
.call(d3.axisLeft(y));
yAxis.append("text")
.classed("standard", true)
.attr("x", -margin.top-80)
.attr("y", 10)
.style("text-anchor", "end")
.style("font-weight", "bold")
.style("transform", "rotate(-90deg)")
.text("Count (i.e. number of eruptions) /");
yAxis.append("text")
.classed("weighted", true)
.attr("x", -margin.top)
.attr("y", 10)
.style("text-anchor", "end")
.style("font-weight", "bold")
.style("transform", "rotate(-90deg)")
.text("Weighted Count");
}
function drawDatGui() {
var controls = new dat.GUI({width: 200});
var standardControler = controls.addFolder("Standard");
standardControler.add(settings.standard, "showBins")
.listen().onChange(function(value) {
d3.select("#standard-histogram").classed("hide", !value);
});
standardControler.add(settings.standard, "showKDE")
.listen().onChange(function(value) {
d3.select("#standard-kde").classed("hide", !value);
});
standardControler.open();
var weightedControler = controls.addFolder("Weighted");
weightedControler.add(settings.weighted, "showWeightedBins")
.listen().onChange(function(value) {
d3.select("#weighted-histogram").classed("hide", !value);
});
weightedControler.add(settings.weighted, "showWeightedKDE")
.listen().onChange(function(value) {
d3.select("#weighted-kde").classed("hide", !value);
});
weightedControler.add(settings.weighted, "showVoronoiMaps")
.listen().onChange(function(value) {
d3.select("#voronoi-map-container").classed("hide", !value);
});
weightedControler.open();
}
</script>
</body>
</html>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment