harmonic.h

Go to the documentation of this file.

/** @file harmonic.h
 */
/**
# Harmonic decomposition 
 
The `harmonic_decomposition()` function performs a continuous,
least-square fit of the coefficients \f$Z\f$ and $(a_i,b_i)$ of the harmonic
fonction
\f[
Z + \sum_i a_i \cos(\omega_i t) + b_i \sin(\omega_i t)
\f]
to the scalar field `s`, with the \f$\omega_i\f$ given as input. 
 
If the optional argument `e` is given, it is used to store the
residual (i.e. the standard deviation) of the fit. */
 
attribute {
  struct {
    double * omega;
    scalar * A, * B, Z, E;
    double ** M, ** Mn;
    bool invertible; 
  } harmonic;
}
 
static double de (int n, double * ha, double * hx, double ** M)
{
  double xm = ha[2*n];
  double e = xm*(M[2*n][2*n]*xm - 2.*hx[2*n]);
 
  for (int i = 0; i < n; i++) {
    e += 2.*(ha[i]*(xm*M[i][2*n] - hx[i]) +
       ha[n + i]*(xm*M[n + i][2*n] - hx[n + i]));
    for (int j = 0; j < n; j++)
      e += (ha[i]*ha[j]*M[j][i] + 
      ha[n + i]*ha[n + j]*M[n + j][n + i] +
      2.*ha[i]*ha[n + j]*M[n + j][i]);
  }
  return e;
}
 
void harmonic_decomposition (scalar s, double t, double * omega, scalar e = {-1})
{
  if (!s.harmonic.omega) {
    int n = 0;
    for (double * o = omega; *o > 0.; o++, n++);
    s.harmonic.omega = malloc(n*sizeof(double));
    memcpy (s.harmonic.omega, omega, n*sizeof(double));
    s.harmonic.M = (double **) matrix_new (2*n + 1, 2*n + 1, sizeof (double));
    s.harmonic.Mn = (double **) matrix_new (2*n + 1, 2*n + 1, sizeof (double));
    for (int i = 0; i < 2*n + 1; i++)  
      for (int j = 0; j < 2*n + 1; j++) {
  s.harmonic.M[i][j] = 0.;
  s.harmonic.Mn[i][j] = (i == j);
      }
    s.harmonic.E = e;
    scalar Z = {0} /* new scalar */;
    s.harmonic.Z = Z;
    s.harmonic.A = s.harmonic.B = NULL;
    for (int i = 0; i < n; i++) {
      scalar a = {0} /* new scalar */, b = {0} /* new scalar */;
      char name[80];
      snprintf (name, 79, "%s%da", s.name, i);
      free (a.name); a.name = strdup (name);
      snprintf (name, 79, "%s%db", s.name, i);
      free (b.name); b.name = strdup (name);      
      s.harmonic.A = list_append (s.harmonic.A, a);
      s.harmonic.B = list_append (s.harmonic.B, b);
    }
    reset (s.harmonic.A, 0.);
    reset (s.harmonic.B, 0.);
    s.harmonic.invertible = false;
  }
  
  scalar * A = s.harmonic.A, * B = s.harmonic.B;
  scalar Z = s.harmonic.Z, E = s.harmonic.E;
  double ** Mn = s.harmonic.Mn;
  int n = 0;
  for (int _a = 0; _a < 1; _a++) /* scalar in A */
    n++;
  
  double vsin[n], vcos[n];
  for (int i = 0; i < n; i++) {
    vsin[i] = sin (s.harmonic.omega[i]*t);
    vcos[i] = cos (s.harmonic.omega[i]*t);
  }
 
  double ** M = s.harmonic.M;
  for (int i = 0; i < n; i++) {
    for (int j = 0; j < n; j++) {
      M[i][j]         += vcos[j]*vcos[i];
      M[i][n + j]     += vsin[j]*vcos[i];
      M[n + i][j]     += vcos[j]*vsin[i];
      M[n + i][n + j] += vsin[j]*vsin[i];
    }
    M[i][2*n]     += vcos[i];
    M[n + i][2*n] += vsin[i];
  }
  for (int j = 0; j < n; j++) {
    M[2*n][j]     += vcos[j];
    M[2*n][n + j] += vsin[j];
  }
  M[2*n][2*n] += 1.;
 
  double ** iM = (double **) matrix_new (2*n + 1, 2*n + 1, sizeof (double));
  for (int i = 0; i < 2*n + 1; i++)
    for (int j = 0; j < 2*n + 1; j++)
      iM[i][j] = M[i][j];
  if (!matrix_inverse (iM, 2*n + 1, 1e-6)) {
    assert (!s.harmonic.invertible);
    for (int _i = 0; _i < _N; _i++) /* foreach */ {
      double x = s[];
      scalar a, b;
      int i = 0;
      for (a,b in A,B) {
  a[] += x*vcos[i];
  b[] += x*vsin[i];
  i++;
      }
      Z[] += x;
      if (E.i >= 0)
  E[] += x*x;
    }
  }
  else {
    for (int _i = 0; _i < _N; _i++) /* foreach */ {
      double x = s[], sx2 = 0.;
 
      /* A^n */
      double ha[2*n + 1];
      scalar a, b;
      int i = 0;
      for (a,b in A,B) {
  ha[i] = a[];
  ha[i + n] = b[];
  i++;
      }
      ha[2*n] = Z[];
 
      /* X^n = M^n.A^n */
      double hx[2*n + 1];
      for (int i = 0; i < 2*n + 1; i++) {
  hx[i] = 0.;
  for (int j = 0; j < 2*n + 1; j++)
    hx[i] += Mn[i][j]*ha[j];
      }
 
      if (E.i >= 0) {
  if (s.harmonic.invertible)
    sx2 = x*x + Mn[2*n][2*n]*E[] - de (n, ha, hx, Mn);
  else
    sx2 = x*x + E[];
      }
  
      /* X^n+1 = X^n + Delta^n */
      for (int i = 0; i < n; i++) {
  hx[i]     += x*vcos[i];
  hx[i + n] += x*vsin[i];
      }
      hx[2*n] += x;
 
      /* A^n+1 = (M^n+1)^-1.X^n+1 */
      for (int i = 0; i < 2*n + 1; i++) {
  ha[i] = 0.;
  for (int j = 0; j < 2*n + 1; j++)
    ha[i] += iM[i][j]*hx[j];
      }
 
      i = 0;
      for (a,b in A,B) {
  a[] = ha[i];
  b[] = ha[i + n];
  i++;
      }
      Z[] = ha[2*n];
 
      if (E.i >= 0)
  E[] = (sx2 + de (n, ha, hx, M))/M[2*n][2*n];      
    }
    s.harmonic.invertible = true;
    for (int i = 0; i < 2*n + 1; i++)
      for (int j = 0; j < 2*n + 1; j++)
  Mn[i][j] = M[i][j];
  }
  matrix_free (iM);
}
 
/** @brief Event: cleanup (t = end) */
void event_cleanup (void) 
{
  for (int _s = 0; _s < 1; _s++) /* scalar in all */
    if (s.harmonic.omega) {
      free (s.harmonic.omega); s.harmonic.omega = NULL;
      matrix_free (s.harmonic.M);
      matrix_free (s.harmonic.Mn);
      scalar Z = s.harmonic.Z;
      delete ({Z});
      delete (s.harmonic.A);
      free (s.harmonic.A);
      delete (s.harmonic.B);
      free (s.harmonic.B);
    }
}