8 files changed, 791 insertions, 0 deletions
diff --git a/FEAST/MIToolbox/demonstration_algorithms/CMIM.m b/FEAST/MIToolbox/demonstration_algorithms/CMIM.m
new file mode 100644
index 0000000..8ae1f7c
--- /dev/null
+++ b/FEAST/MIToolbox/demonstration_algorithms/CMIM.m
@@ -0,0 +1,49 @@
+function selectedFeatures = CMIM(k, featureMatrix, classColumn)
+%function selectedFeatures = CMIM(k, featureMatrix, classColumn)
+%Computes conditional mutual information maximisation algorithm from
+%"Fast Binary Feature Selection with Conditional Mutual Information"
+%by F. Fleuret (2004)
+
+%Computes the top k features from
+%a dataset featureMatrix with n training examples and m features
+%with the classes held in classColumn.
+
+noOfTraining = size(classColumn,1);
+noOfFeatures = size(featureMatrix,2);
+
+partialScore = zeros(noOfFeatures,1);
+m = zeros(noOfFeatures,1);
+score = 0;
+answerFeatures = zeros(k,1);
+highestMI = 0;
+highestMICounter = 0;
+
+for n = 1 : noOfFeatures
+  partialScore(n) = mi(featureMatrix(:,n),classColumn);
+  if partialScore(n) > highestMI
+    highestMI = partialScore(n);
+    highestMICounter = n;
+  end
+end
+
+answerFeatures(1) = highestMICounter;
+
+for i = 2 : k
+  score = 0;
+  limitI = i - 1;
+  for n = 1 : noOfFeatures
+    while ((partialScore(n) >= score) && (m(n) < limitI))
+      m(n) = m(n) + 1;
+      conditionalInfo = cmi(featureMatrix(:,n),classColumn,featureMatrix(:,answerFeatures(m(n))));
+      if partialScore(n) > conditionalInfo
+        partialScore(n) = conditionalInfo;
+      end
+    end
+    if partialScore(n) >= score
+      score = partialScore(n);
+      answerFeatures(i) = n;
+    end
+  end
+end
+
+selectedFeatures = answerFeatures;
diff --git a/FEAST/MIToolbox/demonstration_algorithms/CMIM_Mex.c b/FEAST/MIToolbox/demonstration_algorithms/CMIM_Mex.c
new file mode 100644
index 0000000..daacb34
--- /dev/null
+++ b/FEAST/MIToolbox/demonstration_algorithms/CMIM_Mex.c
@@ -0,0 +1,158 @@
+/*******************************************************************************
+** Demonstration feature selection algorithm - MATLAB r2009a
+**
+** Initial Version - 13/06/2008
+** Updated - 07/07/2010
+** based on CMIM.m
+**
+** Conditional Mutual Information Maximisation
+** in
+** "Fast Binary Feature Selection using Conditional Mutual Information Maximisation
+** F. Fleuret (2004)
+**
+** Author - Adam Pocock
+** Demonstration code for MIToolbox
+*******************************************************************************/
+
+#include "mex.h"
+#include "MutualInformation.h"
+  
+void CMIMCalculation(int k, int noOfSamples, int noOfFeatures,double *featureMatrix, double *classColumn, double *outputFeatures)
+{
+  /*holds the class MI values
+  **the class MI doubles as the partial score from the CMIM paper
+  */
+  double *classMI = (double *)mxCalloc(noOfFeatures,sizeof(double));
+  /*in the CMIM paper, m = lastUsedFeature*/
+  int *lastUsedFeature = (int *)mxCalloc(noOfFeatures,sizeof(int));
+  
+  double score, conditionalInfo;
+  int iMinus, currentFeature;
+  
+  double maxMI = 0.0;
+  int maxMICounter = -1;
+  
+  int j,i;
+
+  double **feature2D = (double**) mxCalloc(noOfFeatures,sizeof(double*));
+
+  for(j = 0; j < noOfFeatures; j++)
+  {
+    feature2D[j] = featureMatrix + (int)j*noOfSamples;
+  }
+  
+  for (i = 0; i < noOfFeatures;i++)
+  {
+    classMI[i] = calculateMutualInformation(feature2D[i], classColumn, noOfSamples);
+    
+    if (classMI[i] > maxMI)
+    {
+      maxMI = classMI[i];
+      maxMICounter = i;
+    }/*if bigger than current maximum*/
+  }/*for noOfFeatures - filling classMI*/
+  
+  outputFeatures[0] = maxMICounter;
+  
+  /*****************************************************************************
+  ** We have populated the classMI array, and selected the highest
+  ** MI feature as the first output feature
+  ** Now we move into the CMIM algorithm
+  *****************************************************************************/
+  
+  for (i = 1; i < k; i++)
+  {
+    score = 0.0;
+    iMinus = i-1;
+    
+    for (j = 0; j < noOfFeatures; j++)
+    {
+      while ((classMI[j] > score) && (lastUsedFeature[j] < i))
+      {
+        /*double calculateConditionalMutualInformation(double *firstVector, double *targetVector, double *conditionVector, int vectorLength);*/
+        currentFeature = (int) outputFeatures[lastUsedFeature[j]];
+        conditionalInfo = calculateConditionalMutualInformation(feature2D[j],classColumn,feature2D[currentFeature],noOfSamples);
+        if (classMI[j] > conditionalInfo)
+        {
+          classMI[j] = conditionalInfo;
+        }/*reset classMI*/
+        /*moved due to C indexing from 0 rather than 1*/
+        lastUsedFeature[j] += 1;
+      }/*while partial score greater than score & not reached last feature*/
+      if (classMI[j] > score)
+      {
+        score = classMI[j];
+        outputFeatures[i] = j;
+      }/*if partial score still greater than score*/
+	  }/*for number of features*/
+  }/*for the number of features to select*/
+  
+  
+  for (i = 0; i < k; i++)
+  {
+    outputFeatures[i] += 1; /*C indexes from 0 not 1*/
+  }/*for number of selected features*/
+  
+}/*CMIMCalculation*/
+
+/*entry point for the mex call
+**nlhs - number of outputs
+**plhs - pointer to array of outputs
+**nrhs - number of inputs
+**prhs - pointer to array of inputs
+*/
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
+{
+  /*************************************************************
+  ** this function takes 3 arguments:
+  ** k = number of features to select,
+  ** featureMatrix[][] = matrix of features
+  ** classColumn[] = targets
+  ** the arguments should all be discrete integers.
+  ** and has one output:
+  ** selectedFeatures[] of size k
+  *************************************************************/
+  
+  int k, numberOfFeatures, numberOfSamples, numberOfTargets;
+  double *featureMatrix, *targets, *output;
+  
+  
+  if (nlhs != 1)
+  {
+    printf("Incorrect number of output arguments");
+  }/*if not 1 output*/
+  if (nrhs != 3)
+  {
+    printf("Incorrect number of input arguments");
+  }/*if not 3 inputs*/
+  
+  /*get the number of features to select, cast out as it is a double*/
+  k = (int) mxGetScalar(prhs[0]);
+
+  numberOfFeatures = mxGetN(prhs[1]);
+  numberOfSamples = mxGetM(prhs[1]);
+  
+  numberOfTargets = mxGetM(prhs[2]);
+  
+  if (numberOfTargets != numberOfSamples)
+  {
+    printf("Number of targets must match number of samples\n");
+    printf("Number of targets = %d, Number of Samples = %d, Number of Features = %d\n",numberOfTargets,numberOfSamples,numberOfFeatures);
+    
+    plhs[0] = mxCreateDoubleMatrix(0,0,mxREAL);
+  }/*if size mismatch*/
+  else
+  {
+    
+    featureMatrix = mxGetPr(prhs[1]);
+    targets = mxGetPr(prhs[2]);
+    
+    plhs[0] = mxCreateDoubleMatrix(k,1,mxREAL);
+    output = (double *)mxGetPr(plhs[0]);
+    
+    /*void CMIMCalculation(int k, int noOfSamples, int noOfFeatures,double *featureMatrix, double *classColumn, double *outputFeatures)*/
+    CMIMCalculation(k,numberOfSamples,numberOfFeatures,featureMatrix,targets,output);
+  }  
+  
+  return;
+}/*mexFunction()*/
diff --git a/FEAST/MIToolbox/demonstration_algorithms/DISR.m b/FEAST/MIToolbox/demonstration_algorithms/DISR.m
new file mode 100644
index 0000000..c8f5669
--- /dev/null
+++ b/FEAST/MIToolbox/demonstration_algorithms/DISR.m
@@ -0,0 +1,73 @@
+function selectedFeatures = DISR(k, featureMatrix, classColumn)
+%function selectedFeatures = DISR(k, featureMatrix, classColumn)
+%
+%Computers optimal features according to DISR algorithm from 
+%On the Use of variable "complementarity for feature selection" 
+%by P Meyer, G Bontempi (2006)
+%
+%Computes the top k features from
+%a dataset featureMatrix with n training examples and m features
+%with the classes held in classColumn.
+%
+%DISR - arg(Xi) max(sum(Xj mem XS)(SimRel(Xij,Y)))
+%where SimRel = MI(Xij,Y) / H(Xij,Y)
+
+totalFeatures = size(featureMatrix,2);
+classMI = zeros(totalFeatures,1);
+unselectedFeatures = ones(totalFeatures,1);
+score = 0;
+currentScore = 0;
+innerScore = 0;
+iMinus = 0;
+answerFeatures = zeros(k,1);
+highestMI = 0;
+highestMICounter = 0;
+currentHighestFeature = 0;
+
+%create a matrix to hold the SRs of a feature pair. 
+%initialised to -1 as you can't get a negative SR.
+featureSRMatrix = -(ones(k,totalFeatures));
+
+for n = 1 : totalFeatures
+	classMI(n) = mi(featureMatrix(:,n),classColumn);
+	if classMI(n) > highestMI
+		highestMI = classMI(n);
+		highestMICounter = n;
+	end
+end
+
+answerFeatures(1) = highestMICounter;
+unselectedFeatures(highestMICounter) = 0;
+
+for i = 2 : k
+	score = 0;
+	currentHighestFeature = 0;
+	iMinus = i-1;
+  for j = 1 : totalFeatures
+		if unselectedFeatures(j) == 1
+			%DISR - arg(Xi) max(sum(Xj mem XS)(SimRel(Xij,Y)))
+			%where SimRel = MI(Xij,Y) / H(Xij,Y)
+			currentScore = 0;
+			for m = 1 : iMinus
+        if featureSRMatrix(m,j) == -1
+          unionedFeatures = joint([featureMatrix(:,answerFeatures(m)),featureMatrix(:,j)]);
+				  tempUnionMI = mi(unionedFeatures,classColumn);
+				  tempTripEntropy = h([unionedFeatures,classColumn]);
+				  featureSRMatrix(m,j) = tempUnionMI/tempTripEntropy;
+        end
+        
+				currentScore =  currentScore + featureSRMatrix(m,j);
+			end
+			if (currentScore > score)
+				score = currentScore;
+				currentHighestFeature = j;
+			end
+		end
+	end
+	%now highest feature is selected in currentHighestFeature
+	%store it
+	unselectedFeatures(currentHighestFeature) = 0;
+	answerFeatures(i) = currentHighestFeature;
+end
+
+selectedFeatures = answerFeatures;
diff --git a/FEAST/MIToolbox/demonstration_algorithms/DISR_Mex.c b/FEAST/MIToolbox/demonstration_algorithms/DISR_Mex.c
new file mode 100644
index 0000000..617ca81
--- /dev/null
+++ b/FEAST/MIToolbox/demonstration_algorithms/DISR_Mex.c
@@ -0,0 +1,199 @@
+/*******************************************************************************
+** Demonstration feature selection algorithm - MATLAB r2009a
+**
+** Initial Version - 13/06/2008
+** Updated - 07/07/2010
+** based on DISR.m
+**
+** Double Input Symmetrical Relevance
+** in
+** "On the Use of Variable Complementarity for Feature Selection in Cancer Classification"
+** P. Meyer and G. Bontempi (2006)
+**
+** Author - Adam Pocock
+** Demonstration code for MIToolbox 
+*******************************************************************************/
+
+#include "mex.h"
+#include "MutualInformation.h"
+#include "Entropy.h"
+#include "ArrayOperations.h"
+
+void DISRCalculation(int k, int noOfSamples, int noOfFeatures,double *featureMatrix, double *classColumn, double *outputFeatures)
+{
+  /*holds the class MI values*/
+  double *classMI = (double *)mxCalloc(noOfFeatures,sizeof(double));
+  
+  char *selectedFeatures = (char *)mxCalloc(noOfFeatures,sizeof(char));
+  
+  /*holds the intra feature MI values*/
+  int sizeOfMatrix = k*noOfFeatures;
+  double *featureMIMatrix = (double *)mxCalloc(sizeOfMatrix,sizeof(double));
+  
+  double maxMI = 0.0;
+  int maxMICounter = -1;
+  
+  double **feature2D = (double**) mxCalloc(noOfFeatures,sizeof(double*));
+  
+  double score, currentScore, totalFeatureMI;
+  int currentHighestFeature;
+  
+  double *mergedVector = (double *) mxCalloc(noOfSamples,sizeof(double));
+  
+  int arrayPosition;
+  double mi, tripEntropy;
+  
+  int i,j,x;
+  
+  for(j = 0; j < noOfFeatures; j++)
+  {
+    feature2D[j] = featureMatrix + (int)j*noOfSamples;
+  }
+  
+  for (i = 0; i < sizeOfMatrix;i++)
+  {
+    featureMIMatrix[i] = -1;
+  }/*for featureMIMatrix - blank to -1*/
+
+
+  for (i = 0; i < noOfFeatures;i++)
+  {    
+    /*calculate mutual info
+    **double calculateMutualInformation(double *firstVector, double *secondVector, int vectorLength);
+    */
+    classMI[i] = calculateMutualInformation(feature2D[i], classColumn, noOfSamples);
+    
+    if (classMI[i] > maxMI)
+    {
+      maxMI = classMI[i];
+      maxMICounter = i;
+    }/*if bigger than current maximum*/
+  }/*for noOfFeatures - filling classMI*/
+  
+  selectedFeatures[maxMICounter] = 1;
+  outputFeatures[0] = maxMICounter;
+  
+  /*****************************************************************************
+  ** We have populated the classMI array, and selected the highest
+  ** MI feature as the first output feature
+  ** Now we move into the DISR algorithm
+  *****************************************************************************/
+  
+  for (i = 1; i < k; i++)
+  {
+    score = 0.0;
+    currentHighestFeature = 0;
+    currentScore = 0.0;
+    totalFeatureMI = 0.0;
+    
+    for (j = 0; j < noOfFeatures; j++)
+    {
+      /*if we haven't selected j*/
+      if (selectedFeatures[j] == 0)
+      {
+        currentScore = 0.0;
+        totalFeatureMI = 0.0;
+        
+        for (x = 0; x < i; x++)
+        {
+          arrayPosition = x*noOfFeatures + j;
+          if (featureMIMatrix[arrayPosition] == -1)
+          {
+            /*
+            **double calculateMutualInformation(double *firstVector, double *secondVector, int vectorLength);
+            **double calculateJointEntropy(double *firstVector, double *secondVector, int vectorLength);
+            */
+            
+            mergeArrays(feature2D[(int) outputFeatures[x]], feature2D[j],mergedVector,noOfSamples);
+            mi = calculateMutualInformation(mergedVector, classColumn, noOfSamples);
+            tripEntropy = calculateJointEntropy(mergedVector, classColumn, noOfSamples);
+            
+            featureMIMatrix[arrayPosition] = mi / tripEntropy;
+          }/*if not already known*/
+          currentScore += featureMIMatrix[arrayPosition];
+        }/*for the number of already selected features*/
+        
+			  if (currentScore > score)
+			  {
+				  score = currentScore;
+				  currentHighestFeature = j;
+			  }
+			}/*if j is unselected*/
+	  }/*for number of features*/
+  
+    selectedFeatures[currentHighestFeature] = 1;
+    outputFeatures[i] = currentHighestFeature;
+  
+  }/*for the number of features to select*/
+  
+  mxFree(mergedVector);
+  mergedVector = NULL;
+            
+  for (i = 0; i < k; i++)
+  {
+    outputFeatures[i] += 1; /*C indexes from 0 not 1*/
+  }/*for number of selected features*/
+  
+}/*DISRCalculation(double[][],double[])*/
+
+/*entry point for the mex call
+**nlhs - number of outputs
+**plhs - pointer to array of outputs
+**nrhs - number of inputs
+**prhs - pointer to array of inputs
+*/
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
+{
+  /*************************************************************
+  ** this function takes 3 arguments:
+  ** k = number of features to select,
+  ** featureMatrix[][] = matrix of features
+  ** classColumn[] = targets
+  ** the arguments should all be discrete integers.
+  ** and has one output:
+  ** selectedFeatures[] of size k
+  *************************************************************/
+  
+  int k, numberOfFeatures, numberOfSamples, numberOfTargets;
+  double *featureMatrix, *targets, *output;
+  
+  
+  if (nlhs != 1)
+  {
+    printf("Incorrect number of output arguments");
+  }/*if not 1 output*/
+  if (nrhs != 3)
+  {
+    printf("Incorrect number of input arguments");
+  }/*if not 3 inputs*/
+  
+  /*get the number of features to select, cast out as it is a double*/
+  k = (int) mxGetScalar(prhs[0]);
+
+  numberOfFeatures = mxGetN(prhs[1]);
+  numberOfSamples = mxGetM(prhs[1]);
+  
+  numberOfTargets = mxGetM(prhs[2]);
+  
+  if (numberOfTargets != numberOfSamples)
+  {
+    printf("Number of targets must match number of samples\n");
+    printf("Number of targets = %d, Number of Samples = %d, Number of Features = %d\n",numberOfTargets,numberOfSamples,numberOfFeatures);
+    
+    plhs[0] = mxCreateDoubleMatrix(0,0,mxREAL);
+  }/*if size mismatch*/
+  else
+  {
+    
+    featureMatrix = mxGetPr(prhs[1]);
+    targets = mxGetPr(prhs[2]);
+    
+    plhs[0] = mxCreateDoubleMatrix(k,1,mxREAL);
+    output = (double *)mxGetPr(plhs[0]);
+    
+    /*void DISRCalculation(int k, int noOfSamples, int noOfFeatures,double *featureMatrix, double *classColumn, double *outputFeatures)*/
+    DISRCalculation(k,numberOfSamples,numberOfFeatures,featureMatrix,targets,output);
+  }
+    
+  return;
+}/*mexFunction()*/
diff --git a/FEAST/MIToolbox/demonstration_algorithms/IAMB.m b/FEAST/MIToolbox/demonstration_algorithms/IAMB.m
new file mode 100644
index 0000000..1011260
--- /dev/null
+++ b/FEAST/MIToolbox/demonstration_algorithms/IAMB.m
@@ -0,0 +1,56 @@
+function [cmb association] = IAMB( data, targetindex, THRESHOLD)
+%function [cmb association] = IAMB( data, targetindex, THRESHOLD)
+%
+%Performs the IAMB algorithm of Tsmardinos et al. (2003)
+%from "Towards principled feature selection: Relevancy, filters and wrappers"
+
+if (nargin == 2)
+  THRESHOLD = 0.02;
+end
+
+numf = size(data,2);
+targets = data(:,targetindex);
+data(:,targetindex) = -10;
+
+cmb = [];
+
+finished = false;
+while ~finished
+    for n = 1:numf
+        cmbVector = joint(data(:,cmb));
+        if isempty(cmb)
+            association(n) = mi( data(:,n), targets );
+        end
+        
+        if ismember(n,cmb)
+            association(n) = -10; %arbtirary large negative constant
+        else
+            association(n) = cmi( data(:,n), targets, cmbVector);
+        end
+    end
+
+    [maxval maxidx] = max(association);
+    if maxval < THRESHOLD
+        finished = true;
+    else
+        cmb = [ cmb maxidx ];
+    end
+end
+
+finished = false;
+while ~finished && ~isempty(cmb)
+    association = [];
+    for n = 1:length(cmb)
+        cmbwithoutn = cmb;
+        cmbwithoutn(n)=[];
+        association(n) = cmi( data(:,cmb(n)), targets, data(:,cmbwithoutn) );
+    end
+    
+    [minval minidx] = min(association);
+    if minval > THRESHOLD
+        finished = true;
+    else        
+        cmb(minidx) = [];
+    end
+end
+
diff --git a/FEAST/MIToolbox/demonstration_algorithms/compile_demos.m b/FEAST/MIToolbox/demonstration_algorithms/compile_demos.m
new file mode 100644
index 0000000..65464b3
--- /dev/null
+++ b/FEAST/MIToolbox/demonstration_algorithms/compile_demos.m
@@ -0,0 +1,3 @@
+mex -I.. CMIM_Mex.c ../MutualInformation.c ../Entropy.c ../CalculateProbability.c ../ArrayOperations.c
+mex -I.. DISR_Mex.c ../MutualInformation.c ../Entropy.c ../CalculateProbability.c ../ArrayOperations.c
+mex -I.. mRMR_D_Mex.c ../MutualInformation.c ../Entropy.c ../CalculateProbability.c ../ArrayOperations.c
+\ No newline at end of file
diff --git a/FEAST/MIToolbox/demonstration_algorithms/mRMR_D.m b/FEAST/MIToolbox/demonstration_algorithms/mRMR_D.m
new file mode 100644
index 0000000..50b14bc
--- /dev/null
+++ b/FEAST/MIToolbox/demonstration_algorithms/mRMR_D.m
@@ -0,0 +1,69 @@
+function selectedFeatures = mRMR_D(k, featureMatrix, classColumn)
+%function selectedFeatures = mRMR_D(k, featureMatrix, classColumn)
+%
+%Selects optimal features according to the mRMR-D algorithm from 
+%"Feature Selection Based on Mutual Information: Criteria of Max-Dependency, Max-Relevance, and Min-Redundancy"
+%by H. Peng et al. (2005)
+%
+%Calculates the top k features
+%a dataset featureMatrix with n training examples and m features
+%with the classes held in classColumn (an n x 1 vector)
+
+noOfTraining = size(classColumn,1);
+noOfFeatures = size(featureMatrix,2);
+unselectedFeatures = ones(noOfFeatures,1);
+
+classMI = zeros(noOfFeatures,1);
+answerFeatures = zeros(k,1);
+highestMI = 0;
+highestMICounter = 0;
+currentHighestFeature = 0;
+
+featureMIMatrix = -(ones(k,noOfFeatures));
+
+%setup the mi against the class
+for n = 1 : noOfFeatures
+  classMI(n) = mi(featureMatrix(:,n),classColumn);
+	if classMI(n) > highestMI
+		highestMI = classMI(n);
+		highestMICounter = n;
+	end
+end
+
+answerFeatures(1) = highestMICounter;
+unselectedFeatures(highestMICounter) = 0;
+
+%iterate over the number of features to select
+for i = 2:k
+  score = -100;
+	currentHighestFeature = 0;
+	iMinus = i-1;
+  for j = 1 : noOfFeatures
+		if unselectedFeatures(j) == 1
+			currentMIScore = 0;
+			for m = 1 : iMinus
+        if featureMIMatrix(m,j) == -1
+          featureMIMatrix(m,j) = mi(featureMatrix(:,j),featureMatrix(:,answerFeatures(m)));
+        end
+				currentMIScore =  currentMIScore + featureMIMatrix(m,j);
+			end
+			currentScore = classMI(j) - (currentMIScore/iMinus);
+			
+			if (currentScore > score)
+				score = currentScore;
+				currentHighestFeature = j;
+			end
+		end
+	end
+	
+  if score < 0 
+    disp(['at selection ' int2str(j) ' mRMRD is negative with value ' num2str(score)]);
+  end
+
+	%now highest feature is selected in currentHighestFeature
+	%store it
+	unselectedFeatures(currentHighestFeature) = 0;
+	answerFeatures(i) = currentHighestFeature;
+end
+
+selectedFeatures = answerFeatures;
diff --git a/FEAST/MIToolbox/demonstration_algorithms/mRMR_D_Mex.c b/FEAST/MIToolbox/demonstration_algorithms/mRMR_D_Mex.c
new file mode 100644
index 0000000..8d7b074
--- /dev/null
+++ b/FEAST/MIToolbox/demonstration_algorithms/mRMR_D_Mex.c
@@ -0,0 +1,184 @@
+/*******************************************************************************
+** Demonstration feature selection algorithm - MATLAB r2009a
+**
+** Initial Version - 13/06/2008
+** Updated - 07/07/2010
+** based on mRMR_D.m
+**
+** Minimum Relevance Maximum Redundancy
+** in
+** "Feature Selection Based on Mutual Information: Criteria of Max-Dependency, Max-Relevance, and Min-Redundancy"
+** H. Peng et al. (2005)
+**
+** Author - Adam Pocock
+** Demonstration code for MIToolbox
+*******************************************************************************/
+
+#include "mex.h"
+#include "MutualInformation.h"
+
+void mRMRCalculation(int k, int noOfSamples, int noOfFeatures,double *featureMatrix, double *classColumn, double *outputFeatures)
+{
+  double **feature2D = (double**) mxCalloc(noOfFeatures,sizeof(double*));
+  /*holds the class MI values*/
+  double *classMI = (double *)mxCalloc(noOfFeatures,sizeof(double));
+  int *selectedFeatures = (int *)mxCalloc(noOfFeatures,sizeof(int));
+  /*holds the intra feature MI values*/
+  int sizeOfMatrix = k*noOfFeatures;
+  double *featureMIMatrix = (double *)mxCalloc(sizeOfMatrix,sizeof(double));
+  
+  double maxMI = 0.0;
+  int maxMICounter = -1;
+  
+  /*init variables*/
+  
+  double score, currentScore, totalFeatureMI;
+  int currentHighestFeature;
+  
+  int arrayPosition, i, j, x;
+  
+  for(j = 0; j < noOfFeatures; j++)
+  {
+    feature2D[j] = featureMatrix + (int)j*noOfSamples;
+  }
+  
+  for (i = 0; i < sizeOfMatrix;i++)
+  {
+    featureMIMatrix[i] = -1;
+  }/*for featureMIMatrix - blank to -1*/
+  
+
+  for (i = 0; i < noOfFeatures;i++)
+  {
+    classMI[i] = calculateMutualInformation(feature2D[i], classColumn, noOfSamples);
+    if (classMI[i] > maxMI)
+    {
+      maxMI = classMI[i];
+      maxMICounter = i;
+    }/*if bigger than current maximum*/
+  }/*for noOfFeatures - filling classMI*/
+  
+  selectedFeatures[maxMICounter] = 1;
+  outputFeatures[0] = maxMICounter;
+  
+  /*************
+  ** Now we have populated the classMI array, and selected the highest
+  ** MI feature as the first output feature
+  ** Now we move into the mRMR-D algorithm
+  *************/
+  
+  for (i = 1; i < k; i++)
+  {
+    /*to ensure it selects some features
+    **if this is zero then it will not pick features where the redundancy is greater than the 
+    **relevance
+    */
+    score = -1000.0;
+    currentHighestFeature = 0;
+    currentScore = 0.0;
+    totalFeatureMI = 0.0;
+    
+    for (j = 0; j < noOfFeatures; j++)
+    {
+      /*if we haven't selected j*/
+      if (selectedFeatures[j] == 0)
+      {
+        currentScore = classMI[j];
+        totalFeatureMI = 0.0;
+        
+        for (x = 0; x < i; x++)
+        {
+          arrayPosition = x*noOfFeatures + j;
+          if (featureMIMatrix[arrayPosition] == -1)
+          {
+            /*work out intra MI*/
+            
+            /*double calculateMutualInformation(double *firstVector, double *secondVector, int vectorLength);*/
+            featureMIMatrix[arrayPosition] = calculateMutualInformation(feature2D[(int) outputFeatures[x]], feature2D[j], noOfSamples);
+          }
+          
+          totalFeatureMI += featureMIMatrix[arrayPosition];
+        }/*for the number of already selected features*/
+        
+        currentScore -= (totalFeatureMI/i);
+			  if (currentScore > score)
+			  {
+				  score = currentScore;
+				  currentHighestFeature = j;
+			  }
+			}/*if j is unselected*/
+	  }/*for number of features*/
+  
+    selectedFeatures[currentHighestFeature] = 1;
+    outputFeatures[i] = currentHighestFeature;
+  
+  }/*for the number of features to select*/
+  
+  for (i = 0; i < k; i++)
+  {
+    outputFeatures[i] += 1; /*C indexes from 0 not 1*/
+  }/*for number of selected features*/
+  
+}/*mRMRCalculation(double[][],double[])*/
+
+/*entry point for the mex call
+**nlhs - number of outputs
+**plhs - pointer to array of outputs
+**nrhs - number of inputs
+**prhs - pointer to array of inputs
+*/
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
+{
+  /*************************************************************
+  ** this function takes 3 arguments:
+  ** k = number of features to select,
+  ** featureMatrix[][] = matrix of features
+  ** classColumn[] = targets
+  ** the arguments should all be discrete integers.
+  ** and has one output:
+  ** selectedFeatures[] of size k
+  *************************************************************/
+  
+  int k, numberOfFeatures, numberOfSamples, numberOfTargets;
+  double *featureMatrix, *targets, *output;
+  
+  
+  if (nlhs != 1)
+  {
+    printf("Incorrect number of output arguments");
+  }/*if not 1 output*/
+  if (nrhs != 3)
+  {
+    printf("Incorrect number of input arguments");
+  }/*if not 3 inputs*/
+  
+  /*get the number of features to select, cast out as it is a double*/
+  k = (int) mxGetScalar(prhs[0]);
+
+  numberOfFeatures = mxGetN(prhs[1]);
+  numberOfSamples = mxGetM(prhs[1]);
+  
+  numberOfTargets = mxGetM(prhs[2]);
+  
+  if (numberOfTargets != numberOfSamples)
+  {
+    printf("Number of targets must match number of samples\n");
+    printf("Number of targets = %d, Number of Samples = %d, Number of Features = %d\n",numberOfTargets,numberOfSamples,numberOfFeatures);
+    
+    plhs[0] = mxCreateDoubleMatrix(0,0,mxREAL);
+  }/*if size mismatch*/
+  else
+  {
+    
+    featureMatrix = mxGetPr(prhs[1]);
+    targets = mxGetPr(prhs[2]);
+    
+    plhs[0] = mxCreateDoubleMatrix(k,1,mxREAL);
+    output = (double *)mxGetPr(plhs[0]);
+    
+    /*void mRMRCalculation(int k, int noOfSamples, int noOfFeatures,double *featureMatrix, double *classColumn, double *outputFeatures)*/
+    mRMRCalculation(k,numberOfSamples,numberOfFeatures,featureMatrix,targets,output);
+  }
+  
+  return;
+}/*mexFunction()*/