#include <gsl/gsl_fit.h>
#include <stdio.h>
#include <math.h>
#include "lib.h"


// Wrapper for the gsl_function structure
//
double function (double * x, size_t dim, void * params)
  {
  return exp(x[0]); 
  }

///////////////////////////////////////////////////////////////////////////////  

// Results printer.
//
void results (size_t calls, double result, double error, double chi)
  {
  if (calls != 0)
    {
    printf ("%6.0e | ", (double)calls);
    }
  printf ("%5f | ", result);
  printf ("%5f | ", error);
  if (chi != 0)
    {
    printf ("%5f", chi);
    }
  } 

///////////////////////////////////////////////////////////////////////////////  

// Perform a fit in order to compare the data with the expected function:
//
//   y = a⋅x^b  →  ln(y) = ln(a) + b⋅ln(x)  →  Y = A + B⋅X
//
// with:
//
//   A = ln(a)  →  a = e^A
//   B = b      →  b = B
//
// using a linear regression. For b, the results is hence compared with the
// expected one:
//
//   b_exp = - 0.5
//
void fit (struct bag full_bag, double* p,
          double* a, double* a_err,
          double* b, double* b_err)
  {
  // Expected value.
  //
  double b_exp = -0.5;

  // Parse arguments.
  //
  double size = full_bag.size;
  double* x = full_bag.pokets.x;
  double* y = full_bag.pokets.y;
  for (size_t i = 0; i < size; i++)
    {
    x[i] = log(x[i]);
    y[i] = log(y[i]);
    }

  // Do fit.
  //
  double A, B, A_err, B_err, AB_cov, sum2;
  gsl_fit_linear(x, 1, y, 1, size, &A, &B, &A_err, &AB_cov, &B_err, &sum2);
  
  // Parse results.
  //
  A_err = sqrt(A_err);
  B_err = sqrt(B_err);
  *a = exp(A);
  *b = B;
  *a_err = *a * A_err;
  *b_err = B_err;

  // Check compatibility with expected values.
  //
  double t = fabs(*b - b_exp)/ *b_err;
  *p = 1 - erf(t/sqrt(2));
  }