findcurve.c

   1 /*
   2  */
   3
   4 #include <math.h>
   5 #include <string.h>
   6
   7 #include <gsl/gsl_errno.h>
   8 #include <gsl/gsl_sf.h>
   9 #include <gsl/gsl_siman.h>
  10 #include <gsl/gsl_randist.h>
  11
  12 #include "symbolic.c"
  13
  14 #define J_END_COL(i) \
  15   for (j=0; j<PN; j++) gsl_matrix_set(J,i,j,J_COL[j]);
  16
  17 static inline _Bool IS_SMALL(double v) {
  18   return v < 1E6;
  19 }
  20
  21 static inline double sinc(double x) {
  22   return gsl_sf_sinc(x / M_PI);
  23 }
  24
  25 static int NP;
  26 static double *INPUT; /* dyanmic array, on main's stack */
  27 static double PREP[NPREP];
  28
  29 static void printcore(const double *X) {
  30   int i, j;
  31   DECLARE_F_G;
  32   CALCULATE_F_G;
  33   printf("[");
  34   for (i=0; i<NP; i++)
  35     for (j=0; j<3; j++)
  36       printf(" %25.18g,", POINT(i)[j]);
  37   printf(" ]\n");
  38 }
  39
  40 static void prepare(double X[] /* startpoint */) {
  41   /* fills in PREP and startpoint */
  42   PREPARE;
  43 }
  44
  45 static double cb_Efunc(void *xp) {
  46   const double *X = xp;
  47   DECLARE_F_G;
  48   CALCULATE_F_G;
  49
  50   printf(" Efunc ");
  51   printcore(X);
  52
  53   double e = 0;
  54   int P;
  55   for (P=0; P<NP-3; P++) {
  56     double P_cost;
  57     CALCULATE_COST;
  58     printf(" %.18g", P_cost);
  59     e += P_cost;
  60   }
  61   printf("\n");
  62   return e;
  63 }
  64
  65 static void cb_step(const gsl_rng *rng, void *xp, double step_size) {
  66   double *x = xp;
  67   int i;
  68   double step[NX];
  69   gsl_ran_dir_nd(rng, NX, step);
  70   for (i=0; i<NX; i++)
  71     x[i] += step_size * step[i];
  72   //printf("\n cb_step %p %10.7f [", xp, step_size);
  73   //for (i=0; i<N; i++) printf(" %10.7f,", step[i]);
  74   //printf("]\n");
  75 }
  76
  77 static double cb_metric(void *xp, void *yp) {
  78   const double *x=xp, *y=yp;
  79   int i;
  80   double s;
  81   for (i=0, s=0; i<NX; i++) {
  82     double d = x[i] - y[i];
  83     s += d*d;
  84   }
  85   return sqrt(s);
  86 }
  87
  88 static void __attribute__((unused)) cb_print(void *xp) {
  89   const double *x = xp;
  90   printf("\n");
  91   printcore(x);
  92 }
  93
  94 static double scan1double(void) {
  95   double v;
  96   int r;
  97
  98   r = scanf("%lf",&v);
  99   if (ferror(stdin)) { perror("stdin"); exit(-1); }
 100   if (feof(stdin)) exit(0);
 101   if (r!=1) { fputs("bad input\n",stderr); exit(-1); }
 102   return v;
 103 }
 104
 105 int main(int argc, const char *const *argv) {
 106   double epsilon;
 107   int i;
 108
 109   NP = atoi(argv[1]);
 110   epsilon = atof(argv[2]);
 111
 112   gsl_rng *rng = gsl_rng_alloc(gsl_rng_ranlxd2);
 113
 114   double input[NINPUT]; INPUT = input;
 115   double startpoint[NX];
 116
 117   for (;;) {
 118     /* NINPUT + 1 doubles: startpoint, epsilon for residual */
 119     for (i=0; i<NINPUT; i++)
 120       INPUT[i] = scan1double();
 121
 122     gsl_rng_set(rng,0);
 123
 124     gsl_siman_params_t siman_params = {
 125       .k = 1.0,
 126       .t_initial = 0.5,
 127       .mu_t = 1.001,
 128       .t_min = epsilon * 1E-3,
 129       .iters_fixed_T = 100,
 130       .n_tries = 10,
 131       .step_size = 0.05,
 132     };
 133
 134     prepare(startpoint);
 135
 136     gsl_siman_solve(rng,
 137                     startpoint,
 138                     cb_Efunc, cb_step, cb_metric,
 139                     0, // cb_print,
 140                     0,0,0, sizeof(startpoint), siman_params);
 141
 142     printcore(startpoint);
 143
 144     printf("[]\n");
 145     fflush(stdout);
 146   }
 147 }