/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/*                                                                        */
/*                 BURST EXPECTATION SEARCH (BES) Method                  */
/*                ---------------------------------------                 */
/*                                                                        */
/*                                                                        */
/*  This demonstration program provides routines to study binned time     */
/*  series data to search for statistically significant fast bursts or    */
/*  spikes. The method used has been called the Burst Expectation Search  */
/*  (BES) method. A detailed description of the method can be found in:   */
/*                                                                        */
/*  "An Efficient Algorithm for the Detection of Infrequent Rapid         */
/*   Bursts in Time Series Data", Giles, A.B., 1997, ApJ, 474, 464        */
/*                                                                        */
/*  That paper contains suficient information to enable this C code       */
/*  program to be understood but I also provide a README file to aid      */
/*  in using the program. This file also describes some of the usefull    */
/*  options that are in the menu but not formally part of the method      */
/*  not hence not covered in the paper.                                   */
/*                                                                        */
/*                                                                        */
/*                                                                        */
/*  Main File:        BESmethod.c (this file)                             */
/*                                                                        */
/*  Include Files:    Poisson.c                                           */
/*                    ArrayCalc.c                                         */
/*                                                                        */
/*  Help File:        README                                              */
/*                                                                        */
/*                                                                        */
/*                                                                        */
/*  Author:   Dr.A.B.Giles    RXTE PCA group, USRA, NASA / GSFC           */
/*                                                                        */
/*  Date:     June 21 1997                                                */
/*                                                                        */
/*  Version:  1.0                                                         */
/*                                                                        */
/*  Contact:  Mail:   Code 662, NASA / GSFC, Greenbelt, MD 20771, USA     */
/*            Tel:    301 286 6628                                        */
/*            FAX:    301 286 1684                                        */
/*            E-mail  barry@pcasun3.gsfc.nasa.gov                         */
/*                                                                        */
/*                                                                        */
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/

#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>
#include <math.h>


/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/*  Use Defined Variables                                                     */
/*  Change these to suit your application                                     */

                                  /*  Definitions for Analysis Table          */
# define AnalysisWidth    256     /*  No. of columns in Analysis table        */
                                  /*  must be a power of 2                    */
# define ResultHeight       7     /*  No. rows in Analysis & Result tables    */
 
                                  /*  Definitions for Poisson table           */
# define MaxXrayHt        101     /*  Maximum bin height in Poisson table     */
# define ResultWidth        9     /*  No. columns in Result & Poisson tables  */

# define binsize          500     /*  Basic bin duration in micro seconds     */

                                  /*  Probability limits for table columns    */
float Prob_Limit_Array[ResultWidth] = {0.04, 0.02, 0.01, 0.005, 0.0025, 0.001, 0.0005, 0.0002, 0.0001 };

# define Sim              'Y'     /*  Sim = Y for simulation, N for user data */

/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/


/*  Definitions for Numerical Recipies routines  */
# define M 714025
# define IA 1366
# define IC 150889
# define PI 3.141592654

/*  Test listing to screen of details on all significant bursts  */
# define Test_Detections     'N'

/*  Test listing to screen of details of all calls to bursts test routine  */
# define Test_Evaluate       'N'

/*  Only ht >= 2 is meaningful in table  */
# define StartOfTable         2


/*  Routines in this file after main  */
void Initialise();
void ShowMenu();
void ClearScreen();
char keypress();
char getch();
void ChangeRepeats();
void ChangeRow();
void ChangeExampleLimit();
void ChangeSimMean();
void GetSomeData();
int GetUserData();

/*  Routines in include file Poisson.c  */
void SetupPoissonTable();
void ShowPoissonTable();
void SavePoissonTable();
void Threshold();
void PoissonProb();
float GetRandom();
float GetNoise();
void TestPoisson();
void HeightFactorial();
float gammln();
float poidev();
float ran2();

/*  Routines in include file ArrayCalc.c  */
void ProcessAnalysis();
void EvaluateBurst();
void ListAnalysisSummary();
void ListAnalysis();
void ListResults();
void ListNormalisation();
void ListScale();
void ListAll();
void ClearResults();
void ClearAnalysis();
void FileWriteResults();
void SaveLightCurve();


double      Comp[MaxXrayHt][ResultWidth];
double     HtFac[MaxXrayHt];
int     Analysis[ResultHeight][AnalysisWidth];
int AnalysisMark[ResultHeight];
long AnalysisSum[ResultHeight];
int  AnalysisBin[ResultHeight];
int  AnalysisCen[ResultHeight];
int  AnalysisTry[ResultHeight];
int       Result[ResultHeight][ResultWidth];
float     Normal[ResultHeight][ResultWidth];
float      Scale[ResultHeight][ResultWidth];

int NoRepeats, NoExamples, BlockNo, SaveRowNo, ExampleLimit;
int AnalysisWidth_1, ResultHeight_1;
int HalfWidth, Complete;
int AnalysisLive, FullArrayNo, NoFrames, NoLoops;

float SuperTime,PoissonMean;
int Rlog,Elog;

FILE *zzcon;
FILE *eecon;

#include "Poisson.c"
#include "ArrayCalc.c"

void main()
{
	int i,j,k,Quit;
	char cmd;

/*  Initialise variables  */
	AnalysisLive  =     0;
	NoRepeats     =     4;
	NoExamples    =     0;
	Complete      =     0;
	SaveRowNo     =     1;
	Quit          =     0;
	ExampleLimit  =     5;
	SuperTime     =   0.0;
	PoissonMean   =   1.0;

	Initialise();
	Rlog = 0;
	Elog = 0;

	cmd = 'Z';
	while (cmd != 'q') {
		ShowMenu();
		cmd = keypress();
		switch (cmd) {
		  case 't' : ShowPoissonTable(); break;
		  case 's' : SavePoissonTable(); break;
		  case 'p' : TestPoisson(20,PoissonMean); break;
		  case 'm' : ChangeSimMean(); break;

		  case '1' : {
				ClearResults();
				NoLoops = 1;
				GetSomeData(NoLoops,1); break;
			     }
		  case '2' : {
				ClearResults();
				NoLoops = NoFrames;
				GetSomeData(NoLoops,1); break;
			     }
		  case '3' : {
				ClearScreen();
			      if (Rlog==1) zzcon = fopen("Result_Log.dat","w");
			      if (Elog==1) {
					eecon = fopen("Bursts_Log.dat","w");
					fprintf(eecon,"Column Limit >= %3d\n",
						ExampleLimit);
					   }
				NoLoops = NoFrames;
				for (i=0; i<NoRepeats; i++) {
					printf("Processing Block %6d\r\n",i); 
					GetSomeData(NoLoops,0); 
					if (Rlog == 1)
					  FileWriteResults(i);
					BlockNo = i;
					ClearScreen();
				}
				if (Rlog == 1) {
			fprintf(zzcon,"Quit at Bin No. =    %12.1f\n",
					SuperTime);
					fclose(zzcon);
				}
				if (Elog == 1) {
			fprintf(eecon,"\nQuit at Bin No. =    %12.1f\n",
					SuperTime);
					fclose(eecon);
				}
				break;
			     }
		  case 'n' : ChangeRepeats(); break;

		  case 'l' : SaveLightCurve(0); break;
		  case 'r' : ChangeRow();  break;

		  case '4' : ListAnalysisSummary(1,1); break;
		  case '5' : ListAnalysis(1,1); break;
		  case '6' : ListResults(1,1); break;
		  case '7' : ListNormalisation(1,1); break;
		  case '8' : ListScale(1,1); break;
		  case '9' : ListAll(); break;

		  case 'y' : ClearResults(); break;
		  case 'a' : if (Rlog == 1) Rlog = 0; else Rlog = 1; break;

		  case 'b' : if (Elog == 1) Elog = 0; else Elog = 1; break;
		  case 'c' : ChangeExampleLimit(); break;

		  case 'q' : break;
		  default: printf("**  invalid key  **\n\n\n");
		}
	}
	ClearScreen();
}



void Initialise()
/*+++++++++++++++++++++++++*/
/*   Initialise variables  */
/*+++++++++++++++++++++++++*/
{
	int i;

	AnalysisWidth_1 = AnalysisWidth - 1;
	ResultHeight_1 = ResultHeight - 1;
	HalfWidth = AnalysisWidth / 2;
	HeightFactorial();
	ClearAnalysis();
	ClearResults();

/*  Set up random number generator  */
	srand(11);

/*  Compute no of basic inputs to reach "loaded " condition  */
	FullArrayNo = 1;
	for (i=1; i<ResultHeight; i++) FullArrayNo = FullArrayNo * 2;
	NoFrames = FullArrayNo;

/*  Set up the poisson statistics table  */
	SetupPoissonTable();
}




void ShowMenu()
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/*  Show menu options                                    */
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
{
ClearScreen();
printf("    MENU OPTIONS\n\n\r");
printf("      List Poisson Table                             ");
printf("t\n");
printf("      Save Poisson Table                             ");
printf("s\n");
printf("      Poisson Statistics - Test Sample               ");
printf("p\n");
printf("      Change Poisson Simulation Mean (%8.3f)      ",PoissonMean);
printf("m\n");
printf("      ------------------------------------------------\n\r");
printf("      Load & Process (%6d) Data Samples           ",AnalysisWidth);
printf("1\n");
printf("      Load & Process ( 1  Full Cycle)                ");
printf("2\n");
printf("      Load & Process ('N' Full Cycles)               ");
printf("3\n");
printf("      Set No. of Full Cycles ('N' =%5d)            ",NoRepeats);
printf("n\n");
printf("      ------------------------------------------------\n\r");
printf("      Save Light Curve - Working Array Row 'R'       ");
printf("l\n");
printf("      Set Row Value ('R' =%3d)                       ",SaveRowNo);
printf("r\n");
printf("      ------------------------------------------------\n\r");
printf("      List Working Array Summary                     ");
printf("4\n");
printf("      List Working Array                             ");
printf("5\n");
printf("      List Result Array                              ");
printf("6\n");
printf("      List Normalisation Array                       ");
printf("7\n");
printf("      List Scale Array                               ");
printf("8\n");
printf("      List All                                       ");
printf("9\n");
printf("      ------------------------------------------------\n\r");
printf("      Clear Result Array                             ");
printf("y\n");
if (Rlog==1) printf("      Result Array Log                     ON        ");
if (Rlog==0) printf("      Result Array Log                    OFF        ");
printf("a\n");
printf("      ------------------------------------------------\n\r");
if (Elog==1) printf("      Significant Bursts Log               ON        ");
if (Elog==0) printf("      Significant Bursts Log              OFF        ");
printf("b\n");
printf("      Set column limit for Bursts Log  (%4d)        ",ExampleLimit);
printf("c\n");
printf("      ------------------------------------------------\n\r");
printf("      Quit                                           ");
printf("q\n\n");
printf("\n\r");
printf("      Input Choice  ?   ");
}



void ClearScreen()
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/*  Utility to reset to top of screen                    */
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
{
	system("clear");
}


char keypress()
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/*  get a character from keyboard for menu               */
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
{
	char ch;

	scanf("%c", &ch);
	getchar();
	return (ch);
}


char getch()
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/*  get a continuation character from keyboard           */
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
{
	char ch;
	int i;

	printf("\n\n\nPress <Return> To Continue  >  ");
	i = getchar();
	ch = 'y';
	return (ch);
}



void ChangeRepeats()
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/*  Change value of number of Full repeat Cycles         */
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
{
	char ch;
	int i;

	ClearScreen();
	printf("Select The Number Of Full Analysis Cycles\n\n");
	printf("Present Value of No. of Cycles = %6d\n\n",NoRepeats);
	printf("Input new value  ? ");
	scanf("%d",&i);
	printf("New Value of No. of Cycles     = %6d\n\n",i);
	printf("\n\n");
	NoRepeats = i;
	ch = getch();
}



void ChangeRow()
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/*  Change value of analysis row number for light curve save   */
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
{
	char ch;
	int i;

	ClearScreen();
	printf("Select the Analysis Row to use for saving a Light Curve\n\n");
	printf("Present Value of Row No. = %6d\n\n",SaveRowNo);
	printf("Min =   1,  Max = %3d\n\n",ResultHeight);
	printf("Input new value  ? ");
	scanf("%d",&i);
	if (i <= 0) i = 1;
	if (i > ResultHeight) i = ResultHeight;
	printf("New Value of Row No.        = %6d\n\n",i);
	printf("\n\n");
	SaveRowNo = i;
	ch = getch();
}



void ChangeExampleLimit()
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/*  Change value of analysis row number for light curve save   */
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
{
	char ch;
	int i;

	ClearScreen();
	printf("Select the Probability Column Limit for Significant ");
	printf("Bursts Log\n\n");
	printf("Present Value of Column Limit = %6d\n\n",ExampleLimit);
	printf("Min =   1,  Max = %3d\n\n",ResultWidth);
	printf("Input new value  ? ");
	scanf("%d",&i);
	if (i <= 0) i = 1;
	if (i > ResultWidth) i = ResultWidth;
	printf("New Value of Cokumn Limit     = %6d\n\n",i);
	printf("\n\n");
	ExampleLimit = i;
	ch = getch();
}



void ChangeSimMean()
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/*  Change value of the mean used in the Poisson simulation   */
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
{
	char ch;
	float f;

	ClearScreen();
	printf("Select the Mean Value used in the Poisson Simulation\n\n");
	printf("Present Value of Mean = %8.3f\n\n",PoissonMean);
	printf("Input new value  ? ");
	scanf("%f",&f);
	if (f <= 0.0) f = 0.0000001;
	printf("New Value of Simulation Mean     = %8.3f\n\n",f);
	printf("\n\n");
	PoissonMean = f;
	ch = getch();
}



void GetSomeData(NoLoops, Info)
int NoLoops, Info;
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/*  Load arrays with real data & do analysis             */
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
{
	int i,j,k;
	int binval,idum;
	float mean;

	mean = PoissonMean;
	for (j=0; j<NoLoops; j++) {
		for (i=0; i<AnalysisWidth; i++) {
			if (Sim =='Y') binval = (int)poidev(mean, &idum);
				else binval = GetUserData();
			if (binval >= MaxXrayHt) binval = MaxXrayHt - 1;
			if (binval <= 0) binval = 0;
			ProcessAnalysis(binval);
			SuperTime = SuperTime + 1.0;
		}
	}
	if (Info == 1) {
		ClearScreen();
		printf("Full set of Real Data Loaded & Processed \n\n");
		ListAll();
	}
}



int GetUserData()
/*++++++++++++++++++++++++++++++++++++++++++*/
/*  User supplied single bin integer count  */
/*++++++++++++++++++++++++++++++++++++++++++*/
{
	int i;

	i = 2;

	return i;
}