- double e1v[D3], e2v[D3], av[D3];
- K {
- e1v[k]= vertices[v1][k] - vertices[v0][k];
- e2v[k]= vertices[v2][k] - vertices[v0][k];
- }
- xprod(av, e1v, e2v);
- total += magnD(av);
- count++;
- }
- areas[v0]= total / count;
- }
-}
-
-/*---------- use of GSL ----------*/
-
- /* We want to do multidimensional minimisation.
- *
- * We don't think there are any local minima. Or at least, if there
- * are, the local minimum which will be found from the starting
- * state is the one we want.
- *
- * We don't want to try to provide a derivative of the cost
- * function. That's too tedious (and anyway the polynomial
- * approximation to our our cost function sometimes has high degree
- * in the inputs which means the quadratic model implied by most of
- * the gradient descent minimisers is not ideal).
- *
- * This eliminates most of the algorithms. Nelder and Mead's
- * simplex algorithm is still available and we will try that.
- *
- * In our application we are searching for the optimal locations of
- * N actualvertices in D3 (3) dimensions - ie, we are searching for
- * the optimal metapoint in an N*D3-dimensional space.
- *
- * So eg with X=Y=100, the simplex will contain 300 metavertices
- * each of which is an array of 300 doubles for the actualvertex
- * coordinates. Hopefully this won't be too slow ...
- */
-
-static gsl_multimin_fminimizer *minimiser;
-
-static const double stop_epsilon= 1e-4;
-
-static double minfunc_f(const gsl_vector *x, void *params) {
- assert(x->size == DIM);
- assert(x->stride == 1);
- return compute_energy((const double(*)[D3])x->data);
-}
-
-int main(int argc, const char *const *argv) {
- gsl_multimin_function multimin_function;
- double size;
- Vertices initial, step_size;
- FILE *initial_f;
- gsl_vector initial_gsl, step_size_gsl;
- int r, v, k;
-
- if (argc!=3 || argv[1][0]=='-' || strncmp(argv[2],"-o",2))
- { fputs("usage: minimise <input> -o<output\n",stderr); exit(8); }
-
- input_file= argv[1];
- output_file= argv[2]+2;
- if (asprintf(&output_file_tmp,"%s.new",output_file) <= 0) diee("asprintf");