/*  MELODIC - Multivariate exploratory linear optimized decomposition into 
              independent components
    
    ggmix.h - class for Gaussian/Gamma Mixture Model

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

/*  CCOPYRIGHT */

#ifndef __GGMIX_h
#define __GGMIX_h

#include "newimage/newimageall.h"
/*#include "utils/log.h"
#include "melodic.h"
#include "utils/options.h"
#include "meloptions.h"
*/

//using namespace Utilities;
using namespace NEWIMAGE;

namespace GGMIX{
  
  class ggmix
    {
    public:
    
      ggmix(){}
      //MelodicOptions &popts, Log &plogger):
      //	opts(popts),
      //	logger(plogger)
      //	{
      //	}

      ~ggmix() { 
	  }

      void save();

      void setup(const RowVector& dat, const string dirname,
		 int here, volume<float> themask, 
		 volume<float> themean, int num_mix = 3, 
		 float eps = 0.0, bool fixdim = false);
      
      void gmmfit();
      void ggmfit();

      inline void fit(string mtype = string("GGM"))
	{
	  mmtype = mtype;
	  if(mmtype==string("GGM")) 
	    this->ggmfit(); 
	  else 
	    this->gmmfit();

	  //re-insert mean and stdev

	  data = data*datastdev + datamean;
	  //threshmaps = threshmaps*datastdev + datamean;
	  means = means*datastdev + datamean;
	  vars = vars*datastdev*datastdev;
	}

      inline Matrix threshold(string levels)
	{return this->threshold(data, levels);}
      inline Matrix threshold(RowVector& levels) 
	{return this->threshold(data, levels);}
      Matrix threshold(const RowVector& dat, Matrix& levels);
      Matrix threshold(const RowVector& dat, string levels);

      void status(const string &txt);
 
      inline RowVector& get_means() {return means;}
      inline void set_means(RowVector& Arg) {means = Arg;}
    
      inline RowVector& get_vars() {return vars;}
      inline void set_vars(RowVector& Arg) {vars = Arg;}
      
      inline RowVector& get_pi() {return props;}
      inline void set_pi(RowVector& Arg) {props = Arg;}
      
      inline RowVector& get_data() {return data;}
      inline void set_data(RowVector& Arg) {data = Arg;}

      inline RowVector& get_prob() {return probmap;}

      inline float get_eps() {return epsilon;}
      inline void set_eps(float Arg) {epsilon = Arg;}

      inline Matrix& get_threshmaps() {return threshmaps;}
      inline void set_threshmaps(Matrix& Arg) {threshmaps = Arg;}

      inline bool isfitted(){return fitted;}

      inline int mixtures(){return nummix;}
     
      inline string get_type() { return mmtype;}
      inline void set_type(string Arg) { mmtype = Arg;}
      
      inline string get_prefix() { return prefix;}
      inline void  set_prefix(string Arg) { prefix = Arg;}

      inline RowVector get_probmap() {return probmap;}

      inline float get_offset() {return offset;}
      inline void set_offset(float Arg) {offset = Arg;}

      inline void flipres(int num){
	means = -means;
	data = -data;
	threshmaps = -threshmaps;
	if(mmtype=="GGM"){
	  float tmp;
	  tmp= means(2);means(2)=means(3);means(3)=tmp;
	  tmp=vars(2);vars(2)=vars(3);vars(3)=tmp;
	  tmp=props(2);props(2)=props(3);props(3)=tmp;
	}
      }      

      inline void add_infstr(string what){
	threshinfo.push_back(what);
      }

      inline string get_infstr(int num){
	if((threshinfo.size()<(unsigned int)(num-1))||(num<1))
	  return string("");
	else
	  return threshinfo[num-1];
      }

      inline int size_infstr(){
      return threshinfo.size();
      }

      inline void clear_infstr(){
	threshinfo.clear();
      }

      inline void smooth_probs(float howmuch){
	volume4D<float> tempVol;
	tempVol.setmatrix(probmap,Mask);
        tempVol[0]= smooth(tempVol[0],howmuch);
        probmap = tempVol.matrix(Mask);
      }



      double datamean;
      double datastdev;

    private:
      //    MelodicOptions &opts;     
      // Log &logger; //global log file

      //Log mainhtml;

      void gmmupdate();
      float gmmevidence();
      void gmmreducemm();
      void add_params(Matrix& mu, Matrix& sig, Matrix& pi, 
		      float logLH, float MDL, float Evi, bool advance = false);
      void get_params(int index, Matrix& mu, Matrix& sig, Matrix& pi, 
		      float logLH, float MDL, float Evi);

      Matrix Params;
      Matrix threshmaps;

      RowVector means;
      RowVector vars;
      RowVector props;
      RowVector data;
      RowVector probmap;

      volume<float> Mean;
      volume<float> Mask;

      float epsilon;
      float logprobY;
      float MDL;
      float Evi;
      float offset;

      int nummix;
      int numdata;
      int cnumber;

      bool fitted;
      bool fixdim;

      string prefix;
      string mmtype;
      string dirname;

      vector<string> threshinfo;

    };
}

#endif