pax_global_header00006660000000000000000000000064126445441770014530gustar00rootroot0000000000000052 comment=efccae11bff5394a77723f6d89d819a8cbf2f554 MIToolbox-2.1.2/000077500000000000000000000000001264454417700134065ustar00rootroot00000000000000MIToolbox-2.1.2/.gitignore000066400000000000000000000003201264454417700153710ustar00rootroot00000000000000# Object files *.o *.ko # Libraries *.lib *.a # Shared objects (inc. Windows DLLs) *.dll *.so *.so.* *.dylib # Executables *.exe *.out *.app # Matlab/Octave executables *.mex* *.oct* # Editor backups *~ MIToolbox-2.1.2/ArrayOperations.c000066400000000000000000000211221264454417700166720ustar00rootroot00000000000000/******************************************************************************* ** ArrayOperations.c ** Part of the mutual information toolbox ** ** Contains functions to floor arrays, and to merge arrays into a joint ** state. ** ** Author: Adam Pocock ** Created 17/2/2010 ** Updated - 22/02/2014 - Added checking on calloc. ** ** Copyright 2010,2014 Adam Pocock, The University Of Manchester ** www.cs.manchester.ac.uk ** ** This file is part of MIToolbox, licensed under the 3-clause BSD license. *******************************************************************************/ #include #include "MIToolbox.h" #include "ArrayOperations.h" void* checkedCalloc(size_t vectorLength, size_t sizeOfType) { void *allocated = CALLOC_FUNC(vectorLength, sizeOfType); if(allocated == NULL) { #ifdef MEX_IMPLEMENTATION /* This call returns control to Matlab, with the associated error message */ mexErrMsgTxt("Failed to allocate memory\n"); #else fprintf(stderr, "Error: %s\n", strerror(errno)); exit(EXIT_FAILURE); #endif } return allocated; } void incrementVector(double* vector, int vectorLength) { int i = 0; for (i = 0; i < vectorLength; i++) { vector[i]++; /*C indexes from 0 not 1*/ }/*for length of array */ }/* incrementVector(double*,int) */ void printDoubleVector(double *vector, int vectorLength) { int i; for (i = 0; i < vectorLength; i++) { if (vector[i] > 0) printf("Val at i=%d, is %f\n",i,vector[i]); }/*for number of items in vector*/ }/*printDoubleVector(double*,int)*/ void printIntVector(int *vector, int vectorLength) { int i; for (i = 0; i < vectorLength; i++) { printf("Val at i=%d, is %d\n",i,vector[i]); }/*for number of items in vector*/ }/*printIntVector(int*,int)*/ int numberOfUniqueValues(double *featureVector, int vectorLength) { int uniqueValues = 0; double *valuesArray = (double *) checkedCalloc(vectorLength,sizeof(double)); int found = 0; int j = 0; int i; for (i = 0; i < vectorLength; i++) { found = 0; j = 0; while ((j < uniqueValues) && (found == 0)) { if (valuesArray[j] == featureVector[i]) { found = 1; featureVector[i] = (double) (j+1); } j++; } if (!found) { valuesArray[uniqueValues] = featureVector[i]; uniqueValues++; featureVector[i] = (double) uniqueValues; } }/*for vectorlength*/ FREE_FUNC(valuesArray); valuesArray = NULL; return uniqueValues; }/*numberOfUniqueValues(double*,int)*/ /******************************************************************************* ** normaliseArray takes an input vector and writes an output vector ** which is a normalised version of the input, and returns the number of states ** A normalised array has min value = 0, max value = number of states ** and all values are integers ** ** length(inputVector) == length(outputVector) == vectorLength otherwise there ** is a memory leak *******************************************************************************/ int normaliseArray(double *inputVector, int *outputVector, int vectorLength) { int minVal = 0; int maxVal = 0; int currentValue; int i; if (vectorLength > 0) { minVal = (int) floor(inputVector[0]); maxVal = (int) floor(inputVector[0]); for (i = 0; i < vectorLength; i++) { currentValue = (int) floor(inputVector[i]); outputVector[i] = currentValue; if (currentValue < minVal) { minVal = currentValue; } else if (currentValue > maxVal) { maxVal = currentValue; } }/*for loop over vector*/ for (i = 0; i < vectorLength; i++) { outputVector[i] = outputVector[i] - minVal; } maxVal = (maxVal - minVal) + 1; } return maxVal; }/*normaliseArray(double*,double*,int)*/ /******************************************************************************* ** mergeArrays takes in two arrays and writes the joint state of those arrays ** to the output vector, returning the number of joint states ** ** the length of the vectors must be the same and equal to vectorLength ** outputVector must be malloc'd before calling this function *******************************************************************************/ int mergeArrays(double *firstVector, double *secondVector, double *outputVector, int vectorLength) { int *firstNormalisedVector; int *secondNormalisedVector; int firstNumStates; int secondNumStates; int i; int *stateMap; int stateCount; int curIndex; firstNormalisedVector = (int *) checkedCalloc(vectorLength,sizeof(int)); secondNormalisedVector = (int *) checkedCalloc(vectorLength,sizeof(int)); firstNumStates = normaliseArray(firstVector,firstNormalisedVector,vectorLength); secondNumStates = normaliseArray(secondVector,secondNormalisedVector,vectorLength); /* ** printVector(firstNormalisedVector,vectorLength); ** printVector(secondNormalisedVector,vectorLength); */ stateMap = (int *) checkedCalloc(firstNumStates*secondNumStates,sizeof(int)); stateCount = 1; for (i = 0; i < vectorLength; i++) { curIndex = firstNormalisedVector[i] + (secondNormalisedVector[i] * firstNumStates); if (stateMap[curIndex] == 0) { stateMap[curIndex] = stateCount; stateCount++; } outputVector[i] = stateMap[curIndex]; } FREE_FUNC(firstNormalisedVector); FREE_FUNC(secondNormalisedVector); FREE_FUNC(stateMap); firstNormalisedVector = NULL; secondNormalisedVector = NULL; stateMap = NULL; /*printVector(outputVector,vectorLength);*/ return stateCount; }/*mergeArrays(double *,double *,double *, int, bool)*/ int mergeArraysArities(double *firstVector, int numFirstStates, double *secondVector, int numSecondStates, double *outputVector, int vectorLength) { int *firstNormalisedVector; int *secondNormalisedVector; int i; int totalStates; int firstStateCheck, secondStateCheck; firstNormalisedVector = (int *) checkedCalloc(vectorLength,sizeof(int)); secondNormalisedVector = (int *) checkedCalloc(vectorLength,sizeof(int)); firstStateCheck = normaliseArray(firstVector,firstNormalisedVector,vectorLength); secondStateCheck = normaliseArray(secondVector,secondNormalisedVector,vectorLength); if ((firstStateCheck <= numFirstStates) && (secondStateCheck <= numSecondStates)) { for (i = 0; i < vectorLength; i++) { outputVector[i] = firstNormalisedVector[i] + (secondNormalisedVector[i] * numFirstStates) + 1; } totalStates = numFirstStates * numSecondStates; } else { totalStates = -1; } FREE_FUNC(firstNormalisedVector); FREE_FUNC(secondNormalisedVector); firstNormalisedVector = NULL; secondNormalisedVector = NULL; return totalStates; }/*mergeArraysArities(double *,int,double *,int,double *,int)*/ int mergeMultipleArrays(double *inputMatrix, double *outputVector, int matrixWidth, int vectorLength) { int i = 0; int currentIndex; int currentNumStates; int *normalisedVector; if (matrixWidth > 1) { currentNumStates = mergeArrays(inputMatrix, (inputMatrix + vectorLength), outputVector,vectorLength); for (i = 2; i < matrixWidth; i++) { currentIndex = i * vectorLength; currentNumStates = mergeArrays(outputVector,(inputMatrix + currentIndex),outputVector,vectorLength); } } else { normalisedVector = (int *) checkedCalloc(vectorLength,sizeof(int)); currentNumStates = normaliseArray(inputMatrix,normalisedVector,vectorLength); for (i = 0; i < vectorLength; i++) { outputVector[i] = inputMatrix[i]; } } return currentNumStates; }/*mergeMultipleArrays(double *, double *, int, int, bool)*/ int mergeMultipleArraysArities(double *inputMatrix, double *outputVector, int matrixWidth, int *arities, int vectorLength) { int i = 0; int currentIndex; int currentNumStates; int *normalisedVector; if (matrixWidth > 1) { currentNumStates = mergeArraysArities(inputMatrix, arities[0], (inputMatrix + vectorLength), arities[1], outputVector,vectorLength); for (i = 2; i < matrixWidth; i++) { currentIndex = i * vectorLength; currentNumStates = mergeArraysArities(outputVector,currentNumStates,(inputMatrix + currentIndex),arities[i],outputVector,vectorLength); if (currentNumStates == -1) break; } } else { normalisedVector = (int *) checkedCalloc(vectorLength,sizeof(int)); currentNumStates = normaliseArray(inputMatrix,normalisedVector,vectorLength); for (i = 0; i < vectorLength; i++) { outputVector[i] = inputMatrix[i]; } } return currentNumStates; }/*mergeMultipleArraysArities(double *, double *, int, int, bool)*/ MIToolbox-2.1.2/ArrayOperations.h000066400000000000000000000074541264454417700167130ustar00rootroot00000000000000/******************************************************************************* ** ArrayOperations.h ** Part of the mutual information toolbox ** ** Contains functions to floor arrays, and to merge arrays into a joint ** state. ** ** Author: Adam Pocock ** Created 17/2/2010 ** Updated - 22/02/2014 - Added checking on calloc, and an increment array function. ** ** Copyright 2010,2014 Adam Pocock, The University Of Manchester ** www.cs.manchester.ac.uk ** ** This file is part of MIToolbox, licensed under the 3-clause BSD license. *******************************************************************************/ #ifndef __ArrayOperations_H #define __ArrayOperations_H #ifdef __cplusplus extern "C" { #endif /******************************************************************************* ** A version of calloc which checks to see if memory was allocated. *******************************************************************************/ void* checkedCalloc(size_t vectorLength, size_t sizeOfType); /******************************************************************************* ** Increments each value in a double array *******************************************************************************/ void incrementVector(double* vector, int vectorLength); /******************************************************************************* ** Simple print functions for debugging *******************************************************************************/ void printDoubleVector(double *vector, int vectorlength); void printIntVector(int *vector, int vectorLength); /******************************************************************************* ** numberOfUniqueValues finds the number of unique values in an array by ** repeatedly iterating through the array and checking if a value has been ** seen previously *******************************************************************************/ int numberOfUniqueValues(double *featureVector, int vectorLength); /******************************************************************************* ** normaliseArray takes an input vector and writes an output vector ** which is a normalised version of the input, and returns the number of states ** A normalised array has min value = 0, max value = number of states ** and all values are integers ** ** length(inputVector) == length(outputVector) == vectorLength otherwise there ** is a memory leak *******************************************************************************/ int normaliseArray(double *inputVector, int *outputVector, int vectorLength); /******************************************************************************* ** mergeArrays takes in two arrays and writes the joint state of those arrays ** to the output vector ** ** the length of the vectors must be the same and equal to vectorLength *******************************************************************************/ int mergeArrays(double *firstVector, double *secondVector, double *outputVector, int vectorLength); int mergeArraysArities(double *firstVector, int numFirstStates, double *secondVector, int numSecondStates, double *outputVector, int vectorLength); /******************************************************************************* ** mergeMultipleArrays takes in a matrix and repeatedly merges the matrix using ** merge arrays and writes the joint state of that matrix ** to the output vector ** ** the length of the vectors must be the same and equal to vectorLength ** matrixWidth = the number of columns in the matrix *******************************************************************************/ int mergeMultipleArrays(double *inputMatrix, double *outputVector, int matrixWidth, int vectorLength); int mergeMultipleArraysArities(double *inputMatrix, double *outputVector, int matrixWidth, int *arities, int vectorLength); #ifdef __cplusplus } #endif #endif MIToolbox-2.1.2/CalculateProbability.c000066400000000000000000000243621264454417700176570ustar00rootroot00000000000000/******************************************************************************* ** CalculateProbability.c ** Part of the mutual information toolbox ** ** Contains functions to calculate the probability of each state in the array ** and to calculate the probability of the joint state of two arrays ** ** Author: Adam Pocock ** Created: 17/02/2010 ** Modified - 04/07/2011 - added weighted probability functions ** Updated - 22/02/2014 - Added checking on calloc. ** ** Copyright 2010-2014 Adam Pocock, The University Of Manchester ** www.cs.manchester.ac.uk ** ** This file is part of MIToolbox, licensed under the 3-clause BSD license. *******************************************************************************/ #include "MIToolbox.h" #include "ArrayOperations.h" #include "CalculateProbability.h" JointProbabilityState calculateJointProbability(double *firstVector, double *secondVector, int vectorLength) { int *firstNormalisedVector; int *secondNormalisedVector; int *firstStateCounts; int *secondStateCounts; int *jointStateCounts; double *firstStateProbs; double *secondStateProbs; double *jointStateProbs; int firstNumStates; int secondNumStates; int jointNumStates; int i; double length = vectorLength; JointProbabilityState state; firstNormalisedVector = (int *) checkedCalloc(vectorLength,sizeof(int)); secondNormalisedVector = (int *) checkedCalloc(vectorLength,sizeof(int)); firstNumStates = normaliseArray(firstVector,firstNormalisedVector,vectorLength); secondNumStates = normaliseArray(secondVector,secondNormalisedVector,vectorLength); jointNumStates = firstNumStates * secondNumStates; firstStateCounts = (int *) checkedCalloc(firstNumStates,sizeof(int)); secondStateCounts = (int *) checkedCalloc(secondNumStates,sizeof(int)); jointStateCounts = (int *) checkedCalloc(jointNumStates,sizeof(int)); firstStateProbs = (double *) checkedCalloc(firstNumStates,sizeof(double)); secondStateProbs = (double *) checkedCalloc(secondNumStates,sizeof(double)); jointStateProbs = (double *) checkedCalloc(jointNumStates,sizeof(double)); /* optimised version, less numerically stable double fractionalState = 1.0 / vectorLength; for (i = 0; i < vectorLength; i++) { firstStateProbs[firstNormalisedVector[i]] += fractionalState; secondStateProbs[secondNormalisedVector[i]] += fractionalState; jointStateProbs[secondNormalisedVector[i] * firstNumStates + firstNormalisedVector[i]] += fractionalState; } */ /* Optimised for number of FP operations now O(states) instead of O(vectorLength) */ for (i = 0; i < vectorLength; i++) { firstStateCounts[firstNormalisedVector[i]] += 1; secondStateCounts[secondNormalisedVector[i]] += 1; jointStateCounts[secondNormalisedVector[i] * firstNumStates + firstNormalisedVector[i]] += 1; } for (i = 0; i < firstNumStates; i++) { firstStateProbs[i] = firstStateCounts[i] / length; } for (i = 0; i < secondNumStates; i++) { secondStateProbs[i] = secondStateCounts[i] / length; } for (i = 0; i < jointNumStates; i++) { jointStateProbs[i] = jointStateCounts[i] / length; } FREE_FUNC(firstNormalisedVector); FREE_FUNC(secondNormalisedVector); FREE_FUNC(firstStateCounts); FREE_FUNC(secondStateCounts); FREE_FUNC(jointStateCounts); firstNormalisedVector = NULL; secondNormalisedVector = NULL; firstStateCounts = NULL; secondStateCounts = NULL; jointStateCounts = NULL; /* **typedef struct **{ ** double *jointProbabilityVector; ** int numJointStates; ** double *firstProbabilityVector; ** int numFirstStates; ** double *secondProbabilityVector; ** int numSecondStates; **} JointProbabilityState; */ state.jointProbabilityVector = jointStateProbs; state.numJointStates = jointNumStates; state.firstProbabilityVector = firstStateProbs; state.numFirstStates = firstNumStates; state.secondProbabilityVector = secondStateProbs; state.numSecondStates = secondNumStates; return state; }/*calculateJointProbability(double *,double *, int)*/ WeightedJointProbState calculateWeightedJointProbability(double *firstVector, double *secondVector, double *weightVector, int vectorLength) { int *firstNormalisedVector; int *secondNormalisedVector; int *firstStateCounts; int *secondStateCounts; int *jointStateCounts; double *firstStateProbs; double *secondStateProbs; double *jointStateProbs; double *firstWeightVec; double *secondWeightVec; double *jointWeightVec; int firstNumStates; int secondNumStates; int jointNumStates; int i; double length = vectorLength; WeightedJointProbState state; firstNormalisedVector = (int *) checkedCalloc(vectorLength,sizeof(int)); secondNormalisedVector = (int *) checkedCalloc(vectorLength,sizeof(int)); firstNumStates = normaliseArray(firstVector,firstNormalisedVector,vectorLength); secondNumStates = normaliseArray(secondVector,secondNormalisedVector,vectorLength); jointNumStates = firstNumStates * secondNumStates; firstStateCounts = (int *) checkedCalloc(firstNumStates,sizeof(int)); secondStateCounts = (int *) checkedCalloc(secondNumStates,sizeof(int)); jointStateCounts = (int *) checkedCalloc(jointNumStates,sizeof(int)); firstStateProbs = (double *) checkedCalloc(firstNumStates,sizeof(double)); secondStateProbs = (double *) checkedCalloc(secondNumStates,sizeof(double)); jointStateProbs = (double *) checkedCalloc(jointNumStates,sizeof(double)); firstWeightVec = (double *) checkedCalloc(firstNumStates,sizeof(double)); secondWeightVec = (double *) checkedCalloc(secondNumStates,sizeof(double)); jointWeightVec = (double *) checkedCalloc(jointNumStates,sizeof(double)); for (i = 0; i < vectorLength; i++) { firstStateCounts[firstNormalisedVector[i]] += 1; secondStateCounts[secondNormalisedVector[i]] += 1; jointStateCounts[secondNormalisedVector[i] * firstNumStates + firstNormalisedVector[i]] += 1; firstWeightVec[firstNormalisedVector[i]] += weightVector[i]; secondWeightVec[secondNormalisedVector[i]] += weightVector[i]; jointWeightVec[secondNormalisedVector[i] * firstNumStates + firstNormalisedVector[i]] += weightVector[i]; } for (i = 0; i < firstNumStates; i++) { if (firstStateCounts[i]) { firstStateProbs[i] = firstStateCounts[i] / length; firstWeightVec[i] /= firstStateCounts[i]; } } for (i = 0; i < secondNumStates; i++) { if (secondStateCounts[i]) { secondStateProbs[i] = secondStateCounts[i] / length; secondWeightVec[i] /= secondStateCounts[i]; } } for (i = 0; i < jointNumStates; i++) { if (jointStateCounts[i]) { jointStateProbs[i] = jointStateCounts[i] / length; jointWeightVec[i] /= jointStateCounts[i]; } } FREE_FUNC(firstNormalisedVector); FREE_FUNC(secondNormalisedVector); FREE_FUNC(firstStateCounts); FREE_FUNC(secondStateCounts); FREE_FUNC(jointStateCounts); firstNormalisedVector = NULL; secondNormalisedVector = NULL; firstStateCounts = NULL; secondStateCounts = NULL; jointStateCounts = NULL; /* **typedef struct **{ ** double *jointProbabilityVector; ** double *jointWeightVector; ** int numJointStates; ** double *firstProbabilityVector; ** double *firstWeightVector; ** int numFirstStates; ** double *secondProbabilityVector; ** double *secondWeightVector; ** int numSecondStates; **} WeightedJointProbState; */ state.jointProbabilityVector = jointStateProbs; state.jointWeightVector = jointWeightVec; state.numJointStates = jointNumStates; state.firstProbabilityVector = firstStateProbs; state.firstWeightVector = firstWeightVec; state.numFirstStates = firstNumStates; state.secondProbabilityVector = secondStateProbs; state.secondWeightVector = secondWeightVec; state.numSecondStates = secondNumStates; return state; }/*calculateJointProbability(double *,double *, int)*/ ProbabilityState calculateProbability(double *dataVector, int vectorLength) { int *normalisedVector; int *stateCounts; double *stateProbs; int numStates; /*double fractionalState;*/ ProbabilityState state; int i; double length = vectorLength; normalisedVector = (int *) checkedCalloc(vectorLength,sizeof(int)); numStates = normaliseArray(dataVector,normalisedVector,vectorLength); stateCounts = (int *) checkedCalloc(numStates,sizeof(int)); stateProbs = (double *) checkedCalloc(numStates,sizeof(double)); /* optimised version, may have floating point problems fractionalState = 1.0 / vectorLength; for (i = 0; i < vectorLength; i++) { stateProbs[normalisedVector[i]] += fractionalState; } */ /* Optimised for number of FP operations now O(states) instead of O(vectorLength) */ for (i = 0; i < vectorLength; i++) { stateCounts[normalisedVector[i]] += 1; } for (i = 0; i < numStates; i++) { stateProbs[i] = stateCounts[i] / length; } FREE_FUNC(stateCounts); FREE_FUNC(normalisedVector); stateCounts = NULL; normalisedVector = NULL; state.probabilityVector = stateProbs; state.numStates = numStates; return state; }/*calculateProbability(double *,int)*/ WeightedProbState calculateWeightedProbability(double *dataVector, double *exampleWeightVector, int vectorLength) { int *normalisedVector; int *stateCounts; double *stateProbs; double *stateWeights; int numStates; WeightedProbState state; int i; double length = vectorLength; normalisedVector = (int *) checkedCalloc(vectorLength,sizeof(int)); numStates = normaliseArray(dataVector,normalisedVector,vectorLength); stateCounts = (int *) checkedCalloc(numStates,sizeof(int)); stateProbs = (double *) checkedCalloc(numStates,sizeof(double)); stateWeights = (double *) checkedCalloc(numStates,sizeof(double)); for (i = 0; i < vectorLength; i++) { stateCounts[normalisedVector[i]] += 1; stateWeights[normalisedVector[i]] += exampleWeightVector[i]; } for (i = 0; i < numStates; i++) { stateProbs[i] = stateCounts[i] / length; stateWeights[i] /= stateCounts[i]; } FREE_FUNC(stateCounts); FREE_FUNC(normalisedVector); stateCounts = NULL; normalisedVector = NULL; state.probabilityVector = stateProbs; state.stateWeightVector = stateWeights; state.numStates = numStates; return state; }/*calculateProbability(double *,int)*/ MIToolbox-2.1.2/CalculateProbability.h000066400000000000000000000067401264454417700176640ustar00rootroot00000000000000/******************************************************************************* ** CalculateProbability.h ** Part of the mutual information toolbox ** ** Contains functions to calculate the probability of each state in the array ** and to calculate the probability of the joint state of two arrays ** ** Author: Adam Pocock ** Created 17/2/2010 ** ** Copyright 2010 Adam Pocock, The University Of Manchester ** www.cs.manchester.ac.uk ** ** This file is part of MIToolbox, licensed under the 3-clause BSD license. *******************************************************************************/ #ifndef __CalculateProbability_H #define __CalculateProbability_H #ifdef __cplusplus extern "C" { #endif typedef struct jpState { double *jointProbabilityVector; int numJointStates; double *firstProbabilityVector; int numFirstStates; double *secondProbabilityVector; int numSecondStates; } JointProbabilityState; typedef struct pState { double *probabilityVector; int numStates; } ProbabilityState; typedef struct wjpState { double *jointProbabilityVector; double *jointWeightVector; int numJointStates; double *firstProbabilityVector; double *firstWeightVector; int numFirstStates; double *secondProbabilityVector; double *secondWeightVector; int numSecondStates; } WeightedJointProbState; typedef struct wpState { double *probabilityVector; double *stateWeightVector; int numStates; } WeightedProbState; /******************************************************************************* ** calculateJointProbability returns the joint probability vector of two vectors ** and the marginal probability vectors in a struct. ** It is the base operation for all information theory calculations involving ** two or more variables. ** ** length(firstVector) == length(secondVector) == vectorLength otherwise there ** will be a segmentation fault *******************************************************************************/ JointProbabilityState calculateJointProbability(double *firstVector, double *secondVector, int vectorLength); /******************************************************************************* ** calculateWeightedJointProbability returns the same joint state object ** as calculateJointProbability, with additional weightVectors giving the ** weight assigned to each state. *******************************************************************************/ WeightedJointProbState calculateWeightedJointProbability(double *firstVector, double *secondVector, double *exampleWeightVector, int vectorLength); /******************************************************************************* ** calculateProbability returns the probability vector from one vector. ** It is the base operation for all information theory calculations involving ** one variable ** ** length(dataVector) == vectorLength otherwise there ** will be a segmentation fault *******************************************************************************/ ProbabilityState calculateProbability(double *dataVector, int vectorLength); /******************************************************************************* ** calculateWeightedProbability returns the same state object ** as calculateProbability, with an additional weightVector giving the ** weight assigned to each state. *******************************************************************************/ WeightedProbState calculateWeightedProbability(double *dataVector, double *exampleWeightVector, int vectorLength); #ifdef __cplusplus } #endif #endif MIToolbox-2.1.2/CompileMIToolbox.m000066400000000000000000000005311264454417700167500ustar00rootroot00000000000000% Compiles the MIToolbox functions mex MIToolboxMex.c MutualInformation.c Entropy.c CalculateProbability.c ArrayOperations.c mex RenyiMIToolboxMex.c RenyiMutualInformation.c RenyiEntropy.c CalculateProbability.c ArrayOperations.c mex WeightedMIToolboxMex.c WeightedMutualInformation.c WeightedEntropy.c CalculateProbability.c ArrayOperations.c MIToolbox-2.1.2/Entropy.c000066400000000000000000000064311264454417700152160ustar00rootroot00000000000000/******************************************************************************* ** Entropy.c ** Part of the mutual information toolbox ** ** Contains functions to calculate the entropy of a single variable H(X), ** the joint entropy of two variables H(X,Y), and the conditional entropy ** H(X|Y) ** ** Author: Adam Pocock ** Created 19/2/2010 ** ** Copyright 2010 Adam Pocock, The University Of Manchester ** www.cs.manchester.ac.uk ** ** This file is part of MIToolbox. ** ** This file is part of MIToolbox, licensed under the 3-clause BSD license. *******************************************************************************/ #include "MIToolbox.h" #include "CalculateProbability.h" #include "Entropy.h" double calculateEntropy(double *dataVector, int vectorLength) { double entropy = 0.0; double tempValue = 0.0; int i; ProbabilityState state = calculateProbability(dataVector,vectorLength); /*H(X) = - sum p(x) log p(x)*/ for (i = 0; i < state.numStates; i++) { tempValue = state.probabilityVector[i]; if (tempValue > 0) { entropy -= tempValue * log(tempValue); } } entropy /= log(2.0); FREE_FUNC(state.probabilityVector); state.probabilityVector = NULL; return entropy; }/*calculateEntropy(double *,int)*/ double calculateJointEntropy(double *firstVector, double *secondVector, int vectorLength) { double jointEntropy = 0.0; double tempValue = 0.0; int i; JointProbabilityState state = calculateJointProbability(firstVector,secondVector,vectorLength); /*H(XY) = - sumx sumy p(xy) log p(xy)*/ for (i = 0; i < state.numJointStates; i++) { tempValue = state.jointProbabilityVector[i]; if (tempValue > 0) { jointEntropy -= tempValue * log(tempValue); } } jointEntropy /= log(2.0); FREE_FUNC(state.firstProbabilityVector); state.firstProbabilityVector = NULL; FREE_FUNC(state.secondProbabilityVector); state.secondProbabilityVector = NULL; FREE_FUNC(state.jointProbabilityVector); state.jointProbabilityVector = NULL; return jointEntropy; }/*calculateJointEntropy(double *, double *, int)*/ double calculateConditionalEntropy(double *dataVector, double *conditionVector, int vectorLength) { /* ** Conditional entropy ** H(X|Y) = - sumx sumy p(xy) log p(xy)/p(y) */ double condEntropy = 0.0; double jointValue = 0.0; double condValue = 0.0; int i; JointProbabilityState state = calculateJointProbability(dataVector,conditionVector,vectorLength); /*H(X|Y) = - sumx sumy p(xy) log p(xy)/p(y)*/ /* to index by numFirstStates use modulus of i ** to index by numSecondStates use integer division of i by numFirstStates */ for (i = 0; i < state.numJointStates; i++) { jointValue = state.jointProbabilityVector[i]; condValue = state.secondProbabilityVector[i / state.numFirstStates]; if ((jointValue > 0) && (condValue > 0)) { condEntropy -= jointValue * log(jointValue / condValue); } } condEntropy /= log(2.0); FREE_FUNC(state.firstProbabilityVector); state.firstProbabilityVector = NULL; FREE_FUNC(state.secondProbabilityVector); state.secondProbabilityVector = NULL; FREE_FUNC(state.jointProbabilityVector); state.jointProbabilityVector = NULL; return condEntropy; }/*calculateConditionalEntropy(double *, double *, int)*/ MIToolbox-2.1.2/Entropy.h000066400000000000000000000040641264454417700152230ustar00rootroot00000000000000/******************************************************************************* ** Entropy.h ** Part of the mutual information toolbox ** ** Contains functions to calculate the entropy of a single variable H(X), ** the joint entropy of two variables H(X,Y), and the conditional entropy ** H(X|Y) ** ** Author: Adam Pocock ** Created 19/2/2010 ** ** Copyright 2010 Adam Pocock, The University Of Manchester ** www.cs.manchester.ac.uk ** ** This file is part of MIToolbox, licensed under the 3-clause BSD license. *******************************************************************************/ #ifndef __Entropy_H #define __Entropy_H #ifdef __cplusplus extern "C" { #endif /******************************************************************************* ** calculateEntropy returns the entropy in log base 2 of dataVector ** H(X) ** ** length(dataVector) == vectorLength otherwise there ** will be a segmentation fault *******************************************************************************/ double calculateEntropy(double *dataVector, int vectorLength); /******************************************************************************* ** calculateJointEntropy returns the entropy in log base 2 of the joint ** variable of firstVector and secondVector H(XY) ** ** length(firstVector) == length(secondVector) == vectorLength otherwise there ** will be a segmentation fault *******************************************************************************/ double calculateJointEntropy(double *firstVector, double *secondVector, int vectorLength); /******************************************************************************* ** calculateConditionalEntropy returns the entropy in log base 2 of dataVector ** conditioned on conditionVector, H(X|Y) ** ** length(dataVector) == length(conditionVector) == vectorLength otherwise there ** will be a segmentation fault *******************************************************************************/ double calculateConditionalEntropy(double *dataVector, double *conditionVector, int vectorLength); #ifdef __cplusplus } #endif #endif MIToolbox-2.1.2/LICENSE000066400000000000000000000031121264454417700144100ustar00rootroot00000000000000This is MIToolbox, a mutual information package for C/C++ and MATLAB. Copyright (c) 2010-2015, Adam Pocock, The University of Manchester All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: - Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. - Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. - Neither the name of MIToolbox nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. MIToolbox-2.1.2/MIToolbox.h000066400000000000000000000021511264454417700154320ustar00rootroot00000000000000/******************************************************************************* ** ** MIToolbox.h ** Provides the header files and #defines to ensure compatibility with MATLAB ** and C/C++. By default it compiles to MATLAB, if COMPILE_C is defined it ** links to the C memory allocation functions. ** ** Author: Adam Pocock ** Created: 17/2/2010 ** Modified: 24/06/2011 - added log base #define ** ** Copyright 2010/2011 Adam Pocock, The University Of Manchester ** www.cs.manchester.ac.uk ** ** This file is part of MIToolbox, licensed under the 3-clause BSD license. *******************************************************************************/ #ifndef __MIToolbox_H #define __MIToolbox_H #include #include #define BASE_TWO 2.0 #define BASE_E M_E #define LOG_BASE BASE_TWO #ifdef COMPILE_C #define C_IMPLEMENTATION #include #include #define CALLOC_FUNC calloc #define FREE_FUNC free #else #define MEX_IMPLEMENTATION #include "mex.h" #define CALLOC_FUNC mxCalloc #define FREE_FUNC mxFree #define printf mexPrintf /*for Octave-3.2*/ #endif #endif MIToolbox-2.1.2/MIToolbox.m000066400000000000000000000052031264454417700154400ustar00rootroot00000000000000function [varargout] = MIToolbox(functionName, varargin) %function [varargout] = MIToolbox(functionName, varargin) % %Provides access to the functions in MIToolboxMex % %Expects column vectors, will not work with row vectors % %Function list %"joint" = joint variable of the matrix %"entropy" or "h" = H(X) %"ConditionalEntropy" or "condh" = H(X|Y) %"mi" = I(X;Y) %"ConditionalMI" or "cmi" = I(X;Y|Z) % %Arguments and returned values %[jointVariable] = joint(matrix) %[entropy] = H(X) = H(vector) %[entropy] = H(X|Y) = H(vector,condition) %[mi] = I(X;Y) = I(vector,target) %[mi] = I(X;Y|Z) = I(vector,target,condition) % %Internal MIToolbox function number %Joint = 3 %Entropy = 4 %Conditional Entropy = 6 %Mutual Information = 7 %Conditional MI = 8 for i = 1:length(varargin) if (~isa(varargin{i},'double')) error('Error, MIToolbox requires inputs to be double vector or matrices') end end if (strcmpi(functionName,'Joint') || strcmpi(functionName,'Merge')) [varargout{1}] = MIToolboxMex(3,varargin{1}); elseif (strcmpi(functionName,'Entropy') || strcmpi(functionName,'h')) %disp('Calculating Entropy'); if (size(varargin{1},2)>1) mergedVector = MIToolboxMex(3,varargin{1}); else mergedVector = varargin{1}; end [varargout{1}] = MIToolboxMex(4,mergedVector); elseif ((strcmpi(functionName,'ConditionalEntropy')) || strcmpi(functionName,'condh')) if (size(varargin{1},2)>1) mergedFirst = MIToolboxMex(3,varargin{1}); else mergedFirst = varargin{1}; end if (size(varargin{2},2)>1) mergedSecond = MIToolboxMex(3,varargin{2}); else mergedSecond = varargin{2}; end [varargout{1}] = MIToolboxMex(6,mergedFirst,mergedSecond); elseif (strcmpi(functionName,'mi')) if (size(varargin{1},2)>1) mergedFirst = MIToolboxMex(3,varargin{1}); else mergedFirst = varargin{1}; end if (size(varargin{2},2)>1) mergedSecond = MIToolboxMex(3,varargin{2}); else mergedSecond = varargin{2}; end [varargout{1}] = MIToolboxMex(7,mergedFirst,mergedSecond); elseif (strcmpi(functionName,'ConditionalMI') || strcmpi(functionName,'cmi')) if (size(varargin{1},2)>1) mergedFirst = MIToolboxMex(3,varargin{1}); else mergedFirst = varargin{1}; end if (size(varargin{2},2)>1) mergedSecond = MIToolboxMex(3,varargin{2}); else mergedSecond = varargin{2}; end if (size(varargin{3},2)>1) mergedThird = MIToolboxMex(3,varargin{3}); else mergedThird = varargin{3}; end [varargout{1}] = MIToolboxMex(8,mergedFirst,mergedSecond,mergedThird); else varargout{1} = 0; disp(['Unrecognised functionName ' functionName]); end MIToolbox-2.1.2/MIToolboxMex.c000066400000000000000000000346051264454417700161100ustar00rootroot00000000000000/******************************************************************************* ** MIToolboxMex.c ** is the MATLAB entry point for the MIToolbox functions when called from ** a MATLAB/OCTAVE script. ** ** Copyright 2010 Adam Pocock, The University Of Manchester ** www.cs.manchester.ac.uk ** ** This file is part of MIToolbox, licensed under the 3-clause BSD license. *******************************************************************************/ #include "MIToolbox.h" #include "ArrayOperations.h" #include "CalculateProbability.h" #include "Entropy.h" #include "MutualInformation.h" /******************************************************************************* ** 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 a flag and a variable number of arguments ** depending on the value of the flag and returns either a construct ** containing probability estimates, a merged vector or a double value ** representing an entropy or mutual information *****************************************************************************/ int flag, i, numberOfSamples, checkSamples, thirdCheckSamples, numberOfFeatures, checkFeatures, thirdCheckFeatures; int numArities, errorTest; double *dataVector, *condVector, *targetVector, *firstVector, *secondVector, *output, *numStates; double *matrix, *mergedVector, *arities; int *outputIntVector, *intArities; double *jointOutput, *numJointStates, *firstOutput, *numFirstStates, *secondOutput, *numSecondStates; ProbabilityState state; JointProbabilityState jointState; /*if (nlhs != 1) { printf("Incorrect number of output arguments\n"); }//if not 1 output */ switch (nrhs) { case 2: { /*printf("Must be H(X), calculateProbability(X), merge(X), normaliseArray(X)\n");*/ break; } case 3: { /*printf("Must be H(XY), H(X|Y), calculateJointProbability(XY), I(X;Y)\n");*/ break; } case 4: { /*printf("Must be I(X;Y|Z)\n");*/ break; } default: { printf("Incorrect number of arguments, format is MIToolbox(\"FLAG\",varargin)\n"); break; } } /* number to function map ** 1 = calculateProbability ** 2 = calculateJointProbability ** 3 = mergeArrays ** 4 = H(X) ** 5 = H(XY) ** 6 = H(X|Y) ** 7 = I(X;Y) ** 8 = I(X;Y|Z) ** 9 = normaliseArray */ flag = *mxGetPr(prhs[0]); switch (flag) { case 1: { /* **calculateProbability */ numberOfSamples = mxGetM(prhs[1]); dataVector = (double *) mxGetPr(prhs[1]); /*ProbabilityState calculateProbability(double *dataVector, int vectorLength);*/ state = calculateProbability(dataVector,numberOfSamples); plhs[0] = mxCreateDoubleMatrix(state.numStates,1,mxREAL); plhs[1] = mxCreateDoubleMatrix(1,1,mxREAL); output = (double *)mxGetPr(plhs[0]); numStates = (double *) mxGetPr(plhs[1]); *numStates = state.numStates; for (i = 0; i < state.numStates; i++) { output[i] = state.probabilityVector[i]; } break; }/*case 1 - calculateProbability*/ case 2: { /* **calculateJointProbability */ numberOfSamples = mxGetM(prhs[1]); firstVector = (double *) mxGetPr(prhs[1]); secondVector = (double *) mxGetPr(prhs[2]); /*JointProbabilityState calculateJointProbability(double *firstVector, double *secondVector int vectorLength);*/ jointState = calculateJointProbability(firstVector,secondVector,numberOfSamples); plhs[0] = mxCreateDoubleMatrix(jointState.numJointStates,1,mxREAL); plhs[1] = mxCreateDoubleMatrix(1,1,mxREAL); plhs[2] = mxCreateDoubleMatrix(jointState.numFirstStates,1,mxREAL); plhs[3] = mxCreateDoubleMatrix(1,1,mxREAL); plhs[4] = mxCreateDoubleMatrix(jointState.numSecondStates,1,mxREAL); plhs[5] = mxCreateDoubleMatrix(1,1,mxREAL); jointOutput = (double *)mxGetPr(plhs[0]); numJointStates = (double *) mxGetPr(plhs[1]); firstOutput = (double *)mxGetPr(plhs[2]); numFirstStates = (double *) mxGetPr(plhs[3]); secondOutput = (double *)mxGetPr(plhs[4]); numSecondStates = (double *) mxGetPr(plhs[5]); *numJointStates = jointState.numJointStates; *numFirstStates = jointState.numFirstStates; *numSecondStates = jointState.numSecondStates; for (i = 0; i < jointState.numJointStates; i++) { jointOutput[i] = jointState.jointProbabilityVector[i]; } for (i = 0; i < jointState.numFirstStates; i++) { firstOutput[i] = jointState.firstProbabilityVector[i]; } for (i = 0; i < jointState.numSecondStates; i++) { secondOutput[i] = jointState.secondProbabilityVector[i]; } break; }/*case 2 - calculateJointProbability */ case 3: { /* **mergeArrays */ numberOfSamples = mxGetM(prhs[1]); numberOfFeatures = mxGetN(prhs[1]); numArities = 0; if (nrhs > 2) { numArities = mxGetN(prhs[2]); /*printf("arities = %d, features = %d, samples = %d\n",numArities,numberOfFeatures,numberOfSamples);*/ } plhs[0] = mxCreateDoubleMatrix(0,0,mxREAL); if (numArities == 0) { /* **no arities therefore compress output */ if ((numberOfFeatures > 0) && (numberOfSamples > 0)) { matrix = (double *) mxGetPr(prhs[1]); mergedVector = (double *) mxCalloc(numberOfSamples,sizeof(double)); plhs[0] = mxCreateDoubleMatrix(numberOfSamples,1,mxREAL); output = (double *)mxGetPr(plhs[0]); /*int mergeMultipleArrays(double *inputMatrix, double *outputVector, int matrixWidth, int vectorLength)*/ mergeMultipleArrays(matrix, mergedVector, numberOfFeatures, numberOfSamples); for (i = 0; i < numberOfSamples; i++) { output[i] = mergedVector[i]; } mxFree(mergedVector); mergedVector = NULL; } } else if (numArities == numberOfFeatures) { if ((numberOfFeatures > 0) && (numberOfSamples > 0)) { matrix = (double *) mxGetPr(prhs[1]); mergedVector = (double *) mxCalloc(numberOfSamples,sizeof(double)); arities = (double *) mxGetPr(prhs[2]); intArities = (int *) mxCalloc(numberOfFeatures,sizeof(int)); for (i = 0; i < numArities; i++) { intArities[i] = (int) floor(arities[i]); } /*int mergeMultipleArrays(double *inputMatrix, double *outputVector, int matrixWidth, int *arities, int vectorLength);*/ errorTest = mergeMultipleArraysArities(matrix, mergedVector, numberOfFeatures, intArities, numberOfSamples); if (errorTest != -1) { plhs[0] = mxCreateDoubleMatrix(numberOfSamples,1,mxREAL); output = (double *)mxGetPr(plhs[0]); for (i = 0; i < numberOfSamples; i++) { output[i] = mergedVector[i]; } } else { printf("Incorrect arities supplied. More states in data than specified\n"); } mxFree(mergedVector); mergedVector = NULL; } } else { printf("Number of arities does not match number of features, arities should be a row vector\n"); } break; }/*case 3 - mergeArrays*/ case 4: { /* **H(X) */ numberOfSamples = mxGetM(prhs[1]); numberOfFeatures = mxGetN(prhs[1]); dataVector = (double *) mxGetPr(prhs[1]); plhs[0] = mxCreateDoubleMatrix(1,1,mxREAL); output = (double *)mxGetPr(plhs[0]); if (numberOfFeatures == 1) { /*double calculateEntropy(double *dataVector, int vectorLength);*/ *output = calculateEntropy(dataVector,numberOfSamples); } else { printf("No columns in input\n"); *output = -1.0; } break; }/*case 4 - H(X)*/ case 5: { /* **H(XY) */ numberOfSamples = mxGetM(prhs[1]); checkSamples = mxGetM(prhs[2]); numberOfFeatures = mxGetN(prhs[1]); checkFeatures = mxGetN(prhs[2]); firstVector = mxGetPr(prhs[1]); secondVector = mxGetPr(prhs[2]); plhs[0] = mxCreateDoubleMatrix(1,1,mxREAL); output = (double *)mxGetPr(plhs[0]); if ((numberOfFeatures == 1) && (checkFeatures == 1)) { if ((numberOfSamples == 0) && (checkSamples == 0)) { *output = 0.0; } else if (numberOfSamples == 0) { *output = calculateEntropy(secondVector,numberOfSamples); } else if (checkSamples == 0) { *output = calculateEntropy(firstVector,numberOfSamples); } else if (numberOfSamples == checkSamples) { /*double calculateJointEntropy(double *firstVector, double *secondVector, int vectorLength);*/ *output = calculateJointEntropy(firstVector,secondVector,numberOfSamples); } else { printf("Vector lengths do not match, they must be the same length\n"); *output = -1.0; } } else { printf("No columns in input\n"); *output = -1.0; } break; }/*case 5 - H(XY)*/ case 6: { /* **H(X|Y) */ numberOfSamples = mxGetM(prhs[1]); checkSamples = mxGetM(prhs[2]); numberOfFeatures = mxGetN(prhs[1]); checkFeatures = mxGetN(prhs[2]); dataVector = mxGetPr(prhs[1]); condVector = mxGetPr(prhs[2]); plhs[0] = mxCreateDoubleMatrix(1,1,mxREAL); output = (double *)mxGetPr(plhs[0]); if ((numberOfFeatures == 1) && (checkFeatures == 1)) { if (numberOfSamples == 0) { *output = 0.0; } else if (checkSamples == 0) { *output = calculateEntropy(dataVector,numberOfSamples); } else if (numberOfSamples == checkSamples) { /*double calculateConditionalEntropy(double *dataVector, double *condVector, int vectorLength);*/ *output = calculateConditionalEntropy(dataVector,condVector,numberOfSamples); } else { printf("Vector lengths do not match, they must be the same length\n"); *output = -1.0; } } else { printf("No columns in input\n"); *output = -1.0; } break; }/*case 6 - H(X|Y)*/ case 7: { /* **I(X;Y) */ numberOfSamples = mxGetM(prhs[1]); checkSamples = mxGetM(prhs[2]); numberOfFeatures = mxGetN(prhs[1]); checkFeatures = mxGetN(prhs[2]); firstVector = mxGetPr(prhs[1]); secondVector = mxGetPr(prhs[2]); plhs[0] = mxCreateDoubleMatrix(1,1,mxREAL); output = (double *)mxGetPr(plhs[0]); if ((numberOfFeatures == 1) && (checkFeatures == 1)) { if ((numberOfSamples == 0) || (checkSamples == 0)) { *output = 0.0; } else if (numberOfSamples == checkSamples) { /*double calculateMutualInformation(double *firstVector, double *secondVector, int vectorLength);*/ *output = calculateMutualInformation(firstVector,secondVector,numberOfSamples); } else { printf("Vector lengths do not match, they must be the same length\n"); *output = -1.0; } } else { printf("No columns in input\n"); *output = -1.0; } break; }/*case 7 - I(X;Y)*/ case 8: { /* **I(X;Y|Z) */ numberOfSamples = mxGetM(prhs[1]); checkSamples = mxGetM(prhs[2]); thirdCheckSamples = mxGetM(prhs[3]); numberOfFeatures = mxGetN(prhs[1]); checkFeatures = mxGetN(prhs[2]); thirdCheckFeatures = mxGetN(prhs[3]); firstVector = mxGetPr(prhs[1]); targetVector = mxGetPr(prhs[2]); condVector = mxGetPr(prhs[3]); plhs[0] = mxCreateDoubleMatrix(1,1,mxREAL); output = (double *)mxGetPr(plhs[0]); if ((numberOfFeatures == 1) && (checkFeatures == 1)) { if ((numberOfSamples == 0) || (checkSamples == 0)) { *output = 0.0; } else if ((thirdCheckSamples == 0) || (thirdCheckFeatures != 1)) { *output = calculateMutualInformation(firstVector,targetVector,numberOfSamples); } else if ((numberOfSamples == checkSamples) && (numberOfSamples == thirdCheckSamples)) { /*double calculateConditionalMutualInformation(double *firstVector, double *targetVector, double *condVector, int vectorLength);*/ *output = calculateConditionalMutualInformation(firstVector,targetVector,condVector,numberOfSamples); } else { printf("Vector lengths do not match, they must be the same length\n"); *output = -1.0; } } else { printf("No columns in input\n"); *output = -1.0; } break; }/*case 8 - I(X;Y|Z)*/ case 9: { /* **normaliseArray */ numberOfSamples = mxGetM(prhs[1]); dataVector = (double *) mxGetPr(prhs[1]); outputIntVector = (int *) mxCalloc(numberOfSamples,sizeof(int)); plhs[0] = mxCreateDoubleMatrix(numberOfSamples,1,mxREAL); plhs[1] = mxCreateDoubleMatrix(1,1,mxREAL); output = (double *)mxGetPr(plhs[0]); numStates = (double *) mxGetPr(plhs[1]); /*int normaliseArray(double *inputVector, int *outputVector, int vectorLength);*/ *numStates = normaliseArray(dataVector, outputIntVector, numberOfSamples); for (i = 0; i < numberOfSamples; i++) { output[i] = outputIntVector[i]; } break; }/*case 9 - normaliseArray*/ default: { printf("Unrecognised flag\n"); break; }/*default*/ }/*switch(flag)*/ return; }/*mexFunction()*/ MIToolbox-2.1.2/Makefile000066400000000000000000000077361264454417700150630ustar00rootroot00000000000000#makefile for MIToolbox #Author: Adam Pocock, apocock@cs.man.ac.uk #Created 11/3/2010 #Updated 22/2/2014 - Added make install # # #Copyright 2010 Adam Pocock, The University Of Manchester #www.cs.manchester.ac.uk # #This file is part of MIToolbox. # #MIToolbox is free software: you can redistribute it and/or modify #it under the terms of the GNU Lesser General Public License as published by #the Free Software Foundation, either version 3 of the License, or #(at your option) any later version. # #MIToolbox is distributed in the hope that it will be useful, #but WITHOUT ANY WARRANTY; without even the implied warranty of #MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #GNU Lesser General Public License for more details. # #You should have received a copy of the GNU Lesser General Public License #along with MIToolbox. If not, see . PREFIX = /usr/local CXXFLAGS = -O3 -fPIC COMPILER = gcc objects = ArrayOperations.o CalculateProbability.o Entropy.o \ MutualInformation.o RenyiEntropy.o RenyiMutualInformation.o \ WeightedEntropy.o WeightedMutualInformation.o libMIToolbox.so : $(objects) $(COMPILER) $(CXXFLAGS) -shared -o libMIToolbox.so $(objects) -lm WeightedMutualInformation.o: WeightedMutualInformation.c MIToolbox.h ArrayOperations.h \ CalculateProbability.h WeightedEntropy.h $(COMPILER) $(CXXFLAGS) -DCOMPILE_C -c WeightedMutualInformation.c WeightedEntropy.o: WeightedEntropy.c MIToolbox.h ArrayOperations.h \ CalculateProbability.h $(COMPILER) $(CXXFLAGS) -DCOMPILE_C -c WeightedEntropy.c RenyiMutualInformation.o: RenyiMutualInformation.c MIToolbox.h ArrayOperations.h \ CalculateProbability.h RenyiEntropy.h $(COMPILER) $(CXXFLAGS) -DCOMPILE_C -c RenyiMutualInformation.c RenyiEntropy.o: RenyiEntropy.c MIToolbox.h ArrayOperations.h \ CalculateProbability.h $(COMPILER) $(CXXFLAGS) -DCOMPILE_C -c RenyiEntropy.c MutualInformation.o: MutualInformation.c MIToolbox.h ArrayOperations.h \ CalculateProbability.h Entropy.h MutualInformation.h $(COMPILER) $(CXXFLAGS) -DCOMPILE_C -c MutualInformation.c Entropy.o: Entropy.c MIToolbox.h ArrayOperations.h CalculateProbability.h \ Entropy.h $(COMPILER) $(CXXFLAGS) -DCOMPILE_C -c Entropy.c CalculateProbability.o: CalculateProbability.c MIToolbox.h ArrayOperations.h \ CalculateProbability.h $(COMPILER) $(CXXFLAGS) -DCOMPILE_C -c CalculateProbability.c ArrayOperations.o: ArrayOperations.c MIToolbox.h ArrayOperations.h $(COMPILER) $(CXXFLAGS) -DCOMPILE_C -c ArrayOperations.c .PHONY : debug debug: $(MAKE) libMIToolbox.so "CXXFLAGS = -g -DDEBUG -fPIC" .PHONY : x86 x86: $(MAKE) libMIToolbox.so "CXXFLAGS = -O3 -fPIC -m32" .PHONY : x64 x64: $(MAKE) libMIToolbox.so "CXXFLAGS = -O3 -fPIC -m64" .PHONY : matlab matlab: mex MIToolboxMex.c MutualInformation.c Entropy.c CalculateProbability.c ArrayOperations.c mex RenyiMIToolboxMex.c RenyiMutualInformation.c RenyiEntropy.c CalculateProbability.c ArrayOperations.c mex WeightedMIToolboxMex.c WeightedMutualInformation.c WeightedEntropy.c CalculateProbability.c ArrayOperations.c .PHONY : matlab-debug matlab-debug: mex -g MIToolboxMex.c MutualInformation.c Entropy.c CalculateProbability.c ArrayOperations.c mex -g RenyiMIToolboxMex.c RenyiMutualInformation.c RenyiEntropy.c CalculateProbability.c ArrayOperations.c mex -g WeightedMIToolboxMex.c WeightedMutualInformation.c WeightedEntropy.c CalculateProbability.c ArrayOperations.c .PHONY : intel intel: $(MAKE) libMIToolbox.so "COMPILER = icc" "CXXFLAGS = -O2 -fPIC -xHost" .PHONY : clean clean: rm *.o rm libMIToolbox.so .PHONY : install install: $(MAKE) @echo "Installing libMIToolbox.so to $(PREFIX)/lib" @cp -v libMIToolbox.so $(PREFIX)/lib @mkdir -p $(PREFIX)/include/MIToolbox @echo "Installing MIToolbox's header files to $(PREFIX)/include/MIToolbox" @cp -v MIToolbox.h MutualInformation.h Entropy.h CalculateProbability.h ArrayOperations.h RenyiMutualInformation.h RenyiEntropy.h WeightedMutualInformation.h WeightedEntropy.h $(PREFIX)/include/MIToolbox/ MIToolbox-2.1.2/MutualInformation.c000066400000000000000000000063221264454417700172320ustar00rootroot00000000000000/******************************************************************************* ** MutualInformation.c ** Part of the mutual information toolbox ** ** Contains functions to calculate the mutual information of ** two variables X and Y, I(X;Y), to calculate the joint mutual information ** of two variables X & Z on the variable Y, I(XZ;Y), and the conditional ** mutual information I(x;Y|Z) ** ** Author: Adam Pocock ** Created 19/2/2010 ** Updated - 22/02/2014 - Added checking on calloc. ** ** Copyright 2010-2014 Adam Pocock, The University Of Manchester ** www.cs.manchester.ac.uk ** ** This file is part of MIToolbox, licensed under the 3-clause BSD license. *******************************************************************************/ #include "MIToolbox.h" #include "ArrayOperations.h" #include "CalculateProbability.h" #include "Entropy.h" #include "MutualInformation.h" double calculateMutualInformation(double *dataVector, double *targetVector, int vectorLength) { double mutualInformation = 0.0; int firstIndex,secondIndex; int i; JointProbabilityState state = calculateJointProbability(dataVector,targetVector,vectorLength); /* ** I(X;Y) = sum sum p(xy) * log (p(xy)/p(x)p(y)) */ for (i = 0; i < state.numJointStates; i++) { firstIndex = i % state.numFirstStates; secondIndex = i / state.numFirstStates; if ((state.jointProbabilityVector[i] > 0) && (state.firstProbabilityVector[firstIndex] > 0) && (state.secondProbabilityVector[secondIndex] > 0)) { /*double division is probably more stable than multiplying two small numbers together ** mutualInformation += state.jointProbabilityVector[i] * log(state.jointProbabilityVector[i] / (state.firstProbabilityVector[firstIndex] * state.secondProbabilityVector[secondIndex])); */ mutualInformation += state.jointProbabilityVector[i] * log(state.jointProbabilityVector[i] / state.firstProbabilityVector[firstIndex] / state.secondProbabilityVector[secondIndex]); } } mutualInformation /= log(2.0); FREE_FUNC(state.firstProbabilityVector); state.firstProbabilityVector = NULL; FREE_FUNC(state.secondProbabilityVector); state.secondProbabilityVector = NULL; FREE_FUNC(state.jointProbabilityVector); state.jointProbabilityVector = NULL; return mutualInformation; }/*calculateMutualInformation(double *,double *,int)*/ double calculateConditionalMutualInformation(double *dataVector, double *targetVector, double *conditionVector, int vectorLength) { double mutualInformation = 0.0; double firstCondition, secondCondition; double *mergedVector = (double *) checkedCalloc(vectorLength,sizeof(double)); mergeArrays(targetVector,conditionVector,mergedVector,vectorLength); /* I(X;Y|Z) = H(X|Z) - H(X|YZ) */ /* double calculateConditionalEntropy(double *dataVector, double *conditionVector, int vectorLength); */ firstCondition = calculateConditionalEntropy(dataVector,conditionVector,vectorLength); secondCondition = calculateConditionalEntropy(dataVector,mergedVector,vectorLength); mutualInformation = firstCondition - secondCondition; FREE_FUNC(mergedVector); mergedVector = NULL; return mutualInformation; }/*calculateConditionalMutualInformation(double *,double *,double *,int)*/ MIToolbox-2.1.2/MutualInformation.h000066400000000000000000000036131264454417700172370ustar00rootroot00000000000000/******************************************************************************* ** MutualInformation.h ** Part of the mutual information toolbox ** ** Contains functions to calculate the mutual information of ** two variables X and Y, I(X;Y), to calculate the joint mutual information ** of two variables X & Z on the variable Y, I(XZ;Y), and the conditional ** mutual information I(x;Y|Z) ** ** Author: Adam Pocock ** Created 19/2/2010 ** ** Copyright 2010 Adam Pocock, The University Of Manchester ** www.cs.manchester.ac.uk ** ** This file is part of MIToolbox, licensed under the 3-clause BSD license. *******************************************************************************/ #ifndef __MutualInformation_H #define __MutualInformation_H #ifdef __cplusplus extern "C" { #endif /******************************************************************************* ** calculateMutualInformation returns the log base 2 mutual information between ** dataVector and targetVector, I(X;Y) ** ** length(dataVector) == length(targetVector) == vectorLength otherwise there ** will be a segmentation fault *******************************************************************************/ double calculateMutualInformation(double *dataVector, double *targetVector, int vectorLength); /******************************************************************************* ** calculateConditionalMutualInformation returns the log base 2 ** mutual information between dataVector and targetVector, conditioned on ** conditionVector, I(X;Y|Z) ** ** length(dataVector) == length(targetVector) == length(condtionVector) == vectorLength ** otherwise it will error with a segmentation fault *******************************************************************************/ double calculateConditionalMutualInformation(double *dataVector, double *targetVector, double *conditionVector, int vectorLength); #ifdef __cplusplus } #endif #endif MIToolbox-2.1.2/README.md000066400000000000000000000062721264454417700146740ustar00rootroot00000000000000MIToolbox ========= v2.1.1 for C/C++ and MATLAB/Octave MIToolbox contains a set of functions to calculate information theoretic quantities from data, such as the entropy and mutual information. The toolbox contains implementations of the most popular Shannon entropies, and also the lesser known Renyi entropy. The toolbox also provides implementations of the weighted entropy and weighted mutual information from "Information Theory with Application", S. Guiasu (1977). The toolbox only supports discrete distributions, as opposed to continuous. All real-valued numbers will be processed by x = floor(x). These functions are targeted for use with feature selection algorithms rather than communication channels and so expect all the data to be available before execution and sample their own probability distributions from the data. All functions expect the inputs to be vectors or matrices of doubles. Functions contained: - Entropy - Conditional Entropy - Mutual Information - Conditional Mutual Information - generating a joint variable - generating a probability distribution from a discrete random variable - Renyi's Entropy - Renyi's Mutual Information - Weighted Entropy - Weighted Mutual Information - Weighted Conditional Mutual Information Note: all functions are calculated in log base 2, so return units of "bits". ====== Examples: ``` >> y = [1 1 1 0 0]'; >> x = [1 0 1 1 0]'; ``` ``` >> mi(x,y) %% mutual information I(X;Y) ans = 0.0200 ``` ``` >> h(x) %% entropy H(X) ans = 0.9710 ``` ``` >> condh(x,y) %% conditional entropy H(X|Y) ans = 0.9510 ``` ``` >> h( [x,y] ) %% joint entropy H(X,Y) ans = 1.9219 ``` ``` >> joint([x,y]) %% joint random variable XY ans = 1 2 1 3 4 ``` ====== To compile the library for use in MATLAB/OCTAVE, execute CompileMIToolbox.m from within MATLAB, or run 'make matlab' from a terminal. To compile the library for use with C programs run 'make x86' for a 32-bit library, or 'make x64' for a 64-bit library. Then run 'sudo make install' to install MIToolbox into /usr/local/lib & /usr/local/include. All code is licensed under the 3-clause BSD license. Update History - 10/01/2016 - v2.1.2 - Relicense from LGPL to BSD. Added checks to ensure input MATLAB types are doubles. - 02/02/2015 - v2.1.1 - Fixed up the Makefile so it installs the headers too. - 22/02/2014 - v2.1 - Fixed a couple of bugs related to memory handling. Added a make install for compatibility with PyFeast. - 30/07/2011 - v2.00 - Added implementations of the weighted entropy and weighted mutual information. More cleanup of Mex entry point to further check the inputs. - 08/11/2011 - v1.03 - Minor documentation changes to accompany the JMLR publication. - 15/10/2010 - v1.02 - Fixed bug where MIToolbox would cause a segmentation fault if a x by 0 empty matrix was passed in. Now prints an error message and returns gracefully. - 02/09/2010 - v1.01 - Fixed a bug in CMIM.m where the last feature would not be selected first if it had the highest MI. - 07/07/2010 - v1.00 - Initial Release. MIToolbox-2.1.2/RenyiEntropy.c000066400000000000000000000123171264454417700162250ustar00rootroot00000000000000/******************************************************************************* ** RenyiEntropy.c ** Part of the mutual information toolbox ** ** Contains functions to calculate the Renyi alpha entropy of a single variable ** H_\alpha(X), the Renyi joint entropy of two variables H_\alpha(X,Y), and the ** conditional Renyi entropy H_\alpha(X|Y) ** ** Author: Adam Pocock ** Created 26/3/2010 ** Updated - 22/02/2014 - Added checking on calloc. ** ** Copyright 2010-2014 Adam Pocock, The University Of Manchester ** www.cs.manchester.ac.uk ** ** This file is part of MIToolbox, licensed under the 3-clause BSD license. *******************************************************************************/ #include "MIToolbox.h" #include "ArrayOperations.h" #include "CalculateProbability.h" #include "Entropy.h" double calculateRenyiEntropy(double alpha, double *dataVector, int vectorLength) { double entropy = 0.0; double tempValue = 0.0; int i; ProbabilityState state = calculateProbability(dataVector,vectorLength); /*H_\alpha(X) = 1/(1-alpha) * log(2)(sum p(x)^alpha)*/ for (i = 0; i < state.numStates; i++) { tempValue = state.probabilityVector[i]; if (tempValue > 0) { entropy += pow(tempValue,alpha); /*printf("Entropy = %f, i = %d\n", entropy,i);*/ } } /*printf("Entropy = %f\n", entropy);*/ entropy = log(entropy); entropy /= log(2.0); entropy /= (1.0-alpha); /*printf("Entropy = %f\n", entropy);*/ FREE_FUNC(state.probabilityVector); state.probabilityVector = NULL; return entropy; }/*calculateRenyiEntropy(double,double*,int)*/ double calculateJointRenyiEntropy(double alpha, double *firstVector, double *secondVector, int vectorLength) { double jointEntropy = 0.0; double tempValue = 0.0; int i; JointProbabilityState state = calculateJointProbability(firstVector,secondVector,vectorLength); /*H_\alpha(XY) = 1/(1-alpha) * log(2)(sum p(xy)^alpha)*/ for (i = 0; i < state.numJointStates; i++) { tempValue = state.jointProbabilityVector[i]; if (tempValue > 0) { jointEntropy += pow(tempValue,alpha); } } jointEntropy = log(jointEntropy); jointEntropy /= log(2.0); jointEntropy /= (1.0-alpha); FREE_FUNC(state.firstProbabilityVector); state.firstProbabilityVector = NULL; FREE_FUNC(state.secondProbabilityVector); state.secondProbabilityVector = NULL; FREE_FUNC(state.jointProbabilityVector); state.jointProbabilityVector = NULL; return jointEntropy; }/*calculateJointRenyiEntropy(double,double*,double*,int)*/ double calcCondRenyiEnt(double alpha, double *dataVector, double *conditionVector, int uniqueInCondVector, int vectorLength) { /*uniqueInCondVector = is the number of unique values in the cond vector.*/ /*condEntropy = sum p(y) * sum p(x|y)^alpha(*/ /* ** first generate the seperate variables */ double *seperateVectors = (double *) checkedCalloc(uniqueInCondVector*vectorLength,sizeof(double)); int *seperateVectorCount = (int *) checkedCalloc(uniqueInCondVector,sizeof(int)); double seperateVectorProb = 0.0; int i,j; double entropy = 0.0; double tempValue = 0.0; int currentValue; double tempEntropy; ProbabilityState state; double **seperateVectors2D = (double **) checkedCalloc(uniqueInCondVector,sizeof(double*)); for(j=0; j < uniqueInCondVector; j++) seperateVectors2D[j] = seperateVectors + (int)j*vectorLength; for (i = 0; i < vectorLength; i++) { currentValue = (int) (conditionVector[i] - 1.0); /*printf("CurrentValue = %d\n",currentValue);*/ seperateVectors2D[currentValue][seperateVectorCount[currentValue]] = dataVector[i]; seperateVectorCount[currentValue]++; } for (j = 0; j < uniqueInCondVector; j++) { tempEntropy = 0.0; seperateVectorProb = ((double)seperateVectorCount[j]) / vectorLength; state = calculateProbability(seperateVectors2D[j],seperateVectorCount[j]); /*H_\alpha(X) = 1/(1-alpha) * log(2)(sum p(x)^alpha)*/ for (i = 0; i < state.numStates; i++) { tempValue = state.probabilityVector[i]; if (tempValue > 0) { tempEntropy += pow(tempValue,alpha); /*printf("Entropy = %f, i = %d\n", entropy,i);*/ } } /*printf("Entropy = %f\n", entropy);*/ tempEntropy = log(tempEntropy); tempEntropy /= log(2.0); tempEntropy /= (1.0-alpha); entropy += tempEntropy; FREE_FUNC(state.probabilityVector); } FREE_FUNC(seperateVectors2D); seperateVectors2D = NULL; FREE_FUNC(seperateVectors); FREE_FUNC(seperateVectorCount); seperateVectors = NULL; seperateVectorCount = NULL; return entropy; }/*calcCondRenyiEnt(double *,double *,int)*/ double calculateConditionalRenyiEntropy(double alpha, double *dataVector, double *conditionVector, int vectorLength) { /*calls this: **double calculateConditionalRenyiEntropy(double alpha, double *firstVector, double *condVector, int uniqueInCondVector, int vectorLength) **after determining uniqueInCondVector */ int numUnique = numberOfUniqueValues(conditionVector, vectorLength); return calcCondRenyiEnt(alpha, dataVector, conditionVector, numUnique, vectorLength); }/*calculateConditionalRenyiEntropy(double,double*,double*,int)*/ MIToolbox-2.1.2/RenyiEntropy.h000066400000000000000000000040141264454417700162250ustar00rootroot00000000000000/******************************************************************************* ** RenyiEntropy.h ** Part of the mutual information toolbox ** ** Contains functions to calculate the Renyi alpha entropy of a single variable ** H_\alpha(X), the Renyi joint entropy of two variables H_\alpha(X,Y), and the ** conditional Renyi entropy H_\alpha(X|Y) ** ** Author: Adam Pocock ** Created 26/3/2010 ** ** Copyright 2010 Adam Pocock, The University Of Manchester ** www.cs.manchester.ac.uk ** ** This file is part of MIToolbox. ** ** This file is part of MIToolbox, licensed under the 3-clause BSD license. *******************************************************************************/ #ifndef __Renyi_Entropy_H #define __Renyi_Entropy_H #ifdef __cplusplus extern "C" { #endif /******************************************************************************* ** calculateRenyiEntropy returns the Renyi entropy in log base 2 of dataVector ** H_{\alpha}(X), for \alpha != 1 ** ** length(dataVector) == vectorLength otherwise there ** will be a segmentation fault *******************************************************************************/ double calculateRenyiEntropy(double alpha, double *dataVector, int vectorLength); /******************************************************************************* ** calculateJointRenyiEntropy returns the Renyi entropy in log base 2 of the ** joint variable of firstVector and secondVector H_{\alpha}(XY), ** for \alpha != 1 ** ** length(firstVector) == length(secondVector) == vectorLength otherwise there ** will be a segmentation fault *******************************************************************************/ double calculateJointRenyiEntropy(double alpha, double *firstVector, double *secondVector, int vectorLength); /* This function does not return a valid conditonal entropy as it has no ** meaning in Renyi's extension of entropy double calculateConditionalRenyiEntropy(double alpha, double *dataVector, double *conditionVector, int vectorLength); */ #ifdef __cplusplus } #endif #endif MIToolbox-2.1.2/RenyiMIToolbox.m000066400000000000000000000031371264454417700164530ustar00rootroot00000000000000function [varargout] = RenyiMIToolbox(functionName, alpha, varargin) %function [varargout] = RenyiMIToolbox(functionName, alpha, varargin) % %Provides access to the functions in RenyiMIToolboxMex % %Expects column vectors, will not work with row vectors % %Function list %"Entropy" = H_{\alpha}(X) = 1 %"MI" = I_{\alpha}(X;Y) = 3 % %Arguments and returned values %[entropy] = H_\alpha(X) = H(alpha,vector) %[mi] = I_\alpha(X;Y) = I(alpha,vector,target) % %Internal RenyiMIToolbox function number %Renyi Entropy = 1; %Renyi MI = 3; for i = 1:length(varargin) if (~isa(varargin{i},'double')) error('Error, MIToolbox requires inputs to be double vector or matrices') end end if (alpha ~= 1) if (strcmpi(functionName,'Entropy') || strcmpi(functionName,'h')) %disp('Calculating Entropy'); if (size(varargin{1},2)>1) mergedVector = MIToolboxMex(3,varargin{1}); else mergedVector = varargin{1}; end [varargout{1}] = RenyiMIToolboxMex(1,alpha,mergedVector); elseif (strcmpi(functionName,'MI')) if (size(varargin{1},2)>1) mergedFirst = MIToolboxMex(3,varargin{1}); else mergedFirst = varargin{1}; end if (size(varargin{2},2)>1) mergedSecond = MIToolboxMex(3,varargin{2}); else mergedSecond = varargin{2}; end [varargout{1}] = RenyiMIToolboxMex(3,alpha,mergedFirst,mergedSecond); else varargout{1} = 0; disp(['Unrecognised functionName ' functionName]); end else disp('For alpha = 1 use functions in MIToolbox.m'); disp('as those functions are the implementation of Shannon''s Information Theory'); end MIToolbox-2.1.2/RenyiMIToolboxMex.c000066400000000000000000000122651264454417700171150ustar00rootroot00000000000000/******************************************************************************* ** RenyiMIToolboxMex.c ** is the MATLAB entry point for the Renyi Entropy and MI MIToolbox functions ** when called from a MATLAB/OCTAVE script. ** ** Copyright 2010 Adam Pocock, The University Of Manchester ** www.cs.manchester.ac.uk ** ** This file is part of MIToolbox, licensed under the 3-clause BSD license. *******************************************************************************/ #include "MIToolbox.h" #include "ArrayOperations.h" #include "RenyiEntropy.h" #include "RenyiMutualInformation.h" /******************************************************************************* **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 a flag and 2 or 3 other arguments ** the first is a scalar alpha value, and the remainder are ** arrays. It returns a Renyi entropy or mutual information using the ** alpha divergence. *****************************************************************************/ int flag, numberOfSamples, checkSamples, numberOfFeatures, checkFeatures; double alpha; double *dataVector, *firstVector, *secondVector, *output; /*if (nlhs != 1) { printf("Incorrect number of output arguments\n"); }//if not 1 output */ switch (nrhs) { case 3: { /*printf("Must be H_\alpha(X)\n");*/ break; } case 4: { /*printf("Must be H_\alpha(XY), I_\alpha(X;Y)\n");*/ break; } default: { printf("Incorrect number of arguments, format is RenyiMIToolbox(\"FLAG\",varargin)\n"); break; } } /* number to function map ** 1 = H(X) ** 2 = H(XY) ** 3 = I(X;Y) */ flag = *mxGetPr(prhs[0]); switch (flag) { case 1: { /* **H_{\alpha}(X) */ alpha = mxGetScalar(prhs[1]); numberOfSamples = mxGetM(prhs[2]); numberOfFeatures = mxGetN(prhs[2]); dataVector = (double *) mxGetPr(prhs[2]); plhs[0] = mxCreateDoubleMatrix(1,1,mxREAL); output = (double *) mxGetPr(plhs[0]); if (numberOfFeatures == 1) { /*double calculateRenyiEntropy(double alpha, double *dataVector, long vectorLength);*/ *output = calculateRenyiEntropy(alpha,dataVector,numberOfSamples); } else { printf("No columns in input\n"); *output = -1.0; } break; }/*case 1 - H_{\alpha}(X)*/ case 2: { /* **H_{\alpha}(XY) */ alpha = mxGetScalar(prhs[1]); numberOfSamples = mxGetM(prhs[2]); checkSamples = mxGetM(prhs[3]); numberOfFeatures = mxGetN(prhs[2]); checkFeatures = mxGetN(prhs[3]); firstVector = mxGetPr(prhs[2]); secondVector = mxGetPr(prhs[3]); plhs[0] = mxCreateDoubleMatrix(1,1,mxREAL); output = (double *)mxGetPr(plhs[0]); if ((numberOfFeatures == 1) && (checkFeatures == 1)) { if ((numberOfSamples == 0) || (checkSamples == 0)) { *output = 0.0; } else if (numberOfSamples == checkSamples) { /*double calculateJointRenyiEntropy(double alpha, double *firstVector, double *secondVector, long vectorLength);*/ *output = calculateJointRenyiEntropy(alpha,firstVector,secondVector,numberOfSamples); } else { printf("Vector lengths do not match, they must be the same length"); *output = -1.0; } } else { printf("No columns in input\n"); *output = -1.0; } break; }/*case 2 - H_{\alpha}(XY)*/ case 3: { /* **I_{\alpha}(X;Y) */ alpha = mxGetScalar(prhs[1]); numberOfSamples = mxGetM(prhs[2]); checkSamples = mxGetM(prhs[3]); numberOfFeatures = mxGetN(prhs[2]); checkFeatures = mxGetN(prhs[3]); firstVector = mxGetPr(prhs[2]); secondVector = mxGetPr(prhs[3]); plhs[0] = mxCreateDoubleMatrix(1,1,mxREAL); output = (double *)mxGetPr(plhs[0]); if ((numberOfFeatures == 1) && (checkFeatures == 1)) { if ((numberOfSamples == 0) || (checkSamples == 0)) { *output = 0.0; } else if (numberOfSamples == checkSamples) { /*double calculateRenyiMIDivergence(double alpha, double *dataVector, double *targetVector, long vectorLength);*/ *output = calculateRenyiMIDivergence(alpha,firstVector,secondVector,numberOfSamples); } else { printf("Vector lengths do not match, they must be the same length"); *output = -1.0; } } else { printf("No columns in input\n"); *output = -1.0; } break; }/*case 3 - I_{\alpha}(X;Y)*/ default: { printf("Unrecognised flag\n"); break; }/*default*/ }/*switch(flag)*/ return; }/*mexFunction()*/ MIToolbox-2.1.2/RenyiMutualInformation.c000066400000000000000000000056351264454417700202470ustar00rootroot00000000000000/******************************************************************************* ** RenyiMutualInformation.c ** Part of the mutual information toolbox ** ** Contains functions to calculate the Renyi mutual information of ** two variables X and Y, I_\alpha(X;Y), using the Renyi alpha divergence and ** the joint entropy difference ** ** Author: Adam Pocock ** Created 26/3/2010 ** ** Copyright 2010 Adam Pocock, The University Of Manchester ** www.cs.manchester.ac.uk ** ** This file is part of MIToolbox, licensed under the 3-clause BSD license. *******************************************************************************/ #include "MIToolbox.h" #include "ArrayOperations.h" #include "CalculateProbability.h" #include "RenyiEntropy.h" #include "RenyiMutualInformation.h" double calculateRenyiMIDivergence(double alpha, double *dataVector, double *targetVector, int vectorLength) { double mutualInformation = 0.0; int firstIndex,secondIndex; int i; double jointTemp = 0.0; double seperateTemp = 0.0; double invAlpha = 1.0 - alpha; JointProbabilityState state = calculateJointProbability(dataVector,targetVector,vectorLength); /* standard MI is D_KL(p(x,y)||p(x)p(y)) ** which expands to ** D_KL(p(x,y)||p(x)p(y)) = sum(p(x,y) * log(p(x,y)/(p(x)p(y)))) ** ** Renyi alpha divergence D_alpha(p(x,y)||p(x)p(y)) ** expands to ** D_alpha(p(x,y)||p(x)p(y)) = 1/(alpha-1) * log(sum((p(x,y)^alpha)*((p(x)p(y))^(1-alpha)))) */ for (i = 0; i < state.numJointStates; i++) { firstIndex = i % state.numFirstStates; secondIndex = i / state.numFirstStates; if ((state.jointProbabilityVector[i] > 0) && (state.firstProbabilityVector[firstIndex] > 0) && (state.secondProbabilityVector[secondIndex] > 0)) { jointTemp = pow(state.jointProbabilityVector[i],alpha); seperateTemp = state.firstProbabilityVector[firstIndex] * state.secondProbabilityVector[secondIndex]; seperateTemp = pow(seperateTemp,invAlpha); mutualInformation += (jointTemp * seperateTemp); } } mutualInformation = log(mutualInformation); mutualInformation /= log(2.0); mutualInformation /= (alpha-1.0); FREE_FUNC(state.firstProbabilityVector); state.firstProbabilityVector = NULL; FREE_FUNC(state.secondProbabilityVector); state.secondProbabilityVector = NULL; FREE_FUNC(state.jointProbabilityVector); state.jointProbabilityVector = NULL; return mutualInformation; }/*calculateRenyiMIDivergence(double, double *, double *, int)*/ double calculateRenyiMIJoint(double alpha, double *dataVector, double *targetVector, int vectorLength) { double hY = calculateRenyiEntropy(alpha, targetVector, vectorLength); double hX = calculateRenyiEntropy(alpha, dataVector, vectorLength); double hXY = calculateJointRenyiEntropy(alpha, dataVector, targetVector, vectorLength); double answer = hX + hY - hXY; return answer; }/*calculateRenyiMIJoint(double, double*, double*, int)*/ MIToolbox-2.1.2/RenyiMutualInformation.h000066400000000000000000000035151264454417700202470ustar00rootroot00000000000000/******************************************************************************* ** RenyiMutualInformation.h ** Part of the mutual information toolbox ** ** Contains functions to calculate the Renyi mutual information of ** two variables X and Y, I_\alpha(X;Y), using the Renyi alpha divergence and ** the joint entropy difference ** ** Author: Adam Pocock ** Created 26/3/2010 ** ** Copyright 2010 Adam Pocock, The University Of Manchester ** www.cs.manchester.ac.uk ** ** This file is part of MIToolbox, licensed under the 3-clause BSD license. *******************************************************************************/ #ifndef __Renyi_MutualInformation_H #define __Renyi_MutualInformation_H #ifdef __cplusplus extern "C" { #endif /******************************************************************************* ** calculateRenyiMIDivergence returns the log base 2 Renyi mutual information ** between dataVector and targetVector, I_{\alpha}(X;Y), for \alpha != 1 ** This uses Renyi's generalised alpha divergence as the difference measure ** instead of the KL-divergence as in Shannon's Mutual Information ** ** length(dataVector) == length(targetVector) == vectorLength otherwise there ** will be a segmentation fault *******************************************************************************/ double calculateRenyiMIDivergence(double alpha, double *dataVector, double *targetVector, int vectorLength); /****************************************************************************** ** This function returns a different value to the alpha divergence mutual ** information, and thus is not a correct mutual information. ******************************************************************************/ double calculateRenyiMIJoint(double alpha, double *dataVector, double *targetVector, int vectorLength); #ifdef __cplusplus } #endif #endif MIToolbox-2.1.2/WeightedEntropy.c000066400000000000000000000101571264454417700166770ustar00rootroot00000000000000/******************************************************************************* ** WeightedEntropy.c ** Part of the mutual information toolbox ** ** Contains functions to calculate the entropy of a single variable H(X), ** the joint entropy of two variables H(X,Y), and the conditional entropy ** H(X|Y), while using a weight vector to modify the calculation. ** ** Author: Adam Pocock ** Created: 20/06/2011 ** ** Copyright 2010/2011 Adam Pocock, The University Of Manchester ** www.cs.manchester.ac.uk ** ** This file is part of MIToolbox, licensed under the 3-clause BSD license. *******************************************************************************/ #include "MIToolbox.h" #include "CalculateProbability.h" #include "WeightedEntropy.h" double calculateWeightedEntropy(double *dataVector, double *weightVector, int vectorLength) { double entropy = 0.0; double tempValue = 0.0; int i; WeightedProbState state = calculateWeightedProbability(dataVector,weightVector,vectorLength); /*H(X) = - sum p(x) log p(x)*/ for (i = 0; i < state.numStates; i++) { tempValue = state.probabilityVector[i]; if (tempValue > 0) { entropy -= state.stateWeightVector[i] * tempValue * log(tempValue); } } entropy /= log(LOG_BASE); FREE_FUNC(state.probabilityVector); state.probabilityVector = NULL; FREE_FUNC(state.stateWeightVector); state.stateWeightVector = NULL; return entropy; }/*calculateWeightedEntropy(double *,double *,int)*/ double calculateWeightedJointEntropy(double *firstVector, double *secondVector, double *weightVector, int vectorLength) { double jointEntropy = 0.0; double tempValue = 0.0; int i; WeightedJointProbState state = calculateWeightedJointProbability(firstVector,secondVector,weightVector,vectorLength); /*H(XY) = - sumx sumy p(xy) log p(xy)*/ for (i = 0; i < state.numJointStates; i++) { tempValue = state.jointProbabilityVector[i]; if (tempValue > 0) { jointEntropy -= state.jointWeightVector[i] * tempValue * log(tempValue); } } jointEntropy /= log(LOG_BASE); FREE_FUNC(state.firstProbabilityVector); state.firstProbabilityVector = NULL; FREE_FUNC(state.secondProbabilityVector); state.secondProbabilityVector = NULL; FREE_FUNC(state.jointProbabilityVector); state.jointProbabilityVector = NULL; FREE_FUNC(state.firstWeightVector); state.firstWeightVector = NULL; FREE_FUNC(state.secondWeightVector); state.secondWeightVector = NULL; FREE_FUNC(state.jointWeightVector); state.jointWeightVector = NULL; return jointEntropy; }/*calculateWeightedJointEntropy(double *,double *,double *,int)*/ double calculateWeightedConditionalEntropy(double *dataVector, double *conditionVector, double *weightVector, int vectorLength) { /* ** Conditional entropy ** H(X|Y) = - sumx sumy p(xy) log p(xy)/p(y) */ double condEntropy = 0.0; double jointValue = 0.0; double condValue = 0.0; int i; WeightedJointProbState state = calculateWeightedJointProbability(dataVector,conditionVector,weightVector,vectorLength); /*H(X|Y) = - sumx sumy p(xy) log p(xy)/p(y)*/ /* to index by numFirstStates use modulus of i ** to index by numSecondStates use integer division of i by numFirstStates */ for (i = 0; i < state.numJointStates; i++) { jointValue = state.jointProbabilityVector[i]; condValue = state.secondProbabilityVector[i / state.numFirstStates]; if ((jointValue > 0) && (condValue > 0)) { condEntropy -= state.jointWeightVector[i] * jointValue * log(jointValue / condValue); } } condEntropy /= log(LOG_BASE); FREE_FUNC(state.firstProbabilityVector); state.firstProbabilityVector = NULL; FREE_FUNC(state.secondProbabilityVector); state.secondProbabilityVector = NULL; FREE_FUNC(state.jointProbabilityVector); state.jointProbabilityVector = NULL; FREE_FUNC(state.firstWeightVector); state.firstWeightVector = NULL; FREE_FUNC(state.secondWeightVector); state.secondWeightVector = NULL; FREE_FUNC(state.jointWeightVector); state.jointWeightVector = NULL; return condEntropy; }/*calculateWeightedConditionalEntropy(double *,double *,double *,int)*/ MIToolbox-2.1.2/WeightedEntropy.h000066400000000000000000000043651264454417700167100ustar00rootroot00000000000000/******************************************************************************* ** WeightedEntropy.h ** Part of the mutual information toolbox ** ** Contains functions to calculate the entropy of a single variable H(X), ** the joint entropy of two variables H(X,Y), and the conditional entropy ** H(X|Y), while using a weight vector to modify the calculation. ** ** Author: Adam Pocock ** Created: 20/6/2011 ** ** Copyright 2010/2011 Adam Pocock, The University Of Manchester ** www.cs.manchester.ac.uk ** ** This file is part of MIToolbox, licensed under the 3-clause BSD license. *******************************************************************************/ #ifndef __WeightedEntropy_H #define __WeightedEntropy_H #ifdef __cplusplus extern "C" { #endif /******************************************************************************* ** calculateWeightedEntropy returns the entropy in log base 2 of dataVector ** H_w(X), weighted by weightVector ** ** length(vectors) == vectorLength otherwise it will segmentation fault *******************************************************************************/ double calculateWeightedEntropy(double *dataVector, double *weightVector, int vectorLength); /******************************************************************************* ** calculateWeightedJointEntropy returns the entropy in log base 2 of the joint ** variable of firstVector and secondVector, H_w(XY), weighted by weightVector ** ** length(vectors) == vectorLength otherwise it will segmentation fault *******************************************************************************/ double calculateWeightedJointEntropy(double *firstVector, double *secondVector, double *weightVector, int vectorLength); /******************************************************************************* ** calculateWeightedConditionalEntropy returns the entropy in log base 2 of ** dataVector conditioned on conditionVector, H_w(X|Y), weighted by weightVector ** ** length(vectors) == vectorLength otherwise it will segmentation fault *******************************************************************************/ double calculateWeightedConditionalEntropy(double *dataVector, double *conditionVector, double *weightVector, int vectorLength); #ifdef __cplusplus } #endif #endif MIToolbox-2.1.2/WeightedMIToolbox.m000066400000000000000000000057321264454417700171300ustar00rootroot00000000000000function [varargout] = WeightedMIToolbox(functionName, weightVector, varargin) %function [varargout] = WeightedMIToolbox(functionName, weightVector, varargin) % %Provides access to the functions in WeightedMIToolboxMex % %Expects column vectors, will not work with row vectors % %Function list %"entropy" or "h" = H(X) %"ConditionalEntropy" or "condh" = H(X|Y) %"mi" = I(X;Y) %"ConditionalMI" or "cmi" = I(X;Y|Z) % %Arguments and returned values %[entropy] = H_w(X) = H_w(vector) %[entropy] = H_w(X|Y) = H_w(vector,condition) %[mi] = I_w(X;Y) = I_w(vector,target) %[mi] = I_w(X;Y|Z) = I_w(vector,target,condition) % %Internal MIToolbox function number %Entropy = 1 %Conditional Entropy = 3 %Mutual Information = 4 %Conditional MI = 5 for i = 1:length(varargin) if (~isa(varargin{i},'double')) error('Error, MIToolbox requires inputs to be double vector or matrices') end end if (strcmpi(functionName,'Entropy') || strcmpi(functionName,'h')) %disp('Calculating Entropy'); if (size(varargin{1},2)>1) mergedVector = MIToolboxMex(3,varargin{1}); else mergedVector = varargin{1}; end [varargout{1}] = WeightedMIToolboxMex(1,weightVector,mergedVector); elseif ((strcmpi(functionName,'JointEntropy')) || strcmpi(functionName,'jointh')) if (size(varargin{1},2)>1) mergedFirst = MIToolboxMex(3,varargin{1}); else mergedFirst = varargin{1}; end if (size(varargin{2},2)>1) mergedSecond = MIToolboxMex(3,varargin{2}); else mergedSecond = varargin{2}; end [varargout{1}] = WeightedMIToolboxMex(2,weightVector,mergedFirst,mergedSecond); elseif ((strcmpi(functionName,'ConditionalEntropy')) || strcmpi(functionName,'condh')) if (size(varargin{1},2)>1) mergedFirst = MIToolboxMex(3,varargin{1}); else mergedFirst = varargin{1}; end if (size(varargin{2},2)>1) mergedSecond = MIToolboxMex(3,varargin{2}); else mergedSecond = varargin{2}; end [varargout{1}] = WeightedMIToolboxMex(3,weightVector,mergedFirst,mergedSecond); elseif (strcmpi(functionName,'mi')) if (size(varargin{1},2)>1) mergedFirst = MIToolboxMex(3,varargin{1}); else mergedFirst = varargin{1}; end if (size(varargin{2},2)>1) mergedSecond = MIToolboxMex(3,varargin{2}); else mergedSecond = varargin{2}; end [varargout{1}] = WeightedMIToolboxMex(4,weightVector,mergedFirst,mergedSecond); elseif (strcmpi(functionName,'ConditionalMI') || strcmpi(functionName,'cmi')) if (size(varargin{1},2)>1) mergedFirst = MIToolboxMex(3,varargin{1}); else mergedFirst = varargin{1}; end if (size(varargin{2},2)>1) mergedSecond = MIToolboxMex(3,varargin{2}); else mergedSecond = varargin{2}; end if (size(varargin{3},2)>1) mergedThird = MIToolboxMex(3,varargin{3}); else mergedThird = varargin{3}; end [varargout{1}] = WeightedMIToolboxMex(5,weightVector,mergedFirst,mergedSecond,mergedThird); else varargout{1} = 0; disp(['Unrecognised functionName ' functionName]); end MIToolbox-2.1.2/WeightedMIToolboxMex.c000066400000000000000000000307741264454417700175740ustar00rootroot00000000000000/******************************************************************************* ** WeightedMIToolboxMex.c ** is the MATLAB entry point for the WeightedMIToolbox functions when called ** from a MATLAB/OCTAVE script. ** ** Copyright 2011 Adam Pocock, The University Of Manchester ** www.cs.manchester.ac.uk ** ** This file is part of MIToolbox, licensed under the 3-clause BSD license. *******************************************************************************/ #include "MIToolbox.h" #include "ArrayOperations.h" #include "CalculateProbability.h" #include "WeightedEntropy.h" #include "WeightedMutualInformation.h" /******************************************************************************* **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 a flag and a variable number of arguments ** depending on the value of the flag and returns either a construct ** containing probability estimates, a merged vector or a double value ** representing an entropy or mutual information *****************************************************************************/ int flag, numberOfSamples, checkSamples, thirdCheckSamples, numberOfFeatures, checkFeatures, thirdCheckFeatures; int vectorMismatch, columnMismatch, weightLength, weightWidth; double *dataVector, *condVector, *targetVector, *firstVector, *secondVector, *output, *weightVector; double *weightOutput, *jointWeight, *firstWeight, *secondWeight, *jointOutput, *firstOutput, *secondOutput; double *numStates, *numJointStates, *numFirstStates, *numSecondStates; WeightedProbState state; WeightedJointProbState jointState; int i; switch (nrhs) { case 3: { /*printf("Must be H_w(X)\n");*/ break; } case 4: { /*printf("Must be H_w(XY), H_w(X|Y), I_w(X;Y)\n");*/ break; } case 5: { /*printf("Must be I_w(X;Y|Z)\n");*/ break; } default: { printf("Incorrect number of arguments, format is WeightedMIToolbox(\"FLAG\",varargin)\n"); break; } } vectorMismatch = 0; columnMismatch = 0; /* number to function map ** 1 = H(X) ** 2 = H(XY) ** 3 = H(X|Y) ** 4 = I(X;Y) ** 5 = I(X;Y|Z) ** 6 = calculateWeightedProbability ** 7 = calculateWeightedJointProbability */ flag = *mxGetPr(prhs[0]); weightLength = mxGetM(prhs[1]); weightWidth = mxGetN(prhs[1]); weightVector = (double *) mxGetPr(prhs[1]); switch (flag) { case 1: { /* **H(X) */ numberOfSamples = mxGetM(prhs[2]); numberOfFeatures = mxGetN(prhs[2]); dataVector = (double *) mxGetPr(prhs[2]); plhs[0] = mxCreateDoubleMatrix(1,1,mxREAL); output = (double *)mxGetPr(plhs[0]); if ((numberOfFeatures == 1) && (weightWidth == 1)) { if (numberOfSamples == weightLength) { /*double calculateWeightedEntropy(double *dataVector, double *weightVector, int vectorLength);*/ *output = calculateWeightedEntropy(dataVector,weightVector,numberOfSamples); } else { vectorMismatch = 1; } } else { columnMismatch = 1; } break; }/*case 1 - H(X)*/ case 2: { /* **H(XY) */ numberOfSamples = mxGetM(prhs[2]); checkSamples = mxGetM(prhs[3]); numberOfFeatures = mxGetN(prhs[2]); checkFeatures = mxGetN(prhs[3]); firstVector = mxGetPr(prhs[2]); secondVector = mxGetPr(prhs[3]); plhs[0] = mxCreateDoubleMatrix(1,1,mxREAL); output = (double *)mxGetPr(plhs[0]); if ((numberOfFeatures == 1) && (checkFeatures == 1) && (weightWidth == 1)) { if ((numberOfSamples == 0) && (checkSamples == 0) && (weightLength == 0)) { *output = 0.0; } else if ((numberOfSamples == 0) && (weightLength == checkSamples)) { *output = calculateWeightedEntropy(secondVector,weightVector,numberOfSamples); } else if ((checkSamples == 0) && (weightLength == numberOfSamples)) { *output = calculateWeightedEntropy(firstVector,weightVector,numberOfSamples); } else if ((numberOfSamples == checkSamples) && (numberOfSamples == weightLength)) { /*double calculateWeightedJointEntropy(double *firstVector, double *secondVector, double *weightVector, int vectorLength);*/ *output = calculateWeightedJointEntropy(firstVector,secondVector,weightVector,numberOfSamples); } else { vectorMismatch = 1; } } else { columnMismatch = 1; } break; }/*case 2 - H(XY)*/ case 3: { /* **H(X|Y) */ numberOfSamples = mxGetM(prhs[2]); checkSamples = mxGetM(prhs[3]); numberOfFeatures = mxGetN(prhs[2]); checkFeatures = mxGetN(prhs[3]); dataVector = mxGetPr(prhs[2]); condVector = mxGetPr(prhs[3]); plhs[0] = mxCreateDoubleMatrix(1,1,mxREAL); output = (double *)mxGetPr(plhs[0]); if ((numberOfFeatures == 1) && (checkFeatures == 1) && (weightWidth == 1)) { if (numberOfSamples == 0) { *output = 0.0; } else if ((checkSamples == 0) && (weightLength == numberOfSamples)) { *output = calculateWeightedEntropy(dataVector,weightVector,numberOfSamples); } else if ((numberOfSamples == checkSamples) && (numberOfSamples == weightLength)) { /*double calculateWeightedConditionalEntropy(double *dataVector, double *condVector, double *weightVector, int vectorLength);*/ *output = calculateWeightedConditionalEntropy(dataVector,condVector,weightVector,numberOfSamples); } else { vectorMismatch = 1; } } else { columnMismatch = 1; } break; }/*case 3 - H(X|Y)*/ case 4: { /* **I(X;Y) */ numberOfSamples = mxGetM(prhs[2]); checkSamples = mxGetM(prhs[3]); numberOfFeatures = mxGetN(prhs[2]); checkFeatures = mxGetN(prhs[3]); firstVector = mxGetPr(prhs[2]); secondVector = mxGetPr(prhs[3]); plhs[0] = mxCreateDoubleMatrix(1,1,mxREAL); output = (double *)mxGetPr(plhs[0]); if ((numberOfFeatures == 1) && (checkFeatures == 1) && (weightWidth == 1)) { if ((numberOfSamples == 0) || (checkSamples == 0)) { *output = 0.0; } else if ((numberOfSamples == checkSamples) && (numberOfSamples == weightLength)) { /*double calculateWeightedMutualInformation(double *firstVector, double *secondVector, double *weightVector, int vectorLength);*/ *output = calculateWeightedMutualInformation(firstVector,secondVector,weightVector,numberOfSamples); } else { vectorMismatch = 1; } } else { columnMismatch = 1; } break; }/*case 4 - I(X;Y)*/ case 5: { /* **I(X;Y|Z) */ numberOfSamples = mxGetM(prhs[2]); checkSamples = mxGetM(prhs[3]); thirdCheckSamples = mxGetM(prhs[4]); numberOfFeatures = mxGetN(prhs[2]); checkFeatures = mxGetN(prhs[3]); thirdCheckFeatures = mxGetN(prhs[4]); firstVector = mxGetPr(prhs[2]); targetVector = mxGetPr(prhs[3]); condVector = mxGetPr(prhs[4]); plhs[0] = mxCreateDoubleMatrix(1,1,mxREAL); output = (double *)mxGetPr(plhs[0]); if ((numberOfFeatures == 1) && (checkFeatures == 1) && (weightWidth == 1)) { if ((numberOfSamples == 0) || (checkSamples == 0)) { *output = 0.0; } else if ((numberOfSamples == checkSamples) && (numberOfSamples == weightLength)) { if ((thirdCheckSamples == 0) || (thirdCheckFeatures != 1)) { *output = calculateWeightedMutualInformation(firstVector,targetVector,weightVector,numberOfSamples); } else if (numberOfSamples == thirdCheckSamples) { /*double calculateWeightedConditionalMutualInformation(double *firstVector, double *targetVector, double *condVector, double *weightVector, int vectorLength);*/ *output = calculateWeightedConditionalMutualInformation(firstVector,targetVector,condVector,weightVector,numberOfSamples); } else { vectorMismatch = 1; } } else { vectorMismatch = 1; } } else { columnMismatch = 1; } break; }/*case 5 - I(X;Y|Z)*/ case 6: { /* **calculateWeightedProbability */ numberOfSamples = mxGetM(prhs[2]); dataVector = (double *) mxGetPr(prhs[2]); /*WeightedProbabilityState calculateWeightedProbability(double *dataVector, double *weightVector, int vectorLength);*/ state = calculateWeightedProbability(dataVector,weightVector,numberOfSamples); plhs[0] = mxCreateDoubleMatrix(state.numStates,1,mxREAL); plhs[1] = mxCreateDoubleMatrix(state.numStates,1,mxREAL); plhs[2] = mxCreateDoubleMatrix(1,1,mxREAL); output = (double *)mxGetPr(plhs[0]); weightOutput = (double *)mxGetPr(plhs[1]); numStates = (double *) mxGetPr(plhs[2]); *numStates = state.numStates; for (i = 0; i < state.numStates; i++) { output[i] = state.probabilityVector[i]; weightOutput[i] = state.stateWeightVector[i]; } break; }/*case 6 - calculateWeightedProbability*/ case 7: { /* **calculateWeightedJointProbability */ numberOfSamples = mxGetM(prhs[2]); firstVector = (double *) mxGetPr(prhs[2]); secondVector = (double *) mxGetPr(prhs[3]); /*JointProbabilityState calculateJointProbability(double *firstVector, double *secondVector, double *weightVector, int vectorLength);*/ jointState = calculateWeightedJointProbability(firstVector,secondVector,weightVector,numberOfSamples); plhs[0] = mxCreateDoubleMatrix(jointState.numJointStates,1,mxREAL); plhs[1] = mxCreateDoubleMatrix(jointState.numJointStates,1,mxREAL); plhs[2] = mxCreateDoubleMatrix(1,1,mxREAL); plhs[3] = mxCreateDoubleMatrix(jointState.numFirstStates,1,mxREAL); plhs[4] = mxCreateDoubleMatrix(jointState.numFirstStates,1,mxREAL); plhs[5] = mxCreateDoubleMatrix(1,1,mxREAL); plhs[6] = mxCreateDoubleMatrix(jointState.numSecondStates,1,mxREAL); plhs[7] = mxCreateDoubleMatrix(jointState.numSecondStates,1,mxREAL); plhs[8] = mxCreateDoubleMatrix(1,1,mxREAL); jointOutput = (double *)mxGetPr(plhs[0]); jointWeight = (double *)mxGetPr(plhs[1]); numJointStates = (double *) mxGetPr(plhs[2]); firstOutput = (double *)mxGetPr(plhs[3]); firstWeight = (double *)mxGetPr(plhs[4]); numFirstStates = (double *) mxGetPr(plhs[5]); secondOutput = (double *)mxGetPr(plhs[6]); secondWeight = (double *)mxGetPr(plhs[7]); numSecondStates = (double *) mxGetPr(plhs[8]); *numJointStates = jointState.numJointStates; *numFirstStates = jointState.numFirstStates; *numSecondStates = jointState.numSecondStates; for (i = 0; i < jointState.numJointStates; i++) { jointOutput[i] = jointState.jointProbabilityVector[i]; jointWeight[i] = jointState.jointWeightVector[i]; } for (i = 0; i < jointState.numFirstStates; i++) { firstOutput[i] = jointState.firstProbabilityVector[i]; firstWeight[i] = jointState.firstWeightVector[i]; } for (i = 0; i < jointState.numSecondStates; i++) { secondOutput[i] = jointState.secondProbabilityVector[i]; secondWeight[i] = jointState.secondWeightVector[i]; } break; }/*case 7 - calculateWeightedJointProbability */ default: { printf("Unrecognised flag\n"); break; }/*default*/ }/*switch(flag)*/ if (columnMismatch) { printf("No columns in one or more inputs\n"); *output = -1.0; } else if (vectorMismatch) { printf("Vector lengths do not match, they must all be the same length\n"); *output = -1.0; } return; }/*mexFunction()*/ MIToolbox-2.1.2/WeightedMutualInformation.c000066400000000000000000000101151264454417700207060ustar00rootroot00000000000000/******************************************************************************* ** WeightedMutualInformation.c ** Part of the mutual information toolbox ** ** Contains functions to calculate the mutual information of ** two variables X and Y, I(X;Y), to calculate the joint mutual information ** of two variables X & Z on the variable Y, I(XZ;Y), and the conditional ** mutual information I(X;Y|Z), while using a weight vector to modify the ** calculation. ** ** Author: Adam Pocock ** Created: 20/06/2011 ** ** Copyright 2010/2011 Adam Pocock, The University Of Manchester ** www.cs.manchester.ac.uk ** ** This file is part of MIToolbox. ** ** MIToolbox is free software: you can redistribute it and/or modify ** it under the terms of the GNU Lesser General Public License as published by ** the Free Software Foundation, either version 3 of the License, or ** (at your option) any later version. ** ** MIToolbox is distributed in the hope that it will be useful, ** but WITHOUT ANY WARRANTY; without even the implied warranty of ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ** GNU Lesser General Public License for more details. ** ** You should have received a copy of the GNU Lesser General Public License ** along with MIToolbox. If not, see . ** *******************************************************************************/ #include "MIToolbox.h" #include "ArrayOperations.h" #include "CalculateProbability.h" #include "WeightedEntropy.h" #include "WeightedMutualInformation.h" double calculateWeightedMutualInformation(double *dataVector, double *targetVector, double *weightVector, int vectorLength) { double mutualInformation = 0.0; int firstIndex,secondIndex; int i; WeightedJointProbState state = calculateWeightedJointProbability(dataVector,targetVector,weightVector,vectorLength); /* ** I(X;Y) = sum sum p(xy) * log (p(xy)/p(x)p(y)) */ for (i = 0; i < state.numJointStates; i++) { firstIndex = i % state.numFirstStates; secondIndex = i / state.numFirstStates; if ((state.jointProbabilityVector[i] > 0) && (state.firstProbabilityVector[firstIndex] > 0) && (state.secondProbabilityVector[secondIndex] > 0)) { /*double division is probably more stable than multiplying two small numbers together ** mutualInformation += state.jointProbabilityVector[i] * log(state.jointProbabilityVector[i] / (state.firstProbabilityVector[firstIndex] * state.secondProbabilityVector[secondIndex])); */ mutualInformation += state.jointWeightVector[i] * state.jointProbabilityVector[i] * log(state.jointProbabilityVector[i] / state.firstProbabilityVector[firstIndex] / state.secondProbabilityVector[secondIndex]); } } mutualInformation /= log(LOG_BASE); FREE_FUNC(state.firstProbabilityVector); state.firstProbabilityVector = NULL; FREE_FUNC(state.secondProbabilityVector); state.secondProbabilityVector = NULL; FREE_FUNC(state.jointProbabilityVector); state.jointProbabilityVector = NULL; return mutualInformation; }/*calculateWeightedMutualInformation(double *,double *,double *,int)*/ double calculateWeightedConditionalMutualInformation(double *dataVector, double *targetVector, double *conditionVector, double *weightVector, int vectorLength) { double mutualInformation = 0.0; double firstCondition, secondCondition; double *mergedVector = (double *) checkedCalloc(vectorLength,sizeof(double)); mergeArrays(targetVector,conditionVector,mergedVector,vectorLength); /* I(X;Y|Z) = H(X|Z) - H(X|YZ) */ /* double calculateWeightedConditionalEntropy(double *dataVector, double *conditionVector, double *weightVector, int vectorLength); */ firstCondition = calculateWeightedConditionalEntropy(dataVector,conditionVector,weightVector,vectorLength); secondCondition = calculateWeightedConditionalEntropy(dataVector,mergedVector,weightVector,vectorLength); mutualInformation = firstCondition - secondCondition; FREE_FUNC(mergedVector); mergedVector = NULL; return mutualInformation; }/*calculateWeightedConditionalMutualInformation(double *,double *,double *,double *,int)*/ MIToolbox-2.1.2/WeightedMutualInformation.h000066400000000000000000000037731264454417700207270ustar00rootroot00000000000000/******************************************************************************* ** WeightedMutualInformation.h ** Part of the mutual information toolbox ** ** Contains functions to calculate the mutual information of ** two variables X and Y, I(X;Y), to calculate the joint mutual information ** of two variables X & Z on the variable Y, I(XZ;Y), and the conditional ** mutual information I(x;Y|Z), while using a weight vector to modify ** the calculation. ** ** Author: Adam Pocock ** Created: 20/6/2011 ** ** Copyright 2010/2011 Adam Pocock, The University Of Manchester ** www.cs.manchester.ac.uk ** ** This file is part of MIToolbox, licensed under the 3-clause BSD license. *******************************************************************************/ #ifndef __WeightedMutualInformation_H #define __WeightedMutualInformation_H #ifdef __cplusplus extern "C" { #endif /******************************************************************************* ** calculateWeightedMutualInformation returns the log base 2 mutual information ** between dataVector and targetVector, I_w(X;Y), weighted by weightVector ** ** length(vectors) == vectorLength otherwise it will segmentation fault *******************************************************************************/ double calculateWeightedMutualInformation(double *dataVector, double *targetVector, double *weightVector, int vectorLength); /******************************************************************************* ** calculateWeightedConditionalMutualInformation returns the log base 2 ** mutual information between dataVector and targetVector, conditioned on ** conditionVector, I(X;Y|Z), weighted by weightVector ** ** length(vectors) == vectorLength otherwise it will segmentation fault *******************************************************************************/ double calculateWeightedConditionalMutualInformation(double *dataVector, double *targetVector, double *conditionVector, double *weightVector, int vectorLength); #ifdef __cplusplus } #endif #endif MIToolbox-2.1.2/cmi.m000066400000000000000000000016121264454417700143340ustar00rootroot00000000000000function output = cmi(X,Y,Z) %function output = cmi(X,Y,Z) %X, Y & Z can be matrices which are converted into a joint variable %before computation % %expects variables to be column-wise % %returns the mutual information between X and Y conditioned on Z, I(X;Y|Z) if nargin == 3 if (~isa(X,'double') || ~isa(Y,'double') || ~isa(Z,'double')) error('Error, inputs must be double vectors or matrices') end if (size(X,2)>1) mergedFirst = MIToolboxMex(3,X); else mergedFirst = X; end if (size(Y,2)>1) mergedSecond = MIToolboxMex(3,Y); else mergedSecond = Y; end if (size(Z,2)>1) mergedThird = MIToolboxMex(3,Z); else mergedThird = Z; end [output] = MIToolboxMex(8,mergedFirst,mergedSecond,mergedThird); elseif nargin == 2 if (~isa(X,'double') || ~isa(Y,'double')) error('Error, inputs must be double vectors or matrices') end output = mi(X,Y); else output = 0; end MIToolbox-2.1.2/condh.m000066400000000000000000000013521264454417700146600ustar00rootroot00000000000000function output = condh(X,Y) %function output = condh(X,Y) %X & Y can be matrices which are converted into a joint variable %before computation % %expects variables to be column-wise % %returns the conditional entropy of X given Y, H(X|Y) if nargin == 2 if (~isa(X,'double') || ~isa(Y,'double')) error('Error, inputs must be double vectors or matrices') end if (size(X,2)>1) mergedFirst = MIToolboxMex(3,X); else mergedFirst = X; end if (size(Y,2)>1) mergedSecond = MIToolboxMex(3,Y); else mergedSecond = Y; end [output] = MIToolboxMex(6,mergedFirst,mergedSecond); elseif nargin == 1 if (~isa(X,'double')) error('Error, inputs must be double vectors or matrices') end output = h(X); else output = 0; end MIToolbox-2.1.2/demonstration_algorithms/000077500000000000000000000000001264454417700205255ustar00rootroot00000000000000MIToolbox-2.1.2/demonstration_algorithms/CMIM.m000066400000000000000000000025341264454417700214340ustar00rootroot00000000000000function 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; MIToolbox-2.1.2/demonstration_algorithms/CMIM_Mex.c000066400000000000000000000117311264454417700222320ustar00rootroot00000000000000/******************************************************************************* ** 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()*/ MIToolbox-2.1.2/demonstration_algorithms/DISR.m000066400000000000000000000042421264454417700214460ustar00rootroot00000000000000function 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; MIToolbox-2.1.2/demonstration_algorithms/DISR_Mex.c000066400000000000000000000141171264454417700222470ustar00rootroot00000000000000/******************************************************************************* ** 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()*/ MIToolbox-2.1.2/demonstration_algorithms/IAMB.m000066400000000000000000000024551264454417700214210ustar00rootroot00000000000000function [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 MIToolbox-2.1.2/demonstration_algorithms/compile_demos.m000066400000000000000000000004661264454417700235300ustar00rootroot00000000000000mex -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.cMIToolbox-2.1.2/demonstration_algorithms/mRMR_D.m000066400000000000000000000036621264454417700217720ustar00rootroot00000000000000function 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; MIToolbox-2.1.2/demonstration_algorithms/mRMR_D_Mex.c000066400000000000000000000130251264454417700225630ustar00rootroot00000000000000/******************************************************************************* ** 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()*/ MIToolbox-2.1.2/h.m000066400000000000000000000006271264454417700140200ustar00rootroot00000000000000function output = h(X) %function output = h(X) %X can be a matrix which is converted into a joint variable before calculation %expects variables to be column-wise % %returns the entropy of X, H(X) if (~isa(X,'double')) error('Error, inputs must be double vectors or matrices') end if (size(X,2)>1) mergedVector = MIToolboxMex(3,X); else mergedVector = X; end [output] = MIToolboxMex(4,mergedVector); MIToolbox-2.1.2/joint.m000066400000000000000000000011461264454417700147110ustar00rootroot00000000000000function output = joint(X,arities) %function output = joint(X,arities) %returns the joint random variable of the matrix X %assuming the variables are in columns % %if passed a vector of the arities then it produces a correct %joint variable, otherwise it may not include all states % %if the joint variable is only compared with variables using the same samples, %then arity information is not required if (~isa(X,'double') || ~isa(arities,'double')) error('Error, inputs must be double vectors or matrices') end if (nargin == 2) [output] = MIToolboxMex(3,X,arities); else [output] = MIToolboxMex(3,X); end MIToolbox-2.1.2/mi.m000066400000000000000000000010551264454417700141720ustar00rootroot00000000000000function output = mi(X,Y) %function output = mi(X,Y) %X & Y can be matrices which are converted into a joint variable %before computation % %expects variables to be column-wise % %returns the mutual information between X and Y, I(X;Y) if (~isa(X,'double') || ~isa(Y,'double')) error('Error, inputs must be double vectors or matrices') end if (size(X,2)>1) mergedFirst = MIToolboxMex(3,X); else mergedFirst = X; end if (size(Y,2)>1) mergedSecond = MIToolboxMex(3,Y); else mergedSecond = Y; end [output] = MIToolboxMex(7,mergedFirst,mergedSecond);