estpa.cpp

计算两个离散时间序列变量间的信息熵和互信息

//=========================================================
//
//This is a prog in the MutualInfo 0.9 package written by
// Hanchuan Peng.
//
//Disclaimer: The author of program is Hanchuan Peng
//      at <[e-mail]> and <[e-mail]>.
//
//The CopyRight is reserved by the author.
//
//Last modification: April/19/2002
//
//========================================================
//
// estpa.cpp
// Calculate the histgram/probability of one vector/image
// By Hanchuan Peng
// Dec/2000
// April/8/2002

#include "miinclude.h"

//return the number of states
template <class T> void copyvecdata(T * srcdata, long len, int * desdata, int& nstate, int &minn, int& maxx);

template <class T> void copyvecdata(T * srcdata, long len, int * desdata, int& nstate, int &minn, int& maxx)
{
  if(!srcdata || !desdata)
  {
    printf("NULL points in copyvecdata()!\n");
    return;
  } 

  long i;

  //note: originally I added 0.5 before rounding, however seems the negative numbers and 
  //      positive numbers are all rounded towarded 0; hence int(-1+0.5)=0 and int(1+0.5)=1;
  //      This is unwanted because I need the above to be -1 and 1.
  // for this reason I just round with 0.5 adjustment for positive and negative differently

  //copy data
  //int minn,maxx;
  if (srcdata[0]>0)
    maxx = minn = int(srcdata[0]+0.5);
  else
    maxx = minn = int(srcdata[0]-0.5);

  int tmp;
  double tmp1;
  for (i=0;i<len;i++)
  {
    tmp1 = double(srcdata[i]);
    tmp = (tmp1>0)?(int)(tmp1+0.5):(int)(tmp1-0.5);//round to integers
    minn = (minn<tmp)?minn:tmp;
    maxx = (maxx>tmp)?maxx:tmp;
    desdata[i] = tmp;
    //    printf("%i ",desdata[i]);
  }
  //printf("\n");

  //make the vector data begin from 0 (i.e. 1st state)
  for (i=0;i<len;i++)
  {
    desdata[i] -= minn;
  }

  //return the #state
  nstate = (maxx-minn+1);

  return;
}


void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{

  // check the arguments

  if(nrhs != 1 && nrhs !=2 && nrhs!=3)
    mexErrMsgTxt("Usage [marginprob,statelist,cumsumlist] = progname(vector1, maxstatenum, b_returnprob). \n(Both vectors can be images). Max range handled: INT type of the OS");
  if(nlhs > 3)
    mexErrMsgTxt("Too many output argument <marginalprob_list, statelist, cumsumlist>.");

  if (!mxIsInt8(prhs[0]) && !mxIsUint8(prhs[0]) && !mxIsDouble(prhs[0]) ) 
    mexErrMsgTxt("The first input argument must be types of INT8 or UINT8 or DOUBLE.");

  //get and check size information

  long i,j;

  void *img1 = (void *)mxGetData(prhs[0]);
  long len1 = mxGetNumberOfElements(prhs[0]);
  mxClassID type1 = mxGetClassID(prhs[0]);

  if (!img1 || !len1)
    mexErrMsgTxt("The input vector is invalid.");

  int b_findstatenum = 1;
  int nstate1 = 0;
  if (nrhs>=2)
  {
    b_findstatenum = 0;
    long MaxGrayLevel = (long) mxGetScalar(prhs[1]);
    nstate1 = MaxGrayLevel;
    if (MaxGrayLevel<=1)
    {
      printf("The argument #state is invalid. This program will decide #state itself.\n");
      b_findstatenum = 1;
    }
  }

  int b_returnprob = 1;
  if (nrhs>=3)
  {
    b_returnprob = (mxGetScalar(prhs[2])!=0);
  }

  //copy data into new INT type array (hence quantization) and then reange them begin from 0 (i.e. state1)

  int * vec1 = new int[len1];
  int nrealstate1=0, minn,maxx;
  switch(type1)
  {
    case mxINT8_CLASS: copyvecdata((char *)img1,len1,vec1,nrealstate1,minn,maxx); break;
    case mxUINT8_CLASS: copyvecdata((unsigned char *)img1,len1,vec1,nrealstate1,minn,maxx); break;
    case mxDOUBLE_CLASS: copyvecdata((double *)img1,len1,vec1,nrealstate1,minn,maxx); break;
  }  

  //update the #state when necessary
  if (nstate1<nrealstate1)
  {
    nstate1 = nrealstate1;
    //printf("First vector #state = %i\n",nrealstate1);
  }

  //generate the marginal-distribution list

  plhs[0] = mxCreateDoubleMatrix(nstate1,1,mxREAL);
  double *ha = (double *) mxGetPr(plhs[0]);

  for (i=0; i<nstate1;i++)
  {
    ha[i] = 0;
  }
  
  for (i=0;i<len1;i++)
  {
    ha[vec1[i]] += 1;
  }
  
  //return the probabilities, otherwise return count numbers
  if(b_returnprob)
  {
    for (i=0; i<nstate1;i++)
    {
      ha[i] /= len1;
    }
  }

  //return more information
  if (nlhs>=2)  //the second return value is the state list
  {
    plhs[1] = mxCreateDoubleMatrix(nstate1,1,mxREAL);
    double * pstate = (double *)mxGetPr(plhs[1]);
    for (i=0;i<nstate1;i++)
      pstate[i] = minn + i;
  }

  if (nlhs>=3) //the third return value is the cumsum prob list
  {
    plhs[2] = mxCreateDoubleMatrix(nstate1,1,mxREAL);
    double * pcumsum = (double *)mxGetPr(plhs[2]);
    pcumsum[0] = ha[0];
    for (i=1;i<nstate1;i++)
      pcumsum[i] = pcumsum[i-1] + ha[i];
  }

  //free memory
  if (vec1) delete []vec1;


  return;
}