#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);
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) {
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); }
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); }
}
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;
+ }
+}