test/lorentzfit.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <math.h>
   4
   5 #include <nlopt.h>
   6
   7 typedef struct {
   8      int N;
   9      double *x, *y; /* length N; */
  10 } lorentzdata;
  11
  12 static double sqr(double x) { return x * x; }
  13
  14 static int count = 0;
  15
  16 static double lorentzerr(int n, const double *p, double *grad, void *data)
  17 {
  18      lorentzdata *d = (lorentzdata *) data;
  19      int N = d->N;
  20      const double *xs = d->x;
  21      const double *ys = d->y;
  22      double val = 0;
  23      int i, j;
  24
  25      for (i = 0; i < N; ++i) {
  26           double x = xs[i], y = ys[i];
  27           double lorsum = 0;
  28
  29           for (j = 0; j < n; j += 3) {
  30                double A = p[j + 0];
  31                double w = p[j + 1];
  32                double G = p[j + 2];
  33                double lor = A / (sqr(x - w) + G*G);
  34
  35                lorsum += lor;
  36           }
  37
  38           val += sqr(y - lorsum);
  39
  40           if (grad)
  41                for (j = 0; j < n; j += 3) {
  42                     double A = p[j + 0];
  43                     double w = p[j + 1];
  44                     double G = p[j + 2];
  45                     double deninv =  1.0 / (sqr(x - w) + G*G);
  46
  47                     grad[j + 0] += -2 * (y - lorsum) * deninv;
  48                     grad[j + 1] += 4*A * (w - x) * (y - lorsum) * sqr(deninv);
  49                     grad[j + 2] += 4*A * G * (y - lorsum) * sqr(deninv);
  50                }
  51      }
  52      ++count;
  53      // printf("%d: f(%g,%g,%g) = %g\n", count, p[0],p[1],p[2], val);
  54      return val;
  55 }
  56
  57 extern double nlopt_urand(double a, double b);
  58
  59 int main(void)
  60 {
  61      lorentzdata d;
  62      int i;
  63      double A = 1, w = 0, G = 1, noise=0.01;
  64      double lb[3] = {-HUGE_VAL,-HUGE_VAL,0};
  65      double ub[3] = {HUGE_VAL,HUGE_VAL,HUGE_VAL};
  66      double p[3] = {0,1,2}, minf;
  67
  68      nlopt_srand_time();
  69
  70      d.N = 200;
  71      d.x = (double *) malloc(sizeof(double) * d.N * 2);
  72      d.y = d.x + d.N;
  73      for (i = 0; i < d.N; ++i) {
  74           d.x[i] = nlopt_urand(-0.5, 0.5) * 8*G + w;
  75           d.y[i] = 2*noise*nlopt_urand(-0.5,0.5) + A / (sqr(d.x[i]-w) + G*G);
  76      }
  77
  78      nlopt_minimize(NLOPT_LN_NEWUOA_BOUND, 3, lorentzerr, &d,
  79                     lb, ub, p, &minf,
  80                     -HUGE_VAL, 0,0, 1e-6,NULL, 0,0);
  81
  82      printf("%d minf=%g at A=%g, w=%g, G=%g\n", count, minf, p[0],p[1],p[2]);
  83
  84      count = 0;
  85      nlopt_minimize(NLOPT_LN_COBYLA, 3, lorentzerr, &d,
  86                     lb, ub, p, &minf,
  87                     -HUGE_VAL, 0,0, 1e-6,NULL, 0,0);
  88
  89      printf("%d minf=%g at A=%g, w=%g, G=%g\n", count, minf, p[0],p[1],p[2]);
  90
  91      count = 0;
  92      nlopt_minimize(NLOPT_LN_NELDERMEAD, 3, lorentzerr, &d,
  93                     lb, ub, p, &minf,
  94                     -HUGE_VAL, 0,0, 1e-6,NULL, 0,0);
  95
  96      printf("%d minf=%g at A=%g, w=%g, G=%g\n", count, minf, p[0],p[1],p[2]);
  97
  98      count = 0;
  99      nlopt_minimize(NLOPT_LN_SBPLX, 3, lorentzerr, &d,
 100                     lb, ub, p, &minf,
 101                     -HUGE_VAL, 0,0, 1e-6,NULL, 0,0);
 102
 103      printf("%d minf=%g at A=%g, w=%g, G=%g\n", count, minf, p[0],p[1],p[2]);
 104
 105      return 0;
 106 }