melgmix.h

/*  MELODIC - Multivariate exploratory linear optimized decomposition into
    independent components

    melgmix.h - class for Gaussian/Gamma Mixture Model

    Christian F. Beckmann, FMRIB Analysis Group

    Copyright (C) 1999-2013 University of Oxford */

/*  CCOPYRIGHT */

#ifndef __MELGMIX_h
#define __MELGMIX_h


#include "armawrap/newmat.h"
#include "miscmaths/miscmaths.h"
#include "newimage/newimageall.h"
#include "utils/log.h"
#include "utils/options.h"
#include "melodic.h"
#include "meloptions.h"


namespace Melodic{

  class MelGMix{
  public:

    MelGMix(MelodicOptions &popts, Utilities::Log &plogger):
      opts(popts),
      logger(plogger){}

    ~MelGMix() {
      //mainhtml << endl << "<hr></CENTER></BODY></HTML>" << endl;
    }

    void save();

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

    void gmmfit();
    void ggmfit();

    inline void fit(std::string mtype = std::string("GGM")){
      mmtype = mtype;
      if(mmtype==std::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 NEWMAT::Matrix threshold(std::string levels){
      return this->threshold(data, levels);
    }
    inline NEWMAT::Matrix threshold(NEWMAT::RowVector& levels){
      return this->threshold(data, levels);
    }
    NEWMAT::Matrix threshold(const NEWMAT::RowVector& dat, NEWMAT::Matrix& levels);
    NEWMAT::Matrix threshold(const NEWMAT::RowVector& dat, std::string levels);

    void status(const std::string &txt);

    inline NEWMAT::RowVector& get_means() {return means;}
    inline void set_means(NEWMAT::RowVector& Arg) {means = Arg;}

    inline NEWMAT::RowVector& get_vars() {return vars;}
    inline void set_vars(NEWMAT::RowVector& Arg) {vars = Arg;}

    inline NEWMAT::RowVector& get_pi() {return props;}
    inline void set_pi(NEWMAT::RowVector& Arg) {props = Arg;}

    inline NEWMAT::RowVector& get_data() {return data;}
    inline void set_data(NEWMAT::RowVector& Arg) {data = Arg;}

    inline NEWMAT::RowVector& get_prob() {return probmap;}

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

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

    inline bool isfitted(){return fitted;}

    inline int mixtures(){return nummix;}

    inline std::string get_type() { return mmtype;}
    inline void set_type(std::string Arg) { mmtype = Arg;}

    inline std::string get_prefix() { return prefix;}
    inline void  set_prefix(std::string Arg) { prefix = Arg;}

    inline NEWMAT::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;
      }
    }

    void create_rep();

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

    inline std::string get_infstr(int num){
      if((threshinfo.size()<(unsigned int)(num-1))||(num<1))
        return std::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){
      NEWIMAGE::volume4D<float> tempVol;
      tempVol.setmatrix(probmap,Mask);
      tempVol[0]= smooth(tempVol[0],howmuch);
      probmap = tempVol.matrix(Mask);
    }

    inline NEWMAT::Matrix get_params(){
      NEWMAT::Matrix tmp = MISCMATHS::zeros(3,means.Ncols());
      tmp.Row(1) = means;
      tmp.Row(2) = vars;
      tmp.Row(3) = props;
      return tmp;
    }

    double datamean;
    double datastdev;

  private:
    __attribute__((unused)) MelodicOptions &opts;
    __attribute__((unused)) Utilities::Log &logger; //global log file

    //Log mainhtml;

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

    NEWMAT::Matrix Params;
    NEWMAT::Matrix threshmaps;

    NEWMAT::RowVector means;
    NEWMAT::RowVector vars;
    NEWMAT::RowVector props;
    NEWMAT::RowVector data;
    NEWMAT::RowVector probmap;

    NEWIMAGE::volume<float> Mean;
    NEWIMAGE::volume<float> Mask;

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

    int nummix;
    int numdata;
    int cnumber;

    bool fitted;
    bool fixdim;

    std::string prefix;
    std::string mmtype;
    std::string dirname;

    std::vector<std::string> threshinfo;

  };
}

#endif