starting to make it compile

diff --git a/Makefile b/Makefile
new file mode 100644 (file)
index 0000000..9e3aa83
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,12 @@
+
+TARGETS=minimise
+
+CFLAGS=        -Wall -Wwrite-strings -Wpointer-arith -Werror
+
+all:           $(TARGETS)
+
+minimise:      energy.o bgl.o common.o mgraph.o
+
+clean:
+               rm -f *.o $(TARGETS) *.new
+               rm -f best initial

diff --git a/bgl.h b/bgl.h
index 68700b98f611708361ea9bcc24ff4db49499d72f..d9c5d6be86d6e26b5cdcdd246c47b0da8606d43a 100644 (file)
--- a/bgl.h
+++ b/bgl.h
@@ -5,6 +5,8 @@
 #ifndef BGL_H
 #define BGL_H
 
+#include "mgraph.h"
+
 double graph_layout_cost(const Vertices v, const double vertex_areas[N]);
 
 #endif /*BGL_H*/

diff --git a/common.c b/common.c
index 146e5d54ef44a3f6fd9fc3f9879d1bc5d1fc7c34..47edc4ddbe0f6be0e5c275aa8192f190ecaf1967 100644 (file)
--- a/common.c
+++ b/common.c
@@ -4,7 +4,7 @@
 
 #include "common.h"
 
-double hypotD1(const double pq[D3]) {
+double magnD(const double pq[D3]) {
   gsl_vector v;
 
   v.size= D3;

diff --git a/common.h b/common.h
index 1a358d1d4033db9036c41b14abf565f66ecea73a..77589f98c45dcbba5770130e0d13fd317a14ab60 100644 (file)
--- a/common.h
+++ b/common.h
@@ -7,12 +7,17 @@
 
 #define _GNU_SOURCE
 #include <math.h>
+#include <float.h>
 #include <limits.h>
 
+#define D3 3
+
 double hypotD(const double p[D3], const double q[D3]);
 double hypotD2(const double p[D3], const double q[D3]);
 double hypotD2plus(const double p[D3], const double q[D3], double add);
 
+double magnD(const double pq[D3]);
+
 #ifdef FP_FAST_FMA
 # define fma_fast fma
 #else

diff --git a/energy.c b/energy.c
index c3b8754aff63dfa661cfd5dcbbbc41e67ae36885..139e7e1304d11937fe5d5261e6dba81d9ad55c3e 100644 (file)
--- a/energy.c
+++ b/energy.c
@@ -10,9 +10,10 @@
 #define INITIAL_F "initial"
 
 static double edgewise_vertex_displacement_cost(const Vertices vertices);
+static double noncircular_rim_cost(Vertices vertices);
 
 static void compute_vertex_areas(const Vertices vertices, double areas[N]);
-static double best_energy= DOUBLE_MAX;
+static double best_energy= DBL_MAX;
 static void flushoutput(void);
 
 static void cost(double *energy, double tweight, double tcost);
@@ -29,7 +30,7 @@ static double compute_energy(Vertices vertices) {
 
   COST(1000.0, edgewise_vertex_displacement_cost(vertices));
   COST(1.0,    graph_layout_cost(vertices,vertex_areas));
-  COST(1e6,    noncircular_edge_cost(vertices));
+  COST(1e3,    noncircular_rim_cost(vertices));
   
   printf("| total %# e |", energy);
   if (energy < best_energy) {
@@ -126,9 +127,9 @@ static double minfunc_f(const gsl_vector *x, void *params) {
 int main(int argc, const char *const *argv) {
   struct gsl_multimin_function multimin_function;
   double size;
-  Vertices initial;
+  Vertices initial, step_size;
   FILE *initial;
-  gsl_vector initial_gsl, *step_size;
+  gsl_vector initial_gsl, step_size_gsl;
   int r;
   
   if (argc>1) { fputs("takes no arguments\n",stderr); exit(8); }
@@ -148,13 +149,29 @@ int main(int argc, const char *const *argv) {
 
   initial_gsl.size= DIM;
   initial_gsl.stride= 1;
-  initial_gsl.data= initial;
   initial_gsl.block= 0;
   initial_gsl.owner= 0;
+  step_size_gsl= initial_gsl;
+
+  initial_gsl.data= initial;
+  step_size_gsl.data= step_size;
+
+  FOR_VERTEX(v)
+    K step_size[v][k]= 1e-3;
+  FOR_RIM_VERTEX(vx,vy,v)
+    step_size[v][3] *= 0.1;
+  
+  for (vy=0; vy<Y; vy+=Y-1)
+  for (vx=0; vx<x
+  for (i=0; i<DIM; i++) step_size[i]= step_size;
+  
 
   step_size= gsl_vector_alloc(DIM);  if (!step_size) gsldie("alloc step");
   gsl_vector_set_all(step_size, 1e-3);
 
+  assert(step_
+  step_
+
   r= gsl_multimin_fminimizer_set(minimiser, &multimin_function,
                                 &initial_gsl, &step_size);
   if (r) { gsldie("fminimizer_set",r); }
@@ -261,3 +278,24 @@ static double edgewise_vertex_displacement_cost(const Vertices vertices) {
   }
   return sigma_bd2_l3;
 }
+
+/*---------- noncircular rim cost ----------*/
+
+static double noncircular_rim_cost(Vertices vertices) {
+  int vy,vx,v;
+  double cost= 0.0;
+  
+  FOR_RIM_VERTEX(vy,vx,v) {
+    double oncircle[3];
+    /* By symmetry, nearest point on circle is the one with
+     * the same angle subtended at the z axis. */
+    oncircle[0]= vertices[v][0];
+    oncircle[1]= vertices[v][1];
+    oncircle[2]= 0;
+    double mult= 1.0/ magnD(oncircle);
+    oncircle[0] *= mult;
+    oncircle[1] *= mult;
+    double d2= hypotD2(vertices[v], oncircle);
+    cost += d2*d2;
+  }
+}

diff --git a/mgraph.h b/mgraph.h
index 0991a7ff98ff0e92da1dbc6061a3b873c69ae8d9..46383677d2ec6c153e799d5936329529d3e75de2 100644 (file)
--- a/mgraph.h
+++ b/mgraph.h
@@ -69,8 +69,6 @@
 
 #define V6 6
 
-#define D3 3
-
 #define FOR_VERTEX(v) \
   for ((v)=0; (v)<N; (v)++)
 
@@ -103,6 +101,10 @@ extern int edge_end2(unsigned v1, int e);
 
 #define K FOR_COORD(k)
 
+#define FOR_RIM_VERTEX(vy,vx,v)                                        \
+  for ((vy)=0; (vy)<Y; (vy)+=Y-1)                              \
+    for ((vx)=0; (v)= (vy)<<YSHIFT | (vx), (vx)<X; (vx)++)
+
 typedef double Vertices[N][D3];
 
 #endif /*MGRAPH_H*/