% This file contains NESL code to find a partition of a D dimensional "irregular" mesh using the geometric method of Miller, Teng and Vavasis. It has been implemented in part by Guy Blelloch, Gary Miller and Shanghua Teng. % % ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; % % SOME STANDARD MATRIX OPERATIONS % % ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; % % Transposes a matrix % function transpose(a) = let rows = {#a: a}[0]; columns = #a; data = permute(flatten(a),distindex(columns, rows)); in partition(data,dist(columns,rows)) $ function dotp(a,b) = sum({a*b:a;b}) $ function norm(a) = sqrt(dotp(a,a)) $ function vsum(p2,p1,a) = {p2 + a*p1:p1;p2} $ function mvmult(A,x) = {dotp(a,x) : a in transpose(A)} $ % Finds the null of a matrix (i.e. a weighted sum of the rows that adds to all zeros) % function null(A) = let l = #a; (p1,A) = head_rest(A); n1 = dotp(p1,p1); in if (n1 < 1.0e-20) then cons(1.0,dist(0.0,l-1)) else if (l == 1) then [0.0] else let b = {negate(dotp(p,p1)/n1) : p in A}; c = null({vsum(a,p1,b): a ; b}); in cons(dotp(b,c), c) $ % ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; % % SOME UTILITIES % % ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; % % Returns a random unit vector in DIMS dimensions % function random_unit_vector(dims) = let point = {rand(v)-1.0: v in dist(2.0,dims)}; norm = norm(point); in if (norm > 1.0) then random_unit_vector(dims) else {p/norm : p in point} $ % makes a stereo projection of p onto a sphere in one higher dimension % function stereo_up(p) = let norm = 2.0/(1.0 + dotp(p,p)); in snoc({p*norm: p}, 1.0 - norm) $ % does a stereo projection from a sphere onto a plane in one lower dimension % function stereo_down(p) = let (p,z) = rest_tail(p); in {p/(1 - z) : p} $ % ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; % % CONFORMAL MAPPINGS ON A SPHERE -- Written by Gary Miller % % ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; % % DESCRIPTION: stereo-dilation takes: v --- a unit vector normal to some plane thur the origin g --- a distance squared and returns a unit normal for the image plane which goes through (0,...,0,sqrt g) plus the distance of the image plane from origin. TYPE: (V.FLOAT, FLOAT) <- (V.FLOAT, FLOAT) % function stereo_dilation(v,g) = let (v,z) = rest_tail(v); norm = sqrt(max(1.0-z^2*g,0.0)); v = {v/norm : v}; z = z*sqrt(max(1.0-g,0.0))/norm; u = g*z; in (snoc(v,z),u) $ % DESCRIPTION: get-sphere takes unit normal v in R^d+1 and an intercept g returns center of sphere in R^d and radi after appling stereo-down TYPE: (v.float, float) <- (v.float, float) % function get_sphere(normal,g) = let (v,z) = rest_tail(normal); denom = g - z; v = {v/denom : v}; radi = negate(sqrt(max(1.0-g^2,0.0))/denom); in (v,radi) $ % DESCRIPTION: reflect-normal reflects n such that c will go to d+1 axis TYPE: v.float <- (v.float,v.float) % function reflect_normal(c,n) = let len = norm(c); (rest,tail) = rest_tail(c); u = snoc(rest,tail-len); norm = dotp(u,u); alpha = if (norm < 1.0e-20) then 0.0 else negate(2.0*dotp(n,u)/norm); in vsum(n,u,alpha) $ % DESCRIPTION: normal-sphere take a normal unit vector n and a center point c both in R^d+1 and returns a center in R^d and a radi. TYPE: (v.float,float) <- (v.float, v.float) % function normal_sphere(n,c) = let c2 = dotp(c,c); (n1,g) = stereo_dilation(n,c2); n2 = reflect_normal(c,n1); in get_sphere(n2,g) $ % ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; % % FINDING A CENTERPOINT % % ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; % % DESCRIPTION: Takes the centerpoint of a set of d+2 points in d dimensions. EXAMPLE: center([[2.0, 1.0], [7.0, 2.0], [3.0, 5.0], [8.0, 2.0]]) -> [7.0, 2.0] TYPE: v.float <- v.v.float % function center(A) = let null = null({snoc(a,1.0): a in A}); znull = {max(x,0.0) : x in null}; sum = max(sum(znull),1e-20); in mvmult(A,{x / sum : x in znull}) $ function centerpoint(points) = if (#points == 1) then points[0] else let dims = 2 + #points[0]; groups = partition(points,dist(dims, #points/dims)); in centerpoint({center(g): g in groups}) $ % ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; % % SELECTING A CIRCLE % % ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; % function distance(a,b) = sum({(a-b)^2: a; b}) $ function getfromedges(edges,v) = {(e1,e2) : e1 in v->{e1: (e1,e2) in edges} ; e2 in v->{e2: (e1,e2) in edges} } $ % This routine picks a circle with a center at point C, that encompasses half the POINTS. It returns the radius of this circle and the number of EDGES cut by the circle. % function getcut(points,edges,c) = let distances = {distance(p,c) : p in points}; median = median(distances); flags = {d < median : d in distances}; cuts = count({e1 xor e2 : (e1,e2) in getfromedges(edges,flags)}); in (sqrt(median),cuts) $ % Tries several great circles. For each it projects the circle into the lower dimension and grows or shrinks the circle so that it encompases half the points. Out of the trials, the one that cuts the minimum number of edges is selected. START: the current trial number END: the last trial number (end - start is the number of trials left) POINTS: locations of the points EDGES: the edges (each is a pair of integers) CENTERPOINT: the centerpoint % function try_several(start,end,points,edges,centerpoint) = if (start == end) then (1000000,[0.0],0.0) else let randvect = random_unit_vector(#centerpoint); (center,radius) = normal_sphere(randvect,centerpoint); (radius,cuts) = getcut(points,edges,center); (rcuts, c, r) = try_several(start+1,end,points,edges,centerpoint); in if (cuts < rcuts) then (cuts,center,radius) else (rcuts, c, r) $ % Used to decide what the sample size for the centerpoint routine will be, given the total number of points N. DIMS is the number of dimensions of the points. The sample size will be an integer power of 2 more than the number of dimensions. % function sample_size(n,dims) = let csize = dims + 2; pow = trunc(log(float(n),float(csize))); in if (pow > 2) then csize^(pow-1) else csize^pow $ % Finds a circular partition of the POINTS into two equal sized sets that minimizes the number of edges cuts. It uses the geometric separator method. OTRIALS: specifies the number of centerpoints to try TRIALS: specifies the number of great circles to try for each centerpoint The routine returns 1) the number of edges cut by the best cut 2) the center of the circle for that cut 3) the radius of the circle for that cut % function find_circle(points,edges,trials,otrials) = if (otrials == 0) then (1000000,[0.0],0.0) else let sample_size = sample_size(#points,#points[0]+1); sample = points->{rand(i) : i in dist(#points,sample_size)}; centerpoint = centerpoint({stereo_up(p): p in sample}); (cut1,c1,r1) = try_several(0,trials,points,edges,centerpoint); (cut2,c2,r2) = find_circle(points,edges,trials,otrials-1); in if (cut1 < cut2) then (cut1,c1,r1) else (cut2,c2,r2) $ % ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; % % SPLITTING A GRAPH % % ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; % function new_numbers(flags) = let down = enumerate({not(flags): flags}); up = enumerate(flags); in {select(flags,up,down): flags; up; down} $ function split_graph(side,points,edges) = let (sizes,sindex) = split_index(side); new_n = new_numbers(side); new_edges = getfromedges(edges,new_n); newside = getfromedges(edges,side); eleft = pack(zip(new_edges, {not(s1 or s2) : (s1,s2) in newside})); eright = pack(zip(new_edges, {s1 and s2 : (s1,s2) in newside})); in (split(points,side),vpair(eleft,eright),sindex) $ % ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; % % THE FINAL CODE % % ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; % % A recursive routine to find a geometric separator tree. POINTS: the set of points for which it is finding a separator EDGES: the edges between the points. The separator tries to minimize the number of edges that are cut. TRIALS, OTRIALS: see the description of find_circle. DEPTH: the depth of the tree. The number of partitions will be 2^depth. COUNT: Used in the recursion to specify the partition number. The function returns an integer for each point, which specifies the partition the point is in. Partition numbers will be in the range 0 <= pnum < 2^depth. Within 1, all partitions will be of equal size. % function geometric_separator(points,edges,trials,otrials,depth,count) = if (depth == 0) then dist(count,#points) else let (cuts,center,radius) = find_circle(points,edges,trials,otrials); str = "Level: " ++ @depth ++ " Partition: " ++ @count ++ " Edges: " ++ @(#edges) ++ " Cuts: " ++ @cuts ++ [newline]; line = print_string(str); side = {norm({p-center:p;center}) < radius: p in points}; (snodes,sedges,sindex) = split_graph(side,points,edges); result = {geometric_separator(n,e,trials,otrials,depth-1,c) : n in snodes; e in sedges; c in vpair(count*2, 1+count*2)}; in flatten(result)->sindex $ % ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; % % TEST CODE % % ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; % function readpoints(filename) = let points = read_object_from_file([(0.0,0.0)],filename); in {vpair(p1,p2) : (p1,p2) in points} $ function readedges(filename) = read_object_from_file([(0,0)],filename) $ data_directory = "/afs/cs.cmu.edu/project/scandal/nesl/.scratch/data/" $ function read_graph_from_data_directory(graphname) = (readpoints(data_directory ++ graphname ++ ".nodes"), readedges(data_directory ++ graphname ++ ".edges")) $ display = "blelloch.pc.cs.cmu.edu:0"; function drawgraph(points,edges,display) = let w = make_window(default_window,bounding_box(points),"geosep",display); p = draw_points(points,w); s = draw_segments(getfromedges(edges,points),1,w); in close_window(w) $ function drawgraphps(points,edges,display) = let w = make_ps_window(((72,72),(504,504)),bounding_box(points),display); s = draw_segments(getfromedges(edges,points),1,w); in close_window(w) $ function drawgraphfile(filename,display) = let (g,e) = read_graph_from_data_directory(filename); nodes = {(p[0],p[1]): p in g}; in drawgraph(nodes,e,display) $ function internal_edges(part,edges) = {e in edges; (e1,e2) in getfromedges(edges,part) | e1==e2} $ function cross_edges(part,edges) = {e in edges; (e1,e2) in getfromedges(edges,part) | e1/=e2} $ function geometric_separator_test(depth,trials,outtrials,graphname,display) = let (points,edges) = read_graph_from_data_directory(graphname); part = geometric_separator(points,edges,trials,outtrials,depth,0); iedges = internal_edges(part,edges); graph = drawgraph({(a[0],a[1]): a in points}, iedges, display); in #iedges $ function geometric_separator_test_noplot(depth,trials,outtrials,graphname) = let (points,edges) = read_graph_from_data_directory(graphname); part = geometric_separator(points,edges,trials,outtrials,depth,0); iedges = internal_edges(part,edges); in #iedges $ % Results -- AIRFOIL old (geometric_separator 6 10 2) -- 10263 new (geometric_separator 6 1 1) -- 9950 (geometric_separator 6 2 1) -- 10192 (geometric_separator 6 3 1) -- 10287 (geometric_separator 6 4 1) -- 10402 (geometric_separator 6 5 1) -- 10392 (geometric_separator 6 10 1) -- 10484 (geometric_separator 6 20 1) -- 10503 (geometric_separator 6 10 2) -- 10504 (geometric_separator 6 50 1) -- 10542 (geometric_separator 6 10 5) -- 10584 (1705)- 18 percent improve testpart (quad_separator 6) -- 10208 (2081) eigenvalue method top level best cut 102 final cut 127 RESULTS -- CRACK (30380 total) (geometric_separator 6 3 2) -- 27157 (geometric_separator 6 5 3) -- 27249 (3131) - 2.5 percent improve (quad_separator 6) -- 27177 (3203) % function drawgraph_n(points,edges,display) = let w = make_window(((200,0),(800,800)),bounding_box(points),"geosep",display); p = draw_points(points,w); e = getfromedges(edges,points); s = draw_segments(e,1,w); in close_window(w) $ function drawgraph_n2(pts,edges,display) = let w = make_window(((200,0),(800,800)),bounding_box(pts),"geosep",display); p = draw_points(pts,w); e = getfromedges(edges,pts); s = draw_segments(e,1,w); in w $ function subselect_points(points,(xmin,xmax),(ymin,ymax)) = let p_flags = {(x >= xmin) and (x <= xmax) and (y >= ymin) and (y <= ymax): (x,y) in points}; new_points = {points; p_flags| p_flags} in (p_flags,new_points) $ function subselect_edges(edges,p_flags) = let data = zip(p_flags, enumerate(p_flags)); newedges = {(i1,i2): ((f1,i1),(f2,i2)) in getfromedges(edges,data) | f1 and f2} in newedges $ function subselect_graph(points,edges,corners) = let (p_flags,new_points) = subselect_points(points,corners); new_edges = subselect_edges(edges,p_flags); in (new_points,new_edges) $ function get_other_corner(w,x1,y1) = let (x2,y2,b2,s) = get_mouse_info(w); xbounds = (min(x1,x2),max(x1,x2)); ybounds = (min(y1,y2),max(y1,y2)); in (xbounds,ybounds) $ function zoom_track(points,edges,w,display) = let (x1,y1,button,shift) = get_mouse_info(w) in if (button /= 1) then kill_window(w) else let corners = get_other_corner(w,x1,y1); (npoints,nedges) = subselect_graph(points,edges,corners); nw = drawgraph_n2(npoints,nedges,display); rec = zoom_track(npoints,nedges,nw,display) in zoom_track(points,edges,w,display) $ function read_graph_new(filename,dimensions) = let p = read_float_seq_from_file(filename++".nodes"); points = group_by(p,dimensions); foo = print_debug("Points: ", #points); e = read_int_seq_from_file(filename++".edges"); edges = zip(even_elts(e),odd_elts(e)); foo = print_debug("Edges: ", #edges); in (points,edges) $ function drawgraph_test(filename,dimensions,display) = let (points,edges) = read_graph_new(filename,dimensions); points = {p[0],p[1]: p in points} in drawgraph_n(points,edges,display) $ function geometric_sep(dims,depth,trials,outtrials,filename,display) = let foo = print_line("Reading Data"); (points,edges) = read_graph_new(filename,dims); foo = print_line("Starting Computation"); part = geometric_separator(points,edges,trials,outtrials, depth,0); iedges = if (dims == 0) then internal_edges(part,edges) else cross_edges(part,edges); pts = if dims == 2 then {(a[0],a[1]): a in points} else {(.86*a[0]-.86*a[1],.5*a[0]+.5*a[1]+a[2]):a in points}; w = drawgraph_n2(pts, iedges, display); x = zoom_track(pts,iedges,w,display) in #iedges $