#include <iostream.h>
#include <iomanip>
#include "libvis/miscplot.h"
#include "libvis/miscpic.h"
#include <iostream>
#include "newmatap.h"
#include "newmatio.h"
#include "newimage/newimageall.h"
#include "miscmaths/miscmaths.h"
#include "miscmaths/miscprob.h"
#include <string>
#include "utils/log.h"
#include "melgmix.h"
#include "meloptions.h"


using namespace std;
using namespace Utilities;
using namespace NEWIMAGE;
using namespace MISCPLOT;
using namespace MISCPIC;
using namespace Melodic;


Matrix calcFFT(const Matrix& Mat)
  {
    Matrix res;
    for(int ctr=1; ctr <= Mat.Ncols(); ctr++)
      {
	ColumnVector tmpCol;
	tmpCol=Mat.Column(ctr);
	ColumnVector FtmpCol_real;
	ColumnVector FtmpCol_imag;
	ColumnVector tmpPow;
	if(tmpCol.Nrows()%2 != 0){
	  Matrix empty(1,1); empty=0;
	  tmpCol &= empty;}
	RealFFT(tmpCol,FtmpCol_real,FtmpCol_imag);
	tmpPow = pow(FtmpCol_real,2)+pow(FtmpCol_imag,2);
	tmpPow = tmpPow.Rows(2,tmpPow.Nrows());
	//if(opts.logPower.value()) tmpPow = log(tmpPow);
	if(res.Storage()==0){res= tmpPow;}else{res|=tmpPow;}
      }
    return res;
  } //Matrix calc_FFT()

string float2str(float f, int width, int prec, int scientif)
	{
	    ostrstream os;
	    int redw = int(std::abs(std::log10(std::abs(f))))+1;
	    if(width>0)
	      os.width(width);
	    if(scientif>0)
	      os.setf(ios::scientific);
	    os.precision(redw+std::abs(prec));
	    os.setf(ios::internal, ios::adjustfield);
	    os << f << '\0';
	    return os.str();
	}

void usage(void)
{
  cout << "Usage: test infile ICfile [mixfile]" << endl;
  exit(1);
}

