Graphviz 13.1.3~dev.20250829.0113
Loading...
Searching...
No Matches
sgd.c
Go to the documentation of this file.
1#include <assert.h>
2#include <limits.h>
3#include <math.h>
4#include <neatogen/dijkstra.h>
5#include <neatogen/neato.h>
8#include <neatogen/sgd.h>
9#include <stdlib.h>
10#include <util/alloc.h>
11#include <util/bitarray.h>
12#include <util/gv_math.h>
13#include <util/unreachable.h>
14
15static double calculate_stress(double *pos, term_sgd *terms, int n_terms) {
16 double stress = 0;
17 for (int ij = 0; ij < n_terms; ij++) {
18 const double dx = pos[2 * terms[ij].i] - pos[2 * terms[ij].j];
19 const double dy = pos[2 * terms[ij].i + 1] - pos[2 * terms[ij].j + 1];
20 const double r = hypot(dx, dy) - terms[ij].d;
21 stress += terms[ij].w * (r * r);
22 }
23 return stress;
24}
25// it is much faster to shuffle term rather than pointers to term, even though
26// the swap is more expensive
27static void fisheryates_shuffle(term_sgd *terms, int n_terms,
28 rk_state *rstate) {
29 for (int i = n_terms - 1; i >= 1; i--) {
30 int j = rk_interval(i, rstate);
31
32 SWAP(&terms[i], &terms[j]);
33 }
34}
35
36// graph_sgd data structure exists only to make dijkstras faster
37static graph_sgd *extract_adjacency(graph_t *G, int model) {
38 size_t n_nodes = 0, n_edges = 0;
39 for (node_t *np = agfstnode(G); np; np = agnxtnode(G, np)) {
40 assert(ND_id(np) == n_nodes);
41 n_nodes++;
42 for (edge_t *ep = agfstedge(G, np); ep; ep = agnxtedge(G, ep, np)) {
43 if (agtail(ep) != aghead(ep)) { // ignore self-loops and double edges
44 n_edges++;
45 }
46 }
47 }
49 graph->sources = gv_calloc(n_nodes + 1, sizeof(size_t));
50 graph->pinneds = bitarray_new(n_nodes);
51 graph->targets = gv_calloc(n_edges, sizeof(size_t));
52 graph->weights = gv_calloc(n_edges, sizeof(float));
53
54 graph->n = n_nodes;
55 assert(n_edges <= INT_MAX);
56 graph->sources[graph->n] = n_edges; // to make looping nice
57
58 n_nodes = 0, n_edges = 0;
59 for (node_t *np = agfstnode(G); np; np = agnxtnode(G, np)) {
60 assert(n_edges <= INT_MAX);
61 graph->sources[n_nodes] = n_edges;
62 bitarray_set(&graph->pinneds, n_nodes, isFixed(np));
63 for (edge_t *ep = agfstedge(G, np); ep; ep = agnxtedge(G, ep, np)) {
64 if (agtail(ep) == aghead(ep)) { // ignore self-loops and double edges
65 continue;
66 }
67 node_t *target = (agtail(ep) == np)
68 ? aghead(ep)
69 : agtail(ep); // in case edge is reversed
70 graph->targets[n_edges] = (size_t)ND_id(target);
71 graph->weights[n_edges] = ED_dist(ep);
72 assert(graph->weights[n_edges] > 0);
73 n_edges++;
74 }
75 n_nodes++;
76 }
77 assert(n_nodes == graph->n);
78 assert(n_edges <= INT_MAX);
79 assert(n_edges == graph->sources[graph->n]);
80 graph->sources[n_nodes] = n_edges;
81
82 if (model == MODEL_SHORTPATH) {
83 // do nothing
84 } else if (model == MODEL_SUBSET) {
85 // i,j,k refer to actual node indices, while x,y refer to edge indices in
86 // graph->targets initialise to no neighbours
87 bitarray_t neighbours_i = bitarray_new(graph->n);
88 bitarray_t neighbours_j = bitarray_new(graph->n);
89 for (size_t i = 0; i < graph->n; i++) {
90 int deg_i = 0;
91 for (size_t x = graph->sources[i]; x < graph->sources[i + 1]; x++) {
92 size_t j = graph->targets[x];
93 if (!bitarray_get(neighbours_i, j)) { // ignore multiedges
94 bitarray_set(&neighbours_i, j, true); // set up sort of hashset
95 deg_i++;
96 }
97 }
98 for (size_t x = graph->sources[i]; x < graph->sources[i + 1]; x++) {
99 size_t j = graph->targets[x];
100 int intersect = 0;
101 int deg_j = 0;
102 for (size_t y = graph->sources[j]; y < graph->sources[j + 1]; y++) {
103 size_t k = graph->targets[y];
104 if (!bitarray_get(neighbours_j, k)) { // ignore multiedges
105 bitarray_set(&neighbours_j, k, true); // set up sort of hashset
106 deg_j++;
107 if (bitarray_get(neighbours_i, k)) {
108 intersect++;
109 }
110 }
111 }
112 graph->weights[x] = deg_i + deg_j - (2 * intersect);
113 assert(graph->weights[x] > 0);
114 for (size_t y = graph->sources[j]; y < graph->sources[j + 1]; y++) {
115 size_t k = graph->targets[y];
116 bitarray_set(&neighbours_j, k, false); // reset sort of hashset
117 }
118 }
119 for (size_t x = graph->sources[i]; x < graph->sources[i + 1]; x++) {
120 size_t j = graph->targets[x];
121 bitarray_set(&neighbours_i, j, false); // reset sort of hashset
122 }
123 }
124 bitarray_reset(&neighbours_i);
125 bitarray_reset(&neighbours_j);
126 } else {
127 // TODO: model == MODEL_MDS and MODEL_CIRCUIT
128 UNREACHABLE(); // mds and circuit model not supported
129 }
130 return graph;
131}
133 free(graph->sources);
134 bitarray_reset(&graph->pinneds);
135 free(graph->targets);
136 free(graph->weights);
137 free(graph);
138}
139
140void sgd(graph_t *G, /* input graph */
141 int model /* distance model */) {
142 if (model == MODEL_CIRCUIT) {
143 agwarningf("circuit model not yet supported in Gmode=sgd, reverting to "
144 "shortpath model\n");
145 model = MODEL_SHORTPATH;
146 }
147 if (model == MODEL_MDS) {
148 agwarningf("mds model not yet supported in Gmode=sgd, reverting to "
149 "shortpath model\n");
150 model = MODEL_SHORTPATH;
151 }
152 int n = agnnodes(G);
153
154 if (Verbose) {
155 fprintf(stderr, "calculating shortest paths and setting up stress terms:");
156 start_timer();
157 }
158 // calculate how many terms will be needed as fixed nodes can be ignored
159 int n_fixed = 0, n_terms = 0;
160 for (int i = 0; i < n; i++) {
161 if (!isFixed(GD_neato_nlist(G)[i])) {
162 n_fixed++;
163 n_terms += n - n_fixed;
164 }
165 }
166 term_sgd *terms = gv_calloc(n_terms, sizeof(term_sgd));
167 // calculate term values through shortest paths
168 int offset = 0;
170 for (int i = 0; i < n; i++) {
171 if (!isFixed(GD_neato_nlist(G)[i])) {
172 offset += dijkstra_sgd(graph, i, terms + offset);
173 }
174 }
175 assert(offset == n_terms);
177 if (Verbose) {
178 fprintf(stderr, " %.2f sec\n", elapsed_sec());
179 }
180
181 // initialise annealing schedule
182 float w_min = terms[0].w, w_max = terms[0].w;
183 for (int ij = 1; ij < n_terms; ij++) {
184 w_min = fminf(w_min, terms[ij].w);
185 w_max = fmaxf(w_max, terms[ij].w);
186 }
187 // note: Epsilon is different from MODE_KK and MODE_MAJOR as it is a minimum
188 // step size rather than energy threshold
189 // MaxIter is also different as it is a fixed number of iterations
190 // rather than a maximum
191 const double eta_max = 1.0 / w_min;
192 const double eta_min = Epsilon / w_max;
193 const double lambda = log(eta_max / eta_min) / (MaxIter - 1);
194
195 // initialise starting positions (from neatoprocs)
197 // copy initial positions and state into temporary space for speed
198 double *const pos = gv_calloc(2 * n, sizeof(double));
199 bool *unfixed = gv_calloc(n, sizeof(bool));
200 for (int i = 0; i < n; i++) {
202 pos[2 * i] = ND_pos(node)[0];
203 pos[2 * i + 1] = ND_pos(node)[1];
204 unfixed[i] = !isFixed(node);
205 }
206
207 // perform optimisation
208 if (Verbose) {
209 fprintf(stderr, "solving model:");
210 start_timer();
211 }
212 rk_state rstate;
213 rk_seed(0, &rstate); // TODO: get seed from graph
214 for (int t = 0; t < MaxIter; t++) {
215 fisheryates_shuffle(terms, n_terms, &rstate);
216 const double eta = eta_max * exp(-lambda * t);
217 for (int ij = 0; ij < n_terms; ij++) {
218 // cap step size
219 const double mu = fmin(eta * terms[ij].w, 1);
220
221 const double dx = pos[2 * terms[ij].i] - pos[2 * terms[ij].j];
222 const double dy = pos[2 * terms[ij].i + 1] - pos[2 * terms[ij].j + 1];
223 const double mag = hypot(dx, dy);
224
225 const double r = (mu * (mag - terms[ij].d)) / (2 * mag);
226 const double r_x = r * dx;
227 const double r_y = r * dy;
228
229 if (unfixed[terms[ij].i]) {
230 pos[2 * terms[ij].i] -= r_x;
231 pos[2 * terms[ij].i + 1] -= r_y;
232 }
233 if (unfixed[terms[ij].j]) {
234 pos[2 * terms[ij].j] += r_x;
235 pos[2 * terms[ij].j + 1] += r_y;
236 }
237 }
238 if (Verbose) {
239 fprintf(stderr, " %.3f", calculate_stress(pos, terms, n_terms));
240 }
241 }
242 if (Verbose) {
243 fprintf(stderr, "\nfinished in %.2f sec\n", elapsed_sec());
244 }
245 free(terms);
246
247 // copy temporary positions back into graph_t
248 for (int i = 0; i < n; i++) {
250 ND_pos(node)[0] = pos[2 * i];
251 ND_pos(node)[1] = pos[2 * i + 1];
252 }
253 free(pos);
254 free(unfixed);
255}
Memory allocation wrappers that exit on failure.
static void * gv_calloc(size_t nmemb, size_t size)
Definition alloc.h:26
static void * gv_alloc(size_t size)
Definition alloc.h:47
#define Epsilon
Definition arcball.h:209
API for compacted arrays of booleans.
static bitarray_t bitarray_new(size_t size_bits)
create an array of the given element length
Definition bitarray.h:47
static bool bitarray_get(bitarray_t self, size_t index)
get the value of the given element
Definition bitarray.h:65
static void bitarray_set(bitarray_t *self, size_t index, bool value)
set or clear the value of the given element
Definition bitarray.h:80
static void bitarray_reset(bitarray_t *self)
free underlying resources and leave a bit array empty
Definition bitarray.h:114
static float dy
Definition draw.c:40
static float dx
Definition draw.c:39
#define G
Definition gdefs.h:7
int MaxIter
Definition globals.h:61
static bool Verbose
Definition gml2gv.c:24
void free(void *)
int agnnodes(Agraph_t *g)
Definition graph.c:155
#define ED_dist(e)
Definition types.h:602
#define agtail(e)
Definition cgraph.h:988
Agedge_t * agnxtedge(Agraph_t *g, Agedge_t *e, Agnode_t *n)
Definition edge.c:96
#define aghead(e)
Definition cgraph.h:989
Agedge_t * agfstedge(Agraph_t *g, Agnode_t *n)
Definition edge.c:87
void agwarningf(const char *fmt,...)
Definition agerror.c:173
#define GD_neato_nlist(g)
Definition types.h:392
Agnode_t * agnxtnode(Agraph_t *g, Agnode_t *n)
Definition node.c:48
Agnode_t * agfstnode(Agraph_t *g)
Definition node.c:41
#define ND_pos(n)
Definition types.h:520
Agraph_t * graph(char *name)
Definition gv.cpp:30
Arithmetic helper functions.
#define SWAP(a, b)
Definition gv_math.h:131
int dijkstra_sgd(graph_sgd *graph, int source, term_sgd *terms)
Definition dijkstra.c:291
#define isFixed(n)
Definition macros.h:19
#define mu
Definition maze.c:137
#define ND_id(n)
Definition mm2gv.c:40
#define MODEL_SUBSET
Definition neato.h:17
#define MODEL_MDS
Definition neato.h:18
#define MODEL_CIRCUIT
Definition neato.h:16
#define MODEL_SHORTPATH
Definition neato.h:15
NEATOPROCS_API void initial_positions(graph_t *, int)
Definition stuff.c:328
unsigned long rk_interval(unsigned long max, rk_state *state)
Definition randomkit.c:138
void rk_seed(unsigned long seed, rk_state *state)
Definition randomkit.c:75
static void free_adjacency(graph_sgd *graph)
Definition sgd.c:132
static void fisheryates_shuffle(term_sgd *terms, int n_terms, rk_state *rstate)
Definition sgd.c:27
static double calculate_stress(double *pos, term_sgd *terms, int n_terms)
Definition sgd.c:15
void sgd(graph_t *G, int model)
Definition sgd.c:140
static graph_sgd * extract_adjacency(graph_t *G, int model)
Definition sgd.c:37
graph or subgraph
Definition cgraph.h:424
Definition sgd.h:16
int j
Definition sgd.h:17
float w
Definition sgd.h:18
int i
Definition sgd.h:17
double elapsed_sec(void)
Definition timing.c:48
void start_timer(void)
Definition timing.c:43
#define UNREACHABLE()
Definition unreachable.h:30
static bool intersect(Ppoint_t a, Ppoint_t b, Ppoint_t c, Ppoint_t d)
Definition visibility.c:78