diff --git a/docs/spec/builtins/file-io.md b/docs/spec/builtins/file-io.md index e25807d6..ae777c6f 100644 --- a/docs/spec/builtins/file-io.md +++ b/docs/spec/builtins/file-io.md @@ -296,6 +296,13 @@ Writes an `Image` to a raster image file, or a `Graphics` object to an image fil plot). Coordinates may be exact (`1/2`, `Pi/4`, `Sqrt[2]`): they are converted the same way the on-screen renderer converts them. Text uses the PDF base-14 Helvetica, so no font is embedded. This is the recommended format for print and for the book. +- In the PDF, `Text[s, pos, {ox, oy}]` aligns as Mathematica does (`{-1, 0}` puts the left + end of `s` at `pos`, `{0, 0}` centres it, using the Helvetica advance widths), and + `Text[Style[s, n | FontSize -> n | colour, ...], ...]` sets that string's size and + colour. `Arrowheads[s]` fixes the arrowhead length at `s` times the plot width, and an + arrow's shaft stops inside its head rather than poking past the point. + `AspectRatio -> Automatic` maps x and y with one scale (the page height follows the + data unless `ImageSize -> {w, h}` fixes both, in which case the picture is centred). - **PNG** and **JPEG** render through the graphics backend into an offscreen buffer, so the file is pixel-identical to the on-screen plot (the same axes, ticks, labels and text). They therefore need graphics support compiled in (`USE_GRAPHICS`) **and** a usable GUI diff --git a/docs/spec/builtins/graphs.md b/docs/spec/builtins/graphs.md index 9f774980..c4bb93d3 100644 --- a/docs/spec/builtins/graphs.md +++ b/docs/spec/builtins/graphs.md @@ -1133,34 +1133,101 @@ Out[5]= TopologicalSort[Graph[<3 vertices, 3 edges>]] ## GraphPlot -- `GraphPlot[g]`: a `Graphics[...]` object drawing `g`. - -**Features**: -- `Protected`. Vertices are laid out on a circle, each drawn as a `Disk`; edges - are `Line`s. The specification calls for one `Text` label per vertex as well, - but the current binary emits no `Text` primitives (see the example below). -- Renders through the standard graphics path (a window when `USE_GRAPHICS=1`, - the text placeholder otherwise). -- MVP limitations: directed edges are drawn as plain lines (no arrowheads yet); - a force-directed layout is a future hook. Mathematica's `GraphPlot` uses a - spring-electrical layout. -- Unevaluated on a non-graph. - -```mathematica -In[1]:= Head[GraphPlot[CycleGraph[8]]] +- `GraphPlot[g]`: a `Graphics[...]` object drawing the graph `g`. +- `GraphPlot[{u -> v, ...}]`: draws the graph of a list of rules. +- `GraphPlot[g, opts]`: with the options below; any other option (`ImageSize`, + `PlotLabel`, `Background`, ...) is passed through to `Graphics`. + +**Features**: +- `Protected`. Implemented in `src/graph/graphplot.c` over the layout engine + `src/graph/glayout.c`. **Deterministic**: no random numbers anywhere, so the + same graph always gives an identical `Graphics` expression. +- **Default layout** (`GraphLayout -> Automatic`): a forest with a branching + vertex is drawn as tidy layered trees (hanging from the tree centre, or from + the source of an arborescence); an all-directed acyclic graph as a layered + drawing; everything else (paths and cycles included) by **stress + majorization** (SMACOF on BFS graph distances, weights `d^-2`), started from + Pivot MDS. Components of up to 60 vertices also try circle starts and Tutte + (barycentric) starts from shortest cycles, and keep, among the drawings within + 10% of the least stress, the one with the fewest edge crossings; a + crossing-free drawing may cost up to 2.5x the stress when it removes at least + four crossings (so the dodecahedron comes out as its Schlegel diagram while + the cube stays the textbook Necker cube). Each drawing is rotated to a + canonical orientation (principal axis horizontal; snapped to the axes when the + edges are four-fold, then straightened so grids come out exactly as grids; + vertex 1 on top when there is no preferred axis). Every connected component is + laid out on its own and the components are shelf-packed, largest first, with + isolated vertices gathered into a square block. +- **Cost**: full stress majorization up to 1000 vertices per component, Pivot + MDS (50 pivots, `O(k (n + m))`) above that. A 500-vertex random graph takes + about 0.15 s, a 25x20 grid 0.03 s. +- `GraphLayout` values: `"StressEmbedding"`, `"SpringElectricalEmbedding"` (Hu's + spring-electrical model, exact repulsion up to 1000 vertices, grid + cut-off above), `"CircularEmbedding"` (`VertexList` order, vertex 1 on top), + `"LayeredEmbedding"` / `"LayeredDigraphEmbedding"` (longest-path layers for a + DAG, BFS layers from the centre otherwise, dummy vertices on long edges, + barycentre crossing reduction, isotonic-regression x placement; wide shallow + drawings get taller layer spacing), `"BipartiteEmbedding"` (the two parts in + two columns, barycentre-ordered; falls back to stress for a non-bipartite + graph), `"GridEmbedding"` (`VertexList` order on a square grid). +- `VertexCoordinates -> {{x1, y1}, ...}` (one pair per vertex, `VertexList` + order) fixes the drawing; `VertexCoordinates -> {v -> {x, y}, ...}` fixes the + given vertices and lays out the rest. +- `VertexLabels -> None` (default) | `"Name"` | `Automatic` | `True` | `All` + labels each vertex with its name; `{v -> lbl, ...}` labels only those. A label + sits beside its vertex on the side with the widest angular gap between the + incident edges (upper right when free), in 10 pt Helvetica, and the frame + grows so no label is clipped. +- **Directed edges** are `Arrow`s with an `Arrowheads` directive sized to the + vertex disks; each arrow starts outside its source disk and its tip stops just + short of the target disk. A mutual pair `u -> v`, `v -> u` is drawn as two + arrows offset to either side. +- `GraphHighlight -> {v, ..., e, ...}`: highlighted vertices and edges are drawn + red (`RGBColor[1, 0, 0]`), vertices 15% larger and edges 2.5x thicker, on top + of the others. Edges may be written `u <-> v`, `UndirectedEdge[u, v]`, + `u -> v` or `DirectedEdge[u, v]`. +- `VertexStyle -> style` or `{v -> style, ...}`; `EdgeStyle -> style` or + `{e -> style, ...}` (a colour, or a list of directives), e.g. a vertex + colouring from `FindVertexColoring`. +- `EdgeLabels -> "EdgeWeight"` writes each weight at its edge midpoint (offset + towards the inside of the drawing); `EdgeLabels -> {e -> lbl, ...}` labels + chosen edges. `VertexSize -> d` sets the disk diameter to `d` edge lengths. +- **Look**: vertices are disks in `RGBColor[0.368417, 0.506779, 0.709798]` + (Mathematica's `ColorData[97]` blue) with a thin darker rim, edges 1.1 pt in + grey-blue `RGBColor[0.571589, 0.586483, 0.699215]`, the disk radius scaled to + the median edge length and the extent of the drawing. The result carries + `PlotRange`, `AspectRatio -> Automatic`, `Axes -> False` and an explicit + `ImageSize -> {w, h}` (300 pt on the long side, growing gently with the vertex + count), so `Export["g.pdf", GraphPlot[g]]` gives a tight, equal-aspect picture. +- Unevaluated on a non-graph, or when an argument after the graph is not a rule. + +```mathematica +In[1]:= Head[GraphPlot[PetersenGraph[]]] Out[1]= Graphics In[2]:= Count[GraphPlot[CompleteGraph[6]], _Line, Infinity] Out[2]= 15 -In[3]:= Count[GraphPlot[CycleGraph[5]], _Disk, Infinity] -Out[3]= 5 +In[3]:= Count[GraphPlot[Graph[{1 -> 2, 2 -> 3, 3 -> 1}]], _Arrow, Infinity] +Out[3]= 3 -In[4]:= Count[GraphPlot[CycleGraph[5]], _Text, Infinity] -Out[4]= 0 +In[4]:= Cases[GraphPlot[PathGraph[{1, 2, 3}], VertexCoordinates -> {{0, 0}, {1, 0}, {2, 1}}], Disk[p_, _] :> p, Infinity] +Out[4]= {{0.0, 0.0}, {1.0, 0.0}, {2.0, 1.0}} + +In[5]:= Count[GraphPlot[CycleGraph[5], VertexLabels -> "Name"], _Text, Infinity] +Out[5]= 5 + +In[6]:= MemberQ[GraphPlot[CycleGraph[3], GraphHighlight -> {1}], RGBColor[1., 0., 0.], Infinity] +Out[6]= True + +In[7]:= {Axes, AspectRatio} /. Rest[List @@ GraphPlot[CycleGraph[4]]] +Out[7]= {False, Automatic} + +In[8]:= Length[Union[Cases[GraphPlot[GridGraph[{4, 4}]], Disk[{x_, _}, _] :> Round[x, 0.001], Infinity]]] +Out[8]= 4 -In[5]:= GraphPlot[5] -Out[5]= GraphPlot[5] +In[9]:= GraphPlot[5] +Out[9]= GraphPlot[5] ``` ## FindVertexColoring diff --git a/docs/spec/builtins/hypergraphs.md b/docs/spec/builtins/hypergraphs.md index 86e7b6b3..899628c5 100644 --- a/docs/spec/builtins/hypergraphs.md +++ b/docs/spec/builtins/hypergraphs.md @@ -354,6 +354,53 @@ In[3]:= InputForm[HypergraphStarExpansion[Hypergraph[{Hyperedge[1], 2},{{Hypered Out[3]= HypergraphStarExpansion[Hypergraph[{Hyperedge[1], 2}, {{Hyperedge[1], 2}}]] ``` +## HypergraphPlot + +- `HypergraphPlot[h]`: a `Graphics[...]` object drawing the hypergraph `h`. +- `HypergraphPlot[{e1, e2, ...}]`: draws the hypergraph of a list of hyperedges. +- `HypergraphPlot[h, opts]`: with `VertexLabels`, `VertexCoordinates`, + `VertexStyle` and `GraphLayout` as in `GraphPlot`; other options pass through + to `Graphics`. + +**Features**: +- `Protected`. Implemented in `src/graph/hyp_plot.c`; Mathematica has no + built-in hypergraph drawing (the Function Repository's `HypergraphPlot` is the + model). Deterministic, like `GraphPlot`. +- **Layout**: the vertices are placed by the stress layout of the star + expansion (one extra node per hyperedge of two or more vertices, joined to its + members), so the members of a hyperedge sit around a common centre and + hyperedges sharing vertices are drawn side by side. `GraphLayout` picks another + embedding of the star expansion; `VertexCoordinates` overrides any subset. +- **Hyperedges**: each is the convex hull of its (distinct) members, inflated + by a margin with rounded corners (sampled every 15 degrees), drawn as a + translucent (`Opacity[0.22]`) filled `Polygon` with a darker outline, in its + own colour of the `ColorData[97]` palette (cycled by hyperedge index). A + hyperedge of size 2 is therefore a stadium (a thick translucent line with round + caps) and one of size 1 a circle around its vertex; an empty one is not drawn. + Larger shapes are drawn first so smaller ones stay visible; a hyperedge sharing + a vertex with smaller ones gets a wider margin, so nested and repeated + hyperedges show as concentric outlines. +- **Vertices** are dark disks drawn on top; labels avoid the directions of the + vertex's hyperedges. +- Unevaluated on a non-hypergraph. + +```mathematica +In[1]:= Head[HypergraphPlot[Hypergraph[{{1,2,3},{3,4},{4,5,6},{7}}]]] +Out[1]= Graphics + +In[2]:= Count[HypergraphPlot[Hypergraph[{{1,2,3},{3,4},{4,5,6},{7}}]], _Polygon, Infinity] +Out[2]= 4 + +In[3]:= Count[HypergraphPlot[{{1, 2, 3}, {3, 4}}, VertexLabels -> "Name"], _Text, Infinity] +Out[3]= 4 + +In[4]:= HypergraphPlot[Hypergraph[{{1,2,3}}]] === HypergraphPlot[Hypergraph[{{1,2,3}}]] +Out[4]= True + +In[5]:= HypergraphPlot[5] +Out[5]= HypergraphPlot[5] +``` + ## HypergraphToGraph - `HypergraphToGraph[h]`: the directed `Graph` obtained by reading `h` as an diff --git a/docs/spec/changelog/2026-09-28.md b/docs/spec/changelog/2026-09-28.md index fdf68599..43aa18a7 100644 --- a/docs/spec/changelog/2026-09-28.md +++ b/docs/spec/changelog/2026-09-28.md @@ -1,5 +1,40 @@ # Changelog: week of 2026-09-28 (Mon) – 2026-10-04 (Sun) +## `GraphPlot` rewritten for publication-quality drawings; new `HypergraphPlot` (v0.237) + +`GraphPlot` was an MVP (every vertex on a circle, directed edges as plain lines, axes on). +It is now a deterministic layout engine (`src/graph/glayout.c`) plus a drawing layer +(`src/graph/graphplot.c`) shared with the new `HypergraphPlot` (`src/graph/hyp_plot.c`). + +- **Layouts.** Default: tidy layered trees for branching forests, layered drawings for DAGs + (longest-path layers, dummy vertices on long edges, barycentre crossing reduction, + isotonic-regression x placement), stress majorization (SMACOF on BFS distances, Pivot MDS + start) for everything else. Small components try circle and Tutte starts and keep the + drawing with the fewest crossings among those within 10% of the least stress (a + crossing-free drawing may cost 2.5x when it removes four or more crossings), then get a + canonical orientation: principal axis horizontal, axis snap and exact straightening for + grids, vertex 1 on top for symmetric graphs. Components are shelf-packed; isolated vertices + form a square block. `GraphLayout -> "StressEmbedding" | "SpringElectricalEmbedding" | + "CircularEmbedding" | "LayeredEmbedding" | "BipartiteEmbedding" | "GridEmbedding"`. + Full stress up to 1000 vertices per component, Pivot MDS above; 500 vertices in ~0.15 s. +- **Options.** `VertexCoordinates` (list or rules), `VertexLabels` (`"Name"`, `Automatic`, + rules; placed in the widest gap between incident edges, never clipped), `GraphHighlight` + (red, larger vertices, 2.5x thicker edges), `VertexStyle`, `EdgeStyle`, + `EdgeLabels -> "EdgeWeight"` or rules, `VertexSize`; other options pass through to + `Graphics`. `GraphPlot[{u -> v, ...}]` draws a rule list. +- **Look.** Mathematica's vertex blue with a darker rim, grey-blue edges, `Arrow`s whose tips + stop at the target disk (mutual pairs offset apart), no axes, equal aspect ratio, a tight + explicit `PlotRange` and `ImageSize`. +- **`HypergraphPlot[h]`** (new, `Protected`): stress layout of the star expansion; each + hyperedge a translucent rounded convex hull in its own `ColorData[97]` colour (a stadium + for size 2, a circle for size 1), largest first, nested hyperedges with widening margins; + vertex disks and labels on top. +- **PDF export** (`src/graphics/graphics_export.c`): `Text` offsets and + `Style[s, n | FontSize -> n | colour]`, the `Arrowheads[s]` directive, arrow shafts that + stop inside the head, and `AspectRatio -> Automatic` (equal x/y scale). +- Tests: `tests/test_graphplot.c` (determinism, options, degenerate and 500-vertex graphs, + hypergraphs, PDF export). `leaks --atExit`: 0 leaks. + ## Merge: `worked-example-fixes` lands on main (v0.233) The `worked-example-fixes` branch (8 commits, Michael Sollami) and main's `flint_field_gcd` diff --git a/src/graph/glayout.c b/src/graph/glayout.c new file mode 100644 index 00000000..ae257f57 --- /dev/null +++ b/src/graph/glayout.c @@ -0,0 +1,1374 @@ +/* glayout.c - deterministic graph layout engine (see glayout.h). + * + * METHODS + * Stress SMACOF stress majorization on BFS graph distances (Gansner, + * Koren & North 2004), with weights d^-2, initialised by Pivot MDS + * (Brandes & Pich 2006). Components above STRESS_EXACT_MAX + * vertices skip the O(n^2)-per-sweep majorization and keep the + * Pivot MDS layout, which is O(k (n + m)) for k pivots. + * Spring Hu's spring-electrical model (attraction d^2/K, repulsion + * C K^2 / d) with adaptive step length, initialised by Pivot MDS. + * Above SPRING_EXACT_MAX vertices repulsion is cut off at a + * radius and evaluated over a uniform grid (Fruchterman-Reingold). + * Layered Trees: a tidy tree (leaves in DFS order, parents centred over + * their children) from the tree centre, or from the source of an + * arborescence. DAGs: longest-path layering, then barycentre + * sweeps to reduce crossings (best ordering kept), then x + * positions by isotonic regression (pool-adjacent-violators) + * towards neighbour barycentres with unit minimum separation. + * Other graphs: BFS layers from the graph centre. + * Circular, Bipartite, Grid closed forms, barycentre-ordered for Bipartite. + * + * Every component is laid out on its own, normalised to unit mean edge length, + * rotated to a canonical orientation (principal axis horizontal; snapped to the + * axes when the edges are mostly axis-parallel, as in grids), and the + * components are then shelf-packed, largest first. + * + * DETERMINISM. No random numbers anywhere: pivots, orders, starts and + * tie-breaks are all fixed functions of the vertex indices, so a graph always + * gets bit-identical coordinates. + */ + +#include "glayout.h" +#include "sym_names.h" +#include +#include +#include + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +#define STRESS_EXACT_MAX 1000 /* full SMACOF up to this component size */ +#define STRESS_MAX_ITERS 400 +#define STRESS_MULTISTART 60 /* components this small try several starts */ +#define SPRING_EXACT_MAX 1000 /* exact O(n^2) repulsion up to this size */ +#define SPRING_MAX_ITERS 500 +#define PIVOTS 50 +#define CENTER_EXACT_MAX 2000 /* exact graph centre (all BFS) up to this size */ + +/* ------------------------------------------------------------ CSR graph --- */ + +typedef struct { + int n, m; /* vertices, edges */ + int *eu, *ev; /* edge endpoints */ + unsigned char* dir; /* directed flags (never NULL inside) */ + int *off, *adj; /* undirected adjacency, CSR */ +} Sub; + +static void sub_free(Sub* s) { + free(s->eu); free(s->ev); free(s->dir); free(s->off); free(s->adj); + memset(s, 0, sizeof(*s)); +} + +/* Builds the symmetric CSR adjacency of s's edge list. 0 on OOM. */ +static int sub_csr(Sub* s) { + s->off = calloc((size_t)s->n + 1, sizeof(int)); + s->adj = malloc(sizeof(int) * (size_t)(2 * s->m + 1)); + if (!s->off || !s->adj) return 0; + for (int k = 0; k < s->m; k++) { s->off[s->eu[k] + 1]++; s->off[s->ev[k] + 1]++; } + for (int i = 0; i < s->n; i++) s->off[i + 1] += s->off[i]; + int* fill = malloc(sizeof(int) * (size_t)(s->n + 1)); + if (!fill) return 0; + memcpy(fill, s->off, sizeof(int) * (size_t)s->n); + for (int k = 0; k < s->m; k++) { + s->adj[fill[s->eu[k]]++] = s->ev[k]; + s->adj[fill[s->ev[k]]++] = s->eu[k]; + } + free(fill); + return 1; +} + +/* BFS from src over s; d[] gets hop distances (-1 unreachable). queue is + * scratch of size n. */ +static void bfs(const Sub* s, int src, int* d, int* queue) { + for (int i = 0; i < s->n; i++) d[i] = -1; + int h = 0, t = 0; + d[src] = 0; queue[t++] = src; + while (h < t) { + int u = queue[h++]; + for (int a = s->off[u]; a < s->off[u + 1]; a++) { + int w = s->adj[a]; + if (d[w] < 0) { d[w] = d[u] + 1; queue[t++] = w; } + } + } +} + +/* ---------------------------------------------------- small linear algebra */ + +/* Cyclic Jacobi eigen-decomposition of the symmetric k x k matrix A (row + * major, destroyed). Eigenvalues to w[], eigenvectors as COLUMNS of V. */ +static void jacobi_eigen(int k, double* A, double* w, double* V) { + for (int i = 0; i < k * k; i++) V[i] = 0.0; + for (int i = 0; i < k; i++) V[i * k + i] = 1.0; + for (int sweep = 0; sweep < 60; sweep++) { + double off = 0.0, diag = 0.0; + for (int p = 0; p < k; p++) + for (int q = 0; q < k; q++) { + if (p == q) diag += A[p * k + q] * A[p * k + q]; + else off += A[p * k + q] * A[p * k + q]; + } + if (off <= 1e-22 * (diag + 1e-300)) break; + for (int p = 0; p < k - 1; p++) + for (int q = p + 1; q < k; q++) { + double apq = A[p * k + q]; + if (fabs(apq) < 1e-300) continue; + double app = A[p * k + p], aqq = A[q * k + q]; + double theta = (aqq - app) / (2.0 * apq); + double t = (theta >= 0 ? 1.0 : -1.0) + / (fabs(theta) + sqrt(theta * theta + 1.0)); + double c = 1.0 / sqrt(t * t + 1.0), sn = t * c; + for (int r = 0; r < k; r++) { /* columns p, q */ + double arp = A[r * k + p], arq = A[r * k + q]; + A[r * k + p] = c * arp - sn * arq; + A[r * k + q] = sn * arp + c * arq; + } + for (int r = 0; r < k; r++) { /* rows p, q */ + double apr = A[p * k + r], aqr = A[q * k + r]; + A[p * k + r] = c * apr - sn * aqr; + A[q * k + r] = sn * apr + c * aqr; + } + for (int r = 0; r < k; r++) { + double vrp = V[r * k + p], vrq = V[r * k + q]; + V[r * k + p] = c * vrp - sn * vrq; + V[r * k + q] = sn * vrp + c * vrq; + } + } + } + for (int i = 0; i < k; i++) w[i] = A[i * k + i]; +} + +/* ------------------------------------------------------------ utilities --- */ + +/* Scales xy so the mean edge length is 1 (no-op without edges or when all + * edges are degenerate). */ +static void normalise_edge_length(const Sub* s, double* xy) { + if (s->m == 0) return; + double tot = 0.0; + for (int k = 0; k < s->m; k++) { + double dx = xy[2 * s->eu[k]] - xy[2 * s->ev[k]]; + double dy = xy[2 * s->eu[k] + 1] - xy[2 * s->ev[k] + 1]; + tot += sqrt(dx * dx + dy * dy); + } + if (tot <= 1e-300) return; + double f = (double)s->m / tot; + for (int i = 0; i < 2 * s->n; i++) xy[i] *= f; +} + +/* Tiny deterministic perturbation that breaks exact symmetries (coincident or + * collinear starts) without any randomness. */ +static void jitter(int n, double* xy, double amp) { + for (int i = 0; i < n; i++) { + xy[2 * i] += amp * sin(1.0 + 2.399963 * (double)i); + xy[2 * i + 1] += amp * cos(1.0 + 2.399963 * (double)i); + } +} + +/* Rotates to the canonical orientation: principal axis horizontal, then + * snapped to the axes when the edge directions are strongly four-fold (grids, + * cubes), then reflected so that vertex 0 sits left of and above the centre. */ +static void canonical_orientation(const Sub* s, double* xy) { + int n = s->n; + if (n < 2) return; + double cx = 0, cy = 0; + for (int i = 0; i < n; i++) { cx += xy[2 * i]; cy += xy[2 * i + 1]; } + cx /= n; cy /= n; + double sxx = 0, syy = 0, sxy = 0; + for (int i = 0; i < n; i++) { + double dx = xy[2 * i] - cx, dy = xy[2 * i + 1] - cy; + sxx += dx * dx; syy += dy * dy; sxy += dx * dy; + } + double th = 0.5 * atan2(2.0 * sxy, sxx - syy); /* principal axis angle */ + /* Only rotate when there is a genuine principal axis. */ + double aniso = sqrt((sxx - syy) * (sxx - syy) + 4 * sxy * sxy) / (sxx + syy + 1e-300); + double rot = aniso > 0.05 ? -th : 0.0; + /* Four-fold snap, measured after the principal rotation. */ + int snapped = 0; + if (s->m > 0) { + double S = 0, C = 0; + for (int k = 0; k < s->m; k++) { + double dx = xy[2 * s->ev[k]] - xy[2 * s->eu[k]]; + double dy = xy[2 * s->ev[k] + 1] - xy[2 * s->eu[k] + 1]; + double ph = atan2(dy, dx) + rot; + S += sin(4 * ph); C += cos(4 * ph); + } + double R = sqrt(S * S + C * C) / s->m; + if (R > 0.6) { rot -= atan2(S, C) / 4.0; snapped = 1; } + } + if (aniso <= 0.05 && !snapped) { + /* No preferred axis (cycles, K_n, vertex-transitive graphs): stand + * vertex 0 straight above the centre, as a textbook polygon does. */ + double dx = xy[0] - cx, dy = xy[1] - cy; + if (dx * dx + dy * dy > 1e-18) rot = M_PI / 2 - atan2(dy, dx); + } + double c = cos(rot), sn = sin(rot); + for (int i = 0; i < n; i++) { + double dx = xy[2 * i] - cx, dy = xy[2 * i + 1] - cy; + xy[2 * i] = c * dx - sn * dy; + xy[2 * i + 1] = sn * dx + c * dy; + } + /* Clean up rounding noise so an axis-aligned layout is exactly aligned. */ + for (int i = 0; i < 2 * n; i++) if (fabs(xy[i]) < 1e-12) xy[i] = 0.0; + if (xy[0] > 1e-9) for (int i = 0; i < n; i++) xy[2 * i] = -xy[2 * i]; + if (xy[1] < -1e-9) for (int i = 0; i < n; i++) xy[2 * i + 1] = -xy[2 * i + 1]; +} + +/* Union-find root with path halving. */ +static int uf_find(int* p, int i) { + while (p[i] != i) { p[i] = p[p[i]]; i = p[i]; } + return i; +} + +/* Orthogonal straightening. When nearly every edge of an oriented drawing is + * within 15 degrees of an axis (grids, ladders, tori drawn flat), stress leaves + * the rows gently bowed, because hop distance is not Euclidean distance along + * a diagonal. Vertices joined by near-horizontal edges then share one y (their + * mean), and by near-vertical edges one x, so a grid comes out as a grid. */ +static void straighten(const Sub* s, double* xy) { + int n = s->n, m = s->m, aligned = 0; + if (m < 4) return; + const double tol = tan(15.0 * M_PI / 180.0); + for (int k = 0; k < m; k++) { + double dx = fabs(xy[2 * s->ev[k]] - xy[2 * s->eu[k]]); + double dy = fabs(xy[2 * s->ev[k] + 1] - xy[2 * s->eu[k] + 1]); + if (dy <= tol * dx || dx <= tol * dy) aligned++; + } + if (aligned < m) return; /* every edge must be axis-like */ + int* p = malloc(sizeof(int) * (size_t)n); + double* sum = malloc(sizeof(double) * (size_t)n); + int* cnt = malloc(sizeof(int) * (size_t)n); + if (p && sum && cnt) { + for (int axis = 0; axis < 2; axis++) { /* 0: rows share y; 1: columns share x */ + for (int i = 0; i < n; i++) { p[i] = i; sum[i] = 0; cnt[i] = 0; } + for (int k = 0; k < m; k++) { + double dx = fabs(xy[2 * s->ev[k]] - xy[2 * s->eu[k]]); + double dy = fabs(xy[2 * s->ev[k] + 1] - xy[2 * s->eu[k] + 1]); + int horiz = dy <= tol * dx; + if (horiz == (axis == 0)) { + int a = uf_find(p, s->eu[k]), b = uf_find(p, s->ev[k]); + if (a != b) p[a < b ? b : a] = a < b ? a : b; + } + } + int co = axis == 0 ? 1 : 0; + for (int i = 0; i < n; i++) { int r = uf_find(p, i); sum[r] += xy[2 * i + co]; cnt[r]++; } + for (int i = 0; i < n; i++) { int r = uf_find(p, i); xy[2 * i + co] = sum[r] / cnt[r]; } + } + } + free(p); free(sum); free(cnt); +} + +/* --------------------------------------------------------------- circle --- */ + +static void layout_circle(int n, const int* order, double* xy) { + if (n == 1) { xy[0] = xy[1] = 0; return; } + double R = (n == 2) ? 0.5 : 0.5 / sin(M_PI / n); /* unit chord */ + for (int r = 0; r < n; r++) { + int i = order ? order[r] : r; + double t = M_PI / 2 + 2.0 * M_PI * r / n; + xy[2 * i] = R * cos(t); xy[2 * i + 1] = R * sin(t); + } +} + +/* ------------------------------------------------------------ Pivot MDS --- */ + +/* Pivot MDS of s into xy. D (n x n, may be NULL) supplies distances when the + * full matrix is already known. 0 on OOM. */ +static int pivot_mds(const Sub* s, const int* D, double* xy) { + int n = s->n; + if (n <= 2) { + xy[0] = 0; xy[1] = 0; + if (n == 2) { xy[2] = 1; xy[3] = 0; } + return 1; + } + int k = n < PIVOTS ? n : PIVOTS; + double* C = malloc(sizeof(double) * (size_t)n * (size_t)k); + int* piv = malloc(sizeof(int) * (size_t)k); + int* mind = malloc(sizeof(int) * (size_t)n); + int* d = malloc(sizeof(int) * (size_t)n); + int* q = malloc(sizeof(int) * (size_t)n); + double* M = malloc(sizeof(double) * (size_t)k * (size_t)k); + double* V = malloc(sizeof(double) * (size_t)k * (size_t)k); + double* w = malloc(sizeof(double) * (size_t)k); + double* cm = calloc((size_t)k, sizeof(double)); + int ok = C && piv && mind && d && q && M && V && w && cm; + if (ok) { + /* First pivot: highest degree (lowest index on ties); then max-min. */ + int best = 0; + for (int i = 1; i < n; i++) + if (s->off[i + 1] - s->off[i] > s->off[best + 1] - s->off[best]) best = i; + for (int i = 0; i < n; i++) mind[i] = 1 << 29; + for (int j = 0; j < k; j++) { + piv[j] = best; + if (D) for (int i = 0; i < n; i++) d[i] = D[(size_t)best * n + i]; + else bfs(s, best, d, q); + for (int i = 0; i < n; i++) { + double dd = d[i] < 0 ? n : d[i]; + C[(size_t)i * k + j] = dd * dd; + if (d[i] >= 0 && d[i] < mind[i]) mind[i] = d[i]; + } + best = 0; + for (int i = 1; i < n; i++) if (mind[i] > mind[best]) best = i; + } + /* Double centring. */ + double gm = 0; + for (int j = 0; j < k; j++) { + for (int i = 0; i < n; i++) cm[j] += C[(size_t)i * k + j]; + cm[j] /= n; gm += cm[j]; + } + gm /= k; + for (int i = 0; i < n; i++) { + double rm = 0; + for (int j = 0; j < k; j++) rm += C[(size_t)i * k + j]; + rm /= k; + for (int j = 0; j < k; j++) + C[(size_t)i * k + j] = -0.5 * (C[(size_t)i * k + j] - rm - cm[j] + gm); + } + for (int a = 0; a < k; a++) + for (int b = a; b < k; b++) { + double acc = 0; + for (int i = 0; i < n; i++) acc += C[(size_t)i * k + a] * C[(size_t)i * k + b]; + M[a * k + b] = M[b * k + a] = acc; + } + jacobi_eigen(k, M, w, V); + int e1 = 0; + for (int i = 1; i < k; i++) if (w[i] > w[e1]) e1 = i; + int e2 = (e1 == 0) ? 1 : 0; + for (int i = 0; i < k; i++) if (i != e1 && w[i] > w[e2]) e2 = i; + for (int i = 0; i < n; i++) { + double x = 0, y = 0; + for (int j = 0; j < k; j++) { + x += C[(size_t)i * k + j] * V[j * k + e1]; + y += C[(size_t)i * k + j] * V[j * k + e2]; + } + xy[2 * i] = x; xy[2 * i + 1] = y; + } + normalise_edge_length(s, xy); + } + free(C); free(piv); free(mind); free(d); free(q); free(M); free(V); free(w); free(cm); + return ok; +} + +/* --------------------------------------------------------------- stress --- */ + +static double stress_value(int n, const int* D, const double* xy) { + double st = 0; + for (int i = 0; i < n; i++) + for (int j = i + 1; j < n; j++) { + double dij = D[(size_t)i * n + j]; + double dx = xy[2 * i] - xy[2 * j], dy = xy[2 * i + 1] - xy[2 * j + 1]; + double e = sqrt(dx * dx + dy * dy) - dij; + st += e * e / (dij * dij); + } + return st; +} + +/* Localised SMACOF (in-place majorization sweeps) from the start in xy. + * Returns the final stress. */ +static void smacof_sweep(int n, const int* D, double* xy) { + for (int i = 0; i < n; i++) { + double sx = 0, sy = 0, sw = 0; + const int* Di = D + (size_t)i * n; + double xi = xy[2 * i], yi = xy[2 * i + 1]; + for (int j = 0; j < n; j++) { + if (j == i) continue; + double dij = Di[j]; + double w = 1.0 / (dij * dij); + double dx = xi - xy[2 * j], dy = yi - xy[2 * j + 1]; + double dist = sqrt(dx * dx + dy * dy); + sw += w; + if (dist > 1e-12) { + sx += w * (xy[2 * j] + dij * dx / dist); + sy += w * (xy[2 * j + 1] + dij * dy / dist); + } else { + sx += w * xy[2 * j]; sy += w * xy[2 * j + 1]; + } + } + if (sw > 0) { xy[2 * i] = sx / sw; xy[2 * i + 1] = sy / sw; } + } +} + +static double smacof(int n, const int* D, double* xy) { + double prev = stress_value(n, D, xy); + for (int it = 0; it < STRESS_MAX_ITERS; it++) { + smacof_sweep(n, D, xy); + if (it % 4 == 3) { + double cur = stress_value(n, D, xy); + if (prev - cur <= 1e-6 * prev + 1e-12) { prev = cur; break; } + prev = cur; + } + } + return stress_value(n, D, xy); +} + +/* All-pairs BFS distances (n x n ints), or NULL on OOM. */ +static int* all_pairs(const Sub* s) { + int n = s->n; + int* D = malloc(sizeof(int) * (size_t)n * (size_t)n); + int* q = malloc(sizeof(int) * (size_t)n); + if (!D || !q) { free(D); free(q); return NULL; } + for (int i = 0; i < n; i++) bfs(s, i, D + (size_t)i * n, q); + free(q); + return D; +} + +/* Proper crossings between non-adjacent edges of a straight-line drawing. */ +static long drawing_crossings(const Sub* s, const double* xy) { + long c = 0; + for (int a = 0; a < s->m; a++) { + int p = s->eu[a], q = s->ev[a]; + double ax = xy[2 * p], ay = xy[2 * p + 1], bx = xy[2 * q], by = xy[2 * q + 1]; + for (int b = a + 1; b < s->m; b++) { + int r = s->eu[b], t = s->ev[b]; + if (r == p || r == q || t == p || t == q) continue; + double cx = xy[2 * r], cy = xy[2 * r + 1], dx = xy[2 * t], dy = xy[2 * t + 1]; + double d1 = (bx - ax) * (cy - ay) - (by - ay) * (cx - ax); + double d2 = (bx - ax) * (dy - ay) - (by - ay) * (dx - ax); + double d3 = (dx - cx) * (ay - cy) - (dy - cy) * (ax - cx); + double d4 = (dx - cx) * (by - cy) - (dy - cy) * (bx - cx); + if (((d1 > 0 && d2 < 0) || (d1 < 0 && d2 > 0)) + && ((d3 > 0 && d4 < 0) || (d3 < 0 && d4 > 0))) c++; + } + } + return c; +} + +/* A shortest cycle through v (BFS with branch labels): writes its vertices in + * order to cyc and returns its length, or 0 if v lies on no cycle. */ +static int shortest_cycle_through(const Sub* s, int v, int* cyc, + int* d, int* par, int* br, int* q) { + int n = s->n; + for (int i = 0; i < n; i++) { d[i] = -1; par[i] = -1; br[i] = -1; } + int h = 0, t = 0, best = 1 << 30, ba = -1, bb = -1; + d[v] = 0; q[t++] = v; + while (h < t) { + int u = q[h++]; + if (2 * d[u] + 1 >= best) break; + for (int a = s->off[u]; a < s->off[u + 1]; a++) { + int w = s->adj[a]; + if (d[w] < 0) { + d[w] = d[u] + 1; par[w] = u; br[w] = (u == v) ? w : br[u]; + q[t++] = w; + } else if (w != par[u] && u != v && w != v && br[w] != br[u]) { + int len = d[u] + d[w] + 1; + if (len < best) { best = len; ba = u; bb = w; } + } else if (w == v && u != v && par[u] != v) { + /* closing edge straight back to v */ + int len = d[u] + 1; + if (len < best) { best = len; ba = u; bb = v; } + } + } + } + if (ba < 0) return 0; + int k = 0; + for (int x = ba; x != v; x = par[x]) cyc[k++] = x; + cyc[k++] = v; + /* reverse so the cycle reads v .. ba, then append bb .. (towards v) */ + for (int i = 0; i < k / 2; i++) { int tmp = cyc[i]; cyc[i] = cyc[k - 1 - i]; cyc[k - 1 - i] = tmp; } + if (bb != v) for (int x = bb; x != v; x = par[x]) cyc[k++] = x; + return k; +} + +/* Tutte-style start: the cycle cyc on a circle, every other vertex at the + * barycentre of its neighbours (Gauss-Seidel). Crossing-free for a + * 3-connected planar graph whose cycle is a face. */ +static void tutte_start(const Sub* s, const int* cyc, int k, double* xy, unsigned char* pinned) { + int n = s->n; + memset(pinned, 0, (size_t)n); + for (int i = 0; i < 2 * n; i++) xy[i] = 0; + double R = (0.5 / sin(M_PI / k)) * (n > k ? sqrt((double)n / k) : 1.0); + for (int i = 0; i < k; i++) { + double t = M_PI / 2 + 2.0 * M_PI * i / k; + xy[2 * cyc[i]] = R * cos(t); xy[2 * cyc[i] + 1] = R * sin(t); pinned[cyc[i]] = 1; + } + for (int it = 0; it < 300; it++) + for (int i = 0; i < n; i++) { + if (pinned[i] || s->off[i + 1] == s->off[i]) continue; + double ax = 0, ay = 0; int na = s->off[i + 1] - s->off[i]; + for (int a = s->off[i]; a < s->off[i + 1]; a++) { ax += xy[2 * s->adj[a]]; ay += xy[2 * s->adj[a] + 1]; } + xy[2 * i] = ax / na; xy[2 * i + 1] = ay / na; + } +} + +/* Rescales xy by the stress-optimal factor and returns the resulting stress, + * so drawings of different scale compare fairly. */ +static double scaled_stress(int n, const int* D, double* xy) { + double num = 0, den = 0; + for (int i = 0; i < n; i++) + for (int j = i + 1; j < n; j++) { + double dij = D[(size_t)i * n + j], w = 1.0 / (dij * dij); + double dx = xy[2 * i] - xy[2 * j], dy = xy[2 * i + 1] - xy[2 * j + 1]; + double e = sqrt(dx * dx + dy * dy); + num += w * dij * e; den += w * e * e; + } + if (den > 0) for (int i = 0; i < 2 * n; i++) xy[i] *= num / den; + return stress_value(n, D, xy); +} + +/* Planarity-preserving refinement of a Tutte drawing: Jacobi SMACOF sweeps + * that move only the unpinned (interior) vertices -- the outer cycle stays a + * regular polygon, and Jacobi updates keep every symmetry of the start -- + * stopped at the last crossing-free iterate. `nxt` is scratch of 2n doubles. + * Returns the stress after optimal rescaling (so it compares with others). */ +static double smacof_planar(const Sub* s, const int* D, double* xy, + const unsigned char* pinned, double* nxt) { + int n = s->n; + double prev = stress_value(n, D, xy); + for (int it = 0; it < STRESS_MAX_ITERS; it++) { + for (int i = 0; i < n; i++) { + nxt[2 * i] = xy[2 * i]; nxt[2 * i + 1] = xy[2 * i + 1]; + if (pinned[i]) continue; + double sx = 0, sy = 0, sw = 0; + for (int j = 0; j < n; j++) { + if (j == i) continue; + double dij = D[(size_t)i * n + j], w = 1.0 / (dij * dij); + double dx = xy[2 * i] - xy[2 * j], dy = xy[2 * i + 1] - xy[2 * j + 1]; + double dist = sqrt(dx * dx + dy * dy); + sw += w; + sx += w * (xy[2 * j] + (dist > 1e-12 ? dij * dx / dist : 0)); + sy += w * (xy[2 * j + 1] + (dist > 1e-12 ? dij * dy / dist : 0)); + } + if (sw > 0) { nxt[2 * i] = sx / sw; nxt[2 * i + 1] = sy / sw; } + } + if (drawing_crossings(s, nxt) > 0) break; + memcpy(xy, nxt, sizeof(double) * 2 * (size_t)n); + double cur = stress_value(n, D, xy); + if (prev - cur <= 1e-6 * prev + 1e-12) break; + prev = cur; + } + return scaled_stress(n, D, xy); +} + +static int layout_stress(const Sub* s, double* xy) { + int n = s->n; + if (n <= 2) return pivot_mds(s, NULL, xy); + if (n > STRESS_EXACT_MAX) return pivot_mds(s, NULL, xy); + int* D = all_pairs(s); + if (!D) return 0; + int ok = pivot_mds(s, D, xy); + if (ok) { + jitter(n, xy, 1e-3); + double best = smacof(n, D, xy); + if (n <= STRESS_MULTISTART && s->m <= 400) { + /* Extra deterministic starts -- a circle in VertexList order, and + * circles in BFS order from several roots -- and keep, among the + * layouts within 10% of the least stress, the one with the fewest + * edge crossings (least stress on ties). */ + int ncirc = 1 + (n < 8 ? n : 8); + int ntutte = 3; + int nstart = ncirc + ntutte; + int nall = nstart + ntutte; /* + the Tutte drawings unrefined */ + double* alt = malloc(sizeof(double) * 2 * (size_t)n); + double* cand = malloc(sizeof(double) * 2 * (size_t)n * (size_t)(nall + 1)); + double* st = malloc(sizeof(double) * (size_t)(nall + 1)); + long* cr = malloc(sizeof(long) * (size_t)(nall + 1)); + int* order = malloc(sizeof(int) * (size_t)n); + int* sc = malloc(sizeof(int) * 5 * (size_t)n); + unsigned char* pin = malloc((size_t)n); + double* pin_xy = malloc(sizeof(double) * 2 * (size_t)n); + if (alt && cand && st && cr && order && sc && pin && pin_xy) { + memcpy(cand, xy, sizeof(double) * 2 * (size_t)n); + st[0] = best; cr[0] = drawing_crossings(s, xy); + for (int k = 0; k < nstart; k++) { + if (k >= ncirc) { + /* Tutte starts from shortest cycles through vertices + * spread over the index range. */ + int v = (int)(((long)(k - ncirc) * n) / ntutte); + int len = shortest_cycle_through(s, v, order, sc, sc + n, sc + 2 * n, sc + 3 * n); + int raw = nstart + 1 + (k - ncirc); + if (len < 3) { + st[k + 1] = st[raw] = 1e300; cr[k + 1] = cr[raw] = 1L << 40; + continue; + } + tutte_start(s, order, len, alt, pin); + { /* keep the unrefined drawing as a candidate too */ + double* c = cand + 2 * (size_t)n * (size_t)raw; + memcpy(c, alt, sizeof(double) * 2 * (size_t)n); + cr[raw] = drawing_crossings(s, c); + st[raw] = cr[raw] == 0 ? smacof_planar(s, D, c, pin, pin_xy) + : scaled_stress(n, D, c); + } + } else if (k == 0) layout_circle(n, NULL, alt); + else { + int root = (int)(((long)(k - 1) * n) / (ncirc - 1)); + const int* Dr = D + (size_t)root * n; + int t = 0; + for (int lev = 0; t < n; lev++) + for (int i = 0; i < n; i++) if (Dr[i] == lev) order[t++] = i; + layout_circle(n, order, alt); + } + jitter(n, alt, 1e-3); + st[k + 1] = smacof(n, D, alt); + cr[k + 1] = drawing_crossings(s, alt); + memcpy(cand + 2 * (size_t)n * (size_t)(k + 1), alt, sizeof(double) * 2 * (size_t)n); + if (st[k + 1] < best) best = st[k + 1]; + } + /* Crossings of the least-stress drawing: a crossing-free + * (Tutte-derived) drawing may cost up to 2.5x the stress, but + * only when it removes many crossings (a dodecahedron's ten, + * not a cube's two -- the Necker cube is the textbook cube). */ + long crbest = 1L << 40; + for (int k = 0; k <= nall; k++) + if (st[k] <= best * (1 + 1e-9) && cr[k] < crbest) crbest = cr[k]; + int planar_ok = crbest >= 4 && crbest * 6 >= s->m; + int pick = -1; + for (int k = 0; k <= nall; k++) { + double win = (cr[k] == 0 && planar_ok) ? 2.5 : 1.10; + if (st[k] > best * win + 1e-12) continue; + if (pick < 0 || cr[k] < cr[pick] + || (cr[k] == cr[pick] && st[k] < st[pick] * (1 - 1e-9))) pick = k; + } + if (pick > 0) memcpy(xy, cand + 2 * (size_t)n * (size_t)pick, sizeof(double) * 2 * (size_t)n); + } + free(alt); free(cand); free(st); free(cr); free(order); free(sc); free(pin); free(pin_xy); + } + } + free(D); + return ok; +} + +/* ----------------------------------------------------- spring-electrical --- */ + +static int layout_spring(const Sub* s, double* xy) { + int n = s->n; + if (!pivot_mds(s, NULL, xy)) return 0; + if (n <= 2) return 1; + jitter(n, xy, 1e-3); + const double K = 1.0, Cr = 0.2, cool = 0.9; + double* f = malloc(sizeof(double) * 2 * (size_t)n); + int* head = NULL; int* next = NULL; + int use_grid = n > SPRING_EXACT_MAX; + if (use_grid) next = malloc(sizeof(int) * (size_t)n); + if (!f || (use_grid && !next)) { free(f); free(next); return 0; } + double step = K, energy0 = 1e300; + int progress = 0; + const double R = 4.0 * K; /* grid repulsion cut-off */ + for (int it = 0; it < SPRING_MAX_ITERS; it++) { + memset(f, 0, sizeof(double) * 2 * (size_t)n); + if (!use_grid) { + for (int i = 0; i < n; i++) + for (int j = i + 1; j < n; j++) { + double dx = xy[2 * i] - xy[2 * j], dy = xy[2 * i + 1] - xy[2 * j + 1]; + double d2 = dx * dx + dy * dy + 1e-12; + double g = Cr * K * K / d2; + f[2 * i] += g * dx; f[2 * i + 1] += g * dy; + f[2 * j] -= g * dx; f[2 * j + 1] -= g * dy; + } + } else { + double x0 = 1e300, y0 = 1e300, x1 = -1e300, y1 = -1e300; + for (int i = 0; i < n; i++) { + if (xy[2 * i] < x0) x0 = xy[2 * i]; + if (xy[2 * i] > x1) x1 = xy[2 * i]; + if (xy[2 * i + 1] < y0) y0 = xy[2 * i + 1]; + if (xy[2 * i + 1] > y1) y1 = xy[2 * i + 1]; + } + int gx = (int)((x1 - x0) / R) + 1, gy = (int)((y1 - y0) / R) + 1; + if ((double)gx * gy > 4.0 * n) { gx = gy = (int)sqrt(4.0 * n) + 1; } + double cwx = (x1 - x0) / gx + 1e-9, cwy = (y1 - y0) / gy + 1e-9; + free(head); + head = malloc(sizeof(int) * (size_t)gx * (size_t)gy); + if (!head) break; + for (int c = 0; c < gx * gy; c++) head[c] = -1; + for (int i = 0; i < n; i++) { + int cx = (int)((xy[2 * i] - x0) / cwx), cy = (int)((xy[2 * i + 1] - y0) / cwy); + if (cx >= gx) cx = gx - 1; + if (cy >= gy) cy = gy - 1; + int c = cy * gx + cx; + next[i] = head[c]; head[c] = i; + } + for (int i = 0; i < n; i++) { + int cx = (int)((xy[2 * i] - x0) / cwx), cy = (int)((xy[2 * i + 1] - y0) / cwy); + if (cx >= gx) cx = gx - 1; + if (cy >= gy) cy = gy - 1; + for (int ox = -1; ox <= 1; ox++) + for (int oy = -1; oy <= 1; oy++) { + int ax = cx + ox, ay = cy + oy; + if (ax < 0 || ay < 0 || ax >= gx || ay >= gy) continue; + for (int j = head[ay * gx + ax]; j >= 0; j = next[j]) { + if (j == i) continue; + double dx = xy[2 * i] - xy[2 * j], dy = xy[2 * i + 1] - xy[2 * j + 1]; + double d2 = dx * dx + dy * dy + 1e-12; + if (d2 > R * R) continue; + double g = Cr * K * K / d2; + f[2 * i] += g * dx; f[2 * i + 1] += g * dy; + } + } + } + } + for (int k = 0; k < s->m; k++) { + int a = s->eu[k], b = s->ev[k]; + double dx = xy[2 * a] - xy[2 * b], dy = xy[2 * a + 1] - xy[2 * b + 1]; + double d = sqrt(dx * dx + dy * dy); + double g = d / K; /* |F| = d^2 / K */ + f[2 * a] -= g * dx; f[2 * a + 1] -= g * dy; + f[2 * b] += g * dx; f[2 * b + 1] += g * dy; + } + double energy = 0, moved = 0; + for (int i = 0; i < n; i++) { + double fx = f[2 * i], fy = f[2 * i + 1]; + double fl = sqrt(fx * fx + fy * fy); + energy += fl * fl; + if (fl > 1e-300) { + xy[2 * i] += step * fx / fl; xy[2 * i + 1] += step * fy / fl; + moved += step; + } + } + if (energy < energy0) { + if (++progress >= 5) { progress = 0; step /= cool; } + } else { progress = 0; step *= cool; } + energy0 = energy; + if (moved < 1e-3 * K * n) break; + } + free(f); free(head); free(next); + normalise_edge_length(s, xy); + return 1; +} + +/* -------------------------------------------------------------- layered --- */ + +/* Pool-adjacent-violators: the least-squares non-decreasing fit of y[0..n-1], + * in place. blk/sum/cnt are scratch of size n. */ +static void pava(int n, double* y, int* start, double* sum, int* cnt) { + int nb = 0; + for (int i = 0; i < n; i++) { + start[nb] = i; sum[nb] = y[i]; cnt[nb] = 1; nb++; + while (nb > 1 && sum[nb - 2] / cnt[nb - 2] > sum[nb - 1] / cnt[nb - 1]) { + sum[nb - 2] += sum[nb - 1]; cnt[nb - 2] += cnt[nb - 1]; nb--; + } + } + for (int b = 0; b < nb; b++) { + double v = sum[b] / cnt[b]; + for (int i = start[b]; i < start[b] + cnt[b]; i++) y[i] = v; + } +} + +/* Longest-path layering over directed edges. Returns 1 and fills layer[] if + * every edge is directed and they form a DAG; 0 otherwise (or on OOM). */ +static int dag_layers(const Sub* s, int* layer) { + int n = s->n; + for (int k = 0; k < s->m; k++) if (!s->dir[k]) return 0; + int* indeg = calloc((size_t)n, sizeof(int)); + int* off = calloc((size_t)n + 1, sizeof(int)); + int* out = malloc(sizeof(int) * (size_t)(s->m + 1)); + int* q = malloc(sizeof(int) * (size_t)n); + int ok = indeg && off && out && q; + if (ok) { + for (int k = 0; k < s->m; k++) { indeg[s->ev[k]]++; off[s->eu[k] + 1]++; } + for (int i = 0; i < n; i++) off[i + 1] += off[i]; + int* fill = malloc(sizeof(int) * (size_t)(n + 1)); + if (!fill) ok = 0; + else { + memcpy(fill, off, sizeof(int) * (size_t)n); + for (int k = 0; k < s->m; k++) out[fill[s->eu[k]]++] = s->ev[k]; + free(fill); + int h = 0, t = 0; + for (int i = 0; i < n; i++) { layer[i] = 0; if (!indeg[i]) q[t++] = i; } + while (h < t) { + int u = q[h++]; + for (int a = off[u]; a < off[u + 1]; a++) { + int w = out[a]; + if (layer[u] + 1 > layer[w]) layer[w] = layer[u] + 1; + if (--indeg[w] == 0) q[t++] = w; + } + } + ok = (t == n); + } + } + free(indeg); free(off); free(out); free(q); + return ok; +} + +/* A root for BFS layering: the graph centre (minimum eccentricity, lowest + * index on ties) when affordable, else the vertex of maximum degree. */ +static int layout_root(const Sub* s) { + int n = s->n, best = 0; + if (n <= CENTER_EXACT_MAX) { + int* d = malloc(sizeof(int) * (size_t)n); + int* q = malloc(sizeof(int) * (size_t)n); + if (d && q) { + int bestecc = 1 << 30; + for (int i = 0; i < n; i++) { + bfs(s, i, d, q); + int ecc = 0; + for (int j = 0; j < n; j++) if (d[j] > ecc) ecc = d[j]; + if (ecc < bestecc) { bestecc = ecc; best = i; } + } + free(d); free(q); + return best; + } + free(d); free(q); + } + for (int i = 1; i < n; i++) + if (s->off[i + 1] - s->off[i] > s->off[best + 1] - s->off[best]) best = i; + return best; +} + +/* Tidy tree from root: leaves at consecutive x in DFS order, each parent + * centred over its first and last child; y = -depth. */ +static int layout_tree(const Sub* s, int root, double* xy) { + int n = s->n; + int* parent = malloc(sizeof(int) * (size_t)n); + int* depth = malloc(sizeof(int) * (size_t)n); + int* pre = malloc(sizeof(int) * (size_t)n); + int* stack = malloc(sizeof(int) * (size_t)(n + 1)); + int* first = malloc(sizeof(int) * (size_t)n); + int* last = malloc(sizeof(int) * (size_t)n); + int ok = parent && depth && pre && stack && first && last; + if (ok) { + for (int i = 0; i < n; i++) { parent[i] = -2; first[i] = last[i] = -1; } + int sp = 0, np = 0; + stack[sp++] = root; parent[root] = -1; depth[root] = 0; + while (sp > 0) { + int u = stack[--sp]; + pre[np++] = u; + /* push children in reverse adjacency order so they pop in order */ + for (int a = s->off[u + 1] - 1; a >= s->off[u]; a--) { + int w = s->adj[a]; + if (parent[w] != -2) continue; + parent[w] = u; depth[w] = depth[u] + 1; + stack[sp++] = w; + } + } + /* first/last child in preorder */ + for (int t = 0; t < np; t++) { + int u = pre[t], p = parent[u]; + if (p < 0) continue; + if (first[p] < 0) first[p] = u; + last[p] = u; + } + double leaf = 0; + for (int t = 0; t < np; t++) { + int u = pre[t]; + if (first[u] < 0) { xy[2 * u] = leaf; leaf += 1.0; } + xy[2 * u + 1] = -(double)depth[u]; + } + for (int t = np - 1; t >= 0; t--) { + int u = pre[t]; + if (first[u] >= 0) xy[2 * u] = 0.5 * (xy[2 * first[u]] + xy[2 * last[u]]); + } + } + free(parent); free(depth); free(pre); free(stack); free(first); free(last); + return ok; +} + +static int cmp_key_idx(const void* a, const void* b) { + const double* x = (const double*)a; const double* y = (const double*)b; + if (x[0] < y[0]) return -1; + if (x[0] > y[0]) return 1; + return (x[1] < y[1]) ? -1 : (x[1] > y[1]); +} + +/* Crossings between edges joining consecutive layers (O(E^2) per layer pair, + * called only for modest edge counts). */ +static long count_crossings(const Sub* s, const int* layer, const int* pos) { + long c = 0; + for (int a = 0; a < s->m; a++) { + int u1 = s->eu[a], v1 = s->ev[a]; + if (layer[u1] > layer[v1]) { int t = u1; u1 = v1; v1 = t; } + if (layer[v1] - layer[u1] != 1) continue; + for (int b = a + 1; b < s->m; b++) { + int u2 = s->eu[b], v2 = s->ev[b]; + if (layer[u2] > layer[v2]) { int t = u2; u2 = v2; v2 = t; } + if (layer[u2] != layer[u1] || layer[v2] != layer[v1]) continue; + if ((long)(pos[u1] - pos[u2]) * (pos[v1] - pos[v2]) < 0) c++; + } + } + return c; +} + +/* Generic layered drawing given layer[]. */ +static int layout_layers(const Sub* s, const int* layer, double* xy) { + int n = s->n, nl = 0; + for (int i = 0; i < n; i++) if (layer[i] + 1 > nl) nl = layer[i] + 1; + int* lcount = calloc((size_t)nl + 1, sizeof(int)); + int* lstart = calloc((size_t)nl + 1, sizeof(int)); + int* order = malloc(sizeof(int) * (size_t)n); /* vertices by layer, ranked */ + int* pos = malloc(sizeof(int) * (size_t)n); + int* bestpos = malloc(sizeof(int) * (size_t)n); + double* key = malloc(sizeof(double) * 2 * (size_t)n); + double* y = malloc(sizeof(double) * (size_t)n); + double* sum = malloc(sizeof(double) * (size_t)n); + int* st = malloc(sizeof(int) * (size_t)n); + int* cnt = malloc(sizeof(int) * (size_t)n); + int ok = lcount && lstart && order && pos && bestpos && key && y && sum && st && cnt; + if (ok) { + for (int i = 0; i < n; i++) lcount[layer[i]]++; + for (int l = 0; l < nl; l++) lstart[l + 1] = lstart[l] + lcount[l]; + { int* fill = calloc((size_t)nl, sizeof(int)); + if (!fill) ok = 0; + else { + for (int i = 0; i < n; i++) { + int l = layer[i]; + order[lstart[l] + fill[l]] = i; pos[i] = fill[l]++; + } + free(fill); + } } + } + if (ok) { + int track = s->m <= 3000; + long best = track ? count_crossings(s, layer, pos) : 0; + memcpy(bestpos, pos, sizeof(int) * (size_t)n); + for (int sweep = 0; sweep < 12 && (!track || best > 0); sweep++) { + int down = (sweep % 2 == 0); + for (int li = 1; li < nl; li++) { + int l = down ? li : nl - 1 - li; + int cntl = lcount[l]; + for (int r = 0; r < cntl; r++) { + int v = order[lstart[l] + r]; + double acc = 0; int na = 0; + for (int a = s->off[v]; a < s->off[v + 1]; a++) { + int w = s->adj[a]; + if (down ? layer[w] < l : layer[w] > l) { acc += pos[w]; na++; } + } + key[2 * r] = na ? acc / na : (double)pos[v]; + key[2 * r + 1] = (double)pos[v]; + (void)v; + } + /* sort (key, current pos); recover vertices via old ranks */ + for (int r = 0; r < cntl; r++) st[r] = order[lstart[l] + r]; + for (int r = 0; r < cntl; r++) key[2 * r + 1] = (double)r; + qsort(key, (size_t)cntl, 2 * sizeof(double), cmp_key_idx); + for (int r = 0; r < cntl; r++) { + int v = st[(int)key[2 * r + 1]]; + order[lstart[l] + r] = v; pos[v] = r; + } + } + if (track) { + long c = count_crossings(s, layer, pos); + if (c < best) { best = c; memcpy(bestpos, pos, sizeof(int) * (size_t)n); } + } + } + if (track) { + memcpy(pos, bestpos, sizeof(int) * (size_t)n); + for (int i = 0; i < n; i++) order[lstart[layer[i]] + pos[i]] = i; + } + /* x coordinates: start at centred ranks, then pull each layer towards + * its neighbours' barycentres subject to unit separation. */ + for (int i = 0; i < n; i++) xy[2 * i] = pos[i] - 0.5 * (lcount[layer[i]] - 1); + for (int pass = 0; pass < 8; pass++) { + int down = (pass % 2 == 0); + for (int li = 0; li < nl; li++) { + int l = down ? li : nl - 1 - li; + int cntl = lcount[l]; + for (int r = 0; r < cntl; r++) { + int v = order[lstart[l] + r]; + double acc = 0; int na = 0; + for (int a = s->off[v]; a < s->off[v + 1]; a++) { + int w = s->adj[a]; + if (layer[w] != l) { acc += xy[2 * w]; na++; } + } + double want = na ? acc / na : xy[2 * v]; + y[r] = want - r; + } + pava(cntl, y, st, sum, cnt); + for (int r = 0; r < cntl; r++) xy[2 * order[lstart[l] + r]] = y[r] + r; + } + } + /* Unit separation within a layer reads cramped against unit layer + * spacing once arrowheads are drawn: widen the layers a little. */ + for (int i = 0; i < n; i++) { xy[2 * i] *= 1.3; xy[2 * i + 1] = -(double)layer[i]; } + } + free(lcount); free(lstart); free(order); free(pos); free(bestpos); + free(key); free(y); free(sum); free(st); free(cnt); + return ok; +} + +/* Sugiyama-style: an edge spanning several layers is routed through one + * dummy vertex per intermediate layer, so crossing reduction and x placement + * see it, and the straight edge then clears the vertices in between. */ +static int layout_layers_dummy(const Sub* s, const int* layer, double* xy) { + long extra = 0; + for (int k = 0; k < s->m; k++) { + int sp = abs(layer[s->eu[k]] - layer[s->ev[k]]); + if (sp > 1) extra += sp - 1; + } + if (extra == 0 || extra > 20000) return layout_layers(s, layer, xy); + Sub a; memset(&a, 0, sizeof(a)); + a.n = s->n + (int)extra; a.m = s->m + (int)extra; + a.eu = malloc(sizeof(int) * (size_t)a.m); + a.ev = malloc(sizeof(int) * (size_t)a.m); + a.dir = calloc((size_t)a.m, 1); + int* L = malloc(sizeof(int) * (size_t)a.n); + double* axy = malloc(sizeof(double) * 2 * (size_t)a.n); + int ok = a.eu && a.ev && a.dir && L && axy; + if (ok) { + memcpy(L, layer, sizeof(int) * (size_t)s->n); + int nv = s->n, ne = 0; + for (int k = 0; k < s->m; k++) { + int u = s->eu[k], v = s->ev[k]; + if (layer[u] > layer[v]) { int t = u; u = v; v = t; } + int prev = u; + for (int l = layer[u] + 1; l < layer[v]; l++) { + L[nv] = l; + a.eu[ne] = prev; a.ev[ne] = nv; ne++; + prev = nv++; + } + a.eu[ne] = prev; a.ev[ne] = v; ne++; + } + ok = sub_csr(&a) && layout_layers(&a, L, axy); + if (ok) memcpy(xy, axy, sizeof(double) * 2 * (size_t)s->n); + } + sub_free(&a); free(L); free(axy); + return ok; +} + +static int layout_layered(const Sub* s, double* xy) { + int n = s->n; + if (n == 1) { xy[0] = xy[1] = 0; return 1; } + int* layer = malloc(sizeof(int) * (size_t)n); + if (!layer) return 0; + int ok; + int dag = dag_layers(s, layer); + if (s->m == n - 1) { + /* A tree: tidy drawing. An arborescence hangs from its source. */ + int root = -1; + if (dag) { + int nsrc = 0; + for (int i = 0; i < n; i++) if (layer[i] == 0) { nsrc++; root = i; } + if (nsrc != 1) root = -1; + } + if (root < 0) root = layout_root(s); + ok = layout_tree(s, root, xy); + } else { + if (!dag) { + int* q = malloc(sizeof(int) * (size_t)n); + if (!q) { free(layer); return 0; } + bfs(s, layout_root(s), layer, q); + free(q); + } + ok = layout_layers_dummy(s, layer, xy); + } + if (ok) { + /* Wide, shallow drawings (big trees) get taller layer spacing so the + * picture is not a flat strip: height at least 0.35 x width. */ + double x0 = 1e300, x1 = -1e300, y0 = 1e300, y1 = -1e300; + for (int i = 0; i < n; i++) { + if (xy[2 * i] < x0) x0 = xy[2 * i]; + if (xy[2 * i] > x1) x1 = xy[2 * i]; + if (xy[2 * i + 1] < y0) y0 = xy[2 * i + 1]; + if (xy[2 * i + 1] > y1) y1 = xy[2 * i + 1]; + } + if (y1 - y0 > 0 && 0.35 * (x1 - x0) > (y1 - y0)) { + double f = 0.35 * (x1 - x0) / (y1 - y0); + for (int i = 0; i < n; i++) xy[2 * i + 1] *= f; + } + } + free(layer); + return ok; +} + +/* ------------------------------------------------------ whole-graph forms */ + +static void layout_grid(int n, double* xy) { + int cols = (int)ceil(sqrt((double)n)); + if (cols < 1) cols = 1; + for (int i = 0; i < n; i++) { xy[2 * i] = i % cols; xy[2 * i + 1] = -(double)(i / cols); } +} + +/* Two columns, parts ordered by barycentre. 0 if s is not bipartite. */ +static int layout_bipartite(const Sub* s, double* xy) { + int n = s->n; + int* col = malloc(sizeof(int) * (size_t)n); + int* q = malloc(sizeof(int) * (size_t)n); + int* layer = calloc((size_t)n + 1, sizeof(int)); + int ok = col && q && layer; + if (ok) { + for (int i = 0; i < n; i++) col[i] = -1; + for (int r = 0; r < n && ok; r++) { + if (col[r] >= 0) continue; + int h = 0, t = 0; + col[r] = 0; q[t++] = r; + while (h < t && ok) { + int u = q[h++]; + for (int a = s->off[u]; a < s->off[u + 1]; a++) { + int w = s->adj[a]; + if (col[w] < 0) { col[w] = 1 - col[u]; q[t++] = w; } + else if (col[w] == col[u]) { ok = 0; break; } + } + } + } + } + if (ok) { + /* Reuse the layered machinery with two layers, then turn it sideways. */ + for (int i = 0; i < n; i++) layer[i] = col[i]; + ok = layout_layers(s, layer, xy); + if (ok) { + int nl = 0, nr = 0; + for (int i = 0; i < n; i++) { if (col[i]) nr++; else nl++; } + double sep = 0.45 * (nl > nr ? nl : nr); + if (sep < 1.5) sep = 1.5; + for (int i = 0; i < n; i++) { + double x = xy[2 * i]; + xy[2 * i] = col[i] ? sep : 0.0; + xy[2 * i + 1] = -x; + } + } + } + free(col); free(q); free(layer); + return ok; +} + +/* -------------------------------------------------------------- packing --- */ + +typedef struct { int comp; int size; int minv; } CompKey; + +static int cmp_comp(const void* a, const void* b) { + const CompKey* x = (const CompKey*)a; const CompKey* y = (const CompKey*)b; + if (x->size != y->size) return y->size - x->size; + return x->minv - y->minv; +} + +/* ----------------------------------------------------------- dispatcher --- */ + +int glayout_parse_method(const Expr* v, GLMethod* out) { + if (!v) return 0; + if (v->type == EXPR_SYMBOL && v->data.symbol.name == SYM_Automatic) { + *out = GL_AUTOMATIC; return 1; + } + if (v->type != EXPR_STRING) return 0; + const char* s = v->data.string; + static const struct { const char* name; GLMethod m; } T[] = { + {"CircularEmbedding", GL_CIRCULAR}, + {"SpringElectricalEmbedding", GL_SPRING}, + {"SpringEmbedding", GL_SPRING}, + {"StressEmbedding", GL_STRESS}, + {"LayeredEmbedding", GL_LAYERED}, + {"LayeredDigraphEmbedding", GL_LAYERED}, + {"TreeEmbedding", GL_LAYERED}, + {"BipartiteEmbedding", GL_BIPARTITE}, + {"GridEmbedding", GL_GRID}, + }; + for (size_t i = 0; i < sizeof(T) / sizeof(T[0]); i++) + if (strcmp(s, T[i].name) == 0) { *out = T[i].m; return 1; } + return 0; +} + +/* Automatic choice: forests with a branch vertex -> layered (tidy trees); + * all-directed DAGs -> layered; everything else (including paths) -> stress. */ +static GLMethod choose_method(const Sub* g) { + if (g->m == 0) return GL_STRESS; + int anydir = 0, alldir = 1; + for (int k = 0; k < g->m; k++) { if (g->dir[k]) anydir = 1; else alldir = 0; } + if (!anydir) { + /* forest <=> m == n - #components (no parallel edges in a Graph) */ + int* c = malloc(sizeof(int) * (size_t)g->n); + int* q = malloc(sizeof(int) * (size_t)g->n); + GLMethod res = GL_STRESS; + if (c && q) { + int nc = 0, maxdeg = 0; + for (int i = 0; i < g->n; i++) c[i] = -1; + for (int r = 0; r < g->n; r++) { + if (c[r] >= 0) continue; + nc++; + int h = 0, t = 0; + c[r] = 1; q[t++] = r; + while (h < t) { + int u = q[h++]; + for (int a = g->off[u]; a < g->off[u + 1]; a++) + if (c[g->adj[a]] < 0) { c[g->adj[a]] = 1; q[t++] = g->adj[a]; } + } + } + for (int i = 0; i < g->n; i++) + if (g->off[i + 1] - g->off[i] > maxdeg) maxdeg = g->off[i + 1] - g->off[i]; + if (g->m == g->n - nc && maxdeg >= 3) res = GL_LAYERED; + } + free(c); free(q); + return res; + } + if (alldir) { + int* layer = malloc(sizeof(int) * (size_t)g->n); + int dag = layer ? dag_layers(g, layer) : 0; + free(layer); + if (dag) return GL_LAYERED; + } + return GL_STRESS; +} + +int glayout_compute(int n, int m, const int* eu, const int* ev, + const unsigned char* dir, GLMethod method, double* xy) { + if (n <= 0) return method == GL_AUTOMATIC ? GL_STRESS : (int)method; + Sub g; memset(&g, 0, sizeof(g)); + g.n = n; g.m = m; + g.eu = malloc(sizeof(int) * (size_t)(m + 1)); + g.ev = malloc(sizeof(int) * (size_t)(m + 1)); + g.dir = calloc((size_t)m + 1, 1); + if (!g.eu || !g.ev || !g.dir) { sub_free(&g); return -1; } + for (int k = 0; k < m; k++) { + g.eu[k] = eu[k]; g.ev[k] = ev[k]; g.dir[k] = dir ? dir[k] : 0; + } + if (!sub_csr(&g)) { sub_free(&g); return -1; } + + if (method == GL_AUTOMATIC) method = choose_method(&g); + + /* Whole-graph embeddings. */ + if (method == GL_CIRCULAR) { layout_circle(n, NULL, xy); sub_free(&g); return method; } + if (method == GL_GRID) { layout_grid(n, xy); sub_free(&g); return method; } + if (method == GL_BIPARTITE) { + if (layout_bipartite(&g, xy)) { sub_free(&g); return method; } + method = GL_STRESS; + } + + /* Per-component embeddings. */ + int* comp = malloc(sizeof(int) * (size_t)n); + int* loc = malloc(sizeof(int) * (size_t)n); + int* q = malloc(sizeof(int) * (size_t)n); + int* members = malloc(sizeof(int) * (size_t)n); + int* cstart = calloc((size_t)n + 1, sizeof(int)); + double* cxy = malloc(sizeof(double) * 2 * (size_t)n); + double* bbox = malloc(sizeof(double) * 4 * (size_t)n); /* per component */ + CompKey* keys = malloc(sizeof(CompKey) * (size_t)n); + int ok = comp && loc && q && members && cstart && cxy && bbox && keys; + int ncomp = 0; + if (ok) { + for (int i = 0; i < n; i++) comp[i] = -1; + for (int r = 0; r < n; r++) { + if (comp[r] >= 0) continue; + int h = 0, t = 0; + comp[r] = ncomp; q[t++] = r; + while (h < t) { + int u = q[h++]; + for (int a = g.off[u]; a < g.off[u + 1]; a++) + if (comp[g.adj[a]] < 0) { comp[g.adj[a]] = ncomp; q[t++] = g.adj[a]; } + } + ncomp++; + } + /* members grouped by component, increasing vertex index within */ + for (int i = 0; i < n; i++) cstart[comp[i] + 1]++; + for (int c = 0; c < ncomp; c++) cstart[c + 1] += cstart[c]; + { int* fill = malloc(sizeof(int) * (size_t)(ncomp + 1)); + if (!fill) ok = 0; + else { + memcpy(fill, cstart, sizeof(int) * (size_t)ncomp); + for (int i = 0; i < n; i++) { loc[i] = fill[comp[i]] - cstart[comp[i]]; members[fill[comp[i]]++] = i; } + free(fill); + } } + } + /* edges grouped by component */ + int* eorder = ok ? malloc(sizeof(int) * (size_t)(m + 1)) : NULL; + int* estart = ok ? calloc((size_t)ncomp + 1, sizeof(int)) : NULL; + if (ok && (!eorder || !estart)) ok = 0; + if (ok) { + for (int k = 0; k < m; k++) estart[comp[eu[k]] + 1]++; + for (int c = 0; c < ncomp; c++) estart[c + 1] += estart[c]; + int* fill = malloc(sizeof(int) * (size_t)(ncomp + 1)); + if (!fill) ok = 0; + else { + memcpy(fill, estart, sizeof(int) * (size_t)ncomp); + for (int k = 0; k < m; k++) eorder[fill[comp[eu[k]]]++] = k; + free(fill); + } + } + for (int c = 0; ok && c < ncomp; c++) { + int nc = cstart[c + 1] - cstart[c]; + int mc = estart[c + 1] - estart[c]; + double* L = cxy + 2 * (size_t)cstart[c]; + if (nc == 1) { L[0] = L[1] = 0; } + else { + Sub s; memset(&s, 0, sizeof(s)); + s.n = nc; s.m = mc; + s.eu = malloc(sizeof(int) * (size_t)(mc + 1)); + s.ev = malloc(sizeof(int) * (size_t)(mc + 1)); + s.dir = calloc((size_t)mc + 1, 1); + if (!s.eu || !s.ev || !s.dir) { sub_free(&s); ok = 0; break; } + for (int t = 0; t < mc; t++) { + int k = eorder[estart[c] + t]; + s.eu[t] = loc[eu[k]]; s.ev[t] = loc[ev[k]]; s.dir[t] = g.dir[k]; + } + if (!sub_csr(&s)) { sub_free(&s); ok = 0; break; } + int r; + if (method == GL_LAYERED) r = layout_layered(&s, L); + else if (method == GL_SPRING) r = layout_spring(&s, L); + else r = layout_stress(&s, L); + if (r && method != GL_LAYERED) { canonical_orientation(&s, L); straighten(&s, L); } + sub_free(&s); + if (!r) { ok = 0; break; } + } + double x0 = 1e300, x1 = -1e300, y0 = 1e300, y1 = -1e300; + for (int i = 0; i < nc; i++) { + if (L[2 * i] < x0) x0 = L[2 * i]; + if (L[2 * i] > x1) x1 = L[2 * i]; + if (L[2 * i + 1] < y0) y0 = L[2 * i + 1]; + if (L[2 * i + 1] > y1) y1 = L[2 * i + 1]; + } + bbox[4 * c] = x0; bbox[4 * c + 1] = x1; bbox[4 * c + 2] = y0; bbox[4 * c + 3] = y1; + keys[c].comp = c; keys[c].size = nc; keys[c].minv = members[cstart[c]]; + } + if (ok) { + /* Shelf packing, largest component first, rows of bounded width. + * Isolated vertices are gathered into one square block at the end + * rather than strung out one per slot. */ + const double gap = 1.0; + qsort(keys, (size_t)ncomp, sizeof(CompKey), cmp_comp); + int nsingle = 0; + for (int t = 0; t < ncomp; t++) if (keys[t].size == 1) nsingle++; + int nitems = ncomp - nsingle + (nsingle > 1 ? 1 : nsingle); + int bcols = (int)ceil(sqrt((double)(nsingle > 0 ? nsingle : 1))); + if (nsingle > 1) { + /* block: singletons (in key order) on a grid; its bbox is the grid */ + int first = ncomp - nsingle; + for (int t = first; t < ncomp; t++) { + int c = keys[t].comp, r = t - first; + double* L = cxy + 2 * (size_t)cstart[c]; + L[0] = r % bcols; L[1] = -(double)(r / bcols); + } + } + /* item u: a component (keys[u].comp) or, last, the singleton block */ + double* iw = malloc(sizeof(double) * (size_t)(nitems + 1)); + double* ih = malloc(sizeof(double) * (size_t)(nitems + 1)); + double* ix0 = malloc(sizeof(double) * (size_t)(nitems + 1)); + double* iy1 = malloc(sizeof(double) * (size_t)(nitems + 1)); + double* ity = malloc(sizeof(double) * (size_t)(nitems + 1)); + double* itx = malloc(sizeof(double) * (size_t)(nitems + 1)); + if (!iw || !ih || !ix0 || !iy1 || !ity || !itx) ok = 0; + int blockitem = (nsingle > 1) ? nitems - 1 : -1; + for (int u = 0; ok && u < nitems; u++) { + if (u == blockitem) { + int rows = (nsingle + bcols - 1) / bcols; + iw[u] = bcols - 1; ih[u] = rows - 1; ix0[u] = 0; iy1[u] = 0; + } else { + int c = keys[u].comp; + iw[u] = bbox[4 * c + 1] - bbox[4 * c]; ih[u] = bbox[4 * c + 3] - bbox[4 * c + 2]; + ix0[u] = bbox[4 * c]; iy1[u] = bbox[4 * c + 3]; + } + } + if (ok) { + double area = 0, maxw = 0; + for (int u = 0; u < nitems; u++) { + area += (iw[u] + gap) * (ih[u] + gap); + if (iw[u] + gap > maxw) maxw = iw[u] + gap; + } + double rowmax = sqrt(area) * 1.4; + if (rowmax < maxw) rowmax = maxw; + double cx = 0, cy = 0, rowh = 0; + int top_align = (method == GL_LAYERED); + int rowbeg = 0; + for (int u = 0; u <= nitems; u++) { + if (u == nitems || (cx > 0 && cx + iw[u] > rowmax)) { + for (int t = rowbeg; t < u; t++) /* centre in the row */ + if (!top_align) ity[t] -= (rowh - ih[t]) / 2.0; + cy -= rowh + gap; cx = 0; rowh = 0; rowbeg = u; + if (u == nitems) break; + } + itx[u] = cx - ix0[u]; ity[u] = cy - iy1[u]; + cx += iw[u] + gap; + if (ih[u] > rowh) rowh = ih[u]; + } + for (int t = 0; t < ncomp; t++) { + int c = keys[t].comp; + int u = (blockitem >= 0 && t >= ncomp - nsingle) ? blockitem : t; + double* L = cxy + 2 * (size_t)cstart[c]; + int nc = cstart[c + 1] - cstart[c]; + for (int q2 = 0; q2 < nc; q2++) { L[2 * q2] += itx[u]; L[2 * q2 + 1] += ity[u]; } + } + } + free(iw); free(ih); free(ix0); free(iy1); free(ity); free(itx); + for (int c = 0; c < ncomp; c++) + for (int t = cstart[c]; t < cstart[c + 1]; t++) { + int v = members[t]; + xy[2 * v] = cxy[2 * t]; xy[2 * v + 1] = cxy[2 * t + 1]; + } + } + free(comp); free(loc); free(q); free(members); free(cstart); free(cxy); + free(bbox); free(keys); free(eorder); free(estart); + sub_free(&g); + return ok ? (int)method : -1; +} diff --git a/src/graph/glayout.h b/src/graph/glayout.h new file mode 100644 index 00000000..79c80e0c --- /dev/null +++ b/src/graph/glayout.h @@ -0,0 +1,149 @@ +#ifndef GLAYOUT_H +#define GLAYOUT_H + +/* glayout.h - graph layout engine and shared drawing helpers for GraphPlot and + * HypergraphPlot (src/graph/glayout.c, graphplot.c, hyp_plot.c). + * + * The layout half is pure numerics over integer edge arrays: it knows nothing + * about Expr, so it can be driven from a Graph (GraphPlot), from the star + * expansion of a Hypergraph (HypergraphPlot), or from a unit test. Every + * algorithm is deterministic -- no random numbers, fixed iteration orders -- + * so the same graph always yields bit-identical coordinates. + * + * Coordinates come out in "edge units": a typical edge has length about 1. + */ + +#include "expr.h" + +typedef enum { + GL_AUTOMATIC = 0, /* trees/DAGs -> layered, paths/cycles/rest -> stress */ + GL_CIRCULAR, /* "CircularEmbedding": VertexList order on a circle */ + GL_SPRING, /* "SpringElectricalEmbedding": Hu's spring-electrical */ + GL_STRESS, /* "StressEmbedding": SMACOF on BFS distances */ + GL_LAYERED, /* "LayeredEmbedding" / "LayeredDigraphEmbedding" */ + GL_BIPARTITE, /* "BipartiteEmbedding": the two parts in two columns */ + GL_GRID /* "GridEmbedding": VertexList order on a square grid */ +} GLMethod; + +/* Lays out n vertices joined by the m edges eu[k] -- ev[k] (0-based vertex + * indices; dir[k] != 0 marks a directed edge eu -> ev, dir may be NULL for all + * undirected). Writes vertex i's position to xy[2i], xy[2i+1]. Every connected + * component is laid out separately and the components are then packed side by + * side. Returns the method actually used (GL_AUTOMATIC resolved; an inapplicable + * choice, e.g. BipartiteEmbedding of an odd cycle, falls back to GL_STRESS), or + * -1 on allocation failure. n == 0 is a no-op success. */ +int glayout_compute(int n, int m, const int* eu, const int* ev, + const unsigned char* dir, GLMethod method, double* xy); + +/* Parses a GraphLayout option value ("StressEmbedding", Automatic, ...). + * Returns 1 and sets *out on a recognised value, 0 otherwise. */ +int glayout_parse_method(const Expr* v, GLMethod* out); + +/* ---- Shared drawing helpers (graphplot.c) -------------------------------- + * Used by both GraphPlot and HypergraphPlot so the two draw vertices, labels, + * margins and options identically. */ + +/* The value of option `name` (Rule or RuleDelayed with a Symbol lhs) among + * res's arguments from position `first` on, or NULL. Later options do NOT + * override earlier ones: the first occurrence wins, as in Mathematica. */ +const Expr* gd_option(const Expr* res, size_t first, const char* name); + +/* Applies a VertexCoordinates spec -- {{x,y}, ...} in VertexList order, or + * {v -> {x,y}, ...} for any subset of vertices -- onto xy. vindex(v) maps a + * vertex to its position (or -1). Returns the number of vertices set; a + * malformed spec sets none. */ +int gd_apply_vertex_coordinates(const Expr* spec, int n, double* xy, + int (*vindex)(const void* ctx, const Expr* v), + const void* ctx); + +/* A growable array of owned primitive Exprs. */ +typedef struct { Expr** p; size_t n, cap; int oom; } GDPrims; +void gd_push(GDPrims* P, Expr* e); /* takes ownership (NULL -> oom) */ +void gd_prims_free(GDPrims* P); + +/* Small constructors. */ +Expr* gd_pt(double x, double y); /* {x, y} */ +Expr* gd_rgb(double r, double g, double b); /* RGBColor[r, g, b] */ +Expr* gd_head1(const char* head, Expr* a); /* head[a] */ +Expr* gd_head2(const char* head, Expr* a, Expr* b); /* head[a, b] */ + +/* The style given for vertex/edge by a VertexStyle/EdgeStyle-like spec: a + * List of rules (or one Rule) looked up with `match`, else the spec itself as a + * global style. NULL when no style applies. Borrowed. */ +const Expr* gd_style_for(const Expr* spec, const void* item, + int (*match)(const Expr* key, const void* item)); + +/* Label text for an arbitrary expression: a String's contents, else its + * printed form. Caller frees. */ +char* gd_label_text(const Expr* e); + +/* Label policy decoded from a VertexLabels value. */ +typedef enum { GD_LBL_NONE, GD_LBL_NAME, GD_LBL_RULES } GDLabelMode; +GDLabelMode gd_label_mode(const Expr* v); + +/* The rhs of the first rule `key -> rhs` in the List `rules` whose lhs is + * SameQ key, or NULL. */ +const Expr* gd_rule_lookup(const Expr* rules, const Expr* key); + +/* Frames the picture. On entry bb (xmin, xmax, ymin, ymax, world units) boxes + * everything drawn except labels; on exit it also boxes the labels plus a small + * margin, and *pw x *ph is the page size in points at which the world maps with + * equal x/y scale. texts[i] (NULL = unlabelled; texts may be NULL) labels the + * vertex at xy[2i], xy[2i+1] with disk radius r; the label goes on the side with + * the widest angular gap between the directions ang[aoff[i] .. aoff[i+1]-1] + * (incident edges), preferring the upper right. Label Text primitives are + * appended to L. The page size honours an ImageSize option in res. */ +void gd_frame(double* bb, const Expr* res, size_t first, int n, + const double* xy, double r, char* const* texts, + const int* aoff, const double* ang, GDPrims* L, + double* pw, double* ph); + +/* Assembles Graphics[prims, PlotRange -> bb, AspectRatio -> Automatic, + * Axes -> False, ImageSize -> {pw, ph}, passthrough...]. Any option of res + * (from `first` on) whose name is not in the NULL-terminated `consumed` list + * is passed through to Graphics (PlotLabel, Background, ...). Consumes P. */ +Expr* gd_finish(GDPrims* P, const double* bb, double pw, double ph, + const Expr* res, size_t first, const char* const* consumed); + +/* The Mathematica ColorData[97] palette, cycled by index. */ +void gd_palette(int k, double* r, double* g, double* b); + +/* Vertex disks (default colour, or vstyle[i] when non-NULL; hl[i] marks a + * highlighted vertex, drawn red and 15% larger), each with a thin darker rim. + * plot_w_pt is the plot width in points (Thickness is relative to it). */ +void gd_emit_vertices(GDPrims* P, int n, const double* xy, double r, + const Expr* const* vstyle, const unsigned char* hl, + double plot_w_pt); + +/* Label texts per a VertexLabels value for the vertices of the List verts + * (NULL when labels are off; entries NULL for unlabelled vertices). */ +char** gd_vertex_texts(const Expr* spec, int n, const Expr* verts); +void gd_free_texts(char** t, int n); + +/* Typical spacing of a drawing: the median edge length, else the median + * nearest-neighbour distance, else 1. */ +double gd_drawing_unit(int n, int m, const int* eu, const int* ev, const double* xy); + +/* Vertex disk radius for a drawing of spacing `unit`; also writes the + * vertices' bounding box (xmin, xmax, ymin, ymax) to bb. With nncap the + * radius is also held under 0.3 x the median nearest-neighbour distance + * (layered drawings, whose leaves sit closer than their edges are long). */ +double gd_vertex_radius(int n, const double* xy, double unit, int nncap, double* bb); + +/* A lone vertex (or a cluster smaller than one unit in both directions) + * would fill the page: widen bb to a unit square about its centre. */ +void gd_min_extent(double* bb, double unit); + +/* expr_eq(key, (const Expr*)item): the vertex matcher for gd_style_for. */ +int gd_vertex_match(const Expr* key, const void* item); + +/* Is item in a GraphHighlight-like spec (a List of items, or one item)? */ +int gd_in_highlight(const Expr* hl, const void* item, + int (*match)(const Expr*, const void*)); + +/* Default vertex colour (Mathematica's first ColorData[97] entry). */ +#define GD_VERTEX_R 0.368417 +#define GD_VERTEX_G 0.506779 +#define GD_VERTEX_B 0.709798 + +#endif /* GLAYOUT_H */ diff --git a/src/graph/graph.c b/src/graph/graph.c index a75c1ba6..90e76201 100644 --- a/src/graph/graph.c +++ b/src/graph/graph.c @@ -296,9 +296,19 @@ void graph_init(void) { symtab_add_builtin("GraphPlot", builtin_graph_plot); symtab_get_def("GraphPlot")->attributes |= ATTR_PROTECTED; symtab_set_docstring("GraphPlot", - "GraphPlot[g] gives a Graphics object drawing the graph g with a " - "circular vertex layout. Vertex labels are off by default; pass " - "VertexLabels -> True to draw them (in black)."); + "GraphPlot[g, opts] gives a Graphics object drawing the graph g (or a " + "list of rules {u -> v, ...}). The default layout is deterministic: tidy " + "layered trees for branching forests, layered drawings for DAGs, stress " + "majorization otherwise, with components packed side by side. Options: " + "GraphLayout -> \"StressEmbedding\" | \"SpringElectricalEmbedding\" | " + "\"CircularEmbedding\" | \"LayeredEmbedding\" | \"BipartiteEmbedding\" | " + "\"GridEmbedding\"; VertexCoordinates -> {{x,y}, ...} or {v -> {x,y}, ...}; " + "VertexLabels -> None | \"Name\" | Automatic | {v -> lbl, ...}; " + "GraphHighlight -> {v, e, ...} (red, thicker); VertexStyle and EdgeStyle -> " + "a colour or {item -> colour, ...}; EdgeLabels -> \"EdgeWeight\" | {e -> lbl}; " + "VertexSize -> d (diameter in edge lengths). Directed edges get arrowheads " + "that stop at the target vertex. Other options (ImageSize, PlotLabel, ...) " + "pass through to Graphics."); /* ---- Editing, transforms, set operations, cycles (gops_*.c) ------- */ graph_ops_init(); diff --git a/src/graph/graph_hyper.h b/src/graph/graph_hyper.h index 2c956065..80d403ca 100644 --- a/src/graph/graph_hyper.h +++ b/src/graph/graph_hyper.h @@ -84,6 +84,7 @@ Expr* builtin_uniform_hypergraph_q(Expr* res); /* UniformHypergraphQ Expr* builtin_hypergraph_dual(Expr* res); /* HypergraphDual */ Expr* builtin_hypergraph_clique_expansion(Expr* res); /* HypergraphCliqueExpansion */ Expr* builtin_hypergraph_star_expansion(Expr* res); /* HypergraphStarExpansion */ +Expr* builtin_hypergraph_plot(Expr* res); /* HypergraphPlot (hyp_plot.c) */ Expr* builtin_hypergraph_to_graph(Expr* res); /* HypergraphToGraph */ Expr* builtin_hypergraph_line_graph(Expr* res); /* HypergraphLineGraph */ Expr* builtin_hypergraph_connected_components(Expr* res); /* HypergraphConnectedComponents */ diff --git a/src/graph/graphplot.c b/src/graph/graphplot.c index e3778ef1..6615dd97 100644 --- a/src/graph/graphplot.c +++ b/src/graph/graphplot.c @@ -1,19 +1,35 @@ -/* graphplot.c - GraphPlot[g]: render a graph as a Graphics[...] expression. +/* graphplot.c - GraphPlot[g, opts]: draw a graph as a Graphics[...] expression, + * plus the drawing helpers GraphPlot shares with HypergraphPlot (glayout.h). * - * Emits the same primitives the plotting engine uses (Line, Disk, Text), so the - * existing renderer draws it with no renderer changes (and the text placeholder - * is used when USE_GRAPHICS=0). Vertices are laid out on a circle (MVP layout; - * a force-directed spring layout is the documented future hook). Each edge is a - * Line between its endpoints, each vertex a Disk plus a Text label. + * PIPELINE + * 1. Layout (glayout.c): GraphLayout -> Automatic (tidy trees for branching + * forests, layered DAGs, stress majorization otherwise), or an explicit + * embedding; VertexCoordinates overrides any subset of the positions. + * 2. Size: vertex disks get a radius scaled to the median edge length and the + * layout extent, so a 10-vertex and a 1000-vertex graph both read. + * 3. Frame (gd_frame): labels are placed beside their vertex on the side + * with the widest angular gap between incident edges, and the world box + * and page size are solved together so labels are never clipped. + * 4. Emit: edges (Line, or Arrow shortened to stop at the target disk, with + * Arrowheads sized to the disks; mutual pairs u->v, v->u offset apart), + * then vertex disks with a thin darker rim, then labels. * - * Directed edges are drawn as plain lines in the MVP (no arrowheads yet). + * The output uses only primitives every renderer draws -- Line, Arrow, Disk, + * Circle, Text, colour and Thickness/Arrowheads directives -- with + * AspectRatio -> Automatic, Axes -> False and an explicit PlotRange/ImageSize. + * Colours: vertices RGBColor[0.368417, 0.506779, 0.709798] (ColorData[97][1]), + * edges a medium grey-blue, GraphHighlight in red and thicker, as Mathematica. * * Memory (SPEC section 4): returns a freshly-allocated Graphics tree; the - * evaluator frees res. + * evaluator frees res. Nothing borrowed from res outlives the call. */ #include "graph.h" +#include "glayout.h" +#include "graphics_export.h" #include "expr.h" +#include "eval.h" +#include "print.h" #include "sym_names.h" #include #include @@ -24,97 +40,762 @@ #define M_PI 3.14159265358979323846 #endif -#define NODE_RADIUS 0.02 +/* Page geometry of the PDF exporter with Axes -> False (see + * graphics_export.c): 8pt left/bottom and 12pt right/top margins. */ +#define PAGE_MARGIN_X 20.0 +#define PAGE_MARGIN_Y 20.0 +#define LABEL_FS 10.0 /* exporter's default Text size, in points */ +#define EDGE_R 0.571589 /* default edge colour: a medium grey-blue */ +#define EDGE_G 0.586483 +#define EDGE_B 0.699215 -static Expr* point2(double x, double y) { - Expr* xy[2] = { expr_new_real(x), expr_new_real(y) }; - return expr_new_function(expr_new_symbol(SYM_List), xy, 2); +/* =========================================================== shared helpers */ + +static bool is_head(const Expr* e, const char* sym) { + return e && e->type == EXPR_FUNCTION && e->data.function.head + && e->data.function.head->type == EXPR_SYMBOL + && e->data.function.head->data.symbol.name == sym; +} + +static bool is_rule(const Expr* e) { + return (is_head(e, SYM_Rule) || is_head(e, SYM_RuleDelayed)) + && e->data.function.arg_count == 2; } -/* Is a VertexLabels option value "on"? True / All / any string (e.g. "Name", - * "Index") enable labels; None / False (and anything else) leave them off. */ -static bool vertex_labels_on(const Expr* rhs) { - if (!rhs) return false; - if (rhs->type == EXPR_STRING) return true; - if (rhs->type == EXPR_SYMBOL) { - const char* n = rhs->data.symbol.name; - return n == SYM_True || n == SYM_All; +static bool num(const Expr* e, double* out) { + if (!e) return false; + if (e->type == EXPR_INTEGER) { *out = (double)e->data.integer; return true; } + if (e->type == EXPR_REAL) { *out = e->data.real; return true; } + if (is_head(e, SYM_Rational) && e->data.function.arg_count == 2) { + double a, b; + if (num(e->data.function.args[0], &a) && num(e->data.function.args[1], &b) && b != 0) { + *out = a / b; return true; + } } return false; } +static bool pair(const Expr* e, double* x, double* y) { + return is_head(e, SYM_List) && e->data.function.arg_count == 2 + && num(e->data.function.args[0], x) && num(e->data.function.args[1], y); +} + +const Expr* gd_option(const Expr* res, size_t first, const char* name) { + for (size_t i = first; i < res->data.function.arg_count; i++) { + const Expr* a = res->data.function.args[i]; + if (!is_rule(a)) continue; + const Expr* lhs = a->data.function.args[0]; + if (lhs->type == EXPR_SYMBOL && strcmp(lhs->data.symbol.name, name) == 0) + return a->data.function.args[1]; + } + return NULL; +} + +int gd_apply_vertex_coordinates(const Expr* spec, int n, double* xy, + int (*vindex)(const void* ctx, const Expr* v), + const void* ctx) { + if (!is_head(spec, SYM_List)) return 0; + size_t k = spec->data.function.arg_count; + if (k == 0) return 0; + double x, y; + if (!is_rule(spec->data.function.args[0])) { + /* {{x,y}, ...}: exactly one pair per vertex, all numeric. */ + if ((int)k != n) return 0; + for (size_t i = 0; i < k; i++) + if (!pair(spec->data.function.args[i], &x, &y)) return 0; + for (size_t i = 0; i < k; i++) { + pair(spec->data.function.args[i], &x, &y); + xy[2 * i] = x; xy[2 * i + 1] = y; + } + return n; + } + int set = 0; + for (size_t i = 0; i < k; i++) { + const Expr* r = spec->data.function.args[i]; + if (!is_rule(r) || !pair(r->data.function.args[1], &x, &y)) continue; + int v = vindex(ctx, r->data.function.args[0]); + if (v < 0 || v >= n) continue; + xy[2 * v] = x; xy[2 * v + 1] = y; set++; + } + return set; +} + +void gd_push(GDPrims* P, Expr* e) { + if (!e) { P->oom = 1; return; } + if (P->n == P->cap) { + size_t nc = P->cap ? 2 * P->cap : 64; + Expr** np = realloc(P->p, nc * sizeof(Expr*)); + if (!np) { expr_free(e); P->oom = 1; return; } + P->p = np; P->cap = nc; + } + P->p[P->n++] = e; +} + +void gd_prims_free(GDPrims* P) { + for (size_t i = 0; i < P->n; i++) expr_free(P->p[i]); + free(P->p); + P->p = NULL; P->n = P->cap = 0; +} + +/* Rounds to 6 significant decimals so the Graphics tree prints compactly and + * identically across runs. */ +static double tidy(double v) { + if (v == 0 || !isfinite(v)) return 0; + double a = fabs(v); + double sc = pow(10.0, 5 - (int)floor(log10(a))); + double r = floor(v * sc + 0.5) / sc; + return r == 0 ? 0 : r; +} + +Expr* gd_pt(double x, double y) { + Expr* a[2] = { expr_new_real(tidy(x)), expr_new_real(tidy(y)) }; + return expr_new_function(expr_new_symbol(SYM_List), a, 2); +} + +Expr* gd_rgb(double r, double g, double b) { + Expr* a[3] = { expr_new_real(tidy(r)), expr_new_real(tidy(g)), expr_new_real(tidy(b)) }; + return expr_new_function(expr_new_symbol(SYM_RGBColor), a, 3); +} + +Expr* gd_head1(const char* head, Expr* a) { + Expr* v[1] = { a }; + return expr_new_function(expr_new_symbol(head), v, 1); +} + +Expr* gd_head2(const char* head, Expr* a, Expr* b) { + Expr* v[2] = { a, b }; + return expr_new_function(expr_new_symbol(head), v, 2); +} + +const Expr* gd_style_for(const Expr* spec, const void* item, + int (*match)(const Expr* key, const void* item)) { + if (!spec) return NULL; + if (is_rule(spec)) + return match(spec->data.function.args[0], item) ? spec->data.function.args[1] : NULL; + if (is_head(spec, SYM_List)) { + for (size_t i = 0; i < spec->data.function.arg_count; i++) { + const Expr* r = spec->data.function.args[i]; + if (is_rule(r) && match(r->data.function.args[0], item)) + return r->data.function.args[1]; + } + /* A List with no rule-like element is a list of directives: a global + * style. Anything else (a rule list that did not match, or a + * malformed one) gives no style. */ + for (size_t i = 0; i < spec->data.function.arg_count; i++) { + const Expr* a = spec->data.function.args[i]; + if (is_rule(a) || is_head(a, SYM_TwoWayRule) || is_head(a, SYM_UndirectedEdge) + || is_head(a, SYM_DirectedEdge)) + return NULL; + } + return spec; + } + if (spec->type == EXPR_SYMBOL + && (spec->data.symbol.name == SYM_Automatic || spec->data.symbol.name == SYM_None)) + return NULL; + return spec; +} + +char* gd_label_text(const Expr* e) { + if (e->type == EXPR_STRING) { + size_t n = strlen(e->data.string) + 1; + char* s = malloc(n); + if (s) memcpy(s, e->data.string, n); + return s; + } + return expr_to_string((Expr*)e); +} + +GDLabelMode gd_label_mode(const Expr* v) { + if (!v) return GD_LBL_NONE; + if (v->type == EXPR_SYMBOL) { + const char* n = v->data.symbol.name; + if (n == SYM_True || n == SYM_All || n == SYM_Automatic) return GD_LBL_NAME; + return GD_LBL_NONE; + } + if (v->type == EXPR_STRING) return GD_LBL_NAME; /* "Name", "Index"... */ + if (is_rule(v) || is_head(v, SYM_List)) return GD_LBL_RULES; + return GD_LBL_NONE; +} + +const Expr* gd_rule_lookup(const Expr* rules, const Expr* key) { + if (is_rule(rules)) + return expr_eq(rules->data.function.args[0], key) ? rules->data.function.args[1] : NULL; + if (!is_head(rules, SYM_List)) return NULL; + for (size_t i = 0; i < rules->data.function.arg_count; i++) { + const Expr* r = rules->data.function.args[i]; + if (is_rule(r) && expr_eq(r->data.function.args[0], key)) + return r->data.function.args[1]; + } + return NULL; +} + +void gd_palette(int k, double* r, double* g, double* b) { + static const double P[10][3] = { + {0.368417, 0.506779, 0.709798}, {0.880722, 0.611041, 0.142051}, + {0.560181, 0.691569, 0.194885}, {0.922526, 0.385626, 0.209179}, + {0.528488, 0.470624, 0.701351}, {0.772079, 0.431554, 0.102387}, + {0.363898, 0.618501, 0.782349}, {1.0, 0.75, 0.0}, + {0.647624, 0.37816, 0.614037}, {0.571589, 0.586483, 0.699215}, + }; + int i = ((k % 10) + 10) % 10; + *r = P[i][0]; *g = P[i][1]; *b = P[i][2]; +} + +/* ---- label placement and framing ---------------------------------------- */ + +/* Candidate label directions in order of preference (upper right first). */ +static const double LBL_DIRS[8] = { + M_PI / 4, 3 * M_PI / 4, -M_PI / 4, -3 * M_PI / 4, M_PI / 2, 0.0, M_PI, -M_PI / 2 +}; + +static double ang_dist(double a, double b) { + double d = fmod(fabs(a - b), 2 * M_PI); + return d > M_PI ? 2 * M_PI - d : d; +} + +/* Placement of one label: anchor point, Text offset, and world box. */ +typedef struct { double ax, ay, ox, oy, box[4]; } LblPlace; + +static void place_label(const char* text, double x, double y, double r, + const double* avoid, int na, double pt2w, LblPlace* L) { + int best = 0; + if (na > 0) { + double bestscore = -1; + for (int c = 0; c < 8; c++) { + double sc = 1e9; + for (int a = 0; a < na; a++) { + double d = ang_dist(LBL_DIRS[c], avoid[a]); + if (d < sc) sc = d; + } + if (sc > bestscore + 1e-6) { bestscore = sc; best = c; } + } + } + double dx = cos(LBL_DIRS[best]), dy = sin(LBL_DIRS[best]); + if (fabs(dx) < 1e-12) dx = 0; + if (fabs(dy) < 1e-12) dy = 0; + double gap = r + 1.5 * pt2w; + L->ax = x + gap * dx; L->ay = y + gap * dy; + double mx = fabs(dx) > fabs(dy) ? fabs(dx) : fabs(dy); + L->ox = -dx / mx; L->oy = -dy / mx; + double tw = graphics_helvetica_width(text) * LABEL_FS * pt2w; + double th = 0.70 * LABEL_FS * pt2w; + /* box centre = anchor - offset * half size */ + double cx = L->ax - L->ox * tw / 2, cy = L->ay - L->oy * th / 2; + L->box[0] = cx - tw / 2; L->box[1] = cx + tw / 2; + L->box[2] = cy - th / 2 - 0.25 * th; /* allow for descenders */ + L->box[3] = cy + th / 2; +} + +/* Page size and points-per-world-unit for world box W x H. */ +static void page_for(double W, double H, double iw, double ih, double D, + double* pw, double* ph, double* k) { + if (W < 1e-9) W = 1e-9; + if (H < 1e-9) H = 1e-9; + if (iw > 0 && ih > 0) { + double kx = (iw - PAGE_MARGIN_X) / W, ky = (ih - PAGE_MARGIN_Y) / H; + *k = kx < ky ? kx : ky; *pw = iw; *ph = ih; + } else if (iw > 0) { + *k = (iw - PAGE_MARGIN_X) / W; *pw = iw; *ph = H * *k + PAGE_MARGIN_Y; + } else if (W >= H) { + *pw = D; *k = (D - PAGE_MARGIN_X) / W; *ph = H * *k + PAGE_MARGIN_Y; + if (*ph < 90) { *ph = 90; } + } else { + *ph = D; *k = (D - PAGE_MARGIN_Y) / H; *pw = W * *k + PAGE_MARGIN_X; + if (*pw < 90) { *pw = 90; } + } +} + +void gd_frame(double* bb, const Expr* res, size_t first, int n, + const double* xy, double r, char* const* texts, + const int* aoff, const double* ang, GDPrims* L, + double* pw, double* ph) { + double iw = 0, ih = 0, v, w2; + const Expr* is = gd_option(res, first, "ImageSize"); + if (is && num(is, &v) && v > 0) iw = v; + else if (is && pair(is, &v, &w2) && v > 0 && w2 > 0) { iw = v; ih = w2; } + /* Default size: 300pt, growing gently with the vertex count. */ + double D = 300.0 * sqrt(n > 40 ? n / 40.0 : 1.0); + if (D > 600) D = 600; + + double base[4] = { bb[0], bb[1], bb[2], bb[3] }; + if (!(base[0] <= base[1])) { base[0] = -1; base[1] = 1; } + if (!(base[2] <= base[3])) { base[2] = -1; base[3] = 1; } + double cur[4]; + memcpy(cur, base, sizeof(cur)); + double k = 1, pwv = 0, phv = 0; + LblPlace lp; + for (int iter = 0; iter < 6; iter++) { + page_for(cur[1] - cur[0], cur[3] - cur[2], iw, ih, D, &pwv, &phv, &k); + double pt2w = 1.0 / k; + double nb[4]; + memcpy(nb, base, sizeof(nb)); + for (int i = 0; texts && i < n; i++) { + if (!texts[i]) continue; + place_label(texts[i], xy[2 * i], xy[2 * i + 1], r, + ang ? ang + aoff[i] : NULL, ang ? aoff[i + 1] - aoff[i] : 0, pt2w, &lp); + if (lp.box[0] < nb[0]) nb[0] = lp.box[0]; + if (lp.box[1] > nb[1]) nb[1] = lp.box[1]; + if (lp.box[2] < nb[2]) nb[2] = lp.box[2]; + if (lp.box[3] > nb[3]) nb[3] = lp.box[3]; + } + /* Margin: 4pt plus 2% of the extent. */ + double ext = (nb[1] - nb[0]) > (nb[3] - nb[2]) ? nb[1] - nb[0] : nb[3] - nb[2]; + double mg = 4.0 * pt2w + 0.02 * ext; + if (ext <= 1e-12) mg = 1.0; + nb[0] -= mg; nb[1] += mg; nb[2] -= mg; nb[3] += mg; + memcpy(cur, nb, sizeof(cur)); + } + page_for(cur[1] - cur[0], cur[3] - cur[2], iw, ih, D, &pwv, &phv, &k); + /* Emit the labels at the final scale. */ + if (texts) { + int any = 0; + for (int i = 0; i < n; i++) { + if (!texts[i]) continue; + if (!any) { gd_push(L, gd_head1(SYM_GrayLevel, expr_new_real(0.0))); any = 1; } + place_label(texts[i], xy[2 * i], xy[2 * i + 1], r, + ang ? ang + aoff[i] : NULL, ang ? aoff[i + 1] - aoff[i] : 0, 1.0 / k, &lp); + Expr* ta[3] = { expr_new_string(texts[i]), gd_pt(lp.ax, lp.ay), gd_pt(lp.ox, lp.oy) }; + gd_push(L, expr_new_function(expr_new_symbol(SYM_Text), ta, 3)); + } + } + memcpy(bb, cur, sizeof(cur)); + *pw = floor(pwv + 0.5); *ph = floor(phv + 0.5); +} + +Expr* gd_finish(GDPrims* P, const double* bb, double pw, double ph, + const Expr* res, size_t first, const char* const* consumed) { + Expr* list = expr_new_function(expr_new_symbol(SYM_List), P->p, P->n); + free(P->p); P->p = NULL; P->n = P->cap = 0; + size_t argc = res->data.function.arg_count; + Expr** ga = malloc(sizeof(Expr*) * (argc + 6)); + if (!ga) { expr_free(list); return NULL; } + size_t k = 0; + ga[k++] = list; + /* Pass-through options first: the exporter honours the first occurrence, + * so a user's PlotRange/AspectRatio/... beats the computed one. */ + for (size_t i = first; i < argc; i++) { + const Expr* a = res->data.function.args[i]; + if (!is_rule(a) || a->data.function.args[0]->type != EXPR_SYMBOL) continue; + const char* nm = a->data.function.args[0]->data.symbol.name; + bool skip = strcmp(nm, "ImageSize") == 0; + for (const char* const* c = consumed; !skip && c && *c; c++) + if (strcmp(nm, *c) == 0) skip = true; + if (!skip) ga[k++] = expr_copy((Expr*)a); + } + ga[k++] = gd_head2(SYM_Rule, expr_new_symbol(SYM_PlotRange), + gd_head2(SYM_List, gd_head2(SYM_List, expr_new_real(tidy(bb[0])), expr_new_real(tidy(bb[1]))), + gd_head2(SYM_List, expr_new_real(tidy(bb[2])), expr_new_real(tidy(bb[3]))))); + ga[k++] = gd_head2(SYM_Rule, expr_new_symbol(SYM_AspectRatio), expr_new_symbol(SYM_Automatic)); + ga[k++] = gd_head2(SYM_Rule, expr_new_symbol(SYM_Axes), expr_new_symbol(SYM_False)); + ga[k++] = gd_head2(SYM_Rule, expr_new_symbol(SYM_ImageSize), + gd_head2(SYM_List, expr_new_integer((int64_t)pw), expr_new_integer((int64_t)ph))); + Expr* g = expr_new_function(expr_new_symbol(SYM_Graphics), ga, k); + free(ga); + return g; +} + +void gd_min_extent(double* bb, double unit) { + if (bb[1] - bb[0] >= unit || bb[3] - bb[2] >= unit) return; + double cx = (bb[0] + bb[1]) / 2, cy = (bb[2] + bb[3]) / 2; + bb[0] = cx - unit / 2; bb[1] = cx + unit / 2; + bb[2] = cy - unit / 2; bb[3] = cy + unit / 2; +} + +/* ================================================================ GraphPlot */ + +typedef struct { const Expr* g; } VCtx; + +static int graph_vindex(const void* ctx, const Expr* v) { + return graph_vertex_position(((const VCtx*)ctx)->g, v); +} + +int gd_vertex_match(const Expr* key, const void* item) { + return expr_eq(key, (const Expr*)item); +} + +/* An edge as seen by the style / highlight matchers. */ +typedef struct { const Expr* u; const Expr* v; int directed; } EdgeRef; + +/* Does the edge spec `key` (UndirectedEdge/TwoWayRule/DirectedEdge/Rule) + * denote edge e? An undirected spec matches an undirected edge either way + * round; a directed spec matches that directed edge, and (leniently) an + * undirected edge in either orientation. */ +static int edge_match(const Expr* key, const void* item) { + const EdgeRef* e = (const EdgeRef*)item; + if (!key || key->type != EXPR_FUNCTION || key->data.function.arg_count != 2) return 0; + bool und = is_head(key, SYM_UndirectedEdge) || is_head(key, SYM_TwoWayRule); + bool dir = is_head(key, SYM_DirectedEdge) || is_head(key, SYM_Rule); + if (!und && !dir) return 0; + const Expr* a = key->data.function.args[0]; + const Expr* b = key->data.function.args[1]; + bool fwd = expr_eq(a, e->u) && expr_eq(b, e->v); + bool rev = expr_eq(a, e->v) && expr_eq(b, e->u); + if (e->directed) return dir ? fwd : 0; + return fwd || rev; +} + +/* Is `item` in the GraphHighlight spec (a List of items, or one item)? */ +int gd_in_highlight(const Expr* hl, const void* item, + int (*match)(const Expr*, const void*)) { + if (!hl) return false; + if (is_head(hl, SYM_List)) { + for (size_t i = 0; i < hl->data.function.arg_count; i++) + if (match(hl->data.function.args[i], item)) return true; + return false; + } + return match(hl, item) != 0; +} + +/* RGB of a colour directive, if it is one. */ +static bool rgb_of(const Expr* c, double* r, double* g, double* b) { + if (is_head(c, SYM_RGBColor) && c->data.function.arg_count >= 3) + return num(c->data.function.args[0], r) && num(c->data.function.args[1], g) + && num(c->data.function.args[2], b); + if (is_head(c, SYM_GrayLevel) && c->data.function.arg_count >= 1 && num(c->data.function.args[0], r)) { + *g = *b = *r; return true; + } + return false; +} + +/* Pushes `style` (a directive, or a List of directives) unless it equals the + * current one. */ +static void push_style(GDPrims* P, const Expr* style, const Expr** cur) { + if (*cur && expr_eq(*cur, style)) return; + if (is_head(style, SYM_List)) + for (size_t i = 0; i < style->data.function.arg_count; i++) + gd_push(P, expr_copy(style->data.function.args[i])); + else gd_push(P, expr_copy((Expr*)style)); + *cur = style; +} + +static int cmp_double(const void* a, const void* b) { + double x = *(const double*)a, y = *(const double*)b; + return (x > y) - (x < y); +} + +/* Draws n vertices at xy with radius r: disk in its style, thin darker rim. + * vstyle[i] may be NULL (default colour); hl[i] marks highlighted vertices. */ +void gd_emit_vertices(GDPrims* P, int n, const double* xy, double r, + const Expr* const* vstyle, const unsigned char* hl, + double pt2page) { + const Expr* cur = NULL; + Expr* defc = gd_rgb(GD_VERTEX_R, GD_VERTEX_G, GD_VERTEX_B); + Expr* red = gd_rgb(1, 0, 0); + double rim = 0.6 / pt2page; /* 0.6pt rim as a Thickness */ + gd_push(P, gd_head1(SYM_Thickness, expr_new_real(tidy(rim)))); + for (int i = 0; i < n; i++) { + const Expr* st = (hl && hl[i]) ? red : (vstyle && vstyle[i]) ? vstyle[i] : defc; + double rr = (hl && hl[i]) ? r * 1.15 : r; + push_style(P, st, &cur); + gd_push(P, gd_head2(SYM_Disk, gd_pt(xy[2 * i], xy[2 * i + 1]), expr_new_real(tidy(rr)))); + double cr, cg, cb; + const Expr* colour = st; + if (is_head(st, SYM_List)) { + colour = NULL; + for (size_t t = 0; t < st->data.function.arg_count; t++) + if (rgb_of(st->data.function.args[t], &cr, &cg, &cb)) colour = st->data.function.args[t]; + } + Expr* rimc = (colour && rgb_of(colour, &cr, &cg, &cb)) + ? gd_rgb(cr * 0.6, cg * 0.6, cb * 0.6) : gd_head1(SYM_GrayLevel, expr_new_real(0.3)); + gd_push(P, rimc); + gd_push(P, gd_head2(SYM_Circle, gd_pt(xy[2 * i], xy[2 * i + 1]), expr_new_real(tidy(rr)))); + cur = NULL; /* the rim colour replaced it */ + } + expr_free(defc); expr_free(red); +} + +/* Typical spacing of a drawing: median edge length, else median nearest- + * neighbour distance, else 1. */ +/* Median nearest-neighbour distance (0 when unknown / too many vertices). */ +static double nn_median(int n, const double* xy) { + if (n < 2 || n > 3000) return 0; + double* L = malloc(sizeof(double) * (size_t)n); + if (!L) return 0; + int t = 0; + for (int i = 0; i < n; i++) { + double best = 1e300; + for (int j = 0; j < n; j++) { + if (j == i) continue; + double dx = xy[2 * i] - xy[2 * j], dy = xy[2 * i + 1] - xy[2 * j + 1]; + double d = dx * dx + dy * dy; + if (d > 1e-24 && d < best) best = d; + } + if (best < 1e300) L[t++] = sqrt(best); + } + double med = 0; + if (t > 0) { qsort(L, (size_t)t, sizeof(double), cmp_double); med = L[t / 2]; } + free(L); + return med; +} + +double gd_drawing_unit(int n, int m, const int* eu, const int* ev, const double* xy) { + double unit = 0; + if (m > 0) { + double* L = malloc(sizeof(double) * (size_t)m); + if (L) { + int t = 0; + for (int k = 0; k < m; k++) { + double dx = xy[2 * eu[k]] - xy[2 * ev[k]], dy = xy[2 * eu[k] + 1] - xy[2 * ev[k] + 1]; + double d = sqrt(dx * dx + dy * dy); + if (d > 1e-12) L[t++] = d; + } + if (t > 0) { qsort(L, (size_t)t, sizeof(double), cmp_double); unit = L[t / 2]; } + free(L); + } + } + if (unit <= 0) unit = nn_median(n, xy); + return unit > 0 ? unit : 1.0; +} + +/* Vertex radius from the unit spacing and the layout extent. */ +double gd_vertex_radius(int n, const double* xy, double unit, int nncap, double* bb) { + bb[0] = bb[2] = 1e300; bb[1] = bb[3] = -1e300; + for (int i = 0; i < n; i++) { + if (xy[2 * i] < bb[0]) bb[0] = xy[2 * i]; + if (xy[2 * i] > bb[1]) bb[1] = xy[2 * i]; + if (xy[2 * i + 1] < bb[2]) bb[2] = xy[2 * i + 1]; + if (xy[2 * i + 1] > bb[3]) bb[3] = xy[2 * i + 1]; + } + double ext = n ? ((bb[1] - bb[0]) > (bb[3] - bb[2]) ? bb[1] - bb[0] : bb[3] - bb[2]) : 0; + if (ext < unit) ext = unit; + double r = 0.075 * unit; + if (r > 0.02 * ext) r = 0.02 * ext; + if (r < 0.004 * ext) r = 0.004 * ext; + /* Vertices may sit closer than edges are long (the leaves of a wide + * tree): with nncap, never let typical neighbours' disks touch. A + * hairball's disks may overlap, as they do in any force layout. */ + double nn = nncap ? nn_median(n, xy) : 0; + if (nn > 0 && r > 0.3 * nn) r = 0.3 * nn; + return r; +} + +/* Builds the per-vertex label texts per VertexLabels (NULL when none). */ +char** gd_vertex_texts(const Expr* spec, int n, const Expr* verts) { + GDLabelMode mode = gd_label_mode(spec); + if (mode == GD_LBL_NONE || n == 0) return NULL; + char** t = calloc((size_t)n, sizeof(char*)); + if (!t) return NULL; + for (int i = 0; i < n; i++) { + const Expr* v = verts->data.function.args[i]; + if (mode == GD_LBL_NAME) { t[i] = gd_label_text(v); continue; } + const Expr* l = gd_rule_lookup(spec, v); + if (!l) continue; + if (l->type == EXPR_STRING && strcmp(l->data.string, "Name") == 0) t[i] = gd_label_text(v); + else if (!(l->type == EXPR_SYMBOL && l->data.symbol.name == SYM_None)) t[i] = gd_label_text(l); + } + return t; +} + +void gd_free_texts(char** t, int n) { + if (!t) return; + for (int i = 0; i < n; i++) free(t[i]); + free(t); +} + +/* Options GraphPlot interprets itself (everything else passes to Graphics). */ +static const char* const GP_CONSUMED[] = { + "GraphLayout", "VertexCoordinates", "VertexLabels", "GraphHighlight", + "VertexStyle", "EdgeStyle", "EdgeLabels", "VertexSize", "Method", NULL +}; + Expr* builtin_graph_plot(Expr* res) { size_t argc = res->data.function.arg_count; if (argc < 1) return NULL; const Expr* g = res->data.function.args[0]; - if (!graph_is_valid(g)) return NULL; - - /* Options (after the graph). Vertex labels are OFF by default, matching - * Mathematica — GraphPlot draws no vertex names unless asked. Pass - * VertexLabels -> True | All | "Name" to draw them. */ - bool show_labels = false; - for (size_t i = 1; i < argc; i++) { - const Expr* a = res->data.function.args[i]; - if (a->type == EXPR_FUNCTION && a->data.function.arg_count == 2 - && a->data.function.head - && a->data.function.head->type == EXPR_SYMBOL - && a->data.function.head->data.symbol.name == SYM_Rule) { - const Expr* lhs = a->data.function.args[0]; - if (lhs->type == EXPR_SYMBOL - && strcmp(lhs->data.symbol.name, "VertexLabels") == 0) - show_labels = vertex_labels_on(a->data.function.args[1]); - } + Expr* built = NULL; + if (is_head(g, SYM_List)) { + /* GraphPlot[{u -> v, ...}]: the Graph of those rules. */ + Expr* call = gd_head1(SYM_Graph, expr_copy((Expr*)g)); + built = evaluate(call); + expr_free(call); + g = built; } + if (!graph_is_valid(g)) { if (built) expr_free(built); return NULL; } + /* Options must be Rules. */ + for (size_t i = 1; i < argc; i++) + if (!is_rule(res->data.function.args[i])) { if (built) expr_free(built); return NULL; } const Expr* verts = g->data.function.args[0]; const Expr* edges = g->data.function.args[1]; int n = (int)verts->data.function.arg_count; - size_t ne = edges->data.function.arg_count; + int m = (int)edges->data.function.arg_count; - /* Circular layout coordinates. */ - double* x = (n > 0) ? calloc((size_t)n, sizeof(double)) : NULL; - double* y = (n > 0) ? calloc((size_t)n, sizeof(double)) : NULL; - for (int i = 0; i < n; i++) { - if (n == 1) { x[i] = 0.0; y[i] = 0.0; } - else { - double t = 2.0 * M_PI * (double)i / (double)n; - x[i] = cos(t); y[i] = sin(t); + const int *beu = NULL, *bev = NULL; const unsigned char* bdir = NULL; + graph_edge_indices(g, &beu, &bev, &bdir); + int* eu = malloc(sizeof(int) * (size_t)(m + 1)); + int* ev = malloc(sizeof(int) * (size_t)(m + 1)); + unsigned char* dir = malloc((size_t)m + 1); + double* xy = calloc(2 * (size_t)n + 2, sizeof(double)); + if (!eu || !ev || !dir || !xy) { + free(eu); free(ev); free(dir); free(xy); if (built) expr_free(built); return NULL; + } + for (int k = 0; k < m; k++) { eu[k] = beu[k]; ev[k] = bev[k]; dir[k] = bdir[k]; } + + /* ---- layout -------------------------------------------------------- */ + GLMethod method = GL_AUTOMATIC; + const Expr* lo = gd_option(res, 1, "GraphLayout"); + if (lo && !glayout_parse_method(lo, &method)) method = GL_AUTOMATIC; + const Expr* vc = gd_option(res, 1, "VertexCoordinates"); + VCtx vctx = { g }; + int given = 0, used = GL_AUTOMATIC; + if (vc) { + /* Probe: a complete spec skips the layout entirely. */ + given = gd_apply_vertex_coordinates(vc, n, xy, graph_vindex, &vctx); + } + if (given < n) { + used = glayout_compute(n, m, eu, ev, dir, method, xy); + if (used < 0) { + free(eu); free(ev); free(dir); free(xy); if (built) expr_free(built); return NULL; } + if (vc) gd_apply_vertex_coordinates(vc, n, xy, graph_vindex, &vctx); } - /* Primitives: one Line per edge, then one Disk per vertex, then (only when - * VertexLabels is on) a black colour directive + one Text label per vertex. */ - size_t total = ne + (size_t)n + (show_labels ? (size_t)n + 1 : 0); - Expr** prims = (total > 0) ? calloc(total, sizeof(Expr*)) : NULL; - size_t p = 0; + /* ---- sizes ---------------------------------------------------------- */ + double unit = gd_drawing_unit(n, m, eu, ev, xy); + double bb[4]; + double r = gd_vertex_radius(n, xy, unit, used == GL_LAYERED || used == GL_GRID, bb); + const Expr* vs_opt = gd_option(res, 1, "VertexSize"); + double vsz; + if (vs_opt && num(vs_opt, &vsz) && vsz > 0) r = vsz * unit / 2; /* diameter, in units */ + if (n > 0) { bb[0] -= 1.15 * r; bb[1] += 1.15 * r; bb[2] -= 1.15 * r; bb[3] += 1.15 * r; } + else { bb[0] = bb[2] = -1; bb[1] = bb[3] = 1; } + gd_min_extent(bb, unit); - for (size_t k = 0; k < ne; k++) { - const Expr* e = edges->data.function.args[k]; - int iu = graph_vertex_index(verts, e->data.function.args[0]); - int iv = graph_vertex_index(verts, e->data.function.args[1]); - Expr* pts[2] = { point2(x[iu], y[iu]), point2(x[iv], y[iv]) }; - Expr* seg = expr_new_function(expr_new_symbol(SYM_List), pts, 2); - Expr* la[1] = { seg }; - prims[p++] = expr_new_function(expr_new_symbol(SYM_Line), la, 1); + /* ---- labels & frame ------------------------------------------------- */ + char** texts = gd_vertex_texts(gd_option(res, 1, "VertexLabels"), n, verts); + int* aoff = NULL; double* ang = NULL; + if (texts) { + aoff = calloc((size_t)n + 1, sizeof(int)); + ang = malloc(sizeof(double) * (size_t)(2 * m + 1)); + if (aoff && ang) { + for (int k = 0; k < m; k++) { aoff[eu[k] + 1]++; aoff[ev[k] + 1]++; } + for (int i = 0; i < n; i++) aoff[i + 1] += aoff[i]; + int* fill = malloc(sizeof(int) * (size_t)(n + 1)); + if (fill) { + memcpy(fill, aoff, sizeof(int) * (size_t)n); + for (int k = 0; k < m; k++) { + double dx = xy[2 * ev[k]] - xy[2 * eu[k]], dy = xy[2 * ev[k] + 1] - xy[2 * eu[k] + 1]; + ang[fill[eu[k]]++] = atan2(dy, dx); + ang[fill[ev[k]]++] = atan2(-dy, -dx); + } + free(fill); + } else { free(aoff); free(ang); aoff = NULL; ang = NULL; } + } else { free(aoff); free(ang); aoff = NULL; ang = NULL; } } - for (int i = 0; i < n; i++) { - Expr* da[2] = { point2(x[i], y[i]), expr_new_real(NODE_RADIUS) }; - prims[p++] = expr_new_function(expr_new_symbol(SYM_Disk), da, 2); - } - if (show_labels) { - /* Labels are black by default so they read against the vertex disks - * (without this they inherit the disk colour and vanish). */ - Expr* bl[1] = { expr_new_real(0.0) }; - prims[p++] = expr_new_function(expr_new_symbol(SYM_GrayLevel), bl, 1); - for (int i = 0; i < n; i++) { - Expr* ta[2] = { expr_copy(verts->data.function.args[i]), point2(x[i], y[i]) }; - prims[p++] = expr_new_function(expr_new_symbol(SYM_Text), ta, 2); + GDPrims labels = {0}; + double pw, ph; + gd_frame(bb, res, 1, n, xy, r, texts, aoff, ang, &labels, &pw, &ph); + double wpt = (pw - PAGE_MARGIN_X) / (bb[1] - bb[0]); /* points per unit */ + double hpt = (ph - PAGE_MARGIN_Y) / (bb[3] - bb[2]); + if (hpt < wpt) wpt = hpt; + double plot_w_pt = (bb[1] - bb[0]) * wpt; /* exporter's plot_w */ + + /* ---- edges ---------------------------------------------------------- */ + GDPrims P = {0}; + const Expr* hl = gd_option(res, 1, "GraphHighlight"); + const Expr* es = gd_option(res, 1, "EdgeStyle"); + const Expr* vsopt = gd_option(res, 1, "VertexStyle"); + Expr* edge_def = gd_rgb(EDGE_R, EDGE_G, EDGE_B); + Expr* red = gd_rgb(1, 0, 0); + double thin = 1.1 / plot_w_pt, thick = 2.8 / plot_w_pt; + double head_len = 2.2 * r; /* world units */ + if (head_len * wpt < 5.0) head_len = 5.0 / wpt; /* >= 5pt */ + if (head_len * wpt > 12.0) head_len = 12.0 / wpt; /* <= 12pt */ + unsigned char* ehl = calloc((size_t)m + 1, 1); + for (int pass = 0; pass < 2 && ehl; pass++) { + /* pass 0: ordinary edges; pass 1: highlighted edges, drawn on top */ + const Expr* cur = NULL; + int started = 0; + for (int k = 0; k < m; k++) { + const Expr* e = edges->data.function.args[k]; + EdgeRef ref = { e->data.function.args[0], e->data.function.args[1], dir[k] }; + if (pass == 0) ehl[k] = gd_in_highlight(hl, &ref, edge_match); + if ((int)ehl[k] != pass) continue; + if (!started) { + gd_push(&P, gd_head1(SYM_Thickness, expr_new_real(tidy(pass ? thick : thin)))); + gd_push(&P, gd_head1(SYM_Arrowheads, + expr_new_real(tidy(head_len * (pass ? 1.25 : 1.0) * wpt / plot_w_pt)))); + started = 1; + } + const Expr* st = pass ? red : gd_style_for(es, &ref, edge_match); + push_style(&P, st ? st : edge_def, &cur); + double x0 = xy[2 * eu[k]], y0 = xy[2 * eu[k] + 1]; + double x1 = xy[2 * ev[k]], y1 = xy[2 * ev[k] + 1]; + if (!dir[k]) { + gd_push(&P, gd_head1(SYM_Line, gd_head2(SYM_List, gd_pt(x0, y0), gd_pt(x1, y1)))); + continue; + } + double dx = x1 - x0, dy = y1 - y0, len = sqrt(dx * dx + dy * dy); + if (len < 1e-12) continue; + double ux = dx / len, uy = dy / len; + /* A mutual pair u->v, v->u: shift each to its own right. */ + if (graph_has_edge(g, e->data.function.args[1], e->data.function.args[0], 1) == 1) { + double off = 0.8 * r; + x0 += uy * off; y0 -= ux * off; x1 += uy * off; y1 -= ux * off; + } + /* Tip stops just short of the target disk (highlighted disks are + * 15% larger, so clear that radius). */ + double re = r * 1.15 + 0.5 / wpt; + gd_push(&P, gd_head1(SYM_Arrow, gd_head2(SYM_List, + gd_pt(x0 + ux * r, y0 + uy * r), gd_pt(x1 - ux * re, y1 - uy * re)))); + } + } + + /* ---- edge labels ---------------------------------------------------- */ + const Expr* el = gd_option(res, 1, "EdgeLabels"); + if (el && m > 0 && !(el->type == EXPR_SYMBOL && (el->data.symbol.name == SYM_None + || el->data.symbol.name == SYM_False))) { + Expr* weights = NULL; + if (el->type == EXPR_STRING && strcmp(el->data.string, "EdgeWeight") == 0) + weights = graph_resolve_edge_weights(g); + double cx = 0, cy = 0; + for (int i = 0; i < n; i++) { cx += xy[2 * i]; cy += xy[2 * i + 1]; } + if (n) { cx /= n; cy /= n; } + int any = 0; + for (int k = 0; k < m; k++) { + const Expr* e = edges->data.function.args[k]; + EdgeRef ref = { e->data.function.args[0], e->data.function.args[1], dir[k] }; + const Expr* lab = NULL; + if (weights && is_head(weights, SYM_List) && (size_t)k < weights->data.function.arg_count) + lab = weights->data.function.args[k]; + else if (!weights) lab = gd_style_for(el, &ref, edge_match); + if (!lab || (el->type == EXPR_STRING && !weights)) continue; + char* s = gd_label_text(lab); + if (!s) continue; + double x0 = xy[2 * eu[k]], y0 = xy[2 * eu[k] + 1]; + double x1 = xy[2 * ev[k]], y1 = xy[2 * ev[k] + 1]; + double mx = (x0 + x1) / 2, my = (y0 + y1) / 2; + double dx = x1 - x0, dy = y1 - y0, len = sqrt(dx * dx + dy * dy); + double px = len > 0 ? -dy / len : 0, py = len > 0 ? dx / len : 1; + if ((mx - cx) * px + (my - cy) * py > 0) { px = -px; py = -py; } /* inward */ + if (!any) { gd_push(&P, gd_head1(SYM_GrayLevel, expr_new_real(0.25))); any = 1; } + double off = 7.0 / wpt; + Expr* ta[3] = { expr_new_string(s), gd_pt(mx + px * off, my + py * off), gd_pt(0, 0) }; + gd_push(&P, expr_new_function(expr_new_symbol(SYM_Text), ta, 3)); + free(s); } + if (weights) expr_free(weights); + } + + /* ---- vertices ------------------------------------------------------- */ + const Expr** vstyle = n ? calloc((size_t)n, sizeof(Expr*)) : NULL; + unsigned char* vhl = n ? calloc((size_t)n, 1) : NULL; + for (int i = 0; i < n && vstyle && vhl; i++) { + const Expr* v = verts->data.function.args[i]; + vstyle[i] = gd_style_for(vsopt, v, gd_vertex_match); + vhl[i] = gd_in_highlight(hl, v, gd_vertex_match); } + if (n) gd_emit_vertices(&P, n, xy, r, vstyle, vhl, plot_w_pt); + for (size_t i = 0; i < labels.n; i++) gd_push(&P, labels.p[i]); + free(labels.p); - free(x); free(y); - Expr* prim_list = expr_new_function(expr_new_symbol(SYM_List), prims, total); - free(prims); - Expr* ga[1] = { prim_list }; - return expr_new_function(expr_new_symbol(SYM_Graphics), ga, 1); + Expr* out = NULL; + if (!P.oom) out = gd_finish(&P, bb, pw, ph, res, 1, GP_CONSUMED); + gd_prims_free(&P); + expr_free(edge_def); expr_free(red); + free((void*)vstyle); free(vhl); free(ehl); + gd_free_texts(texts, n); free(aoff); free(ang); + free(eu); free(ev); free(dir); free(xy); + if (built) expr_free(built); + return out; } diff --git a/src/graph/hyp_init.c b/src/graph/hyp_init.c index 12030ac0..f2cded48 100644 --- a/src/graph/hyp_init.c +++ b/src/graph/hyp_init.c @@ -48,6 +48,13 @@ void graph_hyper_init(void) { "HypergraphStarExpansion[h] gives the incidence (star) expansion of h: the " "bipartite Graph on VertexList[h] and nodes Hyperedge[1], ..., Hyperedge[m], " "with v<->Hyperedge[j] whenever v lies in hyperedge j."); + reg("HypergraphPlot", builtin_hypergraph_plot, + "HypergraphPlot[h, opts] gives a Graphics object drawing the hypergraph h " + "(or a list of hyperedges): vertices placed by the stress layout of the star " + "expansion, each hyperedge of 3 or more vertices a translucent rounded hull " + "in its own palette colour, one of 2 a stadium and one of 1 a circle around " + "its vertex, vertex disks on top. Options: VertexLabels, VertexCoordinates, " + "VertexStyle and GraphLayout as in GraphPlot; others pass through to Graphics."); reg("HypergraphToGraph", builtin_hypergraph_to_graph, "HypergraphToGraph[h] converts h, read as an ordered hypergraph, to the " "directed Graph with v_a->v_b for every a +#include +#include + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +#define HP_FILL_OPACITY 0.22 + +static int is_head(const Expr* e, const char* sym) { + return e && e->type == EXPR_FUNCTION && e->data.function.head + && e->data.function.head->type == EXPR_SYMBOL + && e->data.function.head->data.symbol.name == sym; +} + +static int is_rule(const Expr* e) { + return (is_head(e, SYM_Rule) || is_head(e, SYM_RuleDelayed)) + && e->data.function.arg_count == 2; +} + +static int hyp_vindex(const void* ctx, const Expr* v) { + return graph_vidx_get((const GraphVIdx*)ctx, v); +} + +static int cmp_pt(const void* a, const void* b) { + const double* p = (const double*)a; const double* q = (const double*)b; + if (p[0] != q[0]) return p[0] < q[0] ? -1 : 1; + if (p[1] != q[1]) return p[1] < q[1] ? -1 : 1; + return 0; +} + +static double cross3(const double* o, const double* a, const double* b) { + return (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0]); +} + +/* Convex hull (Andrew's monotone chain) of the k points in pts (sorted in + * place), counter-clockwise, collinear points dropped, into hull (2k + 2 + * doubles). Returns the hull size (1 for coincident points). */ +static int convex_hull(double* pts, int k, double* hull) { + qsort(pts, (size_t)k, 2 * sizeof(double), cmp_pt); + int u = 0; + for (int i = 0; i < k; i++) { /* dedupe */ + if (u > 0 && fabs(pts[2 * i] - pts[2 * (u - 1)]) < 1e-12 + && fabs(pts[2 * i + 1] - pts[2 * (u - 1) + 1]) < 1e-12) continue; + pts[2 * u] = pts[2 * i]; pts[2 * u + 1] = pts[2 * i + 1]; u++; + } + k = u; + if (k <= 2) { memcpy(hull, pts, sizeof(double) * 2 * (size_t)k); return k; } + int h = 0; + for (int i = 0; i < k; i++) { + while (h >= 2 && cross3(hull + 2 * (h - 2), hull + 2 * (h - 1), pts + 2 * i) <= 1e-12) h--; + hull[2 * h] = pts[2 * i]; hull[2 * h + 1] = pts[2 * i + 1]; h++; + } + for (int i = k - 2, lo = h + 1; i >= 0; i--) { + while (h >= lo && cross3(hull + 2 * (h - 2), hull + 2 * (h - 1), pts + 2 * i) <= 1e-12) h--; + hull[2 * h] = pts[2 * i]; hull[2 * h + 1] = pts[2 * i + 1]; h++; + } + return h - 1; /* last == first */ +} + +/* The hull inflated by delta with round corners, as a closed point List; + * extends bb. */ +static Expr* rounded_hull(const double* hull, int h, double delta, double* bb, double* area) { + enum { STEP_DEG = 15 }; + Expr** pts = malloc(sizeof(Expr*) * (size_t)(h * (360 / STEP_DEG + 2) + 4)); + if (!pts) return NULL; + size_t np = 0; + double a2 = 0; + for (int i = 0; i < h; i++) { + const double* p = hull + 2 * i; + double a0, a1; + if (h == 1) { a0 = 0; a1 = 2 * M_PI; } + else { + const double* pp = hull + 2 * ((i + h - 1) % h); + const double* pn = hull + 2 * ((i + 1) % h); + /* outward normals of the incoming and outgoing edges (CCW hull) */ + a0 = atan2(-(p[0] - pp[0]), p[1] - pp[1]); + a1 = atan2(-(pn[0] - p[0]), pn[1] - p[1]); + while (a1 < a0) a1 += 2 * M_PI; + if (h == 2 && a1 - a0 < 1e-9) a1 += M_PI; /* degenerate safety */ + } + int steps = (int)ceil((a1 - a0) / (STEP_DEG * M_PI / 180.0)); + if (steps < 1) steps = 1; + for (int s = 0; s <= steps; s++) { + if (h == 1 && s == steps) break; + double t = a0 + (a1 - a0) * s / steps; + double x = p[0] + delta * cos(t), y = p[1] + delta * sin(t); + pts[np++] = gd_pt(x, y); + if (x < bb[0]) bb[0] = x; + if (x > bb[1]) bb[1] = x; + if (y < bb[2]) bb[2] = y; + if (y > bb[3]) bb[3] = y; + } + } + /* area of the inflated shape: hull area + perimeter*delta + pi delta^2 */ + double per = 0; + for (int i = 0; i < h && h > 1; i++) { + const double* p = hull + 2 * i; const double* q = hull + 2 * ((i + 1) % h); + a2 += p[0] * q[1] - q[0] * p[1]; + per += sqrt((q[0] - p[0]) * (q[0] - p[0]) + (q[1] - p[1]) * (q[1] - p[1])); + } + *area = fabs(a2) / 2 + per * delta + M_PI * delta * delta; + Expr* list = expr_new_function(expr_new_symbol(SYM_List), pts, np); + free(pts); + return list; +} + +typedef struct { int j; double area; } HOrder; + +static int cmp_horder(const void* a, const void* b) { + const HOrder* x = (const HOrder*)a; const HOrder* y = (const HOrder*)b; + if (x->area != y->area) return x->area > y->area ? -1 : 1; + return x->j - y->j; +} + +static const char* const HP_CONSUMED[] = { + "GraphLayout", "VertexCoordinates", "VertexLabels", "VertexStyle", + "VertexSize", "EdgeStyle", NULL +}; + +Expr* builtin_hypergraph_plot(Expr* res) { + size_t argc = res->data.function.arg_count; + if (argc < 1) return NULL; + for (size_t i = 1; i < argc; i++) + if (!is_rule(res->data.function.args[i])) return NULL; + const Expr* h = res->data.function.args[0]; + Expr* built = NULL; + if (is_head(h, SYM_List)) { + Expr* a[1] = { expr_copy((Expr*)h) }; + Expr* call = expr_new_function(expr_new_symbol(hyp_sym_hypergraph()), a, 1); + built = evaluate(call); + expr_free(call); + h = built; + } + HypView V; + if (!hyp_view(h, &V, 0)) { if (built) expr_free(built); return NULL; } + int n = V.n, m = V.m; + + /* ---- star expansion ------------------------------------------------ */ + int* node = malloc(sizeof(int) * (size_t)(m + 1)); /* star node per hyperedge, or -1 */ + int N = n, M = 0; + for (int j = 0; j < m; j++) { + int sz = V.soff[j + 1] - V.soff[j]; + node[j] = sz >= 2 ? N++ : -1; + if (sz >= 2) M += sz; + } + int* eu = malloc(sizeof(int) * (size_t)(M + 1)); + int* ev = malloc(sizeof(int) * (size_t)(M + 1)); + double* xy = calloc(2 * (size_t)N + 2, sizeof(double)); + if (!node || !eu || !ev || !xy) { + free(node); free(eu); free(ev); free(xy); if (built) expr_free(built); return NULL; + } + { int t = 0; + for (int j = 0; j < m; j++) { + if (node[j] < 0) continue; + for (int a = V.soff[j]; a < V.soff[j + 1]; a++) { eu[t] = V.sv[a]; ev[t] = node[j]; t++; } + } } + GLMethod method = GL_STRESS; + const Expr* lo = gd_option(res, 1, "GraphLayout"); + if (lo && (!glayout_parse_method(lo, &method) || method == GL_AUTOMATIC)) method = GL_STRESS; + const Expr* vc = gd_option(res, 1, "VertexCoordinates"); + int given = vc ? gd_apply_vertex_coordinates(vc, n, xy, hyp_vindex, V.ix) : 0; + if (given < n) { + if (glayout_compute(N, M, eu, ev, NULL, method, xy) < 0) { + free(node); free(eu); free(ev); free(xy); if (built) expr_free(built); return NULL; + } + if (vc) gd_apply_vertex_coordinates(vc, n, xy, hyp_vindex, V.ix); + } + + /* ---- sizes: spacing from the vertex positions alone ----------------- */ + double unit = gd_drawing_unit(n, 0, NULL, NULL, xy); + double bb[4]; + double r = gd_vertex_radius(n, xy, unit, 0, bb); + if (n == 0) { bb[0] = bb[2] = -1; bb[1] = bb[3] = 1; } + else { bb[0] -= 1.15 * r; bb[1] += 1.15 * r; bb[2] -= 1.15 * r; bb[3] += 1.15 * r; } + gd_min_extent(bb, 2 * unit); + double delta0 = 0.24 * unit, dstep = 0.10 * unit; + if (delta0 < 2.2 * r) delta0 = 2.2 * r; + + /* Margin levels: process hyperedges by increasing size; each takes one + * more than the most hyperedges already drawn around any of its vertices. */ + int* cnt = calloc((size_t)n + 1, sizeof(int)); + int* level = calloc((size_t)m + 1, sizeof(int)); + HOrder* ord = malloc(sizeof(HOrder) * (size_t)(m + 1)); + double* pts = malloc(sizeof(double) * 2 * (size_t)(n + 1)); + double* hull = malloc(sizeof(double) * (4 * (size_t)n + 8)); /* chain: <= 2k points */ + Expr** shapes = calloc((size_t)m + 1, sizeof(Expr*)); + GDPrims P = {0}; + int ok = cnt && level && ord && pts && hull && shapes; + if (ok) { + for (int j = 0; j < m; j++) { ord[j].j = j; ord[j].area = (double)(V.soff[j + 1] - V.soff[j]); } + /* increasing size, then index */ + for (int a = 1; a < m; a++) { + HOrder t = ord[a]; int b = a - 1; + while (b >= 0 && (ord[b].area > t.area || (ord[b].area == t.area && ord[b].j > t.j))) { + ord[b + 1] = ord[b]; b--; + } + ord[b + 1] = t; + } + for (int t = 0; t < m; t++) { + int j = ord[t].j, lv = 0; + for (int a = V.soff[j]; a < V.soff[j + 1]; a++) if (cnt[V.sv[a]] > lv) lv = cnt[V.sv[a]]; + level[j] = lv; + for (int a = V.soff[j]; a < V.soff[j + 1]; a++) cnt[V.sv[a]]++; + } + /* shapes */ + for (int j = 0; j < m; j++) { + int k = 0; + for (int a = V.soff[j]; a < V.soff[j + 1]; a++) { + pts[2 * k] = xy[2 * V.sv[a]]; pts[2 * k + 1] = xy[2 * V.sv[a] + 1]; k++; + } + ord[j].j = j; ord[j].area = 0; + if (k == 0) continue; /* empty hyperedge */ + int hs = convex_hull(pts, k, hull); + shapes[j] = rounded_hull(hull, hs, delta0 + dstep * level[j], bb, &ord[j].area); + if (!shapes[j]) ok = 0; + } + qsort(ord, (size_t)m, sizeof(HOrder), cmp_horder); + } + if (ok) { + double pw, ph; + /* label avoidance: directions to the centroids of v's hyperedges */ + char** texts = gd_vertex_texts(gd_option(res, 1, "VertexLabels"), n, V.verts); + int* aoff = NULL; double* ang = NULL; + if (texts) { + aoff = calloc((size_t)n + 1, sizeof(int)); + ang = malloc(sizeof(double) * (size_t)(V.soff[m] + 1)); + if (aoff && ang) { + for (int j = 0; j < m; j++) + if (V.soff[j + 1] - V.soff[j] >= 2) + for (int a = V.soff[j]; a < V.soff[j + 1]; a++) aoff[V.sv[a] + 1]++; + for (int i = 0; i < n; i++) aoff[i + 1] += aoff[i]; + int* fill = malloc(sizeof(int) * (size_t)(n + 1)); + if (fill) { + memcpy(fill, aoff, sizeof(int) * (size_t)n); + for (int j = 0; j < m; j++) { + int sz = V.soff[j + 1] - V.soff[j]; + if (sz < 2) continue; + double cx = 0, cy = 0; + for (int a = V.soff[j]; a < V.soff[j + 1]; a++) { cx += xy[2 * V.sv[a]]; cy += xy[2 * V.sv[a] + 1]; } + cx /= sz; cy /= sz; + for (int a = V.soff[j]; a < V.soff[j + 1]; a++) { + int v = V.sv[a]; + ang[fill[v]++] = atan2(cy - xy[2 * v + 1], cx - xy[2 * v]); + } + } + free(fill); + } else { free(aoff); free(ang); aoff = NULL; ang = NULL; } + } else { free(aoff); free(ang); aoff = NULL; ang = NULL; } + } + GDPrims labels = {0}; + gd_frame(bb, res, 1, n, xy, r, texts, aoff, ang, &labels, &pw, &ph); + double wpt = (pw - 20.0) / (bb[1] - bb[0]), hpt = (ph - 20.0) / (bb[3] - bb[2]); + if (hpt < wpt) wpt = hpt; + double plot_w_pt = (bb[1] - bb[0]) * wpt; + + gd_push(&P, gd_head1(SYM_Thickness, expr_new_real(1.0 / plot_w_pt))); + for (int t = 0; t < m; t++) { + int j = ord[t].j; + if (!shapes[j]) continue; + double cr, cg, cb; + gd_palette(j + 1, &cr, &cg, &cb); /* skip the vertex blue */ + gd_push(&P, gd_rgb(cr, cg, cb)); + gd_push(&P, gd_head1(SYM_Opacity, expr_new_real(HP_FILL_OPACITY))); + gd_push(&P, gd_head1(SYM_Polygon, expr_copy(shapes[j]))); + gd_push(&P, gd_head1(SYM_Opacity, expr_new_real(1.0))); + /* outline: the same points, closed */ + Expr* ring = shapes[j]; + size_t k = ring->data.function.arg_count; + Expr** lp = malloc(sizeof(Expr*) * (k + 1)); + if (!lp) { P.oom = 1; break; } + for (size_t q = 0; q < k; q++) lp[q] = expr_copy(ring->data.function.args[q]); + lp[k] = expr_copy(ring->data.function.args[0]); + Expr* line = expr_new_function(expr_new_symbol(SYM_List), lp, k + 1); + free(lp); + gd_push(&P, gd_rgb(cr * 0.8, cg * 0.8, cb * 0.8)); + gd_push(&P, gd_head1(SYM_Line, line)); + } + /* vertices */ + const Expr* vsopt = gd_option(res, 1, "VertexStyle"); + const Expr** vstyle = n ? calloc((size_t)n, sizeof(Expr*)) : NULL; + for (int i = 0; i < n && vstyle; i++) + vstyle[i] = gd_style_for(vsopt, V.verts->data.function.args[i], gd_vertex_match); + Expr* dark = gd_rgb(0.25, 0.30, 0.42); + for (int i = 0; i < n && vstyle; i++) if (!vstyle[i]) vstyle[i] = dark; + if (n) gd_emit_vertices(&P, n, xy, r, vstyle, NULL, plot_w_pt); + expr_free(dark); + free((void*)vstyle); + for (size_t i = 0; i < labels.n; i++) gd_push(&P, labels.p[i]); + free(labels.p); + gd_free_texts(texts, n); free(aoff); free(ang); + Expr* out = P.oom ? NULL : gd_finish(&P, bb, pw, ph, res, 1, HP_CONSUMED); + gd_prims_free(&P); + for (int j = 0; j < m; j++) if (shapes[j]) expr_free(shapes[j]); + free(shapes); free(cnt); free(level); free(ord); free(pts); free(hull); + free(node); free(eu); free(ev); free(xy); + if (built) expr_free(built); + return out; + } + gd_prims_free(&P); + if (shapes) for (int j = 0; j < m; j++) if (shapes[j]) expr_free(shapes[j]); + free(shapes); free(cnt); free(level); free(ord); free(pts); free(hull); + free(node); free(eu); free(ev); free(xy); + if (built) expr_free(built); + return NULL; +} diff --git a/src/graphics/graphics_export.c b/src/graphics/graphics_export.c index 1f102365..a21f9df5 100644 --- a/src/graphics/graphics_export.c +++ b/src/graphics/graphics_export.c @@ -13,10 +13,14 @@ * Point -> filled dot * Rectangle -> filled + stroked box * Arrow -> polyline + solid arrowhead - * Text -> Helvetica string (BT .. Tj .. ET) + * Text -> Helvetica string (BT .. Tj .. ET); Text[s, pos, {ox, oy}] + * aligns as Mathematica does ({-1, 0}: left end at pos), and + * Text[Style[s, n | FontSize -> n | colour ...], ...] sets + * the size and colour of that one string * - * plus the RGBColor/GrayLevel/Hue/CMYKColor/Opacity/Thickness/PointSize - * directives. PDF's coordinate system is y-up with the origin at the lower + * plus the RGBColor/GrayLevel/Hue/CMYKColor/Opacity/Thickness/PointSize/ + * Arrowheads directives. AspectRatio -> Automatic maps world units with equal + * x and y scale (centred on the page) instead of stretching to fill it. PDF's coordinate system is y-up with the origin at the lower * left, which is exactly the mathematical convention, so world y needs no * flip. Colour, opacity and the "nice" tick policy mirror the renderer so the * PDF and the on-screen/PNG output agree. @@ -235,6 +239,8 @@ typedef struct { int axes_set; /* Axes option was present */ double bg_r, bg_g, bg_b; int have_bg; double width, height; /* page size in points (from ImageSize) */ + int equal_aspect; /* AspectRatio -> Automatic: equal x/y scale */ + int size_pair; /* ImageSize -> {w, h} fixed both dimensions */ } Opts; /* Look for option `sym -> value` among the Graphics args (index >= 1). */ @@ -275,6 +281,7 @@ static void parse_range(const Expr* v, Opts* o) { static void parse_opts(const Expr* g, Opts* o) { o->have_range = 0; o->frame = 0; o->axes = 1; o->axes_set = 0; o->have_bg = 0; o->width = 504.0; o->height = 360.0; /* 7in x 5in default */ + o->equal_aspect = 0; o->size_pair = 0; const Expr* v; if ((v = find_option(g, SYM_PlotRange))) parse_range(v, o); @@ -290,9 +297,12 @@ static void parse_opts(const Expr* g, Opts* o) { else if (head_is(v, SYM_List) && v->data.function.arg_count == 2 && to_double(v->data.function.args[0], &w) && to_double(v->data.function.args[1], &h) && w > 0 && h > 0) { - o->width = w; o->height = h; + o->width = w; o->height = h; o->size_pair = 1; } } + if ((v = find_option(g, SYM_AspectRatio)) && v->type == EXPR_SYMBOL + && v->data.symbol.name == SYM_Automatic) + o->equal_aspect = 1; } /* ---------------------------------------------------------- tick policy --- */ @@ -331,6 +341,7 @@ typedef struct { double dash[GFX_MAX_DASH]; int ndash; int dash_emitted; /* a non-solid dash pattern is in force */ + double arrowhead; /* Arrowheads[s]: s (fraction of width), 0 = auto */ /* world->page transform */ double ox, oy, sx, sy; } Emit; @@ -396,6 +407,23 @@ static void emit_pdf_string(Buf* c, const char* s) { buf_cat(c, ")"); } +/* Helvetica advance widths (1/1000 em) for ASCII 32..126, from the AFM. */ +static const short HELV_W[95] = { + 278, 278, 355, 556, 556, 889, 667, 191, 333, 333, 389, 584, 278, 333, 278, 278, + 556, 556, 556, 556, 556, 556, 556, 556, 556, 556, 278, 278, 584, 584, 584, 556, + 1015, 667, 667, 722, 722, 667, 611, 778, 722, 278, 500, 667, 556, 833, 722, 778, + 667, 778, 722, 667, 611, 722, 667, 944, 667, 667, 611, 278, 278, 278, 469, 556, + 333, 556, 556, 500, 556, 556, 278, 556, 556, 222, 222, 500, 222, 833, 556, 556, + 556, 556, 333, 500, 278, 556, 500, 722, 500, 500, 500, 334, 260, 334, 584 +}; + +double graphics_helvetica_width(const char* s) { + double w = 0; + for (const unsigned char* p = (const unsigned char*)s; p && *p; p++) + w += (*p >= 32 && *p <= 126) ? HELV_W[*p - 32] : 556; + return w / 1000.0; +} + /* Text content -> a display string (caller frees), or NULL if unrenderable. */ static char* text_string(const Expr* e) { if (!e) return NULL; @@ -479,6 +507,10 @@ static void draw_prim(Emit* e, const Expr* p, double plot_w) { return; } if (apply_directive(e, p, plot_w)) return; /* Thickness/PointSize/Dashing/Directive */ + if (head_is(p, SYM_Arrowheads) && p->data.function.arg_count >= 1) { + double s; if (to_double(p->data.function.args[0], &s) && s >= 0) e->arrowhead = s; + return; + } /* Line -------------------------------------------------------------- */ if (head_is(p, SYM_Line) && p->data.function.arg_count >= 1) { @@ -566,31 +598,50 @@ static void draw_prim(Emit* e, const Expr* p, double plot_w) { if (head_is(p, SYM_Arrow) && p->data.function.arg_count >= 1) { const Expr* pts = p->data.function.args[0]; if (!head_is(pts, SYM_List) || pts->data.function.arg_count < 2) return; - set_stroke(e); - buf_catf(e->c, "%.3f w\n", stroke_w(e, plot_w)); - double lx = 0, ly = 0, px = 0, py = 0; int first = 1; - for (size_t i = 0; i < pts->data.function.arg_count; i++) { + /* Page-space vertices first, so the shaft can stop at the head's base + * (a butt-capped stroke run to the tip pokes out past the point). */ + size_t np = pts->data.function.arg_count, k = 0; + double* P = (double*)malloc(sizeof(double) * 2 * np); + if (!P) return; + for (size_t i = 0; i < np; i++) { double x, y; if (!get_pt(pts->data.function.args[i], &x, &y)) continue; - px = lx; py = ly; lx = X(e, x); ly = Y(e, y); - buf_catf(e->c, "%.3f %.3f %s\n", lx, ly, first ? "m" : "l"); - first = 0; + P[2 * k] = X(e, x); P[2 * k + 1] = Y(e, y); k++; } - if (!first) buf_cat(e->c, "S\n"); - /* Solid arrowhead on the final segment. Its length scales with the - * final segment (capped), so a short arrow gets a small head — a fixed - * head swamped the many short chevrons of a dashed StreamPlot. Long - * arrows (e.g. VectorPlot) still hit the cap, unchanged. */ + if (k < 2) { free(P); return; } + double lx = P[2 * (k - 1)], ly = P[2 * (k - 1) + 1]; + double px = P[2 * (k - 2)], py = P[2 * (k - 2) + 1]; + /* Solid arrowhead on the final segment. Arrowheads[s] fixes its length + * at s times the plot width; otherwise it scales with the final segment + * (capped), so a short arrow gets a small head -- a fixed head swamped + * the many short chevrons of a dashed StreamPlot. */ double dx = lx - px, dy = ly - py, len = sqrt(dx*dx + dy*dy); + double hl = 0, hw = 0, ux = 0, uy = 0, bx = lx, by = ly; if (len > 1e-6) { - double hl = len * 0.55; - if (hl > 8.0) hl = 8.0; /* cap: long arrows keep the old ~8pt head */ - double ux = dx/len, uy = dy/len, hw = hl * 0.38; /* slimmer than before */ - double bx = lx - ux*hl, by = ly - uy*hl; + if (e->arrowhead > 0) { hl = e->arrowhead * plot_w; hw = hl * 0.42; } + else { + hl = len * 0.55; + if (hl > 8.0) hl = 8.0; /* cap: long arrows keep the ~8pt head */ + hw = hl * 0.38; + } + ux = dx/len; uy = dy/len; + bx = lx - ux*hl; by = ly - uy*hl; + } + set_stroke(e); + buf_catf(e->c, "%.3f w\n", stroke_w(e, plot_w)); + for (size_t i = 0; i + 1 < k; i++) + buf_catf(e->c, "%.3f %.3f %s\n", P[2 * i], P[2 * i + 1], i == 0 ? "m" : "l"); + /* End the shaft inside the head (a little past its base), unless the + * head is longer than the last segment. */ + if (hl > 0 && hl < len) buf_catf(e->c, "%.3f %.3f l S\n", lx - ux*hl*0.8, ly - uy*hl*0.8); + else if (hl > 0) buf_cat(e->c, "S\n"); + else buf_catf(e->c, "%.3f %.3f l S\n", lx, ly); + if (hl > 0) { set_fill(e); buf_catf(e->c, "%.3f %.3f m %.3f %.3f l %.3f %.3f l h f\n", lx, ly, bx - uy*hw, by + ux*hw, bx + uy*hw, by - ux*hw); } + free(P); return; } @@ -598,12 +649,33 @@ static void draw_prim(Emit* e, const Expr* p, double plot_w) { if (head_is(p, SYM_Text) && p->data.function.arg_count >= 2) { double x, y; if (!get_pt(p->data.function.args[1], &x, &y)) return; - char* s = text_string(p->data.function.args[0]); - if (!s) return; + const Expr* body = p->data.function.args[0]; double fs = 10.0; - double tx = X(e, x) - 0.25 * fs * (double)strlen(s); /* rough centring */ - double ty = Y(e, y) - 0.35 * fs; - set_fill(e); + double tr = e->r, tg = e->g, tb = e->b; + /* Style[s, n | FontSize -> n | colour, ...]: size and colour. */ + if (head_is(body, SYM_Style) && body->data.function.arg_count >= 1) { + for (size_t i = 1; i < body->data.function.arg_count; i++) { + const Expr* d = body->data.function.args[i]; + double v, a; + if (to_double(d, &v) && v > 0) fs = v; + else if (head_is(d, SYM_Rule) && d->data.function.arg_count == 2 + && d->data.function.args[0]->type == EXPR_SYMBOL + && d->data.function.args[0]->data.symbol.name == SYM_FontSize + && to_double(d->data.function.args[1], &v) && v > 0) fs = v; + else resolve_color(d, &tr, &tg, &tb, &a); + } + body = body->data.function.args[0]; + } + char* s = text_string(body); + if (!s) return; + /* Offset {ox, oy}: the point of the text box at pos, in box-relative + * coordinates from -1 (left/bottom) to 1; {0, 0} centres. */ + double ox = 0, oy = 0; + if (p->data.function.arg_count >= 3) get_pt(p->data.function.args[2], &ox, &oy); + double tw = graphics_helvetica_width(s) * fs, th = 0.70 * fs; /* cap height */ + double tx = X(e, x) - 0.5 * (1.0 + ox) * tw; + double ty = Y(e, y) - 0.5 * (1.0 + oy) * th; + buf_catf(e->c, "%.4f %.4f %.4f rg\n", tr, tg, tb); buf_cat(e->c, "BT /F1 "); buf_catf(e->c, "%.1f Tf %.3f %.3f Td ", fs, tx, ty); emit_pdf_string(e->c, s); @@ -811,9 +883,22 @@ int graphics_export_pdf(const Expr* g, const char* path) { double mT = 12.0, mB = draw_axes ? 30.0 : 8.0; Deco deco; deco_parse(g, draw_axes, &deco); mT += deco.extra_top; mR += deco.extra_right; + if (o.equal_aspect && !o.size_pair) { + /* Width fixed, height follows the data: no letterboxing. */ + double rw0 = W - mL - mR; + if (rw0 < 20) rw0 = 20; + H = rw0 * (dh / dw) + mT + mB; + if (H > 4.0 * W) H = 4.0 * W; + } double rx = mL, ry = mB, rw = W - mL - mR, rh = H - mT - mB; if (rw < 20) rw = 20; if (rh < 20) rh = 20; + if (o.equal_aspect) { + /* Equal scale: shrink the longer side of the plot region and centre. */ + double sc = rw / dw < rh / dh ? rw / dw : rh / dh; + double nw = dw * sc, nh = dh * sc; + rx += (rw - nw) / 2; ry += (rh - nh) / 2; rw = nw; rh = nh; + } /* Build the content stream. */ Buf content; buf_init(&content); diff --git a/src/graphics/graphics_export.h b/src/graphics/graphics_export.h index 110adb92..1dd49382 100644 --- a/src/graphics/graphics_export.h +++ b/src/graphics/graphics_export.h @@ -21,6 +21,12 @@ #include "expr.h" +/* Width of the string s set in Helvetica at 1pt, from the base-14 AFM + * advance widths (characters outside printable ASCII count as 0.556). The PDF + * writer uses it to align Text[]; layout code (GraphPlot) uses it to size + * labels before emitting them. */ +double graphics_helvetica_width(const char* s); + /* Vector PDF. Dependency-free, headless. `path` is the output file. */ int graphics_export_pdf(const Expr* graphics_expr, const char* path); diff --git a/src/sym_names.c b/src/sym_names.c index 503a7855..a97f69a5 100644 --- a/src/sym_names.c +++ b/src/sym_names.c @@ -968,6 +968,18 @@ const char* SYM_FindVertexColoring = NULL; const char* SYM_EdgeWeight = NULL; const char* SYM_EdgeCapacity = NULL; const char* SYM_WeightedAdjacencyMatrix = NULL; +const char* SYM_VertexCoordinates = NULL; +const char* SYM_VertexLabels = NULL; +const char* SYM_GraphLayout = NULL; +const char* SYM_GraphHighlight = NULL; +const char* SYM_VertexStyle = NULL; +const char* SYM_EdgeStyle = NULL; +const char* SYM_EdgeLabels = NULL; +const char* SYM_VertexSize = NULL; +const char* SYM_HypergraphPlot = NULL; +const char* SYM_Arrowheads = NULL; +const char* SYM_Style = NULL; +const char* SYM_FontSize = NULL; /* NumberForm + Row (numeric-display formatting) and NumberForm's options. */ const char* SYM_NumberForm = NULL; @@ -1930,6 +1942,18 @@ void sym_names_init(void) { SYM_EdgeWeight = intern_symbol("EdgeWeight"); SYM_EdgeCapacity = intern_symbol("EdgeCapacity"); SYM_WeightedAdjacencyMatrix = intern_symbol("WeightedAdjacencyMatrix"); + SYM_VertexCoordinates = intern_symbol("VertexCoordinates"); + SYM_VertexLabels = intern_symbol("VertexLabels"); + SYM_GraphLayout = intern_symbol("GraphLayout"); + SYM_GraphHighlight = intern_symbol("GraphHighlight"); + SYM_VertexStyle = intern_symbol("VertexStyle"); + SYM_EdgeStyle = intern_symbol("EdgeStyle"); + SYM_EdgeLabels = intern_symbol("EdgeLabels"); + SYM_VertexSize = intern_symbol("VertexSize"); + SYM_HypergraphPlot = intern_symbol("HypergraphPlot"); + SYM_Arrowheads = intern_symbol("Arrowheads"); + SYM_Style = intern_symbol("Style"); + SYM_FontSize = intern_symbol("FontSize"); /* System symbols that have no kernel implementation and no cached SYM_* * pointer, but must still be recognized as System` (not qualified into a diff --git a/src/sym_names.h b/src/sym_names.h index 2271e426..38ef875f 100644 --- a/src/sym_names.h +++ b/src/sym_names.h @@ -1027,6 +1027,21 @@ extern const char* SYM_EdgeWeight; extern const char* SYM_EdgeCapacity; extern const char* SYM_WeightedAdjacencyMatrix; +/* GraphPlot / HypergraphPlot options, and the Text/Arrow styling heads the + * PDF exporter honours (Style, FontSize, Arrowheads). */ +extern const char* SYM_VertexCoordinates; +extern const char* SYM_VertexLabels; +extern const char* SYM_GraphLayout; +extern const char* SYM_GraphHighlight; +extern const char* SYM_VertexStyle; +extern const char* SYM_EdgeStyle; +extern const char* SYM_EdgeLabels; +extern const char* SYM_VertexSize; +extern const char* SYM_HypergraphPlot; +extern const char* SYM_Arrowheads; +extern const char* SYM_Style; +extern const char* SYM_FontSize; + /* NumberForm + Row (numeric-display formatting) and NumberForm's option * names. NumberForm is a print wrapper handled in print.c; the option-name * pointers are compared by identity in numberform.c's option parser. */ diff --git a/src/version.h b/src/version.h index d2c27a71..1218bb57 100644 --- a/src/version.h +++ b/src/version.h @@ -13,8 +13,8 @@ * from preprocessor macros; see version.c. mathilda_version() returns it. */ -#define MATHILDA_VERSION_NUMBER 0.236 /* keep in sync with the string below */ -#define MATHILDA_VERSION_STRING "0.236" +#define MATHILDA_VERSION_NUMBER 0.237 /* keep in sync with the string below */ +#define MATHILDA_VERSION_STRING "0.237" /* Full descriptive version string, e.g. * "Mathilda 0.01 (Apple LLVM 17.0.0, GMP 6.3.0, MPFR 4.2.2, FLINT 3.6.0, ...)" diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 969a1fa9..0a932458 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -976,6 +976,7 @@ set(COMMON_SRC ../src/graph/spanningtree.c ../src/graph/connectivity.c ../src/graph/graphplot.c + ../src/graph/glayout.c ../src/graph/vertexcoloring.c ../src/graph/edgeweight.c ../src/graph/wtadjmat.c @@ -1001,6 +1002,7 @@ set(COMMON_SRC ../src/graph/hyp_random.c ../src/graph/hyp_transversal.c ../src/graph/hyp_init.c + ../src/graph/hyp_plot.c ../src/graph/galg_common.c ../src/graph/galg_init.c ../src/graph/galg_flow.c @@ -1079,6 +1081,10 @@ add_executable(graph_metrics_tests test_graph_metrics.c $) +target_link_libraries(graphplot_tests m) +target_include_directories(graphplot_tests PRIVATE ../src ../src/graph ../src/graphics) +add_test(NAME graphplot_tests COMMAND graphplot_tests) add_executable(hypergraph_tests test_hypergraph.c $) target_link_libraries(hypergraph_tests m) target_include_directories(hypergraph_tests PRIVATE ../src ../src/graph) diff --git a/tests/test_graphplot.c b/tests/test_graphplot.c new file mode 100644 index 00000000..60af426f --- /dev/null +++ b/tests/test_graphplot.c @@ -0,0 +1,277 @@ +/* test_graphplot.c - GraphPlot / HypergraphPlot and the layout engine + * (src/graph/graphplot.c, hyp_plot.c, glayout.c). + * + * Determinism (the same call gives an identical Graphics expression, and the + * layout engine gives bit-identical coordinates), every option (arrows for + * directed edges, red + thicker highlighting, VertexCoordinates in both forms, + * labels, styles, edge weights, each GraphLayout), the frame (no axes, equal + * aspect ratio), the degenerate inputs (empty graph, one vertex, disconnected + * and edgeless graphs), and 500-vertex graphs in bounded time. + */ + +#include "expr.h" +#include "eval.h" +#include "core.h" +#include "symtab.h" +#include "parse.h" +#include "print.h" +#include "graph.h" +#include "glayout.h" +#include "test_utils.h" +#include +#include +#include +#include + + +#define PETERSEN "PetersenGraph[]" +#define H1 "Hypergraph[{{1,2,3},{3,4},{4,5,6},{7}}]" + +static void test_basic_shape(void) { + assert_eval_eq("Head[GraphPlot[" PETERSEN "]]", "Graphics", 0); + /* one Line per undirected edge, one Disk per vertex, no labels by default */ + assert_eval_eq("Count[GraphPlot[" PETERSEN "], _Line, Infinity]", "15", 0); + assert_eval_eq("Count[GraphPlot[" PETERSEN "], _Disk, Infinity]", "10", 0); + assert_eval_eq("Count[GraphPlot[" PETERSEN "], _Text, Infinity]", "0", 0); + /* frame: no axes, equal aspect ratio, explicit PlotRange and ImageSize */ + assert_eval_eq("Axes /. Rest[List @@ GraphPlot[" PETERSEN "]]", "False", 0); + assert_eval_eq("AspectRatio /. Rest[List @@ GraphPlot[" PETERSEN "]]", "Automatic", 0); + assert_eval_eq("MatchQ[PlotRange /. Rest[List @@ GraphPlot[" PETERSEN "]], " + "{{_Real, _Real}, {_Real, _Real}}]", "True", 0); + assert_eval_eq("MatchQ[ImageSize /. Rest[List @@ GraphPlot[" PETERSEN "]], " + "{_Integer, _Integer}]", "True", 0); + /* default colours: Mathematica's vertex blue */ + assert_eval_eq("MemberQ[GraphPlot[" PETERSEN "], RGBColor[0.368417, 0.506779, 0.709798], Infinity]", + "True", 0); + /* a user option passes through to Graphics */ + assert_eval_eq("PlotLabel /. Rest[List @@ GraphPlot[CycleGraph[3], PlotLabel -> \"C3\"]]", + "\"C3\"", 0); + /* non-graphs and non-option arguments stay unevaluated */ + assert_eval_eq("Head[GraphPlot[5]]", "GraphPlot", 0); + assert_eval_eq("Head[GraphPlot[CycleGraph[3], 7]]", "GraphPlot", 0); + /* a list of rules is taken as the graph it denotes */ + assert_eval_eq("Count[GraphPlot[{1 -> 2, 2 -> 3}], _Arrow, Infinity]", "2", 0); +} + +static void test_determinism(void) { + assert_eval_eq("GraphPlot[" PETERSEN "] === GraphPlot[" PETERSEN "]", "True", 0); + assert_eval_eq("GraphPlot[GridGraph[{5, 4}]] === GraphPlot[GridGraph[{5, 4}]]", "True", 0); + assert_eval_eq("GraphPlot[CompleteKaryTree[4, 3], VertexLabels -> \"Name\"] === " + "GraphPlot[CompleteKaryTree[4, 3], VertexLabels -> \"Name\"]", "True", 0); + assert_eval_eq("(SeedRandom[7]; g = RandomGraph[{60, 120}]; " + "GraphPlot[g, GraphLayout -> \"SpringElectricalEmbedding\"] === " + "GraphPlot[g, GraphLayout -> \"SpringElectricalEmbedding\"])", "True", 0); + assert_eval_eq("HypergraphPlot[" H1 "] === HypergraphPlot[" H1 "]", "True", 0); + + /* The engine itself: bit-identical coordinates on repeated calls, for + * every method. */ + const int n = 12, m = 18; + int eu[18], ev[18]; + unsigned char dir[18]; + for (int k = 0; k < m; k++) { + eu[k] = k % n; ev[k] = (k * 5 + 1) % n; + if (ev[k] == eu[k]) ev[k] = (ev[k] + 1) % n; + dir[k] = 0; + } + GLMethod ms[] = { GL_AUTOMATIC, GL_CIRCULAR, GL_SPRING, GL_STRESS, GL_LAYERED, + GL_BIPARTITE, GL_GRID }; + for (size_t t = 0; t < sizeof(ms) / sizeof(ms[0]); t++) { + double a[24], b[24]; + int ra = glayout_compute(n, m, eu, ev, dir, ms[t], a); + int rb = glayout_compute(n, m, eu, ev, dir, ms[t], b); + ASSERT_MSG(ra >= 0 && ra == rb, "method %d: %d vs %d", (int)ms[t], ra, rb); + ASSERT_MSG(memcmp(a, b, sizeof(a)) == 0, "method %d not deterministic", (int)ms[t]); + for (int i = 0; i < 2 * n; i++) ASSERT_MSG(isfinite(a[i]), "method %d: non-finite", (int)ms[t]); + } + /* n == 0 is a no-op success */ + ASSERT_MSG(glayout_compute(0, 0, NULL, NULL, NULL, GL_AUTOMATIC, NULL) >= 0, "empty layout"); +} + +static void test_directed(void) { + /* directed edges are Arrows (with an Arrowheads directive), undirected Lines */ + assert_eval_eq("Count[GraphPlot[Graph[{1 -> 2, 2 -> 3, 3 -> 1}]], _Arrow, Infinity]", "3", 0); + assert_eval_eq("Count[GraphPlot[Graph[{1 -> 2, 2 -> 3, 3 -> 1}]], _Line, Infinity]", "0", 0); + assert_eval_eq("Count[GraphPlot[Graph[{1 -> 2, 2 <-> 3}]], _Arrow, Infinity]", "1", 0); + assert_eval_eq("Count[GraphPlot[Graph[{1 -> 2, 2 <-> 3}]], _Line, Infinity]", "1", 0); + assert_eval_eq("Count[GraphPlot[Graph[{1 -> 2}]], _Arrowheads, Infinity] >= 1", "True", 0); + /* the arrow stops outside the target disk and starts outside the source */ + assert_eval_eq("Module[{g = GraphPlot[Graph[{1 -> 2}], VertexCoordinates -> {{0, 0}, {1, 0}}], r, s, t}," + " r = Cases[g, Disk[_, x_] :> x, Infinity][[1]];" + " {s, t} = Cases[g, Arrow[{a_, b_}] :> {a, b}, Infinity][[1]];" + " {s[[1]] >= r, 1 - t[[1]] > r}]", "{True, True}", 0); + /* a mutual pair u->v, v->u is drawn as two separated arrows */ + assert_eval_eq("Module[{g = GraphPlot[Graph[{1 -> 2, 2 -> 1}], VertexCoordinates -> {{0, 0}, {1, 0}}]}," + " Length[Union[Cases[g, Arrow[{{_, y_}, _}] :> y, Infinity]]]]", "2", 0); + /* a DAG defaults to a layered drawing: every edge points downwards */ + assert_eval_eq("Module[{g = GraphPlot[Graph[{1 -> 2, 1 -> 3, 2 -> 4, 3 -> 4, 4 -> 5, 2 -> 5}]]}," + " And @@ Cases[g, Arrow[{{_, y0_}, {_, y1_}}] :> y1 < y0, Infinity]]", "True", 0); +} + +static void test_highlight_and_styles(void) { + /* highlighted vertex: red, drawn larger */ + assert_eval_eq("Count[GraphPlot[CycleGraph[3]], RGBColor[1., 0., 0.], Infinity]", "0", 0); + assert_eval_eq("MatchQ[First[GraphPlot[CycleGraph[3], GraphHighlight -> {1}]]," + " {___, RGBColor[1., 0., 0.], _Disk, ___}]", "True", 0); + assert_eval_eq("Length[Union[Cases[GraphPlot[CycleGraph[3], GraphHighlight -> {1}]," + " Disk[_, r_] :> r, Infinity]]]", "2", 0); + /* highlighted edge: red and thicker than the others */ + assert_eval_eq("MatchQ[First[GraphPlot[CycleGraph[3], GraphHighlight -> {1 <-> 2}]]," + " {___, RGBColor[1., 0., 0.], _Line, ___}]", "True", 0); + assert_eval_eq("Module[{t = Cases[GraphPlot[CycleGraph[3], GraphHighlight -> {1 <-> 2}]," + " Thickness[x_] :> x, Infinity]}, Length[t] >= 2 && t[[2]] > 2 t[[1]]]", "True", 0); + /* a highlighted directed edge (Rule syntax) */ + assert_eval_eq("MatchQ[First[GraphPlot[Graph[{1 -> 2, 2 -> 3}], GraphHighlight -> {2 -> 3}]]," + " {___, RGBColor[1., 0., 0.], _Arrow, ___}]", "True", 0); + /* VertexStyle: per-vertex rules and a global colour */ + assert_eval_eq("MatchQ[First[GraphPlot[CycleGraph[3], VertexStyle -> {2 -> RGBColor[0, 1, 0]}]]," + " {___, RGBColor[0, 1, 0], _Disk, ___}]", "True", 0); + assert_eval_eq("Count[First[GraphPlot[CycleGraph[4], VertexStyle -> RGBColor[0, 0, 1]]]," + " RGBColor[0, 0, 1]]", "4", 0); + /* EdgeStyle rules */ + assert_eval_eq("MatchQ[First[GraphPlot[CycleGraph[4], EdgeStyle -> {(2 <-> 3) -> RGBColor[1, 0, 1]}]]," + " {___, RGBColor[1, 0, 1], _Line, ___}]", "True", 0); +} + +static void test_labels(void) { + assert_eval_eq("Count[GraphPlot[" PETERSEN ", VertexLabels -> \"Name\"], _Text, Infinity]", "10", 0); + assert_eval_eq("Count[GraphPlot[" PETERSEN ", VertexLabels -> Automatic], _Text, Infinity]", "10", 0); + assert_eval_eq("Count[GraphPlot[" PETERSEN ", VertexLabels -> None], _Text, Infinity]", "0", 0); + assert_eval_eq("Cases[GraphPlot[CycleGraph[3], VertexLabels -> {2 -> \"two\"}]," + " Text[s_, __] :> s, Infinity]", "{\"two\"}", 0); + /* labels are Strings with an offset, so the PDF writer can align them */ + assert_eval_eq("MatchQ[Cases[GraphPlot[CycleGraph[3], VertexLabels -> \"Name\"], _Text, Infinity]," + " {Text[_String, {_, _}, {_, _}] ..}]", "True", 0); + /* labels fit inside the PlotRange */ + assert_eval_eq("Module[{g = GraphPlot[CompleteGraph[5], VertexLabels -> \"Name\"], pr, pts}," + " pr = PlotRange /. Rest[List @@ g];" + " pts = Cases[g, Text[_, p_, _] :> p, Infinity];" + " And @@ (pr[[1, 1]] < #[[1]] < pr[[1, 2]] && pr[[2, 1]] < #[[2]] < pr[[2, 2]] & /@ pts)]", + "True", 0); + /* edge weights */ + assert_eval_eq("Cases[GraphPlot[Graph[{1, 2, 3}, {1 <-> 2, 2 <-> 3}, EdgeWeight -> {5, 9}]," + " EdgeLabels -> \"EdgeWeight\"], Text[s_, __] :> s, Infinity]", "{\"5\", \"9\"}", 0); +} + +static void test_coordinates_and_layouts(void) { + /* VertexCoordinates as a list, in VertexList order */ + assert_eval_eq("Cases[GraphPlot[PathGraph[{1, 2, 3}], VertexCoordinates -> {{0, 0}, {1, 0}, {2, 1}}]," + " Disk[p_, _] :> p, Infinity]", "{{0.0, 0.0}, {1.0, 0.0}, {2.0, 1.0}}", 0); + /* ... and as rules for a subset (the rest from the layout) */ + assert_eval_eq("Cases[GraphPlot[CycleGraph[4], VertexCoordinates -> {3 -> {5, 5}}]," + " Disk[p_, _] :> p, Infinity][[3]]", "{5.0, 5.0}", 0); + /* every GraphLayout gives a Graphics with all the vertices */ + const char* layouts[] = { "CircularEmbedding", "SpringElectricalEmbedding", "StressEmbedding", + "LayeredEmbedding", "BipartiteEmbedding", "GridEmbedding", NULL }; + for (int i = 0; layouts[i]; i++) { + char buf[256]; + snprintf(buf, sizeof(buf), "Count[GraphPlot[CycleGraph[6], GraphLayout -> \"%s\"], _Disk, Infinity]", + layouts[i]); + assert_eval_eq(buf, "6", 0); + } + /* CircularEmbedding: every vertex at the same distance from the centre */ + assert_eval_eq("Module[{p = Cases[GraphPlot[" PETERSEN ", GraphLayout -> \"CircularEmbedding\"]," + " Disk[q_, _] :> q, Infinity], c}, c = Mean[p];" + " Max[Norm[# - c] & /@ p] - Min[Norm[# - c] & /@ p] < 10^-4]", "True", 0); + /* BipartiteEmbedding: two columns */ + assert_eval_eq("Length[Union[Cases[GraphPlot[CompleteGraph[{3, 4}], GraphLayout -> \"BipartiteEmbedding\"]," + " Disk[{x_, _}, _] :> x, Infinity]]]", "2", 0); + /* a grid is drawn as a grid: 4 distinct x and 4 distinct y coordinates */ + assert_eval_eq("Module[{p = Cases[GraphPlot[GridGraph[{4, 4}]], Disk[q_, _] :> q, Infinity]}," + " {Length[Union[Round[p[[All, 1]], 0.001]]], Length[Union[Round[p[[All, 2]], 0.001]]]}]", + "{4, 4}", 0); + /* trees default to a layered drawing with the root on top */ + assert_eval_eq("Module[{p = Cases[GraphPlot[CompleteKaryTree[3, 2]], Disk[q_, _] :> q, Infinity]}," + " p[[1, 2]] > Max[Rest[p][[All, 2]]]]", "True", 0); + /* a cycle is a regular polygon: all edges the same length */ + assert_eval_eq("Module[{g = GraphPlot[CycleGraph[7]], l}, l = Cases[g, Line[{a_, b_}] :> Norm[a - b], Infinity];" + " Max[l] - Min[l] < 10^-3 Max[l]]", "True", 0); +} + +static void test_degenerate(void) { + assert_eval_eq("Head[GraphPlot[Graph[{}, {}]]]", "Graphics", 0); + assert_eval_eq("Count[GraphPlot[Graph[{}, {}]], _Disk, Infinity]", "0", 0); + assert_eval_eq("Count[GraphPlot[Graph[{a}, {}], VertexLabels -> \"Name\"], _Disk, Infinity]", "1", 0); + /* an edgeless graph: distinct positions for every vertex */ + assert_eval_eq("Length[Union[Cases[GraphPlot[Graph[Range[10], {}]], Disk[p_, _] :> p, Infinity]]]", "10", 0); + /* disconnected: components do not overlap (distinct positions) */ + assert_eval_eq("Module[{g = GraphDisjointUnion[CycleGraph[5], CompleteGraph[4]], p}," + " p = Cases[GraphPlot[g], Disk[q_, _] :> q, Infinity]; {Length[p], Length[Union[p]]}]", + "{9, 9}", 0); + /* a two-vertex graph and a path */ + assert_eval_eq("Count[GraphPlot[Graph[{1 <-> 2}]], _Line, Infinity]", "1", 0); + assert_eval_eq("Count[GraphPlot[PathGraph[Range[5]]], _Disk, Infinity]", "5", 0); +} + +static double now_s(void) { return (double)clock() / CLOCKS_PER_SEC; } + +static void test_large(void) { + double t0 = now_s(); + assert_eval_eq("(SeedRandom[3]; Count[GraphPlot[RandomGraph[{500, 1500}]], _Disk, Infinity])", "500", 0); + assert_eval_eq("Count[GraphPlot[GridGraph[{25, 20}]], _Disk, Infinity]", "500", 0); + assert_eval_eq("Count[GraphPlot[CompleteKaryTree[9, 2]], _Disk, Infinity]", "511", 0); + assert_eval_eq("Count[GraphPlot[Graph[Range[500], {}]], _Disk, Infinity]", "500", 0); + assert_eval_eq("(SeedRandom[4]; Count[GraphPlot[RandomGraph[{500, 1500}]," + " GraphLayout -> \"SpringElectricalEmbedding\"], _Disk, Infinity])", "500", 0); + assert_eval_eq("(SeedRandom[5]; Count[GraphPlot[RandomGraph[{1500, 3000}]], _Disk, Infinity])", "1500", 0); + double dt = now_s() - t0; + ASSERT_MSG(dt < 20.0, "large graphs took %.2fs", dt); + printf("(%.2fs) ", dt); +} + +static void test_hypergraph_plot(void) { + assert_eval_eq("Head[HypergraphPlot[" H1 "]]", "Graphics", 0); + /* one translucent shape per hyperedge (sizes 3, 2, 3 and 1), one disk per vertex */ + assert_eval_eq("Count[HypergraphPlot[" H1 "], _Polygon, Infinity]", "4", 0); + assert_eval_eq("Count[HypergraphPlot[" H1 "], _Disk, Infinity]", "7", 0); + assert_eval_eq("MemberQ[HypergraphPlot[" H1 "], Opacity[x_ /; x < 1], Infinity]", "True", 0); + /* hyperedges get distinct colours */ + assert_eval_eq("Module[{p = First[HypergraphPlot[" H1 "]]}," + " Length[Union[Cases[p, c_RGBColor /; MemberQ[p, c], 1]]] >= 4]", "True", 0); + assert_eval_eq("Axes /. Rest[List @@ HypergraphPlot[" H1 "]]", "False", 0); + assert_eval_eq("Count[HypergraphPlot[" H1 ", VertexLabels -> \"Name\"], _Text, Infinity]", "7", 0); + /* a plain list of hyperedges */ + assert_eval_eq("Count[HypergraphPlot[{{1, 2, 3}, {3, 4}}], _Polygon, Infinity]", "2", 0); + /* VertexCoordinates: every rounded hull encloses its members */ + assert_eval_eq("Cases[HypergraphPlot[Hypergraph[{{1, 2}}], VertexCoordinates -> {{0, 0}, {3, 0}}]," + " Disk[p_, _] :> p, Infinity]", "{{0.0, 0.0}, {3.0, 0.0}}", 0); + assert_eval_eq("Module[{g = HypergraphPlot[Hypergraph[{{1, 2, 3}}]," + " VertexCoordinates -> {{0, 0}, {2, 0}, {1, 2}}], poly}," + " poly = Cases[g, Polygon[p_] :> p, Infinity][[1]];" + " {Min[poly[[All, 1]]] < 0, Max[poly[[All, 1]]] > 2, Max[poly[[All, 2]]] > 2}]", + "{True, True, True}", 0); + /* degenerate: empty and edgeless */ + assert_eval_eq("Head[HypergraphPlot[Hypergraph[{}, {}]]]", "Graphics", 0); + assert_eval_eq("Count[HypergraphPlot[Hypergraph[{1, 2, 3}, {}]], _Disk, Infinity]", "3", 0); + assert_eval_eq("Head[HypergraphPlot[5]]", "HypergraphPlot", 0); + /* larger hypergraph in bounded time */ + double t0 = now_s(); + assert_eval_eq("(SeedRandom[2]; Count[HypergraphPlot[RandomHypergraph[{300, 200}, 3]], _Disk, Infinity])", + "300", 0); + ASSERT_MSG(now_s() - t0 < 10.0, "hypergraph plot too slow"); +} + +static void test_pdf_export(void) { + assert_eval_eq("Export[\"/tmp/mathilda_test_graphplot.pdf\", GraphPlot[" PETERSEN + ", VertexLabels -> \"Name\"]]", "\"/tmp/mathilda_test_graphplot.pdf\"", 0); + assert_eval_eq("Export[\"/tmp/mathilda_test_hypergraphplot.pdf\", HypergraphPlot[" H1 "]]", + "\"/tmp/mathilda_test_hypergraphplot.pdf\"", 0); + remove("/tmp/mathilda_test_graphplot.pdf"); + remove("/tmp/mathilda_test_hypergraphplot.pdf"); +} + +int main(void) { + symtab_init(); + core_init(); + TEST(test_basic_shape); + TEST(test_determinism); + TEST(test_directed); + TEST(test_highlight_and_styles); + TEST(test_labels); + TEST(test_coordinates_and_layouts); + TEST(test_degenerate); + TEST(test_large); + TEST(test_hypergraph_plot); + TEST(test_pdf_export); + printf("All graphplot tests passed!\n"); + return 0; +}