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