// detrended-fluctuation-analysis.c (C) 2017 Claude Heiland-Allen <claude@mathr.co.uk>
// license: CC-BY and/or BSD3 and/or MIT (at your option)
//
// expository detrended fluctuation analysis calculation, after:
// "Detrended Fluctuation Analysis (DFA)"
// http://www.physionet.org/tutorials/fmnc/node5.html
//
// gcc -std=c99 -Wall -Wextra -pedantic -O3 -fopenmp -lm detrended-fluctuation-analysis.c
// ./a.out

#include <assert.h>
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>

// step function from point to using random in [0,1) for chaos game
typedef double (*step_f)(double, double);

double violet(double z, double r)
{
  (void) z;
  static double lastr = 0.5;
  double v = r - lastr;
  lastr = r;
  return v;
}

double white(double z, double r)
{
  (void) z;
  return (r - 0.5);
}

double brown(double z, double r)
{
  return (r - 0.5) + z;
}

// https://graphics.stanford.edu/~seander/bithacks.html#ZerosOnRightLinear
int CTZ(uint32_t v)
{
int c;  // output: c will count v's trailing zero bits,
        // so if v is 1101000 (base 2), then c will be 3
if (v)
{
  v = (v ^ (v - 1)) >> 1;  // Set v's trailing 0s to 1s and zero rest
  for (c = 0; v; c++)
  {
    v >>= 1;
  }
}
else
{
  c = 8 * sizeof(v);
}
  return c;
}

// http://www.firstpr.com.au/dsp/pink-noise/
// James McCartney  Mon, 6 Sep 1999 12:51:50 -0600
double pink(double z, double r)
{
  (void) z;
  (void) r;
    static uint32_t dice[16], total = 0, counter = 0, seed = 0;
    if (seed == 0)
    {
      for (int i = 0; i < 16; ++i)
        total += (dice[i] = (seed = 1664525 * seed + 1013904223));
    }
    int k = CTZ(counter++);         // count trailing zeroes
    k = k & 15;

    // get previous value of this octave
    uint32_t prevrand = dice[k];

    // generate a new random value
    seed = 1664525 * seed + 1013904223;
    // shift to move into mantissa bitfield divided by 16
    uint32_t newrand = seed >> 13;

    // store new value
    dice[k] = newrand;

    // update total
    total += (newrand - prevrand);

    // generate a new random value for the top octave
    seed = 1664525 * seed + 1013904223;
    // shift to move into mantissa bitfield divided by 16
    newrand = seed >> 13;

    // there is a theoretical chance that total + newrand could overflow
    // the mantissa bitfield, but it is *extremely* unlikely due to total
    // being a a gaussian distributed value.

    // convert to floating point value from 2.0 to 4.0 by masking
    union { float f; uint32_t i; } v;
    v.i = (total + newrand) | 0x40000000;

    // subtract 3 to get in range from -1 to +1 and output
    return (v.f - 3.0f);
}

double blue(double z, double r)
{
  (void) z;
  static double lastr = 0.0;
  r = pink(z, r);
  double v = r - lastr;
  lastr = r;
  return v;
}

double black(double z, double r)
{
  return pink(z, r) + z;
}

// compute a time series
// note: assumes 0 is on attractor, should really iterate a bunch to make sure
double *time_series(step_f f, int n)
{
  double *s = malloc(n * sizeof(*s));
  double z = 0;
  for (int i = 0; i < n; ++i)
    s[i] = z = f(z, rand() / (double) RAND_MAX);
  return s;
}

// detrend and integrate a time series
double *detrend(const double *s, int n)
{
  double *t = malloc(n * sizeof(*t));
  double average = 0;
  for (int i = 0; i < n; ++i)
    average += s[n];
  average /= n;
  double x = 0;
  for (int i = 0; i < n; ++i)
    t[i] = (x += s[i] - average);
  return t;
}

// compute the line of best fit y = a x + b for a box y_n(k)
void line_of_best_fit(double *a, double *b, const double *ys, int n)
{
  double xbar = 0;
  double ybar = 0;
  for (int x = 0; x < n; ++x)
  {
    xbar += x;
    ybar += ys[x];
  }
  xbar /= n;
  ybar /= n;
  double num = 0;
  double den = 0;
  for (int x = 0; x < n; ++x)
  {
    double dx = x - xbar;
    double dy = ys[x] - ybar;
    num += dx * dy;
    den += dx * dx;
  }
  *a = num / den;
  *b = ybar - *a * xbar;
}

// compute fluctuations
double *fluctuations(const double *s, int n, int levels)
{
  assert((1 << levels) <= n);
  double *fs = calloc(1, (levels - 2) * sizeof(*fs));
  #pragma omp parallel for schedule(static, 1)
  for (int l = 2; l < levels; ++l)
  {
    int box_size = 1 << l;
    double f = 0;
    for (int box_start = 0; box_start <= n - box_size; box_start += box_size)
    {
      double a, b;
      line_of_best_fit(&a, &b, s + box_start, box_size);
      for (int x = 0; x < box_size; ++x)
      {
        double y = s[box_start + x];
        double yn = a * x + b;
        double dy = y - yn;
        f += dy * dy;
      }
    }
    fs[l - 2] = log2(sqrt(f / n));
  }
  return fs;
}

// compute the slope of an array
double simple_linear_regression(const double *ys, int n)
{
  double a, b;
  line_of_best_fit(&a, &b, ys, n);
  return a;
}

// compute the DFA alpha for a function
double dfa_alpha(step_f f, int n)
{
  int l = round(log2(n));
  while ((1 << l) > n)
    --l;
  double *s = time_series(f, n);
  double *t = detrend(s, n);
  double *fs = fluctuations(t, n, l);
  double a = simple_linear_regression(fs, l - 2);
  free(s);
  free(t);
  free(fs);
  return a;
}

// run the examples
int main()
{
  int d = 24;
  int n = 1 << d; // note: calculation time is O(n log n), O(2^d d)
  printf("white \t%f\n", dfa_alpha(white,  n));
  printf("brown \t%f\n", dfa_alpha(brown,  n));
  printf("violet\t%f\n", dfa_alpha(violet, n));
  printf("pink  \t%f\n", dfa_alpha(pink,   n));
  printf("blue  \t%f\n", dfa_alpha(blue,   n));
  printf("black \t%f\n", dfa_alpha(black,  n));
  return 0;
}

// EOF