/*
   This program should identify stars with large sigma.
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_errno.h>

#include "vast_limits.h"

#include "stetson_variability_indexes.h"

#include "detailed_error_messages.h"

int main(){
 FILE *sigma_selection_curve_log;
 FILE *dmsf;
 double m,sigma,X,Y,mmax,mean;
 double *data_sigma;
 double *x_sigma;
 double *y;
 double *y_limit_sigma;
 y_limit_sigma=malloc(sizeof(double)*MAX_NUMBER_OF_STARS);
 x_sigma=malloc(sizeof(double)*MAX_NUMBER_OF_STARS);
 y=malloc(sizeof(double)*MAX_NUMBER_OF_STARS);
 data_sigma=malloc(sizeof(double)*MAX_NUMBER_OF_STARS);
 char str[256];
 int n=0;
 int i=0;
 int n_drop_high_sigma_stars;

 double *data_modified_sigma;
 double *data_I;
 double *data_L;
 double modified_sigma_series,skewness,kurtosis,I,J,K,L;
 
 double *y_limit_modified_sigma;
 double *y_limit_I;
 double *y_limit_L;

 double *w;

 double lag1_autocorrelation,RoMS,Chi2,I_sign_only,AGN_v;
 
 data_modified_sigma=malloc(sizeof(double)*MAX_NUMBER_OF_STARS);
 data_I=malloc(sizeof(double)*MAX_NUMBER_OF_STARS);
 data_L=malloc(sizeof(double)*MAX_NUMBER_OF_STARS);
 
 y_limit_modified_sigma=malloc(sizeof(double)*MAX_NUMBER_OF_STARS);
 y_limit_I=malloc(sizeof(double)*MAX_NUMBER_OF_STARS);
 y_limit_L=malloc(sizeof(double)*MAX_NUMBER_OF_STARS);
 
 w=malloc(sizeof(double)*MAX_NUMBER_OF_STARS);

 int n_modified_sigma,n_I,n_L;

 double * m_arr;

 double MAD_scaled_to_sigma;
 int points_in_lightcurve;
 
 int interpolation_status_gsl=0;
 
 double double_tmp;
 
 m_arr=malloc(sizeof(double)*MAX_NUMBER_OF_STARS);

 dmsf=fopen("vast_lightcurve_statistics.log","r");
 if( dmsf==NULL ){
  fprintf(stderr,"ERROR: Can't open file vast_lightcurve_statistics.log !\n");
  report_lightcurve_statistics_computation_problem();
  exit(1);
 }

 
 i=0;n=0;mmax=0.0;
 while(-1<fscanf(dmsf,"%lf %lf %lf %lf %s  %lf  %lf %lf  %lf %lf %lf %lf %d %lf %lf %lf %lf %lf %lf  %lf %lf %lf %lf  %lf  %lf  %lf %lf %lf  %lf",&m,&sigma,&X,&Y,str, &modified_sigma_series,&skewness,&kurtosis,&I,&J,&K,&L,&points_in_lightcurve,&MAD_scaled_to_sigma,&lag1_autocorrelation,&RoMS,&Chi2,&I_sign_only,&AGN_v, &double_tmp, &double_tmp, &double_tmp, &double_tmp, &double_tmp, &double_tmp, &double_tmp, &double_tmp, &double_tmp, &double_tmp)){

  if( mmax==0.0 )mmax=m;

  if( m<=mmax+M_SIGMA_BIN_SIZE_M || n<12 ){
   m_arr[n]=m;
   data_sigma[n]=sigma;
   data_modified_sigma[n]=modified_sigma_series;
   data_I[n]=I;
   data_L[n]=L;
   n++;
  }
  else{
   //
   n_modified_sigma=n;
   gsl_sort2(data_modified_sigma, 1, m_arr, 1, n_modified_sigma);
   n_drop_high_sigma_stars=MAX((int)(0.07*n_modified_sigma),M_SIGMA_BIN_DROP);
   if( n_modified_sigma>2*n_drop_high_sigma_stars+2 )
    n_modified_sigma-=2*n_drop_high_sigma_stars;
   else
    if( n_modified_sigma>n_drop_high_sigma_stars+2 )
     n_modified_sigma-=n_drop_high_sigma_stars;
   mean=gsl_stats_median_from_sorted_data(data_modified_sigma,1,n_modified_sigma);
   
   y_limit_modified_sigma[i]=esimate_sigma_from_MAD_of_sorted_data(data_modified_sigma,n_modified_sigma);

   if( 0!=isnan(y_limit_modified_sigma[i]) ){
    y_limit_modified_sigma[i]=0;
   }
   y_limit_modified_sigma[i]=mean+M_SIGMA_BIN_MAG_SIGMA_DETECT*y_limit_modified_sigma[i];
   //
   //
   n_I=n;
   gsl_sort2(data_I, 1, m_arr, 1, n_I);
   //
   n_drop_high_sigma_stars=MAX((int)(0.07*n_I),M_SIGMA_BIN_DROP);
   if( n_I>2*n_drop_high_sigma_stars+2 )
    n_I-=2*n_drop_high_sigma_stars;
   else
    if( n_I>n_drop_high_sigma_stars+2 )
     n_I-=n_drop_high_sigma_stars;
   mean=gsl_stats_median_from_sorted_data(data_I,1,n_I);
   y_limit_I[i]=esimate_sigma_from_MAD_of_sorted_data(data_I,n_I);
   if( 0!=isnan(y_limit_I[i]) ){
    y_limit_I[i]=0;
   }
   y_limit_I[i]=mean+M_SIGMA_BIN_MAG_SIGMA_DETECT*y_limit_I[i];
   //
   n_L=n;
   gsl_sort2(data_L, 1, m_arr, 1, n_L);
   //
   n_drop_high_sigma_stars=MAX((int)(0.07*n_L),M_SIGMA_BIN_DROP);
   if( n_L>2*n_drop_high_sigma_stars+2 )
    n_L-=2*n_drop_high_sigma_stars;
   else
    if( n_L>n_drop_high_sigma_stars+2 )
     n_L-=n_drop_high_sigma_stars;
   ///
   ///
   mean=gsl_stats_median_from_sorted_data(data_L,1,n_L);
   y_limit_L[i]=esimate_sigma_from_MAD_of_sorted_data(data_L,n_L);
   if( 0!=isnan(y_limit_L[i]) ){
    y_limit_L[i]=0;
   }
   y_limit_L[i]=mean+M_SIGMA_BIN_MAG_SIGMA_DETECT*y_limit_L[i];
   gsl_sort2(data_sigma, 1, m_arr, 1, n);
   n_drop_high_sigma_stars=MAX((int)(0.07*n),M_SIGMA_BIN_DROP);
   // experimental
   if( n>2*n_drop_high_sigma_stars+2 )
    n-=2*n_drop_high_sigma_stars;
   else
    if( n>n_drop_high_sigma_stars+2 )
     n-=n_drop_high_sigma_stars;
   // end of experimental
   mean=gsl_stats_median_from_sorted_data(data_sigma,1,n); //gsl_stats_mean(data,1,n ); // would a median be more appropriate here?
   x_sigma[i]=gsl_stats_mean(m_arr,1,n);
   y[i]=mean;
   y_limit_sigma[i]=gsl_stats_sd(data_sigma,1,n);
   if( 0!=isnan(y_limit_sigma[i]) ){
     y_limit_sigma[i]=0;
   }
   y_limit_sigma[i]=y[i]+M_SIGMA_BIN_MAG_SIGMA_DETECT*y_limit_sigma[i]; // now this is the detection limit

   i++;
   n=0;
   mmax=m;
  }
 }
 fclose(dmsf);
 
 if( i<5 ){
  fprintf(stderr,"ERROR: not enough bins (only %d) to construct the selection curve!\n",i);
  exit(1);
 }
 

 gsl_set_error_handler_off(); // The function call to gsl_set_error_handler_off stops the default error handler from aborting the program.

 // Interpolate
 gsl_interp_accel *acc_sigma = gsl_interp_accel_alloc ();
 gsl_spline *spline_sigma = gsl_spline_alloc (gsl_interp_akima, i);
 interpolation_status_gsl = gsl_spline_init (spline_sigma, x_sigma, y_limit_sigma, i);

 gsl_interp_accel *acc_modified_sigma = gsl_interp_accel_alloc ();
 gsl_spline *spline_modified_sigma = gsl_spline_alloc (gsl_interp_akima, i);
 interpolation_status_gsl = gsl_spline_init (spline_modified_sigma, x_sigma, y_limit_modified_sigma, i);

 gsl_interp_accel *acc_I = gsl_interp_accel_alloc ();
 gsl_spline *spline_I = gsl_spline_alloc (gsl_interp_akima, i);
 interpolation_status_gsl = gsl_spline_init (spline_I, x_sigma, y_limit_I, i);

 gsl_interp_accel *acc_L = gsl_interp_accel_alloc ();
 gsl_spline *spline_L = gsl_spline_alloc (gsl_interp_akima, i);
 interpolation_status_gsl = gsl_spline_init (spline_L, x_sigma, y_limit_L, i);


 if( interpolation_status_gsl!=0 ){
  fprintf(stderr,"Interpolation error!\n");
  exit(1);
 }

 // Compute interpolating function and write it to file
 sigma_selection_curve_log=fopen("vast_sigma_selection_curve.log","w");
 if( sigma_selection_curve_log==NULL ){
  fprintf(stderr,"ERROR: Can't open file vast_sigma_selection_curve.log for writing!\n");
  exit(1);
 }
 double xi,yi;
 double yi_modified_sigma,yi_I,yi_L;
 for (xi = x_sigma[0]; xi < x_sigma[i-1]; xi += 0.01){
  yi = gsl_spline_eval (spline_sigma, xi, acc_sigma);
  yi_modified_sigma = gsl_spline_eval (spline_modified_sigma, xi, acc_modified_sigma);
  yi_I = gsl_spline_eval (spline_I, xi, acc_I);
  yi_L = gsl_spline_eval (spline_L, xi, acc_L);
  fprintf(sigma_selection_curve_log,"%lf %lf %lf %lf %lf\n",xi, yi, yi_modified_sigma, yi_I, yi_L);
 } 
 fclose(sigma_selection_curve_log);



 // ---------------------------------------------------------------------------------------------

 // Write stars which pass the selection criteria to data.m_sigma file
 dmsf=fopen("vast_lightcurve_statistics.log","r");
 if( dmsf==NULL ){
  fprintf(stderr,"ERROR: Can't open file vast_lightcurve_statistics.log for reading\n");
  report_lightcurve_statistics_computation_problem();
  exit(1);
 }
 /*
 candidates_list=fopen("vast_autocandidates.log","w");
 if( candidates_list==NULL ){
  fprintf(stderr,"ERROR: Can't open file vast_autocandidates.log for writing!\n");
  exit(1);
 }
 */
 while(-1<fscanf(dmsf,"%lf %lf %lf %lf %s  %lf  %lf %lf  %lf %lf %lf %lf %d %lf %lf %lf %lf %lf %lf  %lf %lf %lf %lf  %lf  %lf  %lf %lf %lf  %lf",&m,&sigma,&X,&Y,str, &modified_sigma_series,&skewness,&kurtosis,&I,&J,&K,&L,&points_in_lightcurve,&MAD_scaled_to_sigma,&lag1_autocorrelation,&RoMS,&Chi2,&I_sign_only,&AGN_v, &double_tmp, &double_tmp, &double_tmp, &double_tmp, &double_tmp, &double_tmp, &double_tmp, &double_tmp, &double_tmp, &double_tmp)){
  //fprintf(stderr,"%lf %lf\n",m,sigma);
  // This should avoid "interpolation error" crash.
  if( m<x_sigma[0] )m=x_sigma[0];
  if( m>x_sigma[i-1] )m=x_sigma[i-1];
  // Compare the measured value with a limit
  if( sigma>=gsl_spline_eval(spline_sigma, m, acc_sigma) )
   fprintf(stdout,"%10.6lf %.6lf %9.3lf %9.3lf %s\n",m,sigma,X,Y,str);
  // New output to the files
  //if( sigma>gsl_spline_eval(spline_sigma, m, acc_sigma) && modified_sigma_series>gsl_spline_eval(spline_modified_sigma, m, acc_modified_sigma) && I>gsl_spline_eval(spline_I, m, acc_I) && L>gsl_spline_eval(spline_L, m, acc_L) )
  // fprintf(candidates_list,"%s\n",str);
 }
 //fclose(candidates_list);
 fclose(dmsf);

 free(y_limit_sigma);
 free(x_sigma);
 free(y);
 free(data_sigma);
 free(data_modified_sigma);
 free(data_I);
 free(data_L);
 free(w);

// free(x_modified_sigma);
// free(x_I);
// free(x_L);

 return 0;
}