/*  MELODIC - Multivariate exploratory linear optimized decomposition into 
              independent components
    
    meldata.h - data container class

    Christian F. Beckmann, FMRIB Image Analysis Group
    
    Copyright (C) 1999-2008 University of Oxford */

/*  CCOPYRIGHT */


#ifndef __MELODICDATA_h
#define __MELODICDATA_h

#include "newimage/newimageall.h"
#include "utils/log.h"
#include "meloptions.h"
#include "melhlprfns.h"

using namespace Utilities;
using namespace NEWIMAGE;

namespace Melodic{
  
  class MelodicData{
    public:

      //constructor
      MelodicData(MelodicOptions &popts, Log &plogger):  
				opts(popts),logger(plogger)
			{	
	  		after_mm = false;
	  		Resels = 0;
			}  
 
      void save();

      Matrix process_file(string fname, int numfiles = 1);

      inline void save4D(Matrix what, string fname){
	 			volume4D<float> tempVol;
	 			tempVol.setmatrix(what,Mask);
	 			save_volume4D(tempVol,logger.appendDir(fname),tempInfo);
	 			message("  " << logger.appendDir(fname) << endl);
      }
      
      inline void saveascii(Matrix what, string fname){
	 			write_ascii_matrix(logger.appendDir(fname),what);   
	 			message("  " << logger.appendDir(fname) << endl);   
      }
 
      inline void savebinary(Matrix what, string fname){
      	write_binary_matrix(what,logger.appendDir(fname));  
	 			message("  " << logger.appendDir(fname) << endl);    
      }

      int  remove_components();
      void setup();
      void status(const string &txt);

      inline Matrix& get_pcaE() {return pcaE;}
      inline void set_pcaE(Matrix& Arg) {pcaE = Arg;}

      inline RowVector& get_pcaD() {return pcaD;}
      inline void set_pcaD(RowVector& Arg) {pcaD = Arg;}

      inline Matrix& get_data() {return Data;}
      inline void set_data(Matrix& Arg) {Data = Arg;}

      inline Matrix& get_IC() {return IC;}
      inline void set_IC(Matrix& Arg) {IC = Arg;}
      inline void set_IC(int ctr, Matrix& Arg) {IC.Row(ctr) = Arg;}
      
      inline vector<Matrix>& get_Smodes() {return Smodes;}
      inline Matrix& get_Smodes(int what) {return Smodes.at(what);}
      inline void add_Smodes(Matrix& Arg) {Smodes.push_back(Arg);}      
      inline void save_Smodes(){
				Matrix tmp = Smodes.at(0); 
				for(unsigned int ctr = 1; ctr < Smodes.size(); ctr++)
	  			tmp |= Smodes.at(ctr);
				  saveascii(tmp,opts.outputfname.value() + "_Smodes");
      }

      inline vector<Matrix>& get_Tmodes() {return Tmodes;}
      inline Matrix& get_Tmodes(int what) {return Tmodes.at(what);}
      inline void add_Tmodes(Matrix& Arg) {Tmodes.push_back(Arg);}
      inline void save_Tmodes(){
				Matrix tmp = Tmodes.at(0); 
				for(unsigned int ctr = 1; ctr < Tmodes.size(); ctr++)
	  			tmp |= Tmodes.at(ctr);
				saveascii(tmp,opts.outputfname.value() + "_Tmodes");
      }

      void set_TSmode();

      inline Matrix& get_white() {return whiteMatrix;}
      inline void set_white(Matrix& Arg) {whiteMatrix = Arg;}
      
      inline Matrix& get_dewhite() {return dewhiteMatrix;}
      inline void set_dewhite(Matrix& Arg) {dewhiteMatrix = Arg;}
      
      inline Matrix& get_meanC() {return meanC;}
      inline Matrix& get_meanR() {return meanR;}

      inline Matrix& get_stdDevi() {return stdDevi;}
      inline void set_stdDevi(Matrix& Arg) {stdDevi = Arg;}
  
      inline Matrix& get_mix() {return mixMatrix;}

      inline void set_mix(Matrix& Arg) {
				mixMatrix = Arg;
				if (Tmodes.size() < 1)
	  			for (int ctr = 1; ctr <= Arg.Ncols(); ctr++){
	    			Matrix tmp = Arg.Column(ctr);
	    			add_Tmodes(tmp);
	  			}
      }

      Matrix expand_mix(); 
      Matrix expand_dimred(const Matrix& Mat); 
      Matrix reduce_dimred(const Matrix& Mat); 
      
      inline Matrix& get_fmix() {return mixFFT;}
      inline void set_fmix(Matrix& Arg) {mixFFT = Arg;}

      inline Matrix& get_unmix() {return unmixMatrix;}
      inline void set_unmix(Matrix& Arg) {unmixMatrix = Arg;}

      inline volume<float>& get_mask() {return Mask;}
      inline void set_mask(volume<float>& Arg) {Mask = Arg;}
  
      inline volume<float>& get_mean() {return Mean;}
      inline void set_mean(volume<float>& Arg) {Mean = Arg;}
   
      inline volume<float>& get_bg() {
				if(opts.bgimage.value()>"")
					return background;
				else
					return Mean;
			}
      inline void set_bg(volume<float>& Arg) {background = Arg;}
   
      inline Matrix& get_Data() {return Data;}
      inline void set_Data(Matrix& Arg) {Data = Arg;}
    
      inline Matrix& get_RXweight() {return RXweight;}
      inline void set_RXweight(Matrix& Arg) {RXweight = Arg;}
 
      inline Matrix& get_ICstats() {return ICstats;}
      inline void set_ICstats(Matrix& Arg) {ICstats = Arg;}
     
      inline Matrix& get_EVP() {return EVP;}
      inline void set_EVP(Matrix& Arg) {if(EVP.Storage()==0)
																					EVP = Arg;}
      
      inline Matrix& get_EV() {return EV;}
      inline void set_EV(Matrix& Arg) {if(EV.Storage()==0)
																				 EV = Arg;}

      inline Matrix& get_PPCA() {return PPCA;}
      inline void set_PPCA(Matrix& Arg) {if(PPCA.Storage()==0)
																					 PPCA = Arg;}

      inline Matrix& get_stdNoisei() {return stdNoisei;}
      inline void set_stdNoisei(Matrix& Arg) {stdNoisei = Arg;}

      inline int data_dim() {return Data.Nrows();}
      inline int data_samples() {return Data.Ncols();}
     
      inline float get_resels() {return Resels;}
      inline void set_resels(float& Arg) {Resels = Arg;}

      inline int get_numfiles() {return numfiles;}

      inline void set_after_mm(bool val) {after_mm = val;}

      inline void flipres(int num){
				IC.Row(num) = -IC.Row(num);
				mixMatrix.Column(num) = -mixMatrix.Column(num);
				mixFFT=calc_FFT(mixMatrix);
				unmixMatrix = pinv(mixMatrix);
        if(ICstats.Storage()>0&&ICstats.Ncols()>3){
	  			double tmp;
	  			tmp = ICstats(num,3);
	  			ICstats(num,3) = -1.0*ICstats(num,4);
	  			ICstats(num,4) = -1.0*tmp;
				}
      }
      
      void sort();

      volumeinfo tempInfo;
      vector<Matrix> DWM, WM;
			basicGLM glmT, glmS;
			Matrix Tdes, Tcon, TconF, Sdes, Scon, SconF;	
			RowVector explained_var;		

    private:
      MelodicOptions &opts;     
      Log &logger;       

      Matrix pcaE;
      RowVector pcaD;
      Matrix whiteMatrix;
      Matrix dewhiteMatrix;
      Matrix meanC;
      Matrix meanR;
      Matrix stdDev;
      Matrix stdDevi;
      Matrix RXweight;
      Matrix mixMatrix;
      Matrix unmixMatrix;
      Matrix mixFFT;
      Matrix IC;
      Matrix ICstats;
      vector<Matrix> Tmodes;
      vector<Matrix> Smodes;

      Matrix EVP;
      Matrix EV;
      Matrix stdNoisei;
      volume<float> Mask;
      volume<float> Mean;
      volume<float> background;

      Matrix Data;
      Matrix PPCA;
      Matrix jointCC;
      
      bool after_mm;

      float Resels;
      int numfiles;

      char Mean_fname[1000];

      void setup_misc();
      void create_mask(volume<float>& theMask);
      void create_RXweight();
      void est_smoothness();

      unsigned long standardise(volume<float>& mask, 
				volume4D<float>& R);
      float est_resels(volume4D<float> R, volume<float> mask);
  };
}

#endif