-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathvenn.html
1902 lines (1628 loc) · 68 KB
/
venn.html
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
<!doctype html>
<html lang="en">
<head>
<meta charset="utf-8">
<title>Simple venn.js example</title>
<style>
body {
font-family: "Helvetica Neue",Helvetica,Arial,sans-serif;
font-size: 14px;
}
</style>
</head>
<body>
<div id="venn"></div>
</body>
<script src="https://d3js.org/d3.v4.min.js"></script>
<script>
(function (global, factory) {
typeof exports === 'object' && typeof module !== 'undefined' ? factory(exports, require('d3-selection'), require('d3-transition')) :
typeof define === 'function' && define.amd ? define(['exports', 'd3-selection', 'd3-transition'], factory) :
(factory((global.venn = {}),global.d3,global.d3));
}(this, (function (exports,d3Selection,d3Transition) { 'use strict';
var SMALL = 1e-10;
/** Returns the intersection area of a bunch of circles (where each circle
is an object having an x,y and radius property) */
function intersectionArea(circles, stats) {
// get all the intersection points of the circles
var intersectionPoints = getIntersectionPoints(circles);
// filter out points that aren't included in all the circles
var innerPoints = intersectionPoints.filter(function (p) {
return containedInCircles(p, circles);
});
var arcArea = 0, polygonArea = 0, arcs = [], i;
// if we have intersection points that are within all the circles,
// then figure out the area contained by them
if (innerPoints.length > 1) {
// sort the points by angle from the center of the polygon, which lets
// us just iterate over points to get the edges
var center = getCenter(innerPoints);
for (i = 0; i < innerPoints.length; ++i ) {
var p = innerPoints[i];
p.angle = Math.atan2(p.x - center.x, p.y - center.y);
}
innerPoints.sort(function(a,b) { return b.angle - a.angle;});
// iterate over all points, get arc between the points
// and update the areas
var p2 = innerPoints[innerPoints.length - 1];
for (i = 0; i < innerPoints.length; ++i) {
var p1 = innerPoints[i];
// polygon area updates easily ...
polygonArea += (p2.x + p1.x) * (p1.y - p2.y);
// updating the arc area is a little more involved
var midPoint = {x : (p1.x + p2.x) / 2,
y : (p1.y + p2.y) / 2},
arc = null;
for (var j = 0; j < p1.parentIndex.length; ++j) {
if (p2.parentIndex.indexOf(p1.parentIndex[j]) > -1) {
// figure out the angle halfway between the two points
// on the current circle
var circle = circles[p1.parentIndex[j]],
a1 = Math.atan2(p1.x - circle.x, p1.y - circle.y),
a2 = Math.atan2(p2.x - circle.x, p2.y - circle.y);
var angleDiff = (a2 - a1);
if (angleDiff < 0) {
angleDiff += 2*Math.PI;
}
// and use that angle to figure out the width of the
// arc
var a = a2 - angleDiff/2,
width = distance(midPoint, {
x : circle.x + circle.radius * Math.sin(a),
y : circle.y + circle.radius * Math.cos(a)
});
// clamp the width to the largest is can actually be
// (sometimes slightly overflows because of FP errors)
if (width > circle.radius * 2) {
width = circle.radius * 2;
}
// pick the circle whose arc has the smallest width
if ((arc === null) || (arc.width > width)) {
arc = { circle : circle,
width : width,
p1 : p1,
p2 : p2};
}
}
}
if (arc !== null) {
arcs.push(arc);
arcArea += circleArea(arc.circle.radius, arc.width);
p2 = p1;
}
}
} else {
// no intersection points, is either disjoint - or is completely
// overlapped. figure out which by examining the smallest circle
var smallest = circles[0];
for (i = 1; i < circles.length; ++i) {
if (circles[i].radius < smallest.radius) {
smallest = circles[i];
}
}
// make sure the smallest circle is completely contained in all
// the other circles
var disjoint = false;
for (i = 0; i < circles.length; ++i) {
if (distance(circles[i], smallest) > Math.abs(smallest.radius - circles[i].radius)) {
disjoint = true;
break;
}
}
if (disjoint) {
arcArea = polygonArea = 0;
} else {
arcArea = smallest.radius * smallest.radius * Math.PI;
arcs.push({circle : smallest,
p1: { x: smallest.x, y : smallest.y + smallest.radius},
p2: { x: smallest.x - SMALL, y : smallest.y + smallest.radius},
width : smallest.radius * 2 });
}
}
polygonArea /= 2;
if (stats) {
stats.area = arcArea + polygonArea;
stats.arcArea = arcArea;
stats.polygonArea = polygonArea;
stats.arcs = arcs;
stats.innerPoints = innerPoints;
stats.intersectionPoints = intersectionPoints;
}
return arcArea + polygonArea;
}
/** returns whether a point is contained by all of a list of circles */
function containedInCircles(point, circles) {
for (var i = 0; i < circles.length; ++i) {
if (distance(point, circles[i]) > circles[i].radius + SMALL) {
return false;
}
}
return true;
}
/** Gets all intersection points between a bunch of circles */
function getIntersectionPoints(circles) {
var ret = [];
for (var i = 0; i < circles.length; ++i) {
for (var j = i + 1; j < circles.length; ++j) {
var intersect = circleCircleIntersection(circles[i],
circles[j]);
for (var k = 0; k < intersect.length; ++k) {
var p = intersect[k];
p.parentIndex = [i,j];
ret.push(p);
}
}
}
return ret;
}
/** Circular segment area calculation. See http://mathworld.wolfram.com/CircularSegment.html */
function circleArea(r, width) {
return r * r * Math.acos(1 - width/r) - (r - width) * Math.sqrt(width * (2 * r - width));
}
/** euclidean distance between two points */
function distance(p1, p2) {
return Math.sqrt((p1.x - p2.x) * (p1.x - p2.x) +
(p1.y - p2.y) * (p1.y - p2.y));
}
/** Returns the overlap area of two circles of radius r1 and r2 - that
have their centers separated by distance d. Simpler faster
circle intersection for only two circles */
function circleOverlap(r1, r2, d) {
// no overlap
if (d >= r1 + r2) {
return 0;
}
// completely overlapped
if (d <= Math.abs(r1 - r2)) {
return Math.PI * Math.min(r1, r2) * Math.min(r1, r2);
}
var w1 = r1 - (d * d - r2 * r2 + r1 * r1) / (2 * d),
w2 = r2 - (d * d - r1 * r1 + r2 * r2) / (2 * d);
return circleArea(r1, w1) + circleArea(r2, w2);
}
/** Given two circles (containing a x/y/radius attributes),
returns the intersecting points if possible.
note: doesn't handle cases where there are infinitely many
intersection points (circles are equivalent):, or only one intersection point*/
function circleCircleIntersection(p1, p2) {
var d = distance(p1, p2),
r1 = p1.radius,
r2 = p2.radius;
// if to far away, or self contained - can't be done
if ((d >= (r1 + r2)) || (d <= Math.abs(r1 - r2))) {
return [];
}
var a = (r1 * r1 - r2 * r2 + d * d) / (2 * d),
h = Math.sqrt(r1 * r1 - a * a),
x0 = p1.x + a * (p2.x - p1.x) / d,
y0 = p1.y + a * (p2.y - p1.y) / d,
rx = -(p2.y - p1.y) * (h / d),
ry = -(p2.x - p1.x) * (h / d);
return [{x: x0 + rx, y : y0 - ry },
{x: x0 - rx, y : y0 + ry }];
}
/** Returns the center of a bunch of points */
function getCenter(points) {
var center = {x: 0, y: 0};
for (var i =0; i < points.length; ++i ) {
center.x += points[i].x;
center.y += points[i].y;
}
center.x /= points.length;
center.y /= points.length;
return center;
}
/** finds the zeros of a function, given two starting points (which must
* have opposite signs */
function bisect(f, a, b, parameters) {
parameters = parameters || {};
var maxIterations = parameters.maxIterations || 100,
tolerance = parameters.tolerance || 1e-10,
fA = f(a),
fB = f(b),
delta = b - a;
if (fA * fB > 0) {
throw "Initial bisect points must have opposite signs";
}
if (fA === 0) return a;
if (fB === 0) return b;
for (var i = 0; i < maxIterations; ++i) {
delta /= 2;
var mid = a + delta,
fMid = f(mid);
if (fMid * fA >= 0) {
a = mid;
}
if ((Math.abs(delta) < tolerance) || (fMid === 0)) {
return mid;
}
}
return a + delta;
}
// need some basic operations on vectors, rather than adding a dependency,
// just define here
function zeros(x) { var r = new Array(x); for (var i = 0; i < x; ++i) { r[i] = 0; } return r; }
function zerosM(x,y) { return zeros(x).map(function() { return zeros(y); }); }
function dot(a, b) {
var ret = 0;
for (var i = 0; i < a.length; ++i) {
ret += a[i] * b[i];
}
return ret;
}
function norm2(a) {
return Math.sqrt(dot(a, a));
}
function scale(ret, value, c) {
for (var i = 0; i < value.length; ++i) {
ret[i] = value[i] * c;
}
}
function weightedSum(ret, w1, v1, w2, v2) {
for (var j = 0; j < ret.length; ++j) {
ret[j] = w1 * v1[j] + w2 * v2[j];
}
}
/** minimizes a function using the downhill simplex method */
function nelderMead(f, x0, parameters) {
parameters = parameters || {};
var maxIterations = parameters.maxIterations || x0.length * 200,
nonZeroDelta = parameters.nonZeroDelta || 1.05,
zeroDelta = parameters.zeroDelta || 0.001,
minErrorDelta = parameters.minErrorDelta || 1e-6,
minTolerance = parameters.minErrorDelta || 1e-5,
rho = (parameters.rho !== undefined) ? parameters.rho : 1,
chi = (parameters.chi !== undefined) ? parameters.chi : 2,
psi = (parameters.psi !== undefined) ? parameters.psi : -0.5,
sigma = (parameters.sigma !== undefined) ? parameters.sigma : 0.5,
maxDiff;
// initialize simplex.
var N = x0.length,
simplex = new Array(N + 1);
simplex[0] = x0;
simplex[0].fx = f(x0);
simplex[0].id = 0;
for (var i = 0; i < N; ++i) {
var point = x0.slice();
point[i] = point[i] ? point[i] * nonZeroDelta : zeroDelta;
simplex[i+1] = point;
simplex[i+1].fx = f(point);
simplex[i+1].id = i+1;
}
function updateSimplex(value) {
for (var i = 0; i < value.length; i++) {
simplex[N][i] = value[i];
}
simplex[N].fx = value.fx;
}
var sortOrder = function(a, b) { return a.fx - b.fx; };
var centroid = x0.slice(),
reflected = x0.slice(),
contracted = x0.slice(),
expanded = x0.slice();
for (var iteration = 0; iteration < maxIterations; ++iteration) {
simplex.sort(sortOrder);
if (parameters.history) {
// copy the simplex (since later iterations will mutate) and
// sort it to have a consistent order between iterations
var sortedSimplex = simplex.map(function (x) {
var state = x.slice();
state.fx = x.fx;
state.id = x.id;
return state;
});
sortedSimplex.sort(function(a,b) { return a.id - b.id; });
parameters.history.push({x: simplex[0].slice(),
fx: simplex[0].fx,
simplex: sortedSimplex});
}
maxDiff = 0;
for (i = 0; i < N; ++i) {
maxDiff = Math.max(maxDiff, Math.abs(simplex[0][i] - simplex[1][i]));
}
if ((Math.abs(simplex[0].fx - simplex[N].fx) < minErrorDelta) &&
(maxDiff < minTolerance)) {
break;
}
// compute the centroid of all but the worst point in the simplex
for (i = 0; i < N; ++i) {
centroid[i] = 0;
for (var j = 0; j < N; ++j) {
centroid[i] += simplex[j][i];
}
centroid[i] /= N;
}
// reflect the worst point past the centroid and compute loss at reflected
// point
var worst = simplex[N];
weightedSum(reflected, 1+rho, centroid, -rho, worst);
reflected.fx = f(reflected);
// if the reflected point is the best seen, then possibly expand
if (reflected.fx < simplex[0].fx) {
weightedSum(expanded, 1+chi, centroid, -chi, worst);
expanded.fx = f(expanded);
if (expanded.fx < reflected.fx) {
updateSimplex(expanded);
} else {
updateSimplex(reflected);
}
}
// if the reflected point is worse than the second worst, we need to
// contract
else if (reflected.fx >= simplex[N-1].fx) {
var shouldReduce = false;
if (reflected.fx > worst.fx) {
// do an inside contraction
weightedSum(contracted, 1+psi, centroid, -psi, worst);
contracted.fx = f(contracted);
if (contracted.fx < worst.fx) {
updateSimplex(contracted);
} else {
shouldReduce = true;
}
} else {
// do an outside contraction
weightedSum(contracted, 1-psi * rho, centroid, psi*rho, worst);
contracted.fx = f(contracted);
if (contracted.fx < reflected.fx) {
updateSimplex(contracted);
} else {
shouldReduce = true;
}
}
if (shouldReduce) {
// if we don't contract here, we're done
if (sigma >= 1) break;
// do a reduction
for (i = 1; i < simplex.length; ++i) {
weightedSum(simplex[i], 1 - sigma, simplex[0], sigma, simplex[i]);
simplex[i].fx = f(simplex[i]);
}
}
} else {
updateSimplex(reflected);
}
}
simplex.sort(sortOrder);
return {fx : simplex[0].fx,
x : simplex[0]};
}
/// searches along line 'pk' for a point that satifies the wolfe conditions
/// See 'Numerical Optimization' by Nocedal and Wright p59-60
/// f : objective function
/// pk : search direction
/// current: object containing current gradient/loss
/// next: output: contains next gradient/loss
/// returns a: step size taken
function wolfeLineSearch(f, pk, current, next, a, c1, c2) {
var phi0 = current.fx, phiPrime0 = dot(current.fxprime, pk),
phi = phi0, phi_old = phi0,
phiPrime = phiPrime0,
a0 = 0;
a = a || 1;
c1 = c1 || 1e-6;
c2 = c2 || 0.1;
function zoom(a_lo, a_high, phi_lo) {
for (var iteration = 0; iteration < 16; ++iteration) {
a = (a_lo + a_high)/2;
weightedSum(next.x, 1.0, current.x, a, pk);
phi = next.fx = f(next.x, next.fxprime);
phiPrime = dot(next.fxprime, pk);
if ((phi > (phi0 + c1 * a * phiPrime0)) ||
(phi >= phi_lo)) {
a_high = a;
} else {
if (Math.abs(phiPrime) <= -c2 * phiPrime0) {
return a;
}
if (phiPrime * (a_high - a_lo) >=0) {
a_high = a_lo;
}
a_lo = a;
phi_lo = phi;
}
}
return 0;
}
for (var iteration = 0; iteration < 10; ++iteration) {
weightedSum(next.x, 1.0, current.x, a, pk);
phi = next.fx = f(next.x, next.fxprime);
phiPrime = dot(next.fxprime, pk);
if ((phi > (phi0 + c1 * a * phiPrime0)) ||
(iteration && (phi >= phi_old))) {
return zoom(a0, a, phi_old);
}
if (Math.abs(phiPrime) <= -c2 * phiPrime0) {
return a;
}
if (phiPrime >= 0 ) {
return zoom(a, a0, phi);
}
phi_old = phi;
a0 = a;
a *= 2;
}
return a;
}
function conjugateGradient(f, initial, params) {
// allocate all memory up front here, keep out of the loop for perfomance
// reasons
var current = {x: initial.slice(), fx: 0, fxprime: initial.slice()},
next = {x: initial.slice(), fx: 0, fxprime: initial.slice()},
yk = initial.slice(),
pk, temp,
a = 1,
maxIterations;
params = params || {};
maxIterations = params.maxIterations || initial.length * 20;
current.fx = f(current.x, current.fxprime);
pk = current.fxprime.slice();
scale(pk, current.fxprime,-1);
for (var i = 0; i < maxIterations; ++i) {
a = wolfeLineSearch(f, pk, current, next, a);
// todo: history in wrong spot?
if (params.history) {
params.history.push({x: current.x.slice(),
fx: current.fx,
fxprime: current.fxprime.slice(),
alpha: a});
}
if (!a) {
// faiiled to find point that satifies wolfe conditions.
// reset direction for next iteration
scale(pk, current.fxprime, -1);
} else {
// update direction using Polak–Ribiere CG method
weightedSum(yk, 1, next.fxprime, -1, current.fxprime);
var delta_k = dot(current.fxprime, current.fxprime),
beta_k = Math.max(0, dot(yk, next.fxprime) / delta_k);
weightedSum(pk, beta_k, pk, -1, next.fxprime);
temp = current;
current = next;
next = temp;
}
if (norm2(current.fxprime) <= 1e-5) {
break;
}
}
if (params.history) {
params.history.push({x: current.x.slice(),
fx: current.fx,
fxprime: current.fxprime.slice(),
alpha: a});
}
return current;
}
/** given a list of set objects, and their corresponding overlaps.
updates the (x, y, radius) attribute on each set such that their positions
roughly correspond to the desired overlaps */
function venn(areas, parameters) {
parameters = parameters || {};
parameters.maxIterations = parameters.maxIterations || 500;
var initialLayout = parameters.initialLayout || bestInitialLayout;
var loss = parameters.lossFunction || lossFunction;
// add in missing pairwise areas as having 0 size
areas = addMissingAreas(areas);
// initial layout is done greedily
var circles = initialLayout(areas, parameters);
// transform x/y coordinates to a vector to optimize
var initial = [], setids = [], setid;
for (setid in circles) {
if (circles.hasOwnProperty(setid)) {
initial.push(circles[setid].x);
initial.push(circles[setid].y);
setids.push(setid);
}
}
var solution = nelderMead(
function(values) {
var current = {};
for (var i = 0; i < setids.length; ++i) {
var setid = setids[i];
current[setid] = {x: values[2 * i],
y: values[2 * i + 1],
radius : circles[setid].radius,
// size : circles[setid].size
};
}
return loss(current, areas);
},
initial,
parameters);
// transform solution vector back to x/y points
var positions = solution.x;
for (var i = 0; i < setids.length; ++i) {
setid = setids[i];
circles[setid].x = positions[2 * i];
circles[setid].y = positions[2 * i + 1];
}
return circles;
}
var SMALL$1 = 1e-10;
/** Returns the distance necessary for two circles of radius r1 + r2 to
have the overlap area 'overlap' */
function distanceFromIntersectArea(r1, r2, overlap) {
// handle complete overlapped circles
if (Math.min(r1, r2) * Math.min(r1,r2) * Math.PI <= overlap + SMALL$1) {
return Math.abs(r1 - r2);
}
return bisect(function(distance$$1) {
return circleOverlap(r1, r2, distance$$1) - overlap;
}, 0, r1 + r2);
}
/** Missing pair-wise intersection area data can cause problems:
treating as an unknown means that sets will be laid out overlapping,
which isn't what people expect. To reflect that we want disjoint sets
here, set the overlap to 0 for all missing pairwise set intersections */
function addMissingAreas(areas) {
areas = areas.slice();
// two circle intersections that aren't defined
var ids = [], pairs = {}, i, j, a, b;
for (i = 0; i < areas.length; ++i) {
var area = areas[i];
if (area.sets.length == 1) {
ids.push(area.sets[0]);
} else if (area.sets.length == 2) {
a = area.sets[0];
b = area.sets[1];
pairs[[a, b]] = true;
pairs[[b, a]] = true;
}
}
ids.sort(function(a, b) { return a > b; });
for (i = 0; i < ids.length; ++i) {
a = ids[i];
for (j = i + 1; j < ids.length; ++j) {
b = ids[j];
if (!([a, b] in pairs)) {
areas.push({'sets': [a, b],
'size': 0});
}
}
}
return areas;
}
/// Returns two matrices, one of the euclidean distances between the sets
/// and the other indicating if there are subset or disjoint set relationships
function getDistanceMatrices(areas, sets, setids) {
// initialize an empty distance matrix between all the points
var distances = zerosM(sets.length, sets.length),
constraints = zerosM(sets.length, sets.length);
// compute required distances between all the sets such that
// the areas match
areas.filter(function(x) { return x.sets.length == 2; })
.map(function(current) {
var left = setids[current.sets[0]],
right = setids[current.sets[1]],
r1 = Math.sqrt(sets[left].size / Math.PI),
r2 = Math.sqrt(sets[right].size / Math.PI),
distance$$1 = distanceFromIntersectArea(r1, r2, current.size);
distances[left][right] = distances[right][left] = distance$$1;
// also update constraints to indicate if its a subset or disjoint
// relationship
var c = 0;
if (current.size + 1e-10 >= Math.min(sets[left].size,
sets[right].size)) {
c = 1;
} else if (current.size <= 1e-10) {
c = -1;
}
constraints[left][right] = constraints[right][left] = c;
});
return {distances: distances, constraints: constraints};
}
/// computes the gradient and loss simulatenously for our constrained MDS optimizer
function constrainedMDSGradient(x, fxprime, distances, constraints) {
var loss = 0, i;
for (i = 0; i < fxprime.length; ++i) {
fxprime[i] = 0;
}
for (i = 0; i < distances.length; ++i) {
var xi = x[2 * i], yi = x[2 * i + 1];
for (var j = i + 1; j < distances.length; ++j) {
var xj = x[2 * j], yj = x[2 * j + 1],
dij = distances[i][j],
constraint = constraints[i][j];
var squaredDistance = (xj - xi) * (xj - xi) + (yj - yi) * (yj - yi),
distance$$1 = Math.sqrt(squaredDistance),
delta = squaredDistance - dij * dij;
if (((constraint > 0) && (distance$$1 <= dij)) ||
((constraint < 0) && (distance$$1 >= dij))) {
continue;
}
loss += 2 * delta * delta;
fxprime[2*i] += 4 * delta * (xi - xj);
fxprime[2*i + 1] += 4 * delta * (yi - yj);
fxprime[2*j] += 4 * delta * (xj - xi);
fxprime[2*j + 1] += 4 * delta * (yj - yi);
}
}
return loss;
}
/// takes the best working variant of either constrained MDS or greedy
function bestInitialLayout(areas, params) {
var initial = greedyLayout(areas, params);
var loss = params.lossFunction || lossFunction;
// greedylayout is sufficient for all 2/3 circle cases. try out
// constrained MDS for higher order problems, take its output
// if it outperforms. (greedy is aesthetically better on 2/3 circles
// since it axis aligns)
if (areas.length >= 8) {
var constrained = constrainedMDSLayout(areas, params),
constrainedLoss = loss(constrained, areas),
greedyLoss = loss(initial, areas);
if (constrainedLoss + 1e-8 < greedyLoss) {
initial = constrained;
}
}
return initial;
}
/// use the constrained MDS variant to generate an initial layout
function constrainedMDSLayout(areas, params) {
params = params || {};
var restarts = params.restarts || 10;
// bidirectionally map sets to a rowid (so we can create a matrix)
var sets = [], setids = {}, i;
for (i = 0; i < areas.length; ++i ) {
var area = areas[i];
if (area.sets.length == 1) {
setids[area.sets[0]] = sets.length;
sets.push(area);
}
}
var matrices = getDistanceMatrices(areas, sets, setids),
distances = matrices.distances,
constraints = matrices.constraints;
// keep distances bounded, things get messed up otherwise.
// TODO: proper preconditioner?
var norm = norm2(distances.map(norm2))/(distances.length);
distances = distances.map(function (row) {
return row.map(function (value) { return value / norm; });});
var obj = function(x, fxprime) {
return constrainedMDSGradient(x, fxprime, distances, constraints);
};
var best, current;
for (i = 0; i < restarts; ++i) {
var initial = zeros(distances.length*2).map(Math.random);
current = conjugateGradient(obj, initial, params);
if (!best || (current.fx < best.fx)) {
best = current;
}
}
var positions = best.x;
// translate rows back to (x,y,radius) coordinates
var circles = {};
for (i = 0; i < sets.length; ++i) {
var set = sets[i];
circles[set.sets[0]] = {
x: positions[2*i] * norm,
y: positions[2*i + 1] * norm,
radius: Math.sqrt(set.size / Math.PI)
};
}
if (params.history) {
for (i = 0; i < params.history.length; ++i) {
scale(params.history[i].x, norm);
}
}
return circles;
}
/** Lays out a Venn diagram greedily, going from most overlapped sets to
least overlapped, attempting to position each new set such that the
overlapping areas to already positioned sets are basically right */
function greedyLayout(areas, params) {
var loss = params && params.lossFunction ? params.lossFunction : lossFunction;
// define a circle for each set
var circles = {}, setOverlaps = {}, set;
for (var i = 0; i < areas.length; ++i) {
var area = areas[i];
if (area.sets.length == 1) {
set = area.sets[0];
circles[set] = {x: 1e10, y: 1e10,
rowid: circles.length,
size: area.size,
radius: Math.sqrt(area.size / Math.PI)};
setOverlaps[set] = [];
}
}
areas = areas.filter(function(a) { return a.sets.length == 2; });
// map each set to a list of all the other sets that overlap it
for (i = 0; i < areas.length; ++i) {
var current = areas[i];
var weight = current.hasOwnProperty('weight') ? current.weight : 1.0;
var left = current.sets[0], right = current.sets[1];
// completely overlapped circles shouldn't be positioned early here
if (current.size + SMALL$1 >= Math.min(circles[left].size,
circles[right].size)) {
weight = 0;
}
setOverlaps[left].push ({set:right, size:current.size, weight:weight});
setOverlaps[right].push({set:left, size:current.size, weight:weight});
}
// get list of most overlapped sets
var mostOverlapped = [];
for (set in setOverlaps) {
if (setOverlaps.hasOwnProperty(set)) {
var size = 0;
for (i = 0; i < setOverlaps[set].length; ++i) {
size += setOverlaps[set][i].size * setOverlaps[set][i].weight;
}
mostOverlapped.push({set: set, size:size});
}
}
// sort by size desc
function sortOrder(a,b) {
return b.size - a.size;
}
mostOverlapped.sort(sortOrder);
// keep track of what sets have been laid out
var positioned = {};
function isPositioned(element) {
return element.set in positioned;
}
// adds a point to the output
function positionSet(point, index) {
circles[index].x = point.x;
circles[index].y = point.y;
positioned[index] = true;
}
// add most overlapped set at (0,0)
positionSet({x: 0, y: 0}, mostOverlapped[0].set);
// get distances between all points. TODO, necessary?
// answer: probably not
// var distances = venn.getDistanceMatrices(circles, areas).distances;
for (i = 1; i < mostOverlapped.length; ++i) {
var setIndex = mostOverlapped[i].set,
overlap = setOverlaps[setIndex].filter(isPositioned);
set = circles[setIndex];
overlap.sort(sortOrder);
if (overlap.length === 0) {
// this shouldn't happen anymore with addMissingAreas
throw "ERROR: missing pairwise overlap information";
}
var points = [];
for (var j = 0; j < overlap.length; ++j) {
// get appropriate distance from most overlapped already added set
var p1 = circles[overlap[j].set],
d1 = distanceFromIntersectArea(set.radius, p1.radius,
overlap[j].size);
// sample positions at 90 degrees for maximum aesthetics
points.push({x : p1.x + d1, y : p1.y});
points.push({x : p1.x - d1, y : p1.y});
points.push({y : p1.y + d1, x : p1.x});
points.push({y : p1.y - d1, x : p1.x});
// if we have at least 2 overlaps, then figure out where the
// set should be positioned analytically and try those too
for (var k = j + 1; k < overlap.length; ++k) {
var p2 = circles[overlap[k].set],
d2 = distanceFromIntersectArea(set.radius, p2.radius,
overlap[k].size);
var extraPoints = circleCircleIntersection(
{ x: p1.x, y: p1.y, radius: d1},
{ x: p2.x, y: p2.y, radius: d2});
for (var l = 0; l < extraPoints.length; ++l) {
points.push(extraPoints[l]);
}
}
}
// we have some candidate positions for the set, examine loss
// at each position to figure out where to put it at
var bestLoss = 1e50, bestPoint = points[0];
for (j = 0; j < points.length; ++j) {
circles[setIndex].x = points[j].x;
circles[setIndex].y = points[j].y;
var localLoss = loss(circles, areas);
if (localLoss < bestLoss) {
bestLoss = localLoss;
bestPoint = points[j];
}
}
positionSet(bestPoint, setIndex);
}
return circles;
}
/** Given a bunch of sets, and the desired overlaps between these sets - computes
the distance from the actual overlaps to the desired overlaps. Note that
this method ignores overlaps of more than 2 circles */
function lossFunction(sets, overlaps) {
var output = 0;
function getCircles(indices) {
return indices.map(function(i) { return sets[i]; });
}
for (var i = 0; i < overlaps.length; ++i) {
var area = overlaps[i], overlap;
if (area.sets.length == 1) {
continue;
} else if (area.sets.length == 2) {
var left = sets[area.sets[0]],
right = sets[area.sets[1]];
overlap = circleOverlap(left.radius, right.radius,
distance(left, right));
} else {
overlap = intersectionArea(getCircles(area.sets));
}
var weight = area.hasOwnProperty('weight') ? area.weight : 1.0;
output += weight * (overlap - area.size) * (overlap - area.size);
}
return output;
}