int main(int argc, char *argv[])
{

  if (argc<3)
    usage();

  Matrix ICs;
  Matrix mixMatrix;
  Matrix fmixMatrix;
  volumeinfo ICvolInfo;
  volume<float> Mask;
  volume<float> Mean;

  string RAWfname;
  RAWfname = string(argv[1]);
  string ICfname;
  ICfname = string(argv[2]);

  string MIXfname;
  if (argc>3)
    MIXfname = string(argv[3]);

  cerr << argc << "  " << RAWfname << " " <<ICfname << " " << MIXfname << endl;
 
  {
    volume4D<float> RawData;
    cout << " Reading orig. data " << RAWfname << " ... ";
    read_volume4D(RawData,RAWfname,ICvolInfo);
    Mean = meanvol(RawData);

    float howmuch = 0.5*(std::min(std::min(std::abs(Mean.xdim()),std::abs(Mean.ydim())),std::abs(Mean.zdim())));

    cout << " Smoothing by " << howmuch << endl;
    volume<float> tmpvol = smooth(Mean,howmuch);
 
    miscpic newpic;
    char instr[10000];

    sprintf(instr," ");
    strcat(instr,"-s 2");
    strcat(instr," -A 950 ");
    strcat(instr,string("./res/m1.png").c_str());      
    newpic.slicer(Mean, instr, &ICvolInfo); 

    char instr2[10000];

    sprintf(instr2," ");
    strcat(instr2,"-s 2");
    strcat(instr2," -A 950 ");
    strcat(instr2,string("./res/m2.png").c_str());      
    newpic.slicer(tmpvol, instr2, &ICvolInfo); 



    cout << " done! " << endl;
  }

   {
    volume4D<float> RawIC;
    cout << "Reading components " << ICfname << "  ... ";
    read_volume4D(RawIC,ICfname);
    cout << " done" << endl;

    cout << "Creating mask   ... ";
    Mask = binarise(RawIC[0],float(RawIC[0].min()),float(RawIC[0].max()));

    ICs = RawIC.matrix(Mask);
    if(ICs.Nrows()>1){
      Matrix DStDev=stdev(ICs);
      
      volume4D<float> tmpMask;
      tmpMask.setmatrix(DStDev,Mask);
      
      float tMmax;
      volume<float> tmpMask2;
      tmpMask2 = tmpMask[0];
      tMmax = tmpMask2.max();
      double st_mean = DStDev.Sum()/DStDev.Ncols();
      double st_std  = stdev(DStDev.t()).AsScalar();
      
      Mask = binarise(tmpMask2,(float) max((float) st_mean-3*st_std,
					   (float) 0.01*st_mean),tMmax);  
      ICs = RawIC.matrix(Mask);
    }
    else{
      Mask = binarise(RawIC[0],float(0.001),float(RawIC[0].max())) 
	+ binarise(RawIC[0],float(RawIC[0].min()),float(-0.001));
      ICs = RawIC.matrix(Mask);
    }

    //cerr << "ICs : " << ICs.Ncols() << ICs.Nrows() << endl;
    cout << " done" << endl;
  }

  if(MIXfname.length()>0){
    cout << "Reading mixing matrix " << MIXfname << " ... ";
    mixMatrix = read_ascii_matrix(MIXfname);
    if (mixMatrix.Storage()<=0) {
      cerr <<" Please specify the mixing matrix correctly" << endl;
      exit(2);
      cout << " done " << endl;
    }
  }

  cout << " ICs: " << ICs.Nrows() << " x " << ICs.Ncols() << endl;
  if(ICs.Nrows()>0){
    cout << " Plotting histogram for map 1" << endl;
  
    miscplot newplot;
    // newplot.add_label("legend1");
    //newplot.add_label("legend2");
    //newplot.add_ylabel(string("yll1"));
    
    //newplot.add_xlabel(string("xl1"));
    //    newplot.set_xysize(600,600);

    Matrix mu(1,3);
    Matrix pi(1,3);
    Matrix std(1,3);
    mu(1,1)=0;mu(1,2)=2.10;mu(1,3)=-1.99388;
    pi(1,1)=0.911806;pi(1,2)=0.066722;pi(1,3)=0.021472;
    std(1,1)=1;std(1,2)=0.831590;std(1,3)=0.639299;
    newplot.gmmfit(ICs.Row(1),mu,std,pi,string("./res/g1.png"),string("Title"));

    mu(1,1)=0;mu(1,2)=2.209024;mu(1,3)=-2.016637;
    pi(1,1)=0.873335;pi(1,2)=0.086185;pi(1,3)=0.040480;
    std(1,1)=1;std(1,2)=1.012970;std(1,3)=0.616374;
    newplot.gmmfit(ICs.Row(2),mu,std,pi,string("./res/g2.png"),string("Title"));

    mu(1,1)=0;mu(1,2)=2.694458;mu(1,3)=-3.418763;
    pi(1,1)=0.9779397;pi(1,2)=0.017768;pi(1,3)=0.004293;
    std(1,1)=1;std(1,2)=2.848577;std(1,3)=0.385187;
    newplot.gmmfit(ICs.Row(3),mu,std,pi,string("./res/g3.png"),string("Title"));

    // newplot.histogram(ICs.Row(1),string("./res/h1.png"),string("Title"));
  }

  cout << endl << endl;
  /*
  Matrix data;
 
  cerr<< "Reading mixing matrix ";
  data = read_ascii_matrix(string("mixmat"));
  if (data.Storage()<=0) {
    cerr <<" Please specify the mixing matrix correctly" << endl;
  }
  cerr << "done " << endl;
  cerr << "Matrix size : " << data.Nrows() << " x " << data.Ncols() << endl;


  //set up data

  if(data.Ncols()==0){

    cerr << " create new data " << M_PI <<endl;

    ColumnVector X(1000);
   
    for(int i=1; i<=X.Nrows(); i++){
      X(i) = (float((i-1))/(X.Nrows()-1))*M_PI;
    }

    RowVector Y(10);
    for(int i=1; i<=Y.Ncols(); i++)
      Y(i)=3*i;

    data = X*Y;

    for(int i=1; i<=data.Ncols(); i++)
      for(int j=1; j<=data.Nrows(); j++)
	data(j,i) = sin(data(j,i));

   cerr << "Matrix size : " << data.Nrows() << " x " << data.Ncols() << endl;

   Log mmat;
mmat.setDir(string("./res/"),string("mix.txt"));
 mmat << data <<endl;

  }

  Log IChtml;
  IChtml.setDir(string("./res/"),string("index.html"));

  IChtml << "<HTML> " << endl
	     << "<TITLE>MELODIC Component " 
	     << "</TITLE>" << endl
	     << "<BODY BACKGROUND=\"file:" << getenv("FSLDIR") 
	     << "/doc/images/fsl-bg.jpg\">" << endl 
	     << "<hr><CENTER><H1>MELODIC Component "
	     << "</H1>"<< endl;


  float tr = 5.0;

  Matrix data2;
  data2 = calcFFT(data);

  for (int i=1; i<=data.Ncols(); i++)
    {
      
      {//plot time course
	miscplot newplot;
	//newplot.add_label("legend1");
	//newplot.add_label("legend2");
	//        newplot.add_ylabel(string("yll1"));
        //newplot.add_ylabel(string("ylabel2guiygyiasyugfuyigfuigasgfuiguiyasgfuigasuigfagifgiuygasisgfuigasgfgasgfigiagfgaggfigaif"));
	//newplot.add_xlabel(string("xl1"));
	//newplot.add_xlabel(string("xlabel2"));
	//newplot.set_xysize(600,200);
	if(tr>0.0)
	  newplot.timeseries(data.Column(i).t(),
			     string("./res/t")+num2str(i)+".png",
			     string("Timecourse (in seconds); TR = ")+
			     float2str(tr,0,2,0)+" s", 
			     tr);
	else
	  newplot.timeseries(data.Column(i).t(),
			     string("./res/t")+num2str(i)+".png",
			     string("Timecourse (in TRs)"));
	//	write_ascii_matrix(string("./res/t")+num2str(i)+".txt"),
	//	   data.Column(i));
	Log tfile;
	tfile.setDir(string("./res/"),string("t")+num2str(i)+".txt");
	tfile << data.Column(i) << endl;

	IChtml << "<A HREF=\"" << string("t")
	  +num2str(i)+".txt" << "\"> ";
	IChtml << "<img BORDER=0 SRC=\"" 
	  +string("t")+num2str(i)+".png\"></A><p>" << endl;
      }//time series plot
      
      {//plot frequency 
	miscplot newplot;
	int fact = int(std::pow(10.0,int(std::log10(float(data.Nrows())))));

	if(tr>0.0)
	  newplot.timeseries(data2.Column(i).t(),
			     string("./res/f")+num2str(i)+".png",
			     string("FFT of timecourse (in Hz / ") +num2str(fact)+")",
			     fact/(tr*data2.Nrows()));
	else
	  newplot.timeseries(data2.Column(i).t(),
			     string("./res/f")+num2str(i)+".png",
			     string(string("FFT of timecourse (in cycles); ")
				  +"frequency(Hz)=cycles/("
				  +num2str(data2.Nrows())
				  +"* TR); period(s)=("
				  +num2str(data2.Nrows())
				  +"* TR)/cycles"));

	Log ffile;
	ffile.setDir(string("./res/"),string("f")+num2str(i)+".txt");
	ffile << data2.Column(i) << endl;

	IChtml << "<A HREF=\"" << string("f")
	  +num2str(i)+".txt" << "\"> ";
	IChtml << "<img BORDER=0 SRC=\"" 
	  +string("f")+num2str(i)+".png\"></A><p>" << endl;
      }//frequency plot

    }


  IChtml<< "<HR><FONT SIZE=1>This page produced automatically by "	
	<< "FMRIB Software Library</A>.</FONT>" << endl
	<< "</BODY></HTML>" << endl;
	cerr << endl;  */
  return 0;
}