From b975fb74a903f9dd176f8e3bc92b8ad780151a07 Mon Sep 17 00:00:00 2001 From: Ben Frederickson Date: Thu, 24 Nov 2016 16:15:44 -0800 Subject: [PATCH] split fmin code to a separate library --- index.js | 1 - package.json | 3 +- src/diagram.js | 6 +- src/fmin.js | 318 --------------------- src/layout.js | 12 +- tests/test_utils.js | 2 +- tests/unittest.js | 41 --- venn.js | 661 +++++++++++++++++++++++--------------------- 8 files changed, 352 insertions(+), 692 deletions(-) delete mode 100644 src/fmin.js diff --git a/index.js b/index.js index 80aa1c1..18f778e 100644 --- a/index.js +++ b/index.js @@ -1,4 +1,3 @@ -export {fmin, minimizeConjugateGradient, bisect} from "./src/fmin"; export {intersectionArea, circleCircleIntersection, circleOverlap, circleArea, distance, circleIntegral} from "./src/circleintersection"; export {venn, greedyLayout, scaleSolution, normalizeSolution, bestInitialLayout, diff --git a/package.json b/package.json index 76f5c30..57000f8 100644 --- a/package.json +++ b/package.json @@ -1,6 +1,6 @@ { "name": "venn.js", - "version": "0.2.13", + "version": "0.2.14", "author": "Ben Frederickson (http:/www.benfrederickson.com)", "url": "https://github.com/benfred/venn.js/issues", "devDependencies": { @@ -44,6 +44,7 @@ "postpublish": "zip -j build/venn.zip -- LICENSE README.md build/venn.js build/venn.min.js" }, "dependencies": { + "fmin": "0.0.2", "d3-selection": "^1.0.2", "d3-transition": "^1.0.1" } diff --git a/src/diagram.js b/src/diagram.js index eb8d5ff..1afe040 100644 --- a/src/diagram.js +++ b/src/diagram.js @@ -3,7 +3,7 @@ import {transition} from "d3-transition"; import {venn, normalizeSolution, scaleSolution} from "./layout"; import {intersectionArea, distance, getCenter} from "./circleintersection"; -import {fmin} from "./fmin"; +import {nelderMead} from "../node_modules/fmin/index.js"; /*global console:true*/ @@ -363,10 +363,10 @@ export function computeTextCentre(interior, exterior) { } // maximize the margin numerically - var solution = fmin( + var solution = nelderMead( function(p) { return -1 * circleMargin({x: p[0], y: p[1]}, interior, exterior); }, [initial.x, initial.y], - {maxIterations:500, minErrorDelta:1e-10}).solution; + {maxIterations:500, minErrorDelta:1e-10}).x; var ret = {x: solution[0], y: solution[1]}; // check solution, fallback as needed (happens if fully overlapped diff --git a/src/fmin.js b/src/fmin.js deleted file mode 100644 index 5364272..0000000 --- a/src/fmin.js +++ /dev/null @@ -1,318 +0,0 @@ -/** finds the zeros of a function, given two starting points (which must - * have opposite signs */ -export 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 -export function zeros(x) { var r = new Array(x); for (var i = 0; i < x; ++i) { r[i] = 0; } return r; } -export function zerosM(x,y) { return zeros(x).map(function() { return zeros(y); }); } - -export function dot(a, b) { - var ret = 0; - for (var i = 0; i < a.length; ++i) { - ret += a[i] * b[i]; - } - return ret; -} - -export function norm2(a) { - return Math.sqrt(dot(a, a)); -} - -export function multiplyBy(a, c) { - for (var i = 0; i < a.length; ++i) { - a[i] *= c; - } -} - -export 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 */ -export function fmin(f, x0, parameters) { - parameters = parameters || {}; - - var maxIterations = parameters.maxIterations || x0.length * 200, - nonZeroDelta = parameters.nonZeroDelta || 1.1, - zeroDelta = parameters.zeroDelta || 0.001, - minErrorDelta = parameters.minErrorDelta || 1e-6, - minTolerance = parameters.minErrorDelta || 1e-5, - rho = parameters.rho || 1, - chi = parameters.chi || 2, - psi = parameters.psi || -0.5, - sigma = parameters.sigma || 0.5, - callback = parameters.callback, - maxDiff, - temp; - - // initialize simplex. - var N = x0.length, - simplex = new Array(N + 1); - simplex[0] = x0; - simplex[0].fx = f(x0); - 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); - } - - 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 (callback) { - callback(simplex); - } - - 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) { - temp = simplex[N]; - simplex[N] = expanded; - expanded = temp; - } else { - temp = simplex[N]; - simplex[N] = reflected; - reflected = temp; - } - } - - // 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) { - temp = simplex[N]; - simplex[N] = contracted; - contracted = temp; - } 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) { - temp = simplex[N]; - simplex[N] = contracted; - contracted = temp; - } else { - shouldReduce = true; - } - } - - if (shouldReduce) { - // do reduction. doesn't actually happen that often - for (i = 1; i < simplex.length; ++i) { - weightedSum(simplex[i], 1 - sigma, simplex[0], sigma, simplex[i]); - simplex[i].fx = f(simplex[i]); - } - } - } else { - temp = simplex[N]; - simplex[N] = reflected; - reflected = temp; - } - - } - - simplex.sort(sortOrder); - return {f : simplex[0].fx, - solution : simplex[0]}; -} - -export function minimizeConjugateGradient(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 * 5; - - current.fx = f(current.x, current.fxprime); - pk = current.fxprime.slice(); - multiplyBy(pk, -1); - - for (var i = 0; i < maxIterations; ++i) { - if (params.history) { - params.history.push({x: current.x.slice(), - fx: current.fx, - fxprime: current.fxprime.slice()}); - } - - a = wolfeLineSearch(f, pk, current, next, a); - if (!a) { - // faiiled to find point that satifies wolfe conditions. - // reset direction for next iteration - for (var j = 0; j < pk.length; ++j) { - pk[j] = -1 * current.fxprime[j]; - } - } 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()}); - } - - return current; -} - -var c1 = 1e-6, c2 = 0.1; - -/// searches along line 'pk' for a point that satifies the wolfe conditions -/// See 'Numerical Optimization' by Nocedal and Wright p59-60 -export function wolfeLineSearch(f, pk, current, next, a) { - var phi0 = current.fx, phiPrime0 = dot(current.fxprime, pk), - phi = phi0, phi_old = phi0, - phiPrime = phiPrime0, - a0 = 0; - - a = a || 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 0; -} - - diff --git a/src/layout.js b/src/layout.js index 4d241c2..f780025 100644 --- a/src/layout.js +++ b/src/layout.js @@ -1,5 +1,5 @@ -import {fmin, bisect, minimizeConjugateGradient, zeros, zerosM, norm2, - multiplyBy} from './fmin'; +import {nelderMead, bisect, conjugateGradient, zeros, zerosM, norm2, + scale} from '../node_modules/fmin/index.js'; import {intersectionArea, circleOverlap, circleCircleIntersection, distance} from './circleintersection'; /** given a list of set objects, and their corresponding overlaps. @@ -28,7 +28,7 @@ export function venn(areas, parameters) { // optimize initial layout from our loss function var totalFunctionCalls = 0; - var solution = fmin( + var solution = nelderMead( function(values) { totalFunctionCalls += 1; var current = {}; @@ -46,7 +46,7 @@ export function venn(areas, parameters) { parameters); // transform solution vector back to x/y points - var positions = solution.solution; + var positions = solution.x; for (var i = 0; i < setids.length; ++i) { setid = setids[i]; circles[setid].x = positions[2 * i]; @@ -228,7 +228,7 @@ export function constrainedMDSLayout(areas, params) { for (i = 0; i < restarts; ++i) { var initial = zeros(distances.length*2).map(Math.random); - current = minimizeConjugateGradient(obj, initial, params); + current = conjugateGradient(obj, initial, params); if (!best || (current.fx < best.fx)) { best = current; } @@ -248,7 +248,7 @@ export function constrainedMDSLayout(areas, params) { if (params.history) { for (i = 0; i < params.history.length; ++i) { - multiplyBy(params.history[i].x, norm); + scale(params.history[i].x, norm); } } return circles; diff --git a/tests/test_utils.js b/tests/test_utils.js index b529231..b52ae60 100644 --- a/tests/test_utils.js +++ b/tests/test_utils.js @@ -86,7 +86,7 @@ function hilightErrors(div, areas, current, duration) { for (var i = 0; i < failedAreas.length; ++i) { var area = failedAreas[i]; - div.select('[data-venn-sets="' + area.sets[0] + "_" + area.sets[1] + '"] path') + div.selectAll('[data-venn-sets="' + area.sets[0] + "_" + area.sets[1] + '"] path') .style("stroke", "red") .style("stroke-opacity", 1) .style("stroke-width", 2) diff --git a/tests/unittest.js b/tests/unittest.js index e698903..c38ebe1 100644 --- a/tests/unittest.js +++ b/tests/unittest.js @@ -62,47 +62,6 @@ tape("greedyLayout", function(test) { test.end(); }); -tape("fmin", function(test) { - // minimize simple 1 diminesial quadratic - var loss = function(values) { return (values[0] - 10) * (values[0] - 10); }; - var solution = venn.fmin(loss, [0], {minErrorDelta:1e-10}).solution; - nearlyEqual(test, solution[0], 10, 1e-10); - test.end(); -}); - -tape("fmin_himmelblau", function(test) { - // due to a bug, this used to not converge to the minimum - var x = 4.9515014216303825, y = 0.07301421370357275; - function himmelblau(x, y) { - return (x * x + y - 11) * ( x * x + y - 11) + (x + y * y - 7) * (x + y * y - 7); - } - var solution = venn.fmin(function (x) { return himmelblau(x[0], x[1]);}, [x, y]); - nearlyEqual(test, solution.f, 0); - test.end(); -}); - -tape("fmin_banana", function(test) { - var x = 1.6084564160555601, y = -1.5980748860165477; - function banana(x, y) { - return (1 - x) * (1 - x) + 100 * (y - x * x) * ( y - x * x); - } - - var solution = venn.fmin(function (x) { return banana(x[0], x[1]);}, [x, y]); - nearlyEqual(test, solution.f, 0); - test.end(); -}); - -tape("minimizeConjugateGradient", function(test) { - // minimize simple 1 diminesial quadratic - var loss = function(x, xprime) { - xprime[0] = 2 * (x[0] - 10); - return (x[0] - 10) * (x[0] - 10); - }; - var solution = venn.minimizeConjugateGradient(loss, [0]).x; - nearlyEqual(test, solution[0], 10, 1e-10); - test.end(); -}); - tape("circleIntegral", function(test) { nearlyEqual(test, venn.circleIntegral(10, 0), 0, SMALL, "empty circle test"); diff --git a/venn.js b/venn.js index 5c22815..06e1274 100644 --- a/venn.js +++ b/venn.js @@ -1,9 +1,232 @@ (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); + (factory((global.venn = 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) + }); + + // 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; + } + + function circleIntegral(r, x) { + var y = Math.sqrt(r * r - x * x); + return x * y + r * r * Math.atan2(x, y); + } + + /** Returns the area of a circle of radius r - up to width */ + function circleArea(r, width) { + return circleIntegral(r, width - r) - circleIntegral(r, -r); + } + + /** 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) { @@ -54,9 +277,9 @@ return Math.sqrt(dot(a, a)); } - function multiplyBy(a, c) { - for (var i = 0; i < a.length; ++i) { - a[i] *= c; + function scale(ret, value, c) { + for (var i = 0; i < value.length; ++i) { + ret[i] = value[i] * c; } } @@ -67,32 +290,39 @@ } /** minimizes a function using the downhill simplex method */ - function fmin(f, x0, parameters) { + function nelderMead(f, x0, parameters) { parameters = parameters || {}; var maxIterations = parameters.maxIterations || x0.length * 200, - nonZeroDelta = parameters.nonZeroDelta || 1.1, + nonZeroDelta = parameters.nonZeroDelta || 1.05, zeroDelta = parameters.zeroDelta || 0.001, minErrorDelta = parameters.minErrorDelta || 1e-6, minTolerance = parameters.minErrorDelta || 1e-5, - rho = parameters.rho || 1, - chi = parameters.chi || 2, - psi = parameters.psi || -0.5, - sigma = parameters.sigma || 0.5, - callback = parameters.callback, - maxDiff, - temp; + 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; }; @@ -104,8 +334,21 @@ for (var iteration = 0; iteration < maxIterations; ++iteration) { simplex.sort(sortOrder); - if (callback) { - callback(simplex); + + 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; @@ -134,17 +377,13 @@ reflected.fx = f(reflected); // if the reflected point is the best seen, then possibly expand - if (reflected.fx <= simplex[0].fx) { + if (reflected.fx < simplex[0].fx) { weightedSum(expanded, 1+chi, centroid, -chi, worst); expanded.fx = f(expanded); if (expanded.fx < reflected.fx) { - temp = simplex[N]; - simplex[N] = expanded; - expanded = temp; + updateSimplex(expanded); } else { - temp = simplex[N]; - simplex[N] = reflected; - reflected = temp; + updateSimplex(reflected); } } @@ -158,9 +397,7 @@ weightedSum(contracted, 1+psi, centroid, -psi, worst); contracted.fx = f(contracted); if (contracted.fx < worst.fx) { - temp = simplex[N]; - simplex[N] = contracted; - contracted = temp; + updateSimplex(contracted); } else { shouldReduce = true; } @@ -168,105 +405,49 @@ // do an outside contraction weightedSum(contracted, 1-psi * rho, centroid, psi*rho, worst); contracted.fx = f(contracted); - if (contracted.fx <= reflected.fx) { - temp = simplex[N]; - simplex[N] = contracted; - contracted = temp; + if (contracted.fx < reflected.fx) { + updateSimplex(contracted); } else { shouldReduce = true; - } - } - - if (shouldReduce) { - // do reduction. doesn't actually happen that often - for (i = 1; i < simplex.length; ++i) { - weightedSum(simplex[i], 1 - sigma, simplex[0], sigma, simplex[i]); - simplex[i].fx = f(simplex[i]); - } - } - } else { - temp = simplex[N]; - simplex[N] = reflected; - reflected = temp; - } - - } - - simplex.sort(sortOrder); - return {f : simplex[0].fx, - solution : simplex[0]}; - } - - function minimizeConjugateGradient(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 * 5; - - current.fx = f(current.x, current.fxprime); - pk = current.fxprime.slice(); - multiplyBy(pk, -1); - - for (var i = 0; i < maxIterations; ++i) { - if (params.history) { - params.history.push({x: current.x.slice(), - fx: current.fx, - fxprime: current.fxprime.slice()}); - } - - a = wolfeLineSearch(f, pk, current, next, a); - if (!a) { - // faiiled to find point that satifies wolfe conditions. - // reset direction for next iteration - for (var j = 0; j < pk.length; ++j) { - pk[j] = -1 * current.fxprime[j]; - } - } 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 (shouldReduce) { + // if we don't contract here, we're done + if (sigma >= 1) break; - if (params.history) { - params.history.push({x: current.x.slice(), - fx: current.fx, - fxprime: current.fxprime.slice()}); + // 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); + } } - return current; + simplex.sort(sortOrder); + return {fx : simplex[0].fx, + x : simplex[0]}; } - var c1 = 1e-6; - var c2 = 0.1; /// searches along line 'pk' for a point that satifies the wolfe conditions /// See 'Numerical Optimization' by Nocedal and Wright p59-60 - function wolfeLineSearch(f, pk, current, next, a) { + /// 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) { @@ -318,230 +499,69 @@ a *= 2; } - return 0; + return a; } - 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; - } + 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; - // 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) - }); + params = params || {}; + maxIterations = params.maxIterations || initial.length * 20; - // pick the circle whose arc has the smallest width - if ((arc === null) || (arc.width > width)) { - arc = { circle : circle, - width : width, - p1 : p1, - p2 : p2}; - } - } - } + current.fx = f(current.x, current.fxprime); + pk = current.fxprime.slice(); + scale(pk, current.fxprime,-1); - 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]; - } - } + for (var i = 0; i < maxIterations; ++i) { + a = wolfeLineSearch(f, pk, current, next, a); - // 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; - } + // 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 (disjoint) { - arcArea = polygonArea = 0; + if (!a) { + // faiiled to find point that satifies wolfe conditions. + // reset direction for next iteration + scale(pk, current.fxprime, -1); } 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 }); - } - } + // update direction using Polak–Ribiere CG method + weightedSum(yk, 1, next.fxprime, -1, current.fxprime); - polygonArea /= 2; - if (stats) { - stats.area = arcArea + polygonArea; - stats.arcArea = arcArea; - stats.polygonArea = polygonArea; - stats.arcs = arcs; - stats.innerPoints = innerPoints; - stats.intersectionPoints = intersectionPoints; - } + var delta_k = dot(current.fxprime, current.fxprime), + beta_k = Math.max(0, dot(yk, next.fxprime) / delta_k); - return arcArea + polygonArea; - } + weightedSum(pk, beta_k, pk, -1, next.fxprime); - /** 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; + temp = current; + current = next; + next = temp; } - } - 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); - } + if (norm2(current.fxprime) <= 1e-5) { + break; } } - return ret; - } - - function circleIntegral(r, x) { - var y = Math.sqrt(r * r - x * x); - return x * y + r * r * Math.atan2(x, y); - } - - /** Returns the area of a circle of radius r - up to width */ - function circleArea(r, width) { - return circleIntegral(r, width - r) - circleIntegral(r, -r); - } - - /** 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 []; + if (params.history) { + params.history.push({x: current.x.slice(), + fx: current.fx, + fxprime: current.fxprime.slice(), + alpha: a}); } - 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; + return current; } /** given a list of set objects, and their corresponding overlaps. @@ -570,7 +590,7 @@ // optimize initial layout from our loss function var totalFunctionCalls = 0; - var solution = fmin( + var solution = nelderMead( function(values) { totalFunctionCalls += 1; var current = {}; @@ -588,7 +608,7 @@ parameters); // transform solution vector back to x/y points - var positions = solution.solution; + var positions = solution.x; for (var i = 0; i < setids.length; ++i) { setid = setids[i]; circles[setid].x = positions[2 * i]; @@ -770,7 +790,7 @@ for (i = 0; i < restarts; ++i) { var initial = zeros(distances.length*2).map(Math.random); - current = minimizeConjugateGradient(obj, initial, params); + current = conjugateGradient(obj, initial, params); if (!best || (current.fx < best.fx)) { best = current; } @@ -790,7 +810,7 @@ if (params.history) { for (i = 0; i < params.history.length; ++i) { - multiplyBy(params.history[i].x, norm); + scale(params.history[i].x, norm); } } return circles; @@ -1556,10 +1576,10 @@ } // maximize the margin numerically - var solution = fmin( + var solution = nelderMead( function(p) { return -1 * circleMargin({x: p[0], y: p[1]}, interior, exterior); }, [initial.x, initial.y], - {maxIterations:500, minErrorDelta:1e-10}).solution; + {maxIterations:500, minErrorDelta:1e-10}).x; var ret = {x: solution[0], y: solution[1]}; // check solution, fallback as needed (happens if fully overlapped @@ -1760,9 +1780,6 @@ } } - exports.fmin = fmin; - exports.minimizeConjugateGradient = minimizeConjugateGradient; - exports.bisect = bisect; exports.intersectionArea = intersectionArea; exports.circleCircleIntersection = circleCircleIntersection; exports.circleOverlap = circleOverlap; @@ -1786,4 +1803,6 @@ exports.circleFromPath = circleFromPath; exports.intersectionAreaPath = intersectionAreaPath; + Object.defineProperty(exports, '__esModule', { value: true }); + })); \ No newline at end of file