#ifndef __HMM_H__
#define __HMM_H__

#include <iostream>
#include <fstream>
#include <string>
#include "emcluster.h"
#include "matrix.h"
#include <vector>
#include <math.h>
#include "timer.h"

///hidden markov chain
/**
 * does multi sequence training
 * - compile flags:
 *   -# __DEBUG_HMM__ runs extra hmm debug checks
 *   -# __DEBUG_NAN__ runs extra nan/inf checks
 *
 *
 * extra checks involve checking for normalized alpha
 * \warning currently will fail somethimes because of zero
 *          probablity with transition matrix
 */
class hmm
{
 public:
  struct model
  {
    int id;       ///< model id currently not used
    double **A;   ///< hmm transition matrix

    ///initial state probablity
    std::vector <double> pi;

    ///gaussian emission models
    gaussianmodel *gem;

    int states;   ///< number of gaussian emmition models

    long double ll; ///< loglikelyhood
  };

  ///serialized version of hmm for saving
  struct hmmserialized
  {
    char *data;  ///< buffer to store seralized data
    int nstates; ///< number of states in the hmm
    int nparam;  ///< parameters in each gaussian model
    int size;    ///< size of data in bytes
  };

  hmm();
  ~hmm();
  ///clears the training data
  int cleardata();
  ///clears the model
  int clearmodel();
  ///adds a data sequence for training
  /**
   * \param data vector of data points representing the sequence
   * \param s number of dimentions for each point
   */
  int adddata( std::vector< double * > &data, int s );

  ///start training process
  /**
   * \param n number of states
   */
  int train( int n);
  model getmodel();
  int print();

  ///serializes the data
  /**
   * memory will be allocated in d.data to store the model
   * does not save the pi parameters
   */
  int serialize( hmmserialized &d );

  ///deserializes the data
  /**
   * \param d structure with serialized data
   */
  int deserialize( hmmserialized &d );

  ///prints graphviz dot format of model
  /**
   * \param fn the file to place the results in
   */
  int printdot( std::string fn);

  /// prepare the model for threaded probability calculation
  int prepthread();

  ///returns the score for a sequence
  /**
   * \return score for the match
   * \warning does not check for matching data point size
   */
  double match( std::vector< matrix > &data );

  /// threaded match with reduced memory allocs
  /**
   * @param data the time series data
   * @param smatrix size of at least <state vector size>*<time series length>
   * @param calc_buffer state vector size
   */
  double match_t( std::vector< matrix > &data,
                  long double *smatrix,
                  long double *calc_buffer );

  ///finds probablity for sequence
  /**
   * \param data array of 1d matrix
   * \param size number of points in array
   * \return score for match
   */
  double match( matrix *data, int size );

  /// returns the size of the hmm
  /**
   * \param mynstates number of states in the hmm
   * \param myvecsize the vector size describing each state
   */
  int size( int &mynstates, int &myvecsize );
 private:

  ///structure for each training data sequence
  struct dseq
  {
    matrix* data;     ///< data 
    int datasize;     ///< size of data
    double **alpha;   ///< alpha matrix
    double **beta;    ///< beta matrix
    double **gamma;   ///< The gamma function
    double ***xi;     ///< xi
    double *c;        ///< scaling factor
  };

  model hmmodel;      ///< the model
  emcluster em;       ///< emcluster object for initialization

  int cleartrainlog(); ///< clear the training log 

  ///allocate the alpha beta memory
  int allocatememory();

  ///calculate the alpha based on the model
  int calcalpha();
  ///calculates alpha matrix along with the scaling values
  /**
   * \param data data sequence
   * calculates the scaled alpha matrix along with scaling factors
   */
  int calcalpha( dseq &data );

  ///calculate the beta matrix
  int calcbeta();
  ///calculates beta for a data sequence
  /**
   * \param data data sequence
   * \warning scaling values must first be calculated from alpha matrix
   */
  int calcbeta( dseq &data );

  ///calculates gamma values for all data sequences
  int calcgamma();
  ///calculates gamma values
  /**
   * \param data data sequence to calculate
   * calculations based on xi matrix
   * \warning xi must be calculated before gamma 
   */
  int calcgamma( dseq &data );

  ///calculates xi matrix
  int calcxi();
  ///calculates xi values for a given sequence
  /**
   * \param s data sequence
   * \warning calculations based on alpha, beta values
   */ 
  int calcxi( dseq &s);

  ///initializes the hmm
  /**
   * initialization done by running gaussian em modeling
   * and using the results as the emmision models
   */
  int inithmm();

  ///expectation step
  int estep();

  ///mstep
  /**
   * \param updategem flag to update the gaussian emmision models
   */
  int mstep(int updategem);

  ///calculates the probablity for emmision from model k
  /**
   * \param m data point
   * \param k gaussian model id
   */
  double emissionprob(matrix &m, int k);

  /// threaded emission prob with buffer
  double emissionprob( matrix &m, int k, long double *calc_buffer);

  ///same as emissionprob but returns the log
  double logemissionprob( matrix &m, int k);

  int copymodel(model &src, model &des);

  std::vector< dseq > datapts; ///< vector of training data

  std::vector< model > trainsteps;

  int nmodels; ///< number of models for emmision
};

#endif
