#include <stddef.h>

#include "grad_desc_mod.h"

gdm_sys * gdm_alloc(unsigned int n)
{
  if (n == 0) {
    fprintf(stderr, "%s: %i: error: Nelze vytvořit system, ktery ma 0 parametru.", __FILE__, __LINE__);
    exit(1);
  }

  unsigned int i;
  gdm_sys *sys = (gdm_sys *) malloc(sizeof(gdm_sys));
  sys->param = (double **) malloc(sizeof(double *)*n);
  sys->param_count = n;
  sys->gradStep = 0.000001;
  sys->data = NULL;
  sys->error = NULL;
  sys->lerning_rate = 0.1;
  sys->grad_type = both;
  sys->logFile = stdout;
  sys->maxError = 0.000001;
  sys->maxIter = 1000;

  for (i = 0; i<n; i++) {
    sys->param[i] == NULL;
  }

  return sys;
}

void gdm_free(gdm_sys *sys)
{
  free(sys->param);
  free(sys);
}

void gdm_define_param(gdm_sys *sys, unsigned int i, double *pparam)
{
  if (i>=sys->param_count) {
    fprintf(stderr, "%s: %i: error: Nelze definovat parametr s indexem (%i) vetsi nez pocet celkovych parametru (%i).", __FILE__, __LINE__, i, sys->param_count);
    exit(1);
  }

  sys->param[i] = pparam;
}

/*
 * Nastavi I-ty parametr na hodnotu VALUE systemu SYS.
 */
void gdm_set_param(gdm_sys *sys, unsigned int i, double value)
{
  if (i>=sys->param_count) {
    fprintf(stderr, "%s: %i: error: Nelze nastavit hodnotu parametru s indexem (%i) vetsi nez pocet celkovych parametru (%i).", __FILE__, __LINE__, i, sys->param_count);
    exit(1);
  }

  *(sys->param[i]) = value;
}

/*
 * Vrati hodnotu I-teho parametru systemu SYS.
 */
double gdm_get_param(gdm_sys *sys, unsigned int i)
{
  if (i>=sys->param_count) {
    fprintf(stderr, "%s: %i: error: Nelze ziskat hodnotu parametru s indexem (%i) vetsi nez pocet celkovych parametru (%i).", __FILE__, __LINE__, i, sys->param_count);
    exit(1);
  }

  return *(sys->param[i]);
}

double gdm_rand(void)
{
  return ((double) rand())/((double) RAND_MAX);
}

double gdm_partial_derivation (gdm_sys *sys, unsigned int i)
{
  double value, Ea, Eb;

  if (sys->grad_type == left) {
    Eb = sys->E_middle;
    value = gdm_get_param(sys, i);
    gdm_set_param(sys, i, value - sys->gradStep);
    Ea = sys->error(sys->data);
    gdm_set_param(sys, i, value);
    return (Eb-Ea)/sys->gradStep;
  }
  
  if (sys->grad_type == right) {
    Ea = sys->E_middle;
    value = gdm_get_param(sys, i);
    gdm_set_param(sys, i, value + sys->gradStep);
    Eb = sys->error(sys->data);
    gdm_set_param(sys, i, value);
    return (Eb-Ea)/sys->gradStep;
  }

  if (sys->grad_type == both) {
    value = gdm_get_param(sys, i);
    gdm_set_param(sys, i, value + sys->gradStep);
    Eb = sys->error(sys->data);
    gdm_set_param(sys, i, value - sys->gradStep);
    Ea = sys->error(sys->data);
    return (Eb-Ea)/(2*sys->gradStep);
  }

}

void gdm_optim (gdm_sys *sys)
{
  unsigned int k, n, iter = 1;
  double error, dp, p, pav, pstdev, t, t_i;
  double gradient[sys->param_count];
  time_t s, e;
  clock_t st, en;

  fprintf(sys->logFile, "lerning rate : %f\n", sys->lerning_rate);
  fprintf(sys->logFile, "maximal iter : %i\n", sys->maxIter);
  fprintf(sys->logFile, "maximal error : %f\n", sys->maxError);
  fprintf(sys->logFile, "param probability : %f\n", sys->param_prob);
  fflush(sys->logFile);
  error = sys->error(sys->data);
  fseek(sys->logFile, 0, SEEK_END);
  fprintf(sys->logFile, "0. E = %f\n", error);
  fflush(sys->logFile);
  sys->E_middle = error;

  while ((iter <= sys->maxIter) & (error > sys->maxError)) {

    time(&s);
    st = clock();

    // Vypocet gradientu
    n = 0;
    for (k = 0; k < sys->param_count; k++) {

      if (gdm_rand() <= sys->param_prob) {
        gradient[k] = gdm_partial_derivation(sys, k);
        n++;
      } else {
        gradient[k] = 0;
      }

    }

    // Uprava parametru
    dp = 0;
    for (k = 0; k < sys->param_count; k++) {

      dp += pow(gradient[k], 2);
      p = gdm_get_param(sys, k) - gradient[k] * sys->lerning_rate;
      gdm_set_param(sys, k, p);

    }

    time(&e);
    en = clock();
    error = sys->error(sys->data);
    t = (((double) en)-((double) st))/((double) CLOCKS_PER_SEC);
    if (n>0) {
      t_i = t/((double) n);
    } else {
      t_i = 0;
    }

    fseek(sys->logFile, 0, SEEK_END);
    fprintf(sys->logFile, "%i. E = %f  |dp| = %f  n = %i  %.2fs time/iter = %f\n", iter, error, sqrt(dp), n, t, t/((double) n));
    fflush(sys->logFile);
    sys->E_middle = error;
    
    iter++;
    
  }
  
}