From 2c2c1cc2745430048709dc3b11aef7071d802c86 Mon Sep 17 00:00:00 2001
From: Ian Jackson <ian@davenant.relativity.greenend.org.uk>
Date: Sun, 30 Dec 2007 12:08:18 +0000
Subject: [PATCH] starting to make it compile

---
 Makefile | 12 ++++++++++++
 bgl.h    |  2 ++
 common.c |  2 +-
 common.h |  5 +++++
 energy.c | 48 +++++++++++++++++++++++++++++++++++++++++++-----
 mgraph.h |  6 ++++--
 6 files changed, 67 insertions(+), 8 deletions(-)
 create mode 100644 Makefile

diff --git a/Makefile b/Makefile
new file mode 100644
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 68700b9..d9c5d6b 100644
--- 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 146e5d5..47edc4d 100644
--- 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 1a358d1..77589f9 100644
--- 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 c3b8754..139e7e1 100644
--- 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 0991a7f..4638367 100644
--- 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*/
-- 
2.30.2