#include "fileio.h"

#ifdef MATLABSUPPORT

#include <stdio.h>
#include <stdlib.h>
#include <matrix.h>
#include <mat.h>
#include <mclmcr.h>

#include "complexdata.h"

//////////////////////////////////////////////////////////////
typedef std::pair<Protocol,Data<float,4> > PDpair;

class LibMatSingleton
{
  public:
    bool initialized;
    LibMatSingleton() {
      Log<FileIO> odinlog("LibMatSingleton","LibMatSingleton");
      /* Initialize Matlab */
      const char *oplist[]={"-nojvm"};
      if (!mclInitializeApplication((const char **)oplist,1)) {
        ODINLOG(odinlog,errorLog) << "Could not initialize MATLAB interface" << STD_endl;
        initialized=false;
      }
      else {
        ODINLOG(odinlog,infoLog) << "MATLAB interface initialized " << STD_endl;
        initialized=true;
      }
    }
    ~LibMatSingleton() {
      Log<FileIO> odinlog("LibMatSingleton","~LibMatSingleton");
      ODINLOG(odinlog,infoLog) << "Shuting down MATLAB interface" << STD_endl;
      mclTerminateApplication();
    }
};

struct MatlabBinFormat : public FileFormat
{
  MATFile *fp;
  mxArray *root;
  FileWriteOpts w_opts;
  static bool initialized;
  MatlabBinFormat():fp(0){}
  STD_string description() const {return "binary data for Matlab";}
  svector suffix() const
  {
    svector result; result.resize(1);
    result[0]="mat";
    return result;
  }
  svector dialects() const {return svector();}

  bool init() {
    static LibMatSingleton single;
    return single.initialized;
  }
  template<class T> mxArray *dat2array(Data<float,4> dat,mxClassID classID) {
    Log<FileIO> odinlog("MatlabBinFormat","dat2array");
    mxArray *ret;
    const TinyVector<int, 4> shape(dat.shape());
    unsigned short dims;
    if(shape(0)>1)      dims=4;
    else if(shape(1)>1) dims=3;
    else if(shape(2)>1) dims=2;
    else                dims=1;

    mwSize size[dims];
    for(int i=0;i<dims;i++)
      size[i]=shape(i+4-dims);
    if(!(ret=mxCreateNumericArray(dims,size,classID,mxREAL))) {
      ODINLOG(odinlog,errorLog) << "could not create array of shape " << shape << STD_endl;
      return NULL;
    }
    T* ptr=(T*)mxGetData(ret);

    // matlab uses column-major ordering
    dat.transposeSelf(3,2,1,0);
    const unsigned int image_size=product(shape);

    //copy/convert data (, scale only down and only when converted to integers)
    Converter::convert_array(dat.c_array(),ptr,image_size,image_size);

    return ret;
  }

  mxArray *createMetaDataArray(const mwSize xdim,const mwSize ydim=1) {
    Log<FileIO> odinlog("MatlabBinFormat","createMetaDataArray");
    mxArray *ret;
    const mwSize dims[]={ydim,xdim};
    if(!(ret=mxCreateNumericArray(2,dims,mxDOUBLE_CLASS,mxREAL))) {
      ODINLOG(odinlog,errorLog) << "could not create " << xdim << "x" << ydim << " array for Metadata" << STD_endl;
      return NULL;
    }
    else return ret;
  }
  bool saveVar(const mxArray *arr,STD_string name) {
    Log<FileIO> odinlog("MatlabBinFormat","saveVar");
    if(matPutVariable(fp,name.c_str(),arr)!=0) {
      ODINLOG(odinlog,errorLog) << "could not add " << name << STD_endl;
      return false;
    }
    return true;
  }
  mxArray* vector2Array(const dvector vect) {
    const mwSize dim=vect.size();
    mxArray *arr=createMetaDataArray(dim);
    if(!arr)return NULL;

    double* ptr=(double*)mxGetData(arr);
    for(unsigned short i=0;i<dim;i++)
      ptr[i]=vect[i];
    return arr;
  }
  mxArray* triple2Array(const LDRtriple triple) {
    Log<FileIO> odinlog("MatlabBinFormat","triple2Array");
    mxArray *arr=createMetaDataArray(3);
    if(!arr) {
      ODINLOG(odinlog,errorLog) << "array creation failed" << STD_endl;
      return NULL;
    }

    double* ptr=(double*)mxGetData(arr);
    for(unsigned short i=0;i<3;i++)
      ptr[i]=triple[i];
    return arr;
  }

  mxArray* rotMatrix2Array(const RotMatrix rot) {
    Log<FileIO> odinlog("MatlabBinFormat","rotMatrix2Array");
    mxArray *arr=createMetaDataArray(3,3);
    if(!arr) {
      ODINLOG(odinlog,errorLog) << "array creation failed" << STD_endl;
      return NULL;
    }

    double* ptr=(double*)mxGetData(arr);
    for(unsigned short i=0;i<3*3;i++)
      ptr[i]=rot[i%3][i/3];
    return arr;
  }

  mxArray *dat2mat(const PDpair &pdpair, const FileWriteOpts& opts) {
    Log<FileIO> odinlog("MatlabBinFormat","dat2mat");
    const Protocol &prot=pdpair.first;
    const Data<float,4> &dat=pdpair.second;
    mxArray *matArray=NULL;

    STD_string type=select_write_datatype(prot,opts);
    if(type=="s8bit")       matArray=dat2array<s8bit> (dat,mxINT8_CLASS);
    else if(type=="u8bit")  matArray=dat2array<u8bit> (dat,mxUINT8_CLASS);
    else if(type=="s16bit") matArray=dat2array<s16bit>(dat,mxINT16_CLASS);
    else if(type=="u16bit") matArray=dat2array<u16bit>(dat,mxUINT16_CLASS);
    else if(type=="s32bit") matArray=dat2array<s32bit>(dat,mxINT32_CLASS);
    else if(type=="u32bit") matArray=dat2array<u32bit>(dat,mxUINT32_CLASS);
    else if(type=="float")  matArray=dat2array<float> (dat,mxSINGLE_CLASS);
    else if(type=="double") matArray=dat2array<double>(dat,mxDOUBLE_CLASS);
    else
      ODINLOG(odinlog,errorLog) << "datatype >" << type << "< unknown/unhandled" << STD_endl;

    return matArray;
  }

  template<class T> void array2dat(mxArray *matArray,Protocol &prot,Data<float,4> &dat) {
    prot.system.set_data_type(TypeTraits::type2label(T()));

    //copy/convert data (, scale only down and only when converted to integers)
    convert_from_ptr(dat,(T*)mxGetData(matArray),dat.shape());
  }

  const unsigned int mat2dat(Protocol &prot,Data<float,4> &dat) {
    Log<FileIO> odinlog("MatlabBinFormat","mat2dat");
    const char *name;
    mxArray *matArray= matGetNextVariable(fp, &name);
    if(!matArray)return 0;

    prot.seqpars.set_description(name);
    size_t dims=mxGetNumberOfDimensions(matArray);
    const mwSize *dim=mxGetDimensions(matArray);

    if(dims>2) {
      prot.seqpars.set_NumOfRepetitions(dim[2]);
      prot.geometry.set_nSlices(dim[2]);
    }
    if(dims>1)
      prot.seqpars.set_MatrixSize(phaseDirection,dim[1]);
    if(dims >0)
      prot.seqpars.set_MatrixSize(readDirection,dim[0]);

    switch(dims) {
      case 4:dat.resize(dim[3],  dim[2], dim[1], dim[0]);break;
      case 3:dat.resize(dim[2],  dim[1], dim[0], 1);break;
      case 2:dat.resize(dim[1],  dim[0], 1,      1);break;
      case 1:dat.resize(dim[0],  1,      1,      1);break;
      default:ODINLOG(odinlog,errorLog) << dims << " Dimensions not (yet) supportted" << STD_endl;
    }
    switch(mxGetClassID(matArray)) {
      case mxINT8_CLASS:  array2dat<mxInt8>   (matArray, prot, dat);break;
      case mxUINT8_CLASS: array2dat<mxUint8>  (matArray, prot, dat);break;
      case mxINT16_CLASS: array2dat<mxInt16>  (matArray, prot, dat);break;
      case mxUINT16_CLASS:array2dat<mxUint16> (matArray, prot, dat);break;
      case mxINT32_CLASS: array2dat<mxInt32>  (matArray, prot, dat);break;
      case mxUINT32_CLASS:array2dat<mxUint32> (matArray, prot, dat);break;
      /*			case mxINT64_CLASS:	array2dat<mxInt64>	(matArray, prot, dat);break;
            case mxUINT64_CLASS:array2dat<mxUint64>	(matArray, prot, dat);break;*/
      case mxSINGLE_CLASS:array2dat<mxSingle> (matArray, prot, dat);break;
      case mxDOUBLE_CLASS:array2dat<mxDouble> (matArray, prot, dat);break;
      default:
        ODINLOG(odinlog,errorLog) << "Datatype \""<< mxGetClassName(matArray) << "\" not (yet) supportted" << STD_endl;
        dat.resize(0,0,0,0);
        break;
    }
    // matlab uses column-major ordering
    dat.transposeSelf(3,2,1,0);
    prot.geometry.set_Mode(voxel_3d);
    prot.geometry.set_FOV(readDirection,dat.shape()(3));
    prot.geometry.set_FOV(phaseDirection,dat.shape()(2));
    prot.geometry.set_FOV(sliceDirection,dat.shape()(1));

    mxDestroyArray(matArray);
    return 1;
  }

  bool open(STD_string name,const char mode[]) {
    if(init() && (fp=matOpen(name.c_str(),mode)))return true;
    else return false;
  }
  void close() {
    matClose(fp);
  }
  int read(FileIO::ProtocolDataMap& pdmap, const STD_string& filename, const FileReadOpts& opts, const Protocol& protocol_template) {
    Log<FileIO> odinlog("MatlabBinFormat","read");
    if(!open(filename,"r")) {
      ODINLOG(odinlog,errorLog) << "could not open " << filename << " for reading: " << strerror(errno) << STD_endl;
      return -1;
    }
    Data<float,4> dat;
    Protocol prot(protocol_template);
    while(mat2dat(prot,dat)) {
      FileIO::ProtocolDataMap::const_iterator found=pdmap.find(prot);
      if(found!=pdmap.end())
        ODINLOG(odinlog,warningLog) << "equal Protocols found, skipping" << STD_endl;
      else pdmap[prot].reference(dat);
    }
    close();
    return 1;
  }

  void write(const PDpair &pair,unsigned short index, const FileWriteOpts& opts) {
    Log<FileIO> odinlog("MatlabBinFormat","write");
    STD_string seqDesc;
    int seqNum;
    pair.first.study.get_Series(seqDesc,seqNum);

    mxSetField(root, index, "seqDesc", mxCreateString(seqDesc.c_str()));
    mxSetField(root, index, "seqNum", mxCreateDoubleScalar(seqNum));
    mxSetField(root, index, "voxel", dat2mat(pair,opts));
    mxSetField(root, index, "center",vector2Array(pair.first.geometry.get_center()));
    mxSetField(root, index, "rotation",rotMatrix2Array(pair.first.geometry.get_gradrotmatrix()));
    LDRtriple *diff=dynamic_cast<LDRtriple *>(const_cast<Protocol&>(pair.first).methpars.get_parameter("Diffusion_bVector"));
    if(diff)mxSetField(root, index, "bVector",triple2Array(*diff));
  }
  int write(const FileIO::ProtocolDataMap& pdmap, const STD_string& filename, const FileWriteOpts& opts) {
    Log<FileIO> odinlog("MatlabBinFormat","write");
    const char *field_names[] = {"bVector","voxel", "center","rotation","seqDesc","seqNum"};
    if(!open(filename,"w")) {
      ODINLOG(odinlog,errorLog) << "could not open " << filename << " for writing: " << strerror(errno) << STD_endl;
      return -1;
    }
    unsigned short cnt=0;
    const mwSize dim = pdmap.size();
    root=mxCreateStructArray(1, &dim, sizeof(field_names)/sizeof(char*), field_names);

    w_opts=opts;
    for(FileIO::ProtocolDataMap::const_iterator pdit=pdmap.begin(); pdit!=pdmap.end(); ++pdit,cnt++)
      write(*pdit,cnt,opts);

    saveVar(root,"root");
    close();
    return pdmap.size();
  }
};

bool MatlabBinFormat::initialized=0;
#endif                                            //MATLABSUPPORT

//////////////////////////////////////////////////////////////

void register_mat_format() {
#ifdef MATLABSUPPORT
  static MatlabBinFormat mat;
  mat.register_format();
#endif                                          //MATLABSUPPORT
}