ghmmwrapper.i

General Hidden Markov Model Library 一个通用的隐马尔科夫模型的C代码库

    
  /* array of sstate structs */
  sstate *arraysstate(int size) { 
    return (sstate *) malloc(size*sizeof(sstate));
  }	
  
   sstate *get_sstate_ptr(sstate *states, int k) {
    return &(states[k]);
  }
		  	
  void call_smodel_free(smodel *smo ) {smodel_free(&smo);}

  void free_smodel_array(smodel **smo) { if (smo){ m_free(smo);} }
  		
  void smodel_print_stdout(smodel *smo) {
    smodel_print(stdout, smo);
  }
 
  /* extract pointer to sstate  */	
  sstate *get_sstate(smodel *smo, int k) {
    return &(smo->s[k]);
  }


  /* creation and assessor functions for an array of smodels */
  smodel **smodel_array(int size){
  	return (smodel**) malloc(size * sizeof(smodel*));
  } 

  smodel *get_smodel_ptr(smodel **smo, int index) { return smo[index]; }
 
  void set_smodel_ptr(smodel **smo_array ,smodel *smo, int index) { smo_array[index] = smo; }
  
  smodel **cast_smodel_ptr(smodel *smo){
     smodel** res = (smodel**) malloc(sizeof(smodel*));
     res[0] = smo;
     return res;
  }   
  
  // write a smodel to a file
  void call_smodel_print(char *filename, smodel *smo) {
  FILE *fp=fopen(filename, "a");
  if (fp == NULL) {
    fprintf(stderr, "call_smodel_print(0): cannot open file %s\n", filename);    
  } else {
    smodel_print(fp, smo);
    fclose(fp);
  }
}

%}



/* =============================================================================================
   ============================== scluster.c  ================================================== */

/**
   Cluster structure: All models and sequences. */
struct scluster_t{
  /** 
  Vector of SHMMs pointers */
  smodel **smo;
  /** 
    Vector of sequence_d_t pointers; to store the sequences, that  belong to the models */
  sequence_d_t **smo_seq;
  /** 
    Number of models to read in */
  int smo_number;
  /** 
    Number of sequences for each model */
  long *seq_counter;
  /** 
    log(P) for the model */
  double *smo_Z_MD;
  /** a posteriori probability for the models to calculate the objective
      fucntion in case of a MAW classificator. Is calculated using smap_bayes */
  double *smo_Z_MAW;
};
/**
 */
typedef struct scluster_t scluster_t;


/**
   Frees the memory allocated for a scluster_t struct.
   @return 0 for success; -1 for error
   @param scl pointer to scl struct
*/
extern int scluster_t_free(scluster_t *scl);


/**
   Creates random labels for a vector of sequences
   @return 0 for success; -1 for error
   @param sqd vector of sequences
   @param smo_number number of models (needed to determine the interval
   for the random numbers)
*/
extern int scluster_random_labels(sequence_d_t *sqd, int smo_number);

/**
   Calculates the logarithmic probability of sequences for a model.
   @param cs sequences and model 
 */
extern void scluster_prob(smosqd_t *cs);

/** 
    Determines form an already calculated probability matrix, which model 
    fits best to a certain sequence. 
    @return index of best model if success, otherwize -1
    @param cl cluster 
    @param seq_id ID of the sequence in question
    @param all_log_p matrix containing the probability of each sequence
    for each model
    @param log_p the probability of the sequence in question for the 
    best fitting model
*/
extern int scluster_best_model(scluster_t *cl, long seq_id, double **all_log_p,
			double *log_p);

/**
   Updates the cluster with additional sequences.
   @return 0 for success; -1 for error
   @param cl cluster to update
   @param sqd sequences to update the cluster with
 */
extern int scluster_update(scluster_t *cl, sequence_d_t *sqd);

/**
   Prints the input vector for scluster_hmm
 */
void scluster_print_likelihood(FILE *outfile, scluster_t *cl);
	

/**
   Writes out the final models.
   @return 0 for success; -1 for error
   @param cl cluster of models to write
   @param sqd
   @param outfile output file
   @param out_filename name of the output file
 */
