From 2377b968f55ddd5d6a8ff79ee9150316a218afc8 Mon Sep 17 00:00:00 2001 From: Ian Jackson Date: Thu, 3 Jan 2008 20:13:06 +0000 Subject: [PATCH 1/1] new energy calculation but anarres crashes --- Makefile | 2 +- energy.c | 93 ++++++++++++++++++++++++++++++----- graph.c | 116 ++++++++++++++++++++++++++++++++++++++++++++ mgraph.h | 2 +- bgl.h => minimise.h | 8 +-- 5 files changed, 203 insertions(+), 18 deletions(-) create mode 100644 graph.c rename bgl.h => minimise.h (70%) diff --git a/Makefile b/Makefile index 3bb3ef9..e1a54dd 100644 --- a/Makefile +++ b/Makefile @@ -14,7 +14,7 @@ o= >$@.new && mv -f $@.new $@ all: $(TARGETS) -minimise: energy.o bgl.o common.o mgraph.o +minimise: energy.o graph.o common.o mgraph.o $(CXX) $(CXXFLAGS) -o $@ $^ $(LIBGSL) primer: primer.o common.o diff --git a/energy.c b/energy.c index 07a4100..73308ac 100644 --- a/energy.c +++ b/energy.c @@ -3,40 +3,52 @@ */ #include "common.h" -#include "bgl.h" +#include "minimise.h" #include "mgraph.h" #include #include +#include +#include + static const char *input_file, *output_file; static char *output_file_tmp; static void compute_vertex_areas(const Vertices vertices, double areas[N]); static double best_energy= DBL_MAX; -static void addcost(double *energy, double tweight, double tcost); -#define COST(weight, compute) addcost(&energy, (weight), (compute)) +enum printing_instance { pr_cost, pr_size, pr__max }; + +static void addcost(double *energy, double tweight, double tcost, int pr); +#define COST(weight, compute) addcost(&energy, (weight), (compute), printing) +static int printing_check(enum printing_instance); +static void printing_init(void); /*---------- main energy computation and subroutines ----------*/ static double compute_energy(const Vertices vertices) { double vertex_areas[N], energy; + int printing; compute_vertex_areas(vertices,vertex_areas); energy= 0; - printf("cost > energy |"); - COST(1e4, edgewise_vertex_displacement_cost(vertices)); + printing= printing_check(pr_cost); + + if (printing) printf("cost > energy |"); + +// COST(1e4, edgewise_vertex_displacement_cost(vertices)); COST(1e2, graph_layout_cost(vertices,vertex_areas)); // COST(1e4, noncircular_rim_cost(vertices)); - printf("| total %# e |", energy); + if (printing) printf("| total %# e |", energy); + if (energy < best_energy) { FILE *best_f; int r; - printf(" BEST"); + if (printing) printf(" BEST"); best_f= fopen(output_file_tmp,"wb"); if (!best_f) diee("fopen new out"); r= fwrite(vertices,sizeof(Vertices),1,best_f); if (r!=1) diee("fwrite"); @@ -45,15 +57,17 @@ static double compute_energy(const Vertices vertices) { best_energy= energy; } - putchar('\n'); - flushoutput(); + if (printing) { + putchar('\n'); + flushoutput(); + } return energy; } -static void addcost(double *energy, double tweight, double tcost) { +static void addcost(double *energy, double tweight, double tcost, int pr) { double tenergy= tweight * tcost; - printf(" %# e > %# e |", tcost, tenergy); + if (pr) printf(" %# e > %# e |", tcost, tenergy); *energy += tenergy; } @@ -134,6 +148,7 @@ int main(int argc, const char *const *argv) { if (asprintf(&output_file_tmp,"%s.new",output_file) <= 0) diee("asprintf"); graph_layout_prepare(); + printing_init(); minimiser= gsl_multimin_fminimizer_alloc (gsl_multimin_fminimizer_nmsimplex, DIM); @@ -172,7 +187,8 @@ int main(int argc, const char *const *argv) { size= gsl_multimin_fminimizer_size(minimiser); r= gsl_multimin_test_size(size, stop_epsilon); - printf("%*s size %# e, r=%d\n", 135,"", size, r); + if (printing_check(pr_size)) + printf("%*s size %# e, r=%d\n", 135,"", size, r); flushoutput(); if (r==GSL_SUCCESS) break; @@ -267,3 +283,56 @@ double noncircular_rim_cost(const Vertices vertices) { } return cost; } + +/*---------- printing rate limit ----------*/ + +static volatile unsigned print_todo; +static sigset_t print_alarmset; + +static int printing_check(enum printing_instance which) { + static int skipped[pr__max]; + + unsigned bits= 1u << which; + int sk; + + if (!(print_todo & bits)) { + skipped[which]++; + return 0;; + } + + sigprocmask(SIG_BLOCK,&print_alarmset,0); + print_todo &= ~bits; + sigprocmask(SIG_UNBLOCK,&print_alarmset,0); + + sk= skipped[which]; + if (sk) printf("[%4d] ",sk); + else printf(" "); + skipped[which]= 0; + + return 1; +} + +static void alarmhandler(int ignored) { + print_todo= ~0u; +} + +static void printing_init(void) { + struct sigaction sa; + struct itimerval itv; + + sigemptyset(&print_alarmset); + sigaddset(&print_alarmset,SIGALRM); + + sa.sa_handler= alarmhandler; + sa.sa_mask= print_alarmset; + sa.sa_flags= SA_RESTART; + if (sigaction(SIGALRM,&sa,0)) diee("sigaction ALRM"); + + itv.it_interval.tv_sec= 0; + itv.it_interval.tv_usec= 200000; + itv.it_value= itv.it_interval; + + if (setitimer(ITIMER_REAL,&itv,0)) diee("setitimer REAL"); + + raise(SIGALRM); +} diff --git a/graph.c b/graph.c new file mode 100644 index 0000000..c9cdd2c --- /dev/null +++ b/graph.c @@ -0,0 +1,116 @@ +/* + * graph layout energy + */ + +#include "mgraph.h" +#include "minimise.h" + +static int sqdistances[N][N]; + +static double alpha, beta, beta_prime; + +static void breadth_first_search(int start, int sqdistances_r[N]) { + int d[N], buffer[N], *buf_pop=buffer, *buf_push=buffer; + int v,e, current, future, dfuture; + + buf_push= buf_pop= buffer; + FOR_VERTEX(v) d[v]= -1; + + d[start]= 0; + *buf_push++= start; + + while (buf_pop < buf_push) { + current= *buf_pop++; + dfuture= d[current] + 1; + FOR_VEDGE(current,e,future) { + if (d[future] >= 0) continue; /* already found this one */ + d[future]= dfuture; + *buf_push++= future; + } + } + assert(buf_pop==buf_push); + assert(buf_push <= buffer+sizeof(buffer)/sizeof(buffer[0])); + + FOR_VERTEX(v) { + assert(d[v] >= 0); + sqdistances_r[v]= d[v] * d[v]; + } +} + +void graph_layout_prepare() { + int v1; + + FOR_VERTEX(v1) + breadth_first_search(v1, sqdistances[v1]); + + alpha= 2; + beta= -log(10)/log(alpha); + beta_prime= (1-beta)/2; + printf("alpha=%g beta=%g beta'=%g\n", alpha,beta,beta_prime); +} + + +double graph_layout_cost(const Vertices v, const double vertex_areas[N]) { + /* For each (vi,vj) computes shortest path s_ij = |vi..vj| + * along edges, and actual distance d_ij = |vi-vj|. + * + * We will also use the `vertex areas': for each vertex vi the + * vertex area a_vi is the mean area of the incident triangles. + * This is computed elsewhere. + * + * Energy contribution is proportional to + * + * -4 2 + * a a . d . [ (s/d) - 1 ] + * vi vj + * + * (In practice we compute d^2+epsilon and use it for the + * divisions, to avoid division by zero.) + */ + //static const double d2_epsilon= 1e-6; + + // double edge_weights[V6<0); + + double s2= dist2 * meanedgelength2; + + /* energy = (d/s)^(1-beta) where beta is -log\_{alpha}(10) + * energy = ((d/s)^2) ^ (1-beta)/2 + * let beta' = (1-beta)/2 + */ + + double cost= pow(d2/s2, beta_prime); + + //double s2= s*s + d2_epsilon; + //double sd2= s2 / d2; + //double cost_contrib= a1*a2 * (sd2 - 1) / (d2*d2); + //double cost_contrib= sd2; + + //printf("layout %03x..%03x dist^2=%d s^2=%g d^2=%g " + //" cost+=%g\n", v1,v2, dist2, + // s2,d2, cost); + total_cost += cost; + } + } + return total_cost; +} diff --git a/mgraph.h b/mgraph.h index 12f57ab..19cd601 100644 --- a/mgraph.h +++ b/mgraph.h @@ -64,7 +64,7 @@ #define XBITS 4 #define X (1<