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

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

/*  CCOPYRIGHT */


#ifndef __MELODICDATA_h
#define __MELODICDATA_h

#include "newimage/newimageall.h"
#include "utils/log.h"
#include "meloptions.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;
	}  

      //destructor
      ~MelodicData()
	{
	  if(opts.glmcleanupmode)
	    GLMcleanup_postproc();
	}
 
      void save();

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

      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 DiagonalMatrix& get_pcaD() {return pcaD;}
      inline void set_pcaD(DiagonalMatrix& 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 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; unmixMatrix=pinv(mixMatrix);}
      inline void set_mix(int ctr, Matrix& Arg) 
	{mixMatrix.Column(ctr) = Arg; unmixMatrix=pinv(mixMatrix);}
   
      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 Matrix& get_Data() {return Data;}
      inline void set_Data(Matrix& Arg) {Data = Arg;}
    
      inline Matrix& get_DataVN() {return DataVN;}
      inline void set_DataVN(Matrix& Arg) {DataVN = 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) {EVP = Arg;}
      
      inline Matrix& get_EV() {return EV;}
      inline void set_EV(Matrix& Arg) {EV = Arg;}

      inline Matrix& get_PPCA() {return PPCA;}
      inline void set_PPCA(Matrix& Arg) {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 void set_after_mm(bool val) {after_mm = val;}

      Matrix smoothColumns(const Matrix &inp);

      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;
	}
      }

      
      inline void varnorm(){
	Data = SP(Data,ones(Data.Nrows(),1)*stdDevi);
      }

      void sort();

      volumeinfo tempInfo;
      
      Matrix calc_FFT(const Matrix& Mat);

    private:
      MelodicOptions &opts;     
      Log &logger;       
      //      static MelodicData* melodat;

      Matrix pcaE;
      DiagonalMatrix 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;
      Matrix EVP;
      Matrix EV;
      Matrix stdNoisei;
      volume<float> Mask;
      volume<float> Mean;
      Matrix Data;
      Matrix DataVN;
      Matrix PPCA;
      
      bool after_mm;

      float Resels;

      char Mean_fname[1000];

      int glm_pes;

      void create_mask(volume4D<float> &theData, volume<float> &theMask);
      //      void remove_comp(Matrix& Mat, const int numIC);
      void create_RXweight();

      void GLMcleanup_preproc();
      void GLMcleanup_postproc();
         
    };
}

#endif