int scluster_out(scluster_t *cl, sequence_d_t *sqd, FILE *outfile, char *argv[]);

/** Calculates the aposteriori prob. $\log(p(\lambda_best | O[seq\_id]))$, 
    where $\lambda_best$ is the model with highest apost. prob.
    @return 0 for success; -1 for error
    @param cl cluster
    @param sqd the sequence in question
    @param seq_id the ID of the sequence
    @param lob_apo the results
*/
extern int  scluster_log_aposteriori(scluster_t *cl, sequence_d_t *sqd, int seq_id, double *log_apo);

/**
   Avoids empty models going out as outputs by assigning them a random 
   sequence. This may lead to a produce of a new empty model - therefore
   change out sequences until a non empty model is found. (Quit after 100 
   iterations to avoid an infinite loop). 
   @return 0 for success; -1 for error
   @param sqd sequences to generate the model from
   @param cl cluster for the models
 */
extern int scluster_avoid_empty_smodel(sequence_d_t *sqd, scluster_t *cl);

/**
   Makes a cluster and reestimates the HMMs.
   @return 0 for success; -1 for error
   @param argv vector of input files, one with sequences, one with models, 
   one for output and one with labels for the sequences - in this order.
 */
extern int scluster_hmm(char *argv[]);


%inline %{
  void scluster_printl_stdout(scluster_t *scl) {
    scluster_print_likelihood(stdout, scl);
  }
  
%}


/*=============================================================================================
  =============================== sreestimate.c  ============================================== */

/** @name struct smosqd_t
    structure that combines a continuous model (smo) and an integer
    sequence struct. Is used by sreestimate\_baum\_welch for 
    parameter reestimation.
 */
struct smosqd_t {
  /** pointer of continuous model*/
  smodel *smo;
  /** sequence\_d\__t pointer */
  sequence_d_t *sqd;
  /** calculated log likelihood */
  double* logp;
  /** leave reestimation loop if diff. between successive logp values 
      is smaller than eps */
  double eps;
  /** max. no of iterations */
  int max_iter;
}; 
typedef struct smosqd_t smosqd_t;


/**
  Baum-Welch Algorithm for SHMMs.
  Training of model parameter with multiple double sequences (incl. scaling).
  New parameters set directly in hmm (no storage of previous values!).
  @return            0/-1 success/error
  @param cs         initial model and train sequences
  */
extern int sreestimate_baum_welch(smosqd_t *cs);


%inline%{
  
  /************ create and manipulate an array of pointers to smosqd_t structs **************/
  smosqd_t *smosqd_t_array(int size) {
     return (smosqd_t *) malloc(size*sizeof(smosqd_t));
  }

  void set_smosq_t_smo(smosqd_t *cs, smodel *smo, int index){
	  cs[index].smo = smo;
  }	  
  
  smosqd_t get_smosqd_t_ptr(smosqd_t *cs, int i){  return cs[i];}
  
  void free_smosqd_t(smosqd_t *s){
	  if(s){
		m_free(s);
	 }	
  }	  
  		
%}	

/*=============================================================================================
  =============================== miscellanous functions ====================================== */


