#include <pocs.h>
#include <adpw.h>
#include <marv.h>
#include <string.h>
#include <math.h>
#include <fstream.h>
#include <unistd.h>
#include <stdlib.h>
#include <tcl.h>
#include <tk.h>

//////////////////////////////////////////////////////////////////
#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
#define RAND_REMOVE 0.9
#define POCS_ID 10
#define ADPW_ID 20
#define MARV_ID 30

#define YES 1
#define NO -1
///////////////// The global Variables needed ////////////////////
Tcl_Interp *interp;

pocs P;
adpw A;
marv M;

status S;
result pocs_result,adpw_result,marv_result; 
// The results of each iter stored in these result class instances.

float* orig_samples = NULL;// Original regularly sampled signal
  int orig_length;    // - ditto - length

float* irr_samples = NULL; // Irregular samples;
  int* irr_xi = NULL;        // irregular sample locations
  int irr_length;      // Sample train length
  int filt_B;

int h,w; // The width and height of the canvas used for drawing.
int pocs_on=1, adpw_on =1, marv_on =1;
float yscale;
//////////////// The Funcs to be bound to the Tcl buttons ///////////////////////

int StartIter(ClientData clientdata, Tcl_Interp* interp, int argc, char* argv[]);
int InitClasses(ClientData clientdata, Tcl_Interp* interp, int argc, char* argv[]);
int DumpIter(ClientData clientdata, Tcl_Interp* interp, int argc, char* argv[]);
int CompareChangedCmd(ClientData clientdata, Tcl_Interp* interp, int argc, char* argv[]);
int TimeFreqCmd(ClientData clientdata, Tcl_Interp* interp, int argc, char* argv[]);
void draw(result &, Tcl_Interp*);
void irregularize();
void fft(float data[], unsigned   nn, int isign);
int get99bandwidth(float*,int);
void puterror(char*);
int Tcl_AppInit(Tcl_Interp *interp);
//////////////////////////////////////////////////////////////////////
void main(int argc, char* argv[])
{
  Tk_Main(argc,argv,&Tcl_AppInit);
}
//////////////////////////////////////////////////////////////////////
void initializeVars(Tcl_Interp *interp)
{
  // To be called in from the AppInit func
  
  Tcl_LinkVar(interp,"params(fs)",(char*)((void*)&S.Fs),TCL_LINK_DOUBLE);
  Tcl_LinkVar(interp,"params(b)",(char*)((void*)&S.B),TCL_LINK_DOUBLE);
  Tcl_LinkVar(interp,"params(tol)",(char*)((void*)&S.tol),TCL_LINK_DOUBLE);
  Tcl_LinkVar(interp,"params(maxiter)",(char*)((void*)&S.maxiter),TCL_LINK_INT); 
  Tcl_LinkVar(interp,"params(set)",(char*)((void*)&S.params_set),TCL_LINK_INT);
  Tcl_LinkVar(interp,"gtnyquist",(char*)((void*)&S.gtnyquist),TCL_LINK_INT);
  Tcl_LinkVar(interp,"showfreq",(char*)((void*)&S.freq_domain),TCL_LINK_INT);
  Tcl_LinkVar(interp,"ready",(char*)((void*)&S.ready),TCL_LINK_INT);
  Tcl_LinkVar(interp,"pocs_compare",(char*)((void*)&S.pocs_compare),TCL_LINK_INT);
  Tcl_LinkVar(interp,"adpw_compare",(char*)((void*)&S.adpw_compare),TCL_LINK_INT);
  Tcl_LinkVar(interp,"marv_compare",(char*)((void*)&S.marv_compare),TCL_LINK_INT);
  Tcl_LinkVar(interp,"h",(char*)&h,TCL_LINK_INT);
  Tcl_LinkVar(interp,"w",(char*)&w,TCL_LINK_INT);
  Tcl_LinkVar(interp,"pocs_on",(char*)&pocs_on,TCL_LINK_INT);
  Tcl_LinkVar(interp,"adpw_on",(char*)&adpw_on,TCL_LINK_INT);
  Tcl_LinkVar(interp,"marv_on",(char*)&marv_on,TCL_LINK_INT);
  
  printf("Test variable showfreq %s\n",Tcl_GetVar(interp,"showfreq",TCL_GLOBAL_ONLY));
}
//////////////////////////////////////////////////////////////////////
void refreshVars(Tcl_Interp* interp)
{
  char *vars;

  vars = Tcl_GetVar(interp,"fs",TCL_GLOBAL_ONLY);
  sscanf(vars,"%lf",&S.Fs);
  vars = Tcl_GetVar(interp,"b",TCL_GLOBAL_ONLY);
  sscanf(vars,"%lf",&S.B);
  vars = Tcl_GetVar(interp,"tol",TCL_GLOBAL_ONLY);
  sscanf(vars,"%lf",&S.tol);
  vars = Tcl_GetVar(interp,"maxiter",TCL_GLOBAL_ONLY);
  sscanf(vars,"%d",&S.maxiter);
  vars = Tcl_GetVar(interp,"p_set",TCL_GLOBAL_ONLY);
  sscanf(vars,"%d",&S.params_set);
  vars = Tcl_GetVar(interp,"gtnyquist",TCL_GLOBAL_ONLY);
  sscanf(vars,"%d",&S.gtnyquist);
  vars = Tcl_GetVar(interp,"showfreq",TCL_GLOBAL_ONLY);
  sscanf(vars,"%d",&S.freq_domain);
  vars = Tcl_GetVar(interp,"ready",TCL_GLOBAL_ONLY);
  sscanf(vars,"%d",&S.ready);
  vars = Tcl_GetVar(interp,"pocs_compare",TCL_GLOBAL_ONLY);
  sscanf(vars,"%d",&S.pocs_compare);
  vars = Tcl_GetVar(interp,"adpw_compare",TCL_GLOBAL_ONLY);
  sscanf(vars,"%d",&S.adpw_compare);
  vars = Tcl_GetVar(interp,"marv_compare",TCL_GLOBAL_ONLY);
  sscanf(vars,"%d",&S.marv_compare);
  vars = Tcl_GetVar(interp,"h",TCL_GLOBAL_ONLY);
  sscanf(vars,"%d",&h);
  vars = Tcl_GetVar(interp,"w",TCL_GLOBAL_ONLY);
  sscanf(vars,"%d",&w);
  vars = Tcl_GetVar(interp,"pocs_on",TCL_GLOBAL_ONLY);
  sscanf(vars,"%d",&pocs_on);
  vars = Tcl_GetVar(interp,"adpw_on",TCL_GLOBAL_ONLY);
  sscanf(vars,"%d",&adpw_on);
  vars = Tcl_GetVar(interp,"marv_on",TCL_GLOBAL_ONLY);
  sscanf(vars,"%d",&marv_on);
}
//////////////////////////////////////////////////////////////////////
void CreateAllCommands(Tcl_Interp *interp)
{
  Tcl_CreateCommand(interp,"StartIter",StartIter, (ClientData)NULL, 
		    (Tcl_CmdDeleteProc *)NULL); 
  Tcl_CreateCommand(interp,"DumpIter",DumpIter, (ClientData)NULL, 
		    (Tcl_CmdDeleteProc *)NULL); 
  Tcl_CreateCommand(interp,"InitClasses", InitClasses, (ClientData)NULL, 
		    (Tcl_CmdDeleteProc *)NULL);
  Tcl_CreateCommand(interp,"CompareChangedCmd", CompareChangedCmd, (ClientData)NULL, 
		    (Tcl_CmdDeleteProc *)NULL);
  Tcl_CreateCommand(interp,"TimeFreqCmd", TimeFreqCmd, (ClientData)NULL, 
		    (Tcl_CmdDeleteProc *)NULL);
}
//////////////////////////////////////////////////////////////////////
int Tcl_AppInit(Tcl_Interp *interp)
{
  if(Tcl_Init(interp) == TCL_ERROR)
    return TCL_ERROR;
  if(Tk_Init(interp) == TCL_ERROR)
    return TCL_ERROR;

  Tcl_SetVar(interp,"tcl_interactive",0,TCL_GLOBAL_ONLY);
  CreateAllCommands(interp); // Create all the commands bound to the buttons
//  initializeVars(interp);    // This will link the C variables to the Tcl ones
  Tcl_EvalFile(interp,"test.tcl");
  Tk_MainLoop();  
  return TCL_OK;

}
//////////////////////////////////////////////////////////////////////
int StartIter(ClientData clientdata, Tcl_Interp* interp, int argc, char* argv[])
{
  
  refreshVars(interp);
  if(S.ready == NO)
     {
       Tcl_Eval(interp," ErrorCmd \"Not ready yet, parameters/filename\"");
       return 0;
     }
  if(S.params_set == NO)
    {
      Tcl_Eval(interp,"ErrorCmd \" Parameters not set! \"");
      return 0;
    }

  if(S.all_done == YES)
    {
	Tcl_Eval(interp,"ErrorCmd \" All methods finished \"");
	return 0;
    }
 
  if(S.all_done == NO)
    {
      if((S.pocs_done == NO) && (pocs_on == YES))
	{
	  S.iter++;
	  pocs_result = P.next_iter();
	  pocs_result.id = POCS_ID;
	  S.pocs_done = P.Done;
	  draw(pocs_result,interp);
	  if(S.pocs_done == YES)
	    cerr<<"POCS done \n";
	}
      if((S.adpw_done == NO)&&(adpw_on == YES))
	{
	  S.iter++;
	  adpw_result = A.next_iter();
	 adpw_result.id = ADPW_ID;
	 S.adpw_done = A.Done;
	 if(S.adpw_done == YES)
	   cerr<<"ADPW done \n";
	 draw(adpw_result,interp);
       }
     if((S.marv_done == NO)&&(marv_on==YES))
	 {
	   S.iter++;
	   marv_result = M.next_iter();
	   marv_result.id = MARV_ID;
	   S.marv_done = M.Done;
	   if(S.marv_done == YES)
	     cerr<<"MARV done \n";
	   draw(marv_result,interp);
	 }
   }
 if( (S.marv_done == YES)  && (S.pocs_done == YES) && (S.adpw_done == YES)) S.all_done = 1;
 
  return TCL_OK;  
}
//////////////////////////////////////////////////////////////////////
int DumpIter(ClientData clientdata, Tcl_Interp* interp, int argc, char* argv[])
{
  char outfile[200];
  refreshVars(interp);
  if( S.iter <=0 || S.ready == -1)
    {
      Tcl_Eval(interp,"ErrorCmd \" No Iterations in memory \"");
      return TCL_ERROR;
    }

  sprintf(outfile,"%sDump%0.4d.dat",argv[1],P.iter);
  ofstream out(outfile);
  if( strcmp(argv[1],"pocs") == 0)
    {
      P.DebugOn();
      P.SetDebugLevel(PRINT_RECONSTRUCTED_SIGNAL);
      out<<P;
      P.DebugOff();
    }
  else if ( strcmp(argv[1],"adpw") == 0)
    {
      out<<A;
    }
  else if ( strcmp(argv[1],"marv") == 0)
    {
      out<<M;
    }
  else
    {
      cerr<<" Incomprehensible argument"<<argv[1]<<"  passed to DumpIter\n";
    }
  out.close();
  
  return TCL_OK;
}
//////////////////////////////////////////////////////////////////////
int InitClasses(ClientData clientdata, Tcl_Interp* interp, int argc, char* argv[])
{
  char message[200];
  char name[400];
  float max = -10e10,min = 10e10;

  pocs_result.clear();
  adpw_result.clear();
  marv_result.clear();
  S.clear();
   
  refreshVars(interp);
  if (S.params_set == NO)
    {
      Tcl_Eval(interp,"ErrorCmd \" Parameters not set, Use the parameters menu button\"");
      return TCL_ERROR;
    }
  strcpy(name,argv[1]);
  cerr<<" NUmber of arguments in Init Classes "<<argc<<'\n';
  for(int l=0; l<argc;l++)
	cout<<argv[l]<<'\n';
  ifstream in(name);
  if( in == NULL) 
    {
      cerr<<"File"<<name<<"  not opened\n";
      exit(1);
    }
  orig_length = 0;
  in>>orig_length;
  sprintf(message,"ErrorCmd \"%d samples in %s \"",orig_length,name);
  Tcl_Eval(interp,message);

  if( orig_length <=0)
    {
      Tcl_Eval(interp,"ErrorCmd \" Empty file \"");
      return TCL_ERROR;
    }
  if( orig_samples != NULL) delete orig_samples;
  orig_samples = new float[orig_length];

  for(  int i=0; i<orig_length; i++)
    in>>orig_samples[i];
  in.close();
  // Irregularly resample the dataset
  irregularize();
  // reload all the class containers
  
  cerr<< "TOlerance is "<<S.tol<<'\n';
  
  P.reset(orig_samples, orig_length,
	    irr_samples, irr_xi, irr_length,
	    S.tol, S.maxiter, S.B, S.Fs);
  A.reset(orig_samples, orig_length,
	  irr_samples, irr_xi, irr_length,
	  S.tol, S.maxiter, S.B, S.Fs);
  M.reset(orig_samples, orig_length,
	  irr_samples, irr_xi, irr_length,
	  S.tol, S.maxiter, S.B, S.Fs);
  char pass[2];
  sprintf(pass,"%d",1);
  Tcl_SetVar(interp,"ready",pass,TCL_GLOBAL_ONLY);
  // Now calculate the x and y scales used for drawing.

  for( i=0; i< orig_length; i++)
    {
      if (orig_samples[i] > max) max = orig_samples[i];
      if (orig_samples[i] < min) min = orig_samples[i];
    }
  max = fabs(max) > fabs(min)?fabs(max):fabs(min);
  yscale = (float)(0.9*h)/(2*max);
  cout<<"finished initclass\n";
  return TCL_OK;
}
//////////////////////////////////////////////////////////////////////
void irregularize()
{
  int* rand_cleared = NULL, nyquist_gap,kk=0;
  int removed = 0,index,li,ui,times=0;
  
  // This func irregularly samples the given regular signal
  // after first bandlimiting it to B creating maximal gaps as
  // decided by the gtnyquist variable.
  // now randomly remove samples from the bandpassed signal 
  nyquist_gap = S.Fs/(2*S.B) -1 ;
  if( S.gtnyquist ==1 ) nyquist_gap *= 1.2;
  
  rand_cleared =  new int[orig_length];
  for(  int i=0; i<orig_length;i++)
    rand_cleared[i] = -1;
  
  while ( removed < RAND_REMOVE*orig_length && times < orig_length)
    {
      times++;
      index = fabs((orig_length-1)*drand48());
       if(rand_cleared[index] == -1)
	{
	  ui = li =index;
	  do {  if(li ==0) break;
		li--;
	      } while (rand_cleared[li] == 1);
	  do {  if(ui == orig_length-1) break;
		ui++;} while (rand_cleared[ui] == 1);
	  
	  if( (ui -li) <= nyquist_gap ) 
	    {
	      rand_cleared[index] = 1;
	      removed++;
	      times = 0;
	    }
	}
    }
  
  irr_length = orig_length - removed;
  cerr<<" Irr length = "<<irr_length<<'\n';
  if( irr_samples != NULL)
    { delete irr_samples;
    irr_samples = new float[irr_length];
    }
  else irr_samples = new float[irr_length];
  if( irr_xi != NULL) 
    { 
      delete irr_xi;
      irr_xi = new int[irr_length];
    }
  else irr_xi = new int[irr_length];
  
  for(i=0; i<orig_length; i++)
     if(rand_cleared[i] == -1)
       {
	 irr_xi[kk] = i;
	 irr_samples[kk] = orig_samples[i];
	 kk++;
	 cerr<<kk<<'\n';
       }
  cout<<" irr_samples indexed upto  "<<kk-1<<'\n'; 
  delete rand_cleared;
  cerr<<"Removed "<<removed<<" Samples randomly\n";  
}
//////////////////////////////////////////////////////////////////////
void fft(float data[], unsigned   nn, int isign)
{
	unsigned   n,mmax,m,j,istep,i;
	double wtemp,wr,wpr,wpi,wi,theta;
	float tempr,tempi;

	n=nn << 1;
	j=1;
	for (i=1;i<n;i+=2) {
		if (j > i) {
			SWAP(data[j],data[i]);
			SWAP(data[j+1],data[i+1]);
		}
		m=n >> 1;
		while (m >= 2 && j > m) {
			j -= m;
			m >>= 1;
		}
		j += m;
	}
	mmax=2;
	while (n > mmax) {
		istep=mmax << 1;
		theta=isign*(6.28318530717959/mmax);
		wtemp=sin(0.5*theta);
		wpr = -2.0*wtemp*wtemp;
		wpi=sin(theta);
		wr=1.0;
		wi=0.0;
		for (m=1;m<mmax;m+=2) {
			for (i=m;i<=n;i+=istep) {
				j=i+mmax;
				tempr=wr*data[j]-wi*data[j+1];
				tempi=wr*data[j+1]+wi*data[j];
				data[j]=data[i]-tempr;
				data[j+1]=data[i+1]-tempi;
				data[i] += tempr;
				data[i+1] += tempi;
			}
			wr=(wtemp=wr)*wpr-wi*wpi+wr;
			wi=wi*wpr+wtemp*wpi+wi;
		}
		mmax=istep;
	}
}
//////////////////////////////////////////////////////////////////////
void draw(result& R, Tcl_Interp* interp)
{
  char outname_o[30],canname_o[30],outname_i[30],canname_i[30];
  char message[300];
  float oy;
  float max = -10e10;
  unsigned   len;
  float xscale;
  int compare;

  refreshVars(interp);
  xscale = (float)w/float(R.final_length);
  len = R.final_length;
  switch (R.id)
    {
    case POCS_ID :
      { 
	strcpy(outname_o,".pocs.main.out.orig");
	strcpy(outname_i,".pocs.main.out.iter");
	strcpy(canname_o,".pocs.main.can.orig");
	strcpy(canname_i,".pocs.main.can.iter");
	compare = S.pocs_compare;
	break;
      }
    case ADPW_ID :
      { 
	strcpy(outname_o,".adpw.main.out.orig");
	strcpy(outname_i,".adpw.main.out.iter");
	strcpy(canname_o,".adpw.main.can.orig");
	strcpy(canname_i,".adpw.main.can.iter");
	compare = S.adpw_compare;
	break;
      } 
    case MARV_ID :
      { 
	strcpy(outname_o,".marv.main.out.orig");
	strcpy(outname_i,".marv.main.out.iter");
	strcpy(canname_o,".marv.main.can.orig");
	strcpy(canname_i,".marv.main.can.iter");
	compare = S.marv_compare;
	break;
      } 
    default : cerr<<" ubknown tag\n";break;
    };

  // first draw the main signal scale to fit the window properly.
  // the zoom in and out buttons will zoom and scale the x coords only.
  // draw the reconstructed signal first and the error signal next
  oy = h/2;
  if(S.iter >0)
    {
      sprintf(message,"%s delete all",canname_o);
      Tcl_Eval(interp,message);
      sprintf(message,"%s delete all",canname_i);
      Tcl_Eval(interp,message);
    }
  
  
  sprintf(message,"%s create line %d %d %d %d ",canname_o,0,h/2,w,h/2);
  Tcl_Eval(interp,message);
  sprintf(message,"%s create line %d %d %d %d ",canname_i,0,h/2,w,h/2);
  Tcl_Eval(interp,message);
  
  //Now the actual drawing of the curves;
  for(  int i =0; i<R.final_length-1 ;i++)
    {
      sprintf(message,"%s create line %f %f %f %f -width 1.5",canname_i,xscale*i, 
	      oy - yscale*R.final[i], xscale*(i+1), oy -yscale*R.final[i+1]);
      Tcl_Eval(interp,message);
      sprintf(message,"%s create line %f %f %f %f -width 1.5",canname_o,xscale*i, 
	      oy -yscale*R.err[i], xscale*(i+1), oy -yscale*R.err[i+1]);
      Tcl_Eval(interp,message);
    }
  if(compare == YES)
    {
      for(i =0; i<orig_length-1 ;i++)
	{
	  sprintf(message,"%s create line %d %d %d %d -fill red",canname_i,int(xscale*i), 
		  int(oy - yscale*orig_samples[i]), int(xscale*(i+1)), 
		  int(oy - yscale*orig_samples[i+1]));
	  Tcl_Eval(interp,message);
	} 
      for( i=0; i<irr_length;i++)
	{
	  sprintf(message,"%s create text %d %d -text \"*\" -fill blue -anchor center",canname_i,
		  (int)(xscale*irr_xi[i]),(int)( oy - yscale*orig_samples[irr_xi[i]])); 
	  Tcl_Eval(interp,message);
	}
    }
  
  sprintf(message,"%s config -text \"Samples: %d\nT Bandwidth: %0.2lf\"",outname_o,orig_length,S.B);
  Tcl_Eval(interp,message);
  sprintf(message,"%s config -text \"Iter: %d\nT incr: %0.4e\nTot T : %0.4e\nRMS %%error: %0.4f\"",
	  outname_i,R.iter, R.iter_time/1e6,R.tot_time/1e6,R.mean_err);
  Tcl_Eval(interp,message);
  
  printf("error %f\n",R.mean_err);
  
}
//////////////////////////////////////////////////////////////////////
void puterror(char *error)
{
  char *message;
  int length;
  
  length = 150;
  message = new char[length+40];
  sprintf(message,"ErrorCmd \"%s\"",error);
  Tcl_Eval(interp,message);
  delete message;
}
//////////////////////////////////////////////////////////////////////
int CompareChangedCmd(ClientData clientdata, Tcl_Interp* interp, int argc, char* argv[])
  {
    if(pocs_result.iter >0 && (strcmp("pocs",argv[1]) ==0))
      draw(pocs_result,interp);
    if(adpw_result.iter >0 && (strcmp("adpw",argv[1]) ==0))
      draw(adpw_result,interp);
    if(marv_result.iter >0 && (strcmp("marv",argv[1]) ==0))
      draw(marv_result,interp);
    return TCL_OK;
  }
//////////////////////////////////////////////////////////////////////
int TimeFreqCmd(ClientData clientdata, Tcl_Interp* interp, int argc, char* argv[])
{
  if(pocs_result.iter>0)
    draw(pocs_result,interp);
  return TCL_OK;
}
//////////////////////////////////////////////////////////////////////
int get99bandwidth(float* sig, int len)
{
  float* data,tot=0,sofar=0;
  int i;
  
  data = new float [2*len];
  for( i =0; i<2*len; i++)
    data[i] = 0;
  for(i = 0; i<len; i++)
    data[2*i] = sig[1];
  fft(data,len,1);
  
  for(i = 1; i<len/2; i++)
    tot+=data[2*i]*data[2*i];
  for(i=1;i<len/2;i++)
    {
      sofar+=data[2*i]*data[2*i];
      if(sofar >= 0.99*tot)
	break;
    }
  delete data;
  return i;
}
//////////////////////////////////////////////////////////////////////