starting to make it compile

[moebius2.git] / energy.c

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;
+  }
+}