%inline %{
   
  // Some array helpers	
  /* Create any sort of int[size] array */
  int *int_array(int size) {
     return (int *) malloc(size*sizeof(int));
  }

  void set_arrayint(int  *ary, int index, int value) {
    ary[index] = value;
  }

  int  get_arrayint(int  *ary, int index) { return ary[index]; }

  void free_arrayi(int *pt ) { 
          m_free(pt); 
          (pt) = NULL;
  }

  /************ Create and access double[size] arrays ************/
 
   double *double_array(int size) {
     return (double *) malloc(size*sizeof(double));
  }

  void set_arrayd(double *ary, int index, double value) {
    ary[index] = value;
  }

  double get_arrayd(double *ary, int index) { return ary[index]; }
  
  void free_arrayd(double *pt) { m_free(pt); (pt) = NULL;}
  
  /************  Create and access sort of long[size] arrays **************/
    
  long *long_array(int size) {
     return (long *) malloc(size*sizeof(long));
  }

  void set_arrayl(long *ary, int index, long value) {
    ary[index] = value;
  }
  
  double get_arrayl(long *ary, int index) { return ary[index]; }
  
  void free_arrayl(long *ary) {m_free(ary);(ary) = NULL;}
  
  /*********** Create and access char** arrays ***********/
  char **char_array(int size) {
     return (char **) malloc(size*sizeof(char*));
  }

  void set_arraychar(char** ary, int index, char* value) {
    ary[index] = value;
  }

  char *get_arraychar(char** ary, int index) { return ary[index]; }

   
  /************  Create and access double[size1][size2] arrays ************/
  
 
  double **double_2d_array(int rows, int cols) {
    return matrix_d_alloc(rows,cols);
  }
  
  double **double_2d_array_nocols(int rows){
	return (double **) malloc(rows*sizeof(double*));
  }	  
  
  void set_2d_arrayd_col(double **ary, int index, double *col){
	  ary[index] = col;
  }	  
  
  void set_2d_arrayd(double **ary, int index1,int index2, double value) {
    ary[index1][index2] = value;
  }
  
  double get_2d_arrayd(double **ary, int index1, int index2) { return ary[index1][index2]; }


  double *get_col_pointer_d(double **ary, int index) {
	  return ary[index];
  }
  
  void double_2d_print(double **ary, int row, int col){
	int i,j;
	printf("(%dx%d) Matrix", row, col);
	for(i =0;i<row;i++){		  
      printf("\n");
      for(j =0;j<col;j++){		  		 
        printf("%f ",ary[i][j]);
	  }
	}
	printf("\n");  	
  }

  double** cast_ptr_d(double* array){
     double ** res = (double **) malloc(sizeof(double*));
	 res[0] = array;
	 return res;
  }	  
  
  void free_2darrayd(double **pt,int row) { matrix_d_free(&pt,row); }
   
  /************  Create and access int[size1][size2] arrays ************/
  
   int **int_2d_array_nocols(int rows){
	return (int **) malloc(rows*sizeof(int*));
  }	  
  
  void set_2d_arrayint_col(int **ary, int index, int *col){
	  ary[index] = col;
  }	  
  
  int *get_col_pointer_int(int **ary, int index) {
	  return ary[index];
  }
  
  
  void set_2d_arrayint(int **ary, int index1,int index2, int value) {
    ary[index1][index2] = value;
  }
  
  /* Get two dimensional array entry */
  int  get_2d_arrayint (int **ary, int index1, int index2) {return ary[index1][index2]; }
   
   int** cast_ptr_int(int* array){

	 int ** res = (int **) malloc(sizeof(int*));
	 res[0] = array;
	 return res; 
  }	 
  
  void free_2darrayint(int **pt, int rows,int cols) {  matrix_i_free(&pt, rows); }


  /**************** generalized deallocation *******************/  
  void freearray(void *pt)  { m_free(pt); }
   
%}


/******************* typemap ********************************/

%typemap(memberin) sequence_d_t * {
  sequence_d_t *src_pt;
  sequence_d_t *dest_pt;

  dest_pt = (sequence_d_t *) malloc(sizeof(sequence_d_t));
  src_pt  = $input;

  dest_pt.seq        = src_pt.seq;
  dest_pt.seq_len    = src_pt.seq_len;
  dest_pt.seq_label  = src_pt.seq_label;
  dest_pt.seq_id     = src_pt.seq_id;
  dest_pt.seq_w      = src_pt.seq_w;
  dest_pt.seq_number = src_pt.seq_number;
  dest_pt.total_w    = src_pt.total_w;

  $1 = dest_pt;
}

%typemap(check) double **{
    if ($1 == 0) {
      PyErr_SetString(PyExc_TypeError,"NULL Pointer not allowed");
      return NULL;
    }
}


%typemap(check) sequence_d_t *{
    if ($1 == 0) {
      PyErr_SetString(PyExc_TypeError,"NULL Pointer not allowed");
      return NULL;
    }
}


%typemap(check) sequence_t *{
    if ($1 == 0) {
      PyErr_SetString(PyExc_TypeError,"NULL Pointer not allowed");
      return NULL;
    }
}