pax_global_header00006660000000000000000000000064130342707140014513gustar00rootroot0000000000000052 comment=4a2ebfb43e08f783d31cd2ba373bf20ce9ac01dd MIToolbox-3.0.0/000077500000000000000000000000001303427071400133675ustar00rootroot00000000000000MIToolbox-3.0.0/.gitignore000066400000000000000000000003471303427071400153630ustar00rootroot00000000000000# 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 *~ # Build folder /build/MIToolbox-3.0.0/LICENSE000066400000000000000000000031121303427071400143710ustar00rootroot00000000000000This is MIToolbox, a mutual information package for C/C++ and MATLAB. Copyright (c) 2010-2017, 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-3.0.0/Makefile000066400000000000000000000030761303427071400150350ustar00rootroot00000000000000#makefile for MIToolbox #Author: Adam Pocock, apocock@cs.man.ac.uk #Created 11/03/2010 #Updated 22/02/2014 - Added make install #Updated 05/09/2016 - Repackaging # # #Copyright 2010-2017 Adam Pocock, The University Of Manchester #www.cs.manchester.ac.uk # #This file is part of MIToolbox, licensed under the 3-clause BSD license. PREFIX = /usr/local CFLAGS = -O3 -fPIC -std=c89 -Wall -Werror INCLUDES = -Iinclude CC = gcc objects = build/ArrayOperations.o build/CalculateProbability.o \ build/Entropy.o build/MutualInformation.o \ build/RenyiEntropy.o build/RenyiMutualInformation.o \ build/WeightedEntropy.o build/WeightedMutualInformation.o libMIToolbox.so : $(objects) $(CC) $(CFLAGS) -shared -o libMIToolbox.so $(objects) -lm build/%.o: src/%.c @mkdir -p build $(CC) $(CFLAGS) $(INCLUDES) -DCOMPILE_C -o build/$*.o -c $< .PHONY : debug x86 x64 matlab matlab-debug intel install test debug: $(MAKE) libMIToolbox.so "CFLAGS = -g -DDEBUG -fPIC" x86: $(MAKE) libMIToolbox.so "CFLAGS = -O3 -fPIC -m32" x64: $(MAKE) libMIToolbox.so "CFLAGS = -O3 -fPIC -m64" intel: $(MAKE) libMIToolbox.so "CC = icc" "CFLAGS = -O2 -fPIC -xHost" clean: -rm -fr build -rm -f matlab/*.o matlab/*.mex* -rm -f libMIToolbox.so install: libMIToolbox.so @echo "Installing libMIToolbox.so to $(PREFIX)/lib" @cp -v libMIToolbox.so $(PREFIX)/lib @mkdir -p $(PREFIX)/include @echo "Installing MIToolbox's header files to $(PREFIX)/include/MIToolbox" @cp -rv include/MIToolbox $(PREFIX)/include/ test: $(CC) -std=c89 $(INCLUDES) -L. -o test.out test/testMIToolbox.c -lMIToolbox MIToolbox-3.0.0/README.md000066400000000000000000000075601303427071400146560ustar00rootroot00000000000000MIToolbox ========= v3.0.0 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". MIToolbox works on discrete inputs, and all continuous values **must** be discretised before use with MIToolbox. Real-valued inputs will be discretised with x = floor(x) to ensure compatibility. MIToolbox produces unreliable results when used with continuous inputs, runs slowly and uses much more memory than usual. The discrete inputs should have small cardinality, MIToolbox will treat values {1,10,100} the same way it treats {1,2,3} and the latter will be both faster and use less memory. This limitation is due to the difficulties in estimating information theoretic functions of continuous variables. ====== 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 in the "matlab" folder. 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 - 07/01/2017 - v3.0.0 - Refactored internals to expose integer information theoretic calculations. - 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-3.0.0/include/000077500000000000000000000000001303427071400150125ustar00rootroot00000000000000MIToolbox-3.0.0/include/MIToolbox/000077500000000000000000000000001303427071400166665ustar00rootroot00000000000000MIToolbox-3.0.0/include/MIToolbox/ArrayOperations.h000066400000000000000000000114161303427071400221640ustar00rootroot00000000000000/******************************************************************************* ** 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-2017 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 #include "MIToolbox/MIToolbox.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); void printUintVector(uint *vector, int vectorLength); /******************************************************************************* ** Generates a 2D representation of a contiguous fortran style matrix *******************************************************************************/ uint **generateIntIndices(uint *featureMatrix, uint noOfSamples, uint noOfFeatures); double **generateDoubleIndices(double *featureMatrix, uint noOfSamples, uint noOfFeatures); /******************************************************************************* ** Finds the maximum state of an int array. *******************************************************************************/ int maxState(uint* 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, uint *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(uint *firstVector, uint *secondVector, uint *outputVector, int vectorLength); int discAndMergeArrays(double *firstVector, double *secondVector, uint *outputVector, int vectorLength); int mergeArraysArities(uint *firstVector, int numFirstStates, uint *secondVector, int numSecondStates, uint *outputVector, int vectorLength); int discAndMergeArraysArities(double *firstVector, int numFirstStates, double *secondVector, int numSecondStates, uint *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, uint *outputVector, int matrixWidth, int vectorLength); int mergeMultipleArraysArities(double *inputMatrix, uint *outputVector, int matrixWidth, int *arities, int vectorLength); #ifdef __cplusplus } #endif #endif MIToolbox-3.0.0/include/MIToolbox/CalculateProbability.h000066400000000000000000000134001303427071400231330ustar00rootroot00000000000000/******************************************************************************* ** 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-2017 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 #include "MIToolbox/MIToolbox.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 it will crash *******************************************************************************/ JointProbabilityState calculateJointProbability(uint *firstVector, uint *secondVector, int vectorLength); /******************************************************************************* ** discAndCalcJointProbability discretises the double vectors into int vectors, ** then generates a JointProbabilityState. ** ** length(firstVector) == length(secondVector) == vectorLength otherwise there ** will be a segmentation fault *******************************************************************************/ JointProbabilityState discAndCalcJointProbability(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(uint *firstVector, uint *secondVector, double *weightVector, int vectorLength); /******************************************************************************* ** discAndCalcWeightedJointProbability returns the same joint state object ** as discAndCalcJointProbability, with additional weightVectors giving the ** weight assigned to each state. *******************************************************************************/ WeightedJointProbState discAndCalcWeightedJointProbability(double *firstVector, double *secondVector, double *weightVector, 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(uint *dataVector, int vectorLength); /******************************************************************************* ** discAndCalcProbability discretises the double vector into an int vector, ** then generates a ProbabilityState. ** ** length(dataVector) == vectorLength otherwise there ** will be a segmentation fault *******************************************************************************/ ProbabilityState discAndCalcProbability(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(uint *dataVector, double *weightVector, int vectorLength); /******************************************************************************* ** discAndCalcWeightedProbability returns the same state object ** as discAndCalcProbability, with an additional weightVector giving the ** weight assigned to each state. *******************************************************************************/ WeightedProbState discAndCalcWeightedProbability(double *dataVector, double *weightVector, int vectorLength); /******************************************************************************* ** Frees the struct members and sets all pointers to NULL. *******************************************************************************/ void freeProbabilityState(ProbabilityState state); void freeJointProbabilityState(JointProbabilityState state); void freeWeightedProbState(WeightedProbState state); void freeWeightedJointProbState(WeightedJointProbState state); #ifdef __cplusplus } #endif #endif MIToolbox-3.0.0/include/MIToolbox/Entropy.h000066400000000000000000000053501303427071400205020ustar00rootroot00000000000000/******************************************************************************* ** 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-2017 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 #include "MIToolbox/MIToolbox.h" #include "MIToolbox/CalculateProbability.h" #ifdef __cplusplus extern "C" { #endif /******************************************************************************* ** calculateEntropy returns the entropy in log base LOG_BASE of dataVector ** H(X) ** ** length(dataVector) == vectorLength otherwise there ** will be a segmentation fault *******************************************************************************/ double discAndCalcEntropy(double *dataVector, int vectorLength); double calcEntropy(uint *dataVector, int vectorLength); /******************************************************************************* ** calculateJointEntropy returns the entropy in log base LOG_BASE of the joint ** variable of firstVector and secondVector H(X,Y) ** ** length(firstVector) == length(secondVector) == vectorLength otherwise there ** will be a segmentation fault *******************************************************************************/ double discAndCalcJointEntropy(double *firstVector, double *secondVector, int vectorLength); double calcJointEntropy(uint *firstVector, uint *secondVector, int vectorLength); /******************************************************************************* ** calculateConditionalEntropy returns the entropy in log base LOG_BASE of dataVector ** conditioned on conditionVector, H(X|Y) ** ** length(dataVector) == length(conditionVector) == vectorLength otherwise there ** will be a segmentation fault *******************************************************************************/ double discAndCalcConditionalEntropy(double *dataVector, double *conditionVector, int vectorLength); double calcConditionalEntropy(uint *dataVector, uint *conditionVector, int vectorLength); /******************************************************************************* ** Inner functions which operate on state structs. *******************************************************************************/ double entropy(ProbabilityState state); double jointEntropy(JointProbabilityState state); double condEntropy(JointProbabilityState state); #ifdef __cplusplus } #endif #endif MIToolbox-3.0.0/include/MIToolbox/MIToolbox.h000066400000000000000000000025121303427071400207130ustar00rootroot00000000000000/******************************************************************************* ** 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-2017 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 typedef unsigned int uint; #ifdef COMPILE_C #define C_IMPLEMENTATION #include #include #define CALLOC_FUNC(a,b) calloc(a,b) #define FREE_FUNC(a) free(a) #elif defined(COMPILE_R) #define R_IMPLEMENTATION #include "R.h" #define CALLOC_FUNC(a,b) Calloc((a)*(b),char) #define FREE_FUNC(a) Free((a)) #define printf Rprintf #else #define MEX_IMPLEMENTATION #include "mex.h" #define CALLOC_FUNC(a,b) mxCalloc(a,b) #define FREE_FUNC(a) mxFree(a) #define printf mexPrintf /*for Octave-3.2*/ #endif #endif MIToolbox-3.0.0/include/MIToolbox/MutualInformation.h000066400000000000000000000046741303427071400225270ustar00rootroot00000000000000/******************************************************************************* ** 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-2017 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 #include "MIToolbox/MIToolbox.h" #include "MIToolbox/CalculateProbability.h" #ifdef __cplusplus extern "C" { #endif /******************************************************************************* ** calculateMutualInformation returns the log base LOG_BASE mutual information between ** dataVector and targetVector, I(X;Y) ** ** length(dataVector) == length(targetVector) == vectorLength otherwise there ** will be a segmentation fault *******************************************************************************/ double calcMutualInformation(uint *dataVector, uint *targetVector, int vectorLength); double discAndCalcMutualInformation(double *dataVector, double *targetVector, int vectorLength); /******************************************************************************* ** calculateConditionalMutualInformation returns the log base LOG_BASE ** 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 calcConditionalMutualInformation(uint *dataVector, uint *targetVector, uint *conditionVector, int vectorLength); double discAndCalcConditionalMutualInformation(double *dataVector, double *targetVector, double *conditionVector, int vectorLength); /******************************************************************************* ** Inner functions which operate on state structs. *******************************************************************************/ double mi(JointProbabilityState state); #ifdef __cplusplus } #endif #endif MIToolbox-3.0.0/include/MIToolbox/RenyiEntropy.h000066400000000000000000000045371303427071400215170ustar00rootroot00000000000000/******************************************************************************* ** 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-2017 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_Entropy_H #define __Renyi_Entropy_H #include "MIToolbox/MIToolbox.h" #include "MIToolbox/CalculateProbability.h" #ifdef __cplusplus extern "C" { #endif /******************************************************************************* ** calculateRenyiEntropy returns the Renyi entropy in log base LOG_BASE of dataVector ** H_{\alpha}(X), for \alpha != 1 ** ** length(dataVector) == vectorLength otherwise there ** will be a segmentation fault *******************************************************************************/ double calcRenyiEntropy(double alpha, uint *dataVector, int vectorLength); double discAndCalcRenyiEntropy(double alpha, double *dataVector, int vectorLength); /******************************************************************************* ** calculateJointRenyiEntropy returns the Renyi entropy in log base LOG_BASE 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 calcJointRenyiEntropy(double alpha, uint *firstVector, uint *secondVector, int vectorLength); double discAndCalcJointRenyiEntropy(double alpha, double *firstVector, double *secondVector, int vectorLength); /******************************************************************************* ** Inner functions which operate on state structs. *******************************************************************************/ double renyiEntropy(ProbabilityState state, double alpha); double jointRenyiEntropy(JointProbabilityState state, double alpha); #ifdef __cplusplus } #endif #endif MIToolbox-3.0.0/include/MIToolbox/RenyiMutualInformation.h000066400000000000000000000047101303427071400235250ustar00rootroot00000000000000/******************************************************************************* ** 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-2017 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 #include "MIToolbox/MIToolbox.h" #include "MIToolbox/CalculateProbability.h" #ifdef __cplusplus extern "C" { #endif /******************************************************************************* ** calculateRenyiMIDivergence returns the log base LOG_BASE 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 calcRenyiMIDivergence(double alpha, uint *dataVector, uint *targetVector, int vectorLength); double discAndCalcRenyiMIDivergence(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. ** It is maintained to show how different Renyi's Information Theory is. ******************************************************************************/ double calcRenyiMIJoint(double alpha, uint *dataVector, uint *targetVector, int vectorLength); double discAndCalcRenyiMIJoint(double alpha, double *dataVector, double *targetVector, int vectorLength); /******************************************************************************* ** Inner functions which operate on state structs. *******************************************************************************/ double renyiMI(JointProbabilityState state, double alpha); #ifdef __cplusplus } #endif #endif MIToolbox-3.0.0/include/MIToolbox/WeightedEntropy.h000066400000000000000000000060031303427071400221570ustar00rootroot00000000000000/******************************************************************************* ** 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-2017 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 #include "MIToolbox/MIToolbox.h" #include "MIToolbox/CalculateProbability.h" #ifdef __cplusplus extern "C" { #endif /******************************************************************************* ** calculateWeightedEntropy returns the entropy in log base LOG_BASE of dataVector ** H_w(X), weighted by weightVector ** ** length(vectors) == vectorLength otherwise it will segmentation fault *******************************************************************************/ double calcWeightedEntropy(uint *dataVector, double *weightVector, int vectorLength); double discAndCalcWeightedEntropy(double *dataVector, double *weightVector, int vectorLength); /******************************************************************************* ** calculateWeightedJointEntropy returns the entropy in log base LOG_BASE of the joint ** variable of firstVector and secondVector, H_w(XY), weighted by weightVector ** ** length(vectors) == vectorLength otherwise it will segmentation fault *******************************************************************************/ double calcWeightedJointEntropy(uint *firstVector, uint *secondVector, double *weightVector, int vectorLength); double discAndCalcWeightedJointEntropy(double *firstVector, double *secondVector, double *weightVector, int vectorLength); /******************************************************************************* ** calculateWeightedConditionalEntropy returns the entropy in log base LOG_BASE of ** dataVector conditioned on conditionVector, H_w(X|Y), weighted by weightVector ** ** length(vectors) == vectorLength otherwise it will segmentation fault *******************************************************************************/ double calcWeightedConditionalEntropy(uint *dataVector, uint *conditionVector, double *weightVector, int vectorLength); double discAndCalcWeightedConditionalEntropy(double *dataVector, double *conditionVector, double *weightVector, int vectorLength); /******************************************************************************* ** Inner functions which operate on state structs. *******************************************************************************/ double wEntropy(WeightedProbState state); double wJointEntropy(WeightedJointProbState state); double wCondEntropy(WeightedJointProbState state); #ifdef __cplusplus } #endif #endif MIToolbox-3.0.0/include/MIToolbox/WeightedMutualInformation.h000066400000000000000000000051231303427071400241760ustar00rootroot00000000000000/******************************************************************************* ** 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-2017 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 #include "MIToolbox/MIToolbox.h" #include "MIToolbox/CalculateProbability.h" #ifdef __cplusplus extern "C" { #endif /******************************************************************************* ** calculateWeightedMutualInformation returns the LOG_BASE mutual information ** between dataVector and targetVector, I_w(X;Y), weighted by weightVector ** ** length(vectors) == vectorLength otherwise it will segmentation fault *******************************************************************************/ double calcWeightedMutualInformation(uint *dataVector, uint *targetVector, double *weightVector, int vectorLength); double discAndCalcWeightedMutualInformation(double *dataVector, double *targetVector, double *weightVector, int vectorLength); /******************************************************************************* ** calculateWeightedConditionalMutualInformation returns the LOG_BASE ** 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 calcWeightedConditionalMutualInformation(uint *dataVector, uint *targetVector, uint *conditionVector, double *weightVector, int vectorLength); double discAndCalcWeightedConditionalMutualInformation(double *dataVector, double *targetVector, double *conditionVector, double *weightVector, int vectorLength); /******************************************************************************* ** Inner functions which operate on state structs. *******************************************************************************/ double wmi(WeightedJointProbState state); #ifdef __cplusplus } #endif #endif MIToolbox-3.0.0/matlab/000077500000000000000000000000001303427071400146275ustar00rootroot00000000000000MIToolbox-3.0.0/matlab/CompileMIToolbox.m000066400000000000000000000007241303427071400201750ustar00rootroot00000000000000% Compiles the MIToolbox functions mex -I../include MIToolboxMex.c ../src/MutualInformation.c ../src/Entropy.c ../src/CalculateProbability.c ../src/ArrayOperations.c mex -I../include RenyiMIToolboxMex.c ../src/RenyiMutualInformation.c ../src/RenyiEntropy.c ../src/CalculateProbability.c ../src/ArrayOperations.c mex -I../include WeightedMIToolboxMex.c ../src/WeightedMutualInformation.c ../src/WeightedEntropy.c ../src/CalculateProbability.c ../src/ArrayOperations.c MIToolbox-3.0.0/matlab/MIToolbox.m000066400000000000000000000052031303427071400166610ustar00rootroot00000000000000function [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-3.0.0/matlab/MIToolboxMex.c000066400000000000000000000352061303427071400173270ustar00rootroot00000000000000/******************************************************************************* ** MIToolboxMex.c ** is the MATLAB entry point for the MIToolbox functions when called from ** a MATLAB/OCTAVE script. ** ** Copyright 2010-2017 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/MIToolbox.h" #include "MIToolbox/ArrayOperations.h" #include "MIToolbox/CalculateProbability.h" #include "MIToolbox/Entropy.h" #include "MIToolbox/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, *arities; uint *mergedVector; uint *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), discAndCalcProbability(X), merge(X), normaliseArray(X)\n");*/ break; } case 3: { /*printf("Must be H(XY), H(X|Y), discAndCalcJointProbability(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 = discAndCalcProbability ** 2 = discAndCalcJointProbability ** 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: { /* **discAndCalcProbability */ numberOfSamples = mxGetM(prhs[1]); dataVector = (double *) mxGetPr(prhs[1]); /*ProbabilityState discAndCalcProbability(double *dataVector, int vectorLength);*/ state = discAndCalcProbability(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]; } freeProbabilityState(state); break; }/*case 1 - discAndCalcProbability*/ case 2: { /* **discAndCalcJointProbability */ numberOfSamples = mxGetM(prhs[1]); firstVector = (double *) mxGetPr(prhs[1]); secondVector = (double *) mxGetPr(prhs[2]); /*JointProbabilityState discAndCalcJointProbability(double *firstVector, double *secondVector int vectorLength);*/ jointState = discAndCalcJointProbability(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]; } freeJointProbabilityState(jointState); break; }/*case 2 - discAndCalcJointProbability */ 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 = (uint *) mxCalloc(numberOfSamples,sizeof(uint)); 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 = (uint *) mxCalloc(numberOfSamples,sizeof(uint)); arities = (double *) mxGetPr(prhs[2]); intArities = (uint *) mxCalloc(numberOfFeatures,sizeof(uint)); for (i = 0; i < numArities; i++) { intArities[i] = (uint) 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 discAndCalcEntropy(double *dataVector, int vectorLength);*/ *output = discAndCalcEntropy(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 = discAndCalcEntropy(secondVector,numberOfSamples); } else if (checkSamples == 0) { *output = discAndCalcEntropy(firstVector,numberOfSamples); } else if (numberOfSamples == checkSamples) { /*double discAndCalcJointEntropy(double *firstVector, double *secondVector, int vectorLength);*/ *output = discAndCalcJointEntropy(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 = discAndCalcEntropy(dataVector,numberOfSamples); } else if (numberOfSamples == checkSamples) { /*double discAndCalcConditionalEntropy(double *dataVector, double *condVector, int vectorLength);*/ *output = discAndCalcConditionalEntropy(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 discAndCalcMutualInformation(double *firstVector, double *secondVector, int vectorLength);*/ *output = discAndCalcMutualInformation(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 = discAndCalcMutualInformation(firstVector,targetVector,numberOfSamples); } else if ((numberOfSamples == checkSamples) && (numberOfSamples == thirdCheckSamples)) { /*double discAndCalcConditionalMutualInformation(double *firstVector, double *targetVector, double *condVector, int vectorLength);*/ *output = discAndCalcConditionalMutualInformation(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 = (uint *) mxCalloc(numberOfSamples,sizeof(uint)); 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]; } mxFree(outputIntVector); outputIntVector = NULL; break; }/*case 9 - normaliseArray*/ default: { printf("Unrecognised flag\n"); break; }/*default*/ }/*switch(flag)*/ return; }/*mexFunction()*/ MIToolbox-3.0.0/matlab/RenyiMIToolbox.m000066400000000000000000000031371303427071400176740ustar00rootroot00000000000000function [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-3.0.0/matlab/RenyiMIToolboxMex.c000066400000000000000000000123561303427071400203370ustar00rootroot00000000000000/******************************************************************************* ** RenyiMIToolboxMex.c ** is the MATLAB entry point for the Renyi Entropy and MI MIToolbox functions ** when called from a MATLAB/OCTAVE script. ** ** Copyright 2010-2017 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/MIToolbox.h" #include "MIToolbox/ArrayOperations.h" #include "MIToolbox/RenyiEntropy.h" #include "MIToolbox/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 discAndCalcRenyiEntropy(double alpha, double *dataVector, long vectorLength);*/ *output = discAndCalcRenyiEntropy(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 discAndCalcJointRenyiEntropy(double alpha, double *firstVector, double *secondVector, long vectorLength);*/ *output = discAndCalcJointRenyiEntropy(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 discAndCalcRenyiMIDivergence(double alpha, double *dataVector, double *targetVector, long vectorLength);*/ *output = discAndCalcRenyiMIDivergence(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-3.0.0/matlab/WeightedMIToolbox.m000066400000000000000000000057321303427071400203510ustar00rootroot00000000000000function [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-3.0.0/matlab/WeightedMIToolboxMex.c000066400000000000000000000311431303427071400210040ustar00rootroot00000000000000/******************************************************************************* ** WeightedMIToolboxMex.c ** is the MATLAB entry point for the WeightedMIToolbox functions when called ** from a MATLAB/OCTAVE script. ** ** Copyright 2011-2017 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/MIToolbox.h" #include "MIToolbox/ArrayOperations.h" #include "MIToolbox/CalculateProbability.h" #include "MIToolbox/WeightedEntropy.h" #include "MIToolbox/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 = discAndCalcWeightedProbability ** 7 = discAndCalcWeightedJointProbability */ 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 discAndCalcWeightedEntropy(double *dataVector, double *weightVector, int vectorLength);*/ *output = discAndCalcWeightedEntropy(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 = discAndCalcWeightedEntropy(secondVector,weightVector,numberOfSamples); } else if ((checkSamples == 0) && (weightLength == numberOfSamples)) { *output = discAndCalcWeightedEntropy(firstVector,weightVector,numberOfSamples); } else if ((numberOfSamples == checkSamples) && (numberOfSamples == weightLength)) { /*double discAndCalcWeightedJointEntropy(double *firstVector, double *secondVector, double *weightVector, int vectorLength);*/ *output = discAndCalcWeightedJointEntropy(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 = discAndCalcWeightedEntropy(dataVector,weightVector,numberOfSamples); } else if ((numberOfSamples == checkSamples) && (numberOfSamples == weightLength)) { /*double discAndCalcWeightedConditionalEntropy(double *dataVector, double *condVector, double *weightVector, int vectorLength);*/ *output = discAndCalcWeightedConditionalEntropy(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 discAndCalcWeightedMutualInformation(double *firstVector, double *secondVector, double *weightVector, int vectorLength);*/ *output = discAndCalcWeightedMutualInformation(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 = discAndCalcWeightedMutualInformation(firstVector,targetVector,weightVector,numberOfSamples); } else if (numberOfSamples == thirdCheckSamples) { /*double discAndCalcWeightedConditionalMutualInformation(double *firstVector, double *targetVector, double *condVector, double *weightVector, int vectorLength);*/ *output = discAndCalcWeightedConditionalMutualInformation(firstVector,targetVector,condVector,weightVector,numberOfSamples); } else { vectorMismatch = 1; } } else { vectorMismatch = 1; } } else { columnMismatch = 1; } break; }/*case 5 - I(X;Y|Z)*/ case 6: { /* **discAndCalcWeightedProbability */ numberOfSamples = mxGetM(prhs[2]); dataVector = (double *) mxGetPr(prhs[2]); /*WeightedProbabilityState discAndCalcWeightedProbability(double *dataVector, double *weightVector, int vectorLength);*/ state = discAndCalcWeightedProbability(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 - discAndCalcWeightedProbability*/ case 7: { /* **discAndCalcWeightedJointProbability */ numberOfSamples = mxGetM(prhs[2]); firstVector = (double *) mxGetPr(prhs[2]); secondVector = (double *) mxGetPr(prhs[3]); /*JointProbabilityState discAndCalcJointProbability(double *firstVector, double *secondVector, double *weightVector, int vectorLength);*/ jointState = discAndCalcWeightedJointProbability(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 - discAndCalcWeightedJointProbability */ 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-3.0.0/matlab/cmi.m000066400000000000000000000016121303427071400155550ustar00rootroot00000000000000function 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-3.0.0/matlab/condh.m000066400000000000000000000013521303427071400161010ustar00rootroot00000000000000function 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-3.0.0/matlab/demonstration_algorithms/000077500000000000000000000000001303427071400217465ustar00rootroot00000000000000MIToolbox-3.0.0/matlab/demonstration_algorithms/CMIM.m000066400000000000000000000025341303427071400226550ustar00rootroot00000000000000function 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-3.0.0/matlab/demonstration_algorithms/CMIM_Mex.c000066400000000000000000000120121303427071400234440ustar00rootroot00000000000000/******************************************************************************* ** 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 "MIToolbox/MIToolbox.h" #include "MIToolbox/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] = discAndCalcMutualInformation(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 discAndCalcConditionalMutualInformation(double *firstVector, double *targetVector, double *conditionVector, int vectorLength);*/ currentFeature = (int) outputFeatures[lastUsedFeature[j]]; conditionalInfo = discAndCalcConditionalMutualInformation(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-3.0.0/matlab/demonstration_algorithms/CompileDemos.m000066400000000000000000000006321303427071400245050ustar00rootroot00000000000000 mex -I../../ mRMR_D_Mex.c ../../src/MutualInformation.c ../../src/ArrayOperations.c ../../src/CalculateProbability.c ../../src/Entropy.c mex -I../../ DISR_Mex.c ../../src/MutualInformation.c ../../src/ArrayOperations.c ../../src/CalculateProbability.c ../../src/Entropy.c mex -I../../ CMIM_Mex.c ../../src/MutualInformation.c ../../src/ArrayOperations.c ../../src/CalculateProbability.c ../../src/Entropy.c MIToolbox-3.0.0/matlab/demonstration_algorithms/DISR.m000066400000000000000000000042421303427071400226670ustar00rootroot00000000000000function 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-3.0.0/matlab/demonstration_algorithms/DISR_Mex.c000066400000000000000000000143771303427071400235000ustar00rootroot00000000000000/******************************************************************************* ** 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 "MIToolbox/MIToolbox.h" #include "MIToolbox/MutualInformation.h" #include "MIToolbox/Entropy.h" #include "MIToolbox/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; int *mergedVector = (int *) mxCalloc(noOfSamples,sizeof(int)); int *classColumnInt = (int *) mxCalloc(noOfSamples,sizeof(int)); int arrayPosition; double mi, tripEntropy; int i,j,x; normaliseArray(classColumn,classColumnInt,noOfSamples); 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 discAndCalcMutualInformation(double *firstVector, double *secondVector, int vectorLength); */ classMI[i] = discAndCalcMutualInformation(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 calcMutualInformation(double *firstVector, double *secondVector, int vectorLength); **double calcJointEntropy(double *firstVector, double *secondVector, int vectorLength); */ discAndMergeArrays(feature2D[(int) outputFeatures[x]], feature2D[j],mergedVector,noOfSamples); mi = calcMutualInformation(mergedVector, classColumnInt, noOfSamples); tripEntropy = calcJointEntropy(mergedVector, classColumnInt, 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-3.0.0/matlab/demonstration_algorithms/IAMB.m000066400000000000000000000024551303427071400226420ustar00rootroot00000000000000function [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-3.0.0/matlab/demonstration_algorithms/mRMR_D.m000066400000000000000000000036621303427071400232130ustar00rootroot00000000000000function 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-3.0.0/matlab/demonstration_algorithms/mRMR_D_Mex.c000066400000000000000000000131101303427071400237770ustar00rootroot00000000000000/******************************************************************************* ** 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 "MIToolbox/MIToolbox.h" #include "MIToolbox/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] = discAndCalcMutualInformation(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 discAndCalcMutualInformation(double *firstVector, double *secondVector, int vectorLength);*/ featureMIMatrix[arrayPosition] = discAndCalcMutualInformation(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-3.0.0/matlab/h.m000066400000000000000000000006271303427071400152410ustar00rootroot00000000000000function 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-3.0.0/matlab/joint.m000066400000000000000000000011461303427071400161320ustar00rootroot00000000000000function 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-3.0.0/matlab/mi.m000066400000000000000000000010551303427071400154130ustar00rootroot00000000000000function 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); MIToolbox-3.0.0/src/000077500000000000000000000000001303427071400141565ustar00rootroot00000000000000MIToolbox-3.0.0/src/ArrayOperations.c000066400000000000000000000276101303427071400174520ustar00rootroot00000000000000/******************************************************************************* ** 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-2017 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/MIToolbox.h" #include "MIToolbox/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"); #elif defined(C_IMPLEMENTATION) fprintf(stderr, "Error: %s\nAttempted to allocate %lu length of size %lu\n", strerror(errno), vectorLength, sizeOfType); exit(EXIT_FAILURE); #endif } return allocated; } void incrementVector(double* vector, int vectorLength) { /*This is used to map from C indices to MATLAB indices*/ int i = 0; for (i = 0; i < vectorLength; i++) { vector[i]++; }/*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("Value 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("Value at i=%d, is %d\n",i,vector[i]); }/*for number of items in vector*/ }/*printIntVector(int*,int)*/ void printUintVector(uint *vector, int vectorLength) { int i; for (i = 0; i < vectorLength; i++) { printf("Value at i=%d, is %d\n",i,vector[i]); }/*for number of items in vector*/ }/*printUintVector(int*,int)*/ uint **generateIntIndices(uint *featureMatrix, uint noOfSamples, uint noOfFeatures) { int j; uint **feature2D = (uint **) checkedCalloc(noOfFeatures,sizeof(uint *)); for (j = 0; j < noOfFeatures; j++) { feature2D[j] = featureMatrix + j*noOfSamples; } return feature2D; } double **generateDoubleIndices(double *featureMatrix, uint noOfSamples, uint noOfFeatures) { int j; double **feature2D = (double **) checkedCalloc(noOfFeatures,sizeof(double *)); for (j = 0; j < noOfFeatures; j++) { feature2D[j] = featureMatrix + j*noOfSamples; } return feature2D; } int maxState(uint *vector, int vectorLength) { int i, max; max = 0; for (i = 0; i < vectorLength; i++) { if (vector[i] > max) { max = vector[i]; } } return max + 1; } 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, uint *outputVector, int vectorLength) { int minVal = 0; int maxVal = 0; int currentValue; int i; if (vectorLength > 0) { int* tempVector = (int*) checkedCalloc(vectorLength,sizeof(int)); minVal = (int) floor(inputVector[0]); maxVal = (int) floor(inputVector[0]); for (i = 0; i < vectorLength; i++) { currentValue = (int) floor(inputVector[i]); tempVector[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] = tempVector[i] - minVal; } maxVal = (maxVal - minVal) + 1; FREE_FUNC(tempVector); tempVector = NULL; } 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(uint *firstVector, uint *secondVector, uint *outputVector, int vectorLength) { int firstNumStates = maxState(firstVector,vectorLength); int secondNumStates = maxState(secondVector,vectorLength); uint *stateMap = (uint *) checkedCalloc(firstNumStates*secondNumStates,sizeof(uint)); int stateCount = 1; int i, curIndex; for (i = 0; i < vectorLength; i++) { curIndex = firstVector[i] + (secondVector[i] * firstNumStates); if (stateMap[curIndex] == 0) { stateMap[curIndex] = stateCount; stateCount++; } outputVector[i] = stateMap[curIndex]; } FREE_FUNC(stateMap); stateMap = NULL; return stateCount; }/*mergeArrays(double *,double *,int *, int)*/ /******************************************************************************* ** discAndMergeArrays takes in two arrays, discretises them 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 discAndMergeArrays(double *firstVector, double *secondVector, uint *outputVector, int vectorLength) { uint *firstNormalisedVector; uint *secondNormalisedVector; int stateCount; firstNormalisedVector = (uint *) checkedCalloc(vectorLength,sizeof(uint)); secondNormalisedVector = (uint *) checkedCalloc(vectorLength,sizeof(uint)); normaliseArray(firstVector,firstNormalisedVector,vectorLength); normaliseArray(secondVector,secondNormalisedVector,vectorLength); stateCount = mergeArrays(firstNormalisedVector,secondNormalisedVector,outputVector,vectorLength); FREE_FUNC(firstNormalisedVector); FREE_FUNC(secondNormalisedVector); firstNormalisedVector = NULL; secondNormalisedVector = NULL; return stateCount; }/*discAndMergeArrays(double *,double *,int *, int)*/ int mergeArraysArities(uint *firstVector, int numFirstStates, uint *secondVector, int numSecondStates, uint *outputVector, int vectorLength) { int i; int totalStates; int firstStateCheck, secondStateCheck; firstStateCheck = maxState(firstVector,vectorLength); secondStateCheck = maxState(secondVector,vectorLength); if ((firstStateCheck <= numFirstStates) && (secondStateCheck <= numSecondStates)) { for (i = 0; i < vectorLength; i++) { outputVector[i] = firstVector[i] + (secondVector[i] * numFirstStates) + 1; } totalStates = numFirstStates * numSecondStates; } else { totalStates = -1; } return totalStates; }/*mergeArraysArities(int *,int,int *,int,int *,int)*/ int discAndMergeArraysArities(double *firstVector, int numFirstStates, double *secondVector, int numSecondStates, uint *outputVector, int vectorLength) { uint *firstNormalisedVector; uint *secondNormalisedVector; int i; int totalStates; int firstStateCheck, secondStateCheck; firstNormalisedVector = (uint *) checkedCalloc(vectorLength,sizeof(uint)); secondNormalisedVector = (uint *) checkedCalloc(vectorLength,sizeof(uint)); 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, uint *outputVector, int matrixWidth, int vectorLength) { int i = 0; int currentIndex; int currentNumStates; uint *normalisedVector = (uint *) checkedCalloc(vectorLength, sizeof(uint)); if (matrixWidth > 1) { currentNumStates = discAndMergeArrays(inputMatrix, (inputMatrix + vectorLength), outputVector, vectorLength); for (i = 2; i < matrixWidth; i++) { currentIndex = i * vectorLength; normaliseArray(inputMatrix+currentIndex, normalisedVector, vectorLength); currentNumStates = mergeArrays(outputVector, normalisedVector, outputVector, vectorLength); } } else { currentNumStates = normaliseArray(inputMatrix, normalisedVector, vectorLength); for (i = 0; i < vectorLength; i++) { outputVector[i] = normalisedVector[i]; } } FREE_FUNC(normalisedVector); normalisedVector = NULL; return currentNumStates; }/*mergeMultipleArrays(double *, double *, int, int)*/ int mergeMultipleArraysArities(double *inputMatrix, uint *outputVector, int matrixWidth, int *arities, int vectorLength) { int i = 0; int currentIndex; int currentNumStates; uint *normalisedVector = (uint *) checkedCalloc(vectorLength, sizeof(uint)); if (matrixWidth > 1) { currentNumStates = discAndMergeArraysArities(inputMatrix, arities[0], (inputMatrix + vectorLength), arities[1], outputVector, vectorLength); for (i = 2; i < matrixWidth; i++) { currentIndex = i * vectorLength; normaliseArray(inputMatrix+currentIndex, normalisedVector, vectorLength); currentNumStates = mergeArraysArities(outputVector, currentNumStates, normalisedVector, arities[i], outputVector, vectorLength); if (currentNumStates == -1) { break; } } } else { currentNumStates = normaliseArray(inputMatrix, normalisedVector, vectorLength); for (i = 0; i < vectorLength; i++) { outputVector[i] = normalisedVector[i]; } } FREE_FUNC(normalisedVector); normalisedVector = NULL; return currentNumStates; }/*mergeMultipleArraysArities(double *, double *, int, int)*/ MIToolbox-3.0.0/src/CalculateProbability.c000066400000000000000000000275221303427071400204300ustar00rootroot00000000000000/******************************************************************************* ** 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-2017 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/MIToolbox.h" #include "MIToolbox/ArrayOperations.h" #include "MIToolbox/CalculateProbability.h" JointProbabilityState calculateJointProbability(uint *firstVector, uint *secondVector, int vectorLength) { 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; firstNumStates = maxState(firstVector,vectorLength); secondNumStates = maxState(secondVector,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 for number of FP operations now O(states) instead of O(vectorLength) */ for (i = 0; i < vectorLength; i++) { firstStateCounts[firstVector[i]] += 1; secondStateCounts[secondVector[i]] += 1; jointStateCounts[secondVector[i] * firstNumStates + firstVector[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(firstStateCounts); FREE_FUNC(secondStateCounts); FREE_FUNC(jointStateCounts); 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; }/*calcJointProbability(uint *,uint *, int)*/ JointProbabilityState discAndCalcJointProbability(double *firstVector, double *secondVector, int vectorLength) { uint *firstNormalisedVector = (uint *) checkedCalloc(vectorLength,sizeof(uint)); uint *secondNormalisedVector = (uint *) checkedCalloc(vectorLength,sizeof(uint)); JointProbabilityState state; normaliseArray(firstVector,firstNormalisedVector,vectorLength); normaliseArray(secondVector,secondNormalisedVector,vectorLength); state = calculateJointProbability(firstNormalisedVector,secondNormalisedVector,vectorLength); FREE_FUNC(firstNormalisedVector); FREE_FUNC(secondNormalisedVector); firstNormalisedVector = NULL; secondNormalisedVector = NULL; return state; }/*discAndCalcJointProbability(double *,double *, int)*/ WeightedJointProbState calculateWeightedJointProbability(uint *firstVector, uint *secondVector, double *weightVector, int vectorLength) { 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; firstNumStates = maxState(firstVector,vectorLength); secondNumStates = maxState(secondVector,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[firstVector[i]] += 1; secondStateCounts[secondVector[i]] += 1; jointStateCounts[secondVector[i] * firstNumStates + firstVector[i]] += 1; firstWeightVec[firstVector[i]] += weightVector[i]; secondWeightVec[secondVector[i]] += weightVector[i]; jointWeightVec[secondVector[i] * firstNumStates + firstVector[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(firstStateCounts); FREE_FUNC(secondStateCounts); FREE_FUNC(jointStateCounts); 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; } WeightedJointProbState discAndCalcWeightedJointProbability(double *firstVector, double *secondVector, double *weightVector, int vectorLength) { uint *firstNormalisedVector = (uint *) checkedCalloc(vectorLength,sizeof(uint)); uint *secondNormalisedVector = (uint *) checkedCalloc(vectorLength,sizeof(uint)); WeightedJointProbState state; normaliseArray(firstVector,firstNormalisedVector,vectorLength); normaliseArray(secondVector,secondNormalisedVector,vectorLength); state = calculateWeightedJointProbability(firstNormalisedVector,secondNormalisedVector,weightVector,vectorLength); FREE_FUNC(firstNormalisedVector); FREE_FUNC(secondNormalisedVector); firstNormalisedVector = NULL; secondNormalisedVector = NULL; return state; }/*discAndCalcWeightedJointProbability(double *,double *,double *,int)*/ ProbabilityState calculateProbability(uint* dataVector, int vectorLength) { int numStates; int *stateCounts; double *stateProbs; ProbabilityState state; int i; double length = vectorLength; numStates = maxState(dataVector,vectorLength); stateCounts = (int *) checkedCalloc(numStates,sizeof(int)); stateProbs = (double *) checkedCalloc(numStates,sizeof(double)); /* Optimised for number of FP operations now O(states) instead of O(vectorLength) */ for (i = 0; i < vectorLength; i++) { stateCounts[dataVector[i]] += 1; } for (i = 0; i < numStates; i++) { stateProbs[i] = stateCounts[i] / length; } FREE_FUNC(stateCounts); stateCounts = NULL; state.probabilityVector = stateProbs; state.numStates = numStates; return state; } ProbabilityState discAndCalcProbability(double *dataVector, int vectorLength) { uint *normalisedVector = (uint *) checkedCalloc(vectorLength,sizeof(uint)); normaliseArray(dataVector,normalisedVector,vectorLength); ProbabilityState state = calculateProbability(normalisedVector,vectorLength); FREE_FUNC(normalisedVector); normalisedVector = NULL; return state; }/*discAndCalcProbability(double *,int)*/ WeightedProbState calculateWeightedProbability(uint *dataVector, double *weightVector, int vectorLength) { int *stateCounts; double *stateProbs; double *stateWeights; int numStates, i; WeightedProbState state; double length = vectorLength; numStates = maxState(dataVector,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[dataVector[i]] += 1; stateWeights[dataVector[i]] += weightVector[i]; } for (i = 0; i < numStates; i++) { stateProbs[i] = stateCounts[i] / length; stateWeights[i] /= stateCounts[i]; } FREE_FUNC(stateCounts); stateCounts = NULL; state.probabilityVector = stateProbs; state.stateWeightVector = stateWeights; state.numStates = numStates; return state; }/*calculateWeightedProbability(int *, double *, int)*/ WeightedProbState discAndCalcWeightedProbability(double *dataVector, double *weightVector, int vectorLength) { WeightedProbState state; uint *normalisedVector = (uint *) checkedCalloc(vectorLength,sizeof(uint)); normaliseArray(dataVector,normalisedVector,vectorLength); state = calculateWeightedProbability(normalisedVector,weightVector,vectorLength); FREE_FUNC(normalisedVector); normalisedVector = NULL; return state; }/*discAndCalcWeightedProbability(double *, double *, int)*/ void freeProbabilityState(ProbabilityState state) { FREE_FUNC(state.probabilityVector); state.probabilityVector = NULL; } void freeJointProbabilityState(JointProbabilityState state) { FREE_FUNC(state.firstProbabilityVector); state.firstProbabilityVector = NULL; FREE_FUNC(state.secondProbabilityVector); state.secondProbabilityVector = NULL; FREE_FUNC(state.jointProbabilityVector); state.jointProbabilityVector = NULL; } void freeWeightedProbState(WeightedProbState state) { FREE_FUNC(state.probabilityVector); state.probabilityVector = NULL; FREE_FUNC(state.stateWeightVector); state.stateWeightVector = NULL; } void freeWeightedJointProbState(WeightedJointProbState state) { 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; } MIToolbox-3.0.0/src/Entropy.c000066400000000000000000000100551303427071400157630ustar00rootroot00000000000000/******************************************************************************* ** 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-2017 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/MIToolbox.h" #include "MIToolbox/CalculateProbability.h" #include "MIToolbox/Entropy.h" double entropy(ProbabilityState state) { double entropy = 0.0; double tempValue = 0.0; int i; /*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(LOG_BASE); return entropy; } double discAndCalcEntropy(double* dataVector, int vectorLength) { ProbabilityState state = discAndCalcProbability(dataVector, vectorLength); double h = entropy(state); freeProbabilityState(state); return h; }/*discAndCalcEntropy(double* ,int)*/ double calcEntropy(uint* dataVector, int vectorLength) { ProbabilityState state = calculateProbability(dataVector, vectorLength); double h = entropy(state); freeProbabilityState(state); return h; }/*calcEntropy(uint* ,int)*/ double jointEntropy(JointProbabilityState state) { double jointEntropy = 0.0; double tempValue = 0.0; int i; /*H(XY) = - \sum_x \sum_y 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(LOG_BASE); return jointEntropy; } double discAndCalcJointEntropy(double *firstVector, double *secondVector, int vectorLength) { JointProbabilityState state = discAndCalcJointProbability(firstVector, secondVector, vectorLength); double h = jointEntropy(state); freeJointProbabilityState(state); return h; }/*discAndCalcJointEntropy(double *, double *, int)*/ double calcJointEntropy(uint *firstVector, uint *secondVector, int vectorLength) { JointProbabilityState state = calculateJointProbability(firstVector, secondVector, vectorLength); double h = jointEntropy(state); freeJointProbabilityState(state); return h; }/*calcJointEntropy(uint *, uint *, int)*/ double condEntropy(JointProbabilityState state) { double condEntropy = 0.0; double jointValue = 0.0; double condValue = 0.0; int i; /*H(X|Y) = - \sum_x \sum_y p(x,y) \log p(x,y)/p(y)*/ /** Indexing by numFirstStates use modulus of i ** Indexing 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(LOG_BASE); return condEntropy; } double discAndCalcConditionalEntropy(double *dataVector, double *conditionVector, int vectorLength) { JointProbabilityState state = discAndCalcJointProbability(dataVector, conditionVector, vectorLength); double h = condEntropy(state); freeJointProbabilityState(state); return h; }/*discAndCalcConditionalEntropy(double *, double *, int)*/ double calcConditionalEntropy(uint *dataVector, uint *conditionVector, int vectorLength) { JointProbabilityState state = calculateJointProbability(dataVector, conditionVector, vectorLength); double h = condEntropy(state); freeJointProbabilityState(state); return h; }/*calcConditionalEntropy(uint *, uint *, int)*/ MIToolbox-3.0.0/src/MutualInformation.c000066400000000000000000000117711303427071400200060ustar00rootroot00000000000000/******************************************************************************* ** 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-2017 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/MIToolbox.h" #include "MIToolbox/ArrayOperations.h" #include "MIToolbox/CalculateProbability.h" #include "MIToolbox/Entropy.h" #include "MIToolbox/MutualInformation.h" double mi(JointProbabilityState state) { double mutualInformation = 0.0; int firstIndex,secondIndex; int i; /* ** I(X;Y) = \sum_x \sum_y p(x,y) * \log (p(x,y)/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(LOG_BASE); return mutualInformation; }/*mi(JointProbabilityState)*/ double calcMutualInformation(uint *dataVector, uint *targetVector, int vectorLength) { JointProbabilityState state = calculateJointProbability(dataVector,targetVector,vectorLength); double mutualInformation = mi(state); freeJointProbabilityState(state); return mutualInformation; }/*calculateMutualInformation(uint *,uint *,int)*/ double discAndCalcMutualInformation(double *dataVector, double *targetVector, int vectorLength) { JointProbabilityState state = discAndCalcJointProbability(dataVector,targetVector,vectorLength); double mutualInformation = mi(state); freeJointProbabilityState(state); return mutualInformation; }/*discAndCalcMutualInformation(double *,double *,int)*/ double calcConditionalMutualInformation(uint *dataVector, uint *targetVector, uint *conditionVector, int vectorLength) { double mutualInformation = 0.0; double firstCondition, secondCondition; uint *mergedVector = (uint *) checkedCalloc(vectorLength,sizeof(uint)); mergeArrays(targetVector,conditionVector,mergedVector,vectorLength); /* I(X;Y|Z) = H(X|Z) - H(X|YZ) */ /* double calculateConditionalEntropy(double *dataVector, double *conditionVector, int vectorLength); */ firstCondition = calcConditionalEntropy(dataVector,conditionVector,vectorLength); secondCondition = calcConditionalEntropy(dataVector,mergedVector,vectorLength); mutualInformation = firstCondition - secondCondition; FREE_FUNC(mergedVector); mergedVector = NULL; return mutualInformation; }/*calculateConditionalMutualInformation(double *,double *,double *,int)*/ double discAndCalcConditionalMutualInformation(double *dataVector, double *targetVector, double *conditionVector, int vectorLength) { double mutualInformation = 0.0; double firstCondition, secondCondition; uint *dataNormVector = (uint *) checkedCalloc(vectorLength,sizeof(uint)); uint *targetNormVector = (uint *) checkedCalloc(vectorLength,sizeof(uint)); uint *conditionNormVector = (uint *) checkedCalloc(vectorLength,sizeof(uint)); uint *mergedVector = (uint *) checkedCalloc(vectorLength,sizeof(uint)); normaliseArray(dataVector,dataNormVector,vectorLength); normaliseArray(targetVector,targetNormVector,vectorLength); normaliseArray(conditionVector,conditionNormVector,vectorLength); mergeArrays(targetNormVector,conditionNormVector,mergedVector,vectorLength); /* I(X;Y|Z) = H(X|Z) - H(X|YZ) */ /* double calculateConditionalEntropy(double *dataVector, double *conditionVector, int vectorLength); */ firstCondition = calcConditionalEntropy(dataNormVector,conditionNormVector,vectorLength); secondCondition = calcConditionalEntropy(dataNormVector,mergedVector,vectorLength); mutualInformation = firstCondition - secondCondition; FREE_FUNC(dataNormVector); FREE_FUNC(targetNormVector); FREE_FUNC(conditionNormVector); FREE_FUNC(mergedVector); dataNormVector = NULL; targetNormVector = NULL; conditionNormVector = NULL; mergedVector = NULL; return mutualInformation; }/*calculateConditionalMutualInformation(double *,double *,double *,int)*/ MIToolbox-3.0.0/src/RenyiEntropy.c000066400000000000000000000060441303427071400167750ustar00rootroot00000000000000/******************************************************************************* ** 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-2017 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/MIToolbox.h" #include "MIToolbox/ArrayOperations.h" #include "MIToolbox/CalculateProbability.h" #include "MIToolbox/Entropy.h" double renyiEntropy(ProbabilityState state, double alpha) { double entropy = 0.0; double tempValue = 0.0; int i; /*H_\alpha(X) = 1/(1-alpha) * \log(\sum_x p(x)^alpha)*/ for (i = 0; i < state.numStates; i++) { tempValue = state.probabilityVector[i]; if (tempValue > 0) { entropy += pow(tempValue,alpha); } } entropy = log(entropy); entropy /= log(LOG_BASE); entropy /= (1.0-alpha); return entropy; } double calcRenyiEntropy(double alpha, uint *dataVector, int vectorLength) { ProbabilityState state = calculateProbability(dataVector,vectorLength); double h = renyiEntropy(state,alpha); freeProbabilityState(state); return h; }/*calcRenyiEntropy(double,uint*,int)*/ double discAndCalcRenyiEntropy(double alpha, double *dataVector, int vectorLength) { ProbabilityState state = discAndCalcProbability(dataVector,vectorLength); double h = renyiEntropy(state,alpha); freeProbabilityState(state); return h; }/*discAndCalcRenyiEntropy(double,double*,int)*/ double jointRenyiEntropy(JointProbabilityState state, double alpha) { double jointEntropy = 0.0; double tempValue = 0.0; int i; /*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(LOG_BASE); jointEntropy /= (1.0-alpha); return jointEntropy; } double calcJointRenyiEntropy(double alpha, uint *firstVector, uint *secondVector, int vectorLength) { JointProbabilityState state = calculateJointProbability(firstVector,secondVector,vectorLength); double h = jointRenyiEntropy(state,alpha); freeJointProbabilityState(state); return h; }/*calcJointRenyiEntropy(double,uint*,uint*,int)*/ double discAndCalcJointRenyiEntropy(double alpha, double *firstVector, double *secondVector, int vectorLength) { JointProbabilityState state = discAndCalcJointProbability(firstVector,secondVector,vectorLength); double h = jointRenyiEntropy(state,alpha); freeJointProbabilityState(state); return h; }/*discAndCalcJointRenyiEntropy(double,double*,double*,int)*/ MIToolbox-3.0.0/src/RenyiMutualInformation.c000066400000000000000000000075431303427071400210170ustar00rootroot00000000000000/******************************************************************************* ** 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-2017 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/MIToolbox.h" #include "MIToolbox/ArrayOperations.h" #include "MIToolbox/CalculateProbability.h" #include "MIToolbox/RenyiEntropy.h" #include "MIToolbox/RenyiMutualInformation.h" double renyiMI(JointProbabilityState state, double alpha) { int firstIndex,secondIndex; int i; double jointTemp, marginalTemp; double invAlpha = 1.0 - alpha; double mutualInformation = 0.0; /* 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); marginalTemp = state.firstProbabilityVector[firstIndex] * state.secondProbabilityVector[secondIndex]; marginalTemp = pow(marginalTemp,invAlpha); mutualInformation += (jointTemp * marginalTemp); } } mutualInformation = log(mutualInformation); mutualInformation /= log(LOG_BASE); mutualInformation /= (alpha-1.0); return mutualInformation; } double calcRenyiMIDivergence(double alpha, uint *dataVector, uint *targetVector, int vectorLength) { JointProbabilityState state = calculateJointProbability(dataVector,targetVector,vectorLength); double mutualInformation = renyiMI(state,alpha); freeJointProbabilityState(state); return mutualInformation; }/*calcRenyiMIDivergence(double, uint *, uint *, int)*/ double discAndCalcRenyiMIDivergence(double alpha, double *dataVector, double *targetVector, int vectorLength) { JointProbabilityState state = discAndCalcJointProbability(dataVector,targetVector,vectorLength); double mutualInformation = renyiMI(state,alpha); freeJointProbabilityState(state); return mutualInformation; }/*discAndCalcRenyiMIDivergence(double, double *, double *, int)*/ double calcRenyiMIJoint(double alpha, uint *dataVector, uint *targetVector, int vectorLength) { double hY = calcRenyiEntropy(alpha, targetVector, vectorLength); double hX = calcRenyiEntropy(alpha, dataVector, vectorLength); double hXY = calcJointRenyiEntropy(alpha, dataVector, targetVector, vectorLength); double answer = hX + hY - hXY; return answer; }/*calcRenyiMIJoint(double, uint*, uint*, int)*/ double discAndCalcRenyiMIJoint(double alpha, double *dataVector, double *targetVector, int vectorLength) { double mi; uint *dataNormVector = (uint *) checkedCalloc(vectorLength, sizeof(uint)); uint *targetNormVector = (uint *) checkedCalloc(vectorLength, sizeof(uint)); normaliseArray(dataVector,dataNormVector,vectorLength); normaliseArray(targetVector,targetNormVector,vectorLength); mi = calcRenyiMIJoint(alpha,dataNormVector,targetNormVector,vectorLength); FREE_FUNC(dataNormVector); FREE_FUNC(targetNormVector); dataNormVector = NULL; targetNormVector = NULL; return mi; }/*discAndCalcRenyiMIJoint(double, double*, double*, int)*/ MIToolbox-3.0.0/src/WeightedEntropy.c000066400000000000000000000111471303427071400174470ustar00rootroot00000000000000/******************************************************************************* ** 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-2017 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/MIToolbox.h" #include "MIToolbox/CalculateProbability.h" #include "MIToolbox/WeightedEntropy.h" double wEntropy(WeightedProbState state) { double entropy = 0.0; double tempValue = 0.0; int i; /*H_w(X) = - \sum_x w(x)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); return entropy; } double calcWeightedEntropy(uint *dataVector, double *weightVector, int vectorLength) { WeightedProbState state = calculateWeightedProbability(dataVector,weightVector,vectorLength); double entropy = wEntropy(state); freeWeightedProbState(state); return entropy; }/*calcWeightedEntropy(uint *,double *,int)*/ double discAndCalcWeightedEntropy(double *dataVector, double *weightVector, int vectorLength) { WeightedProbState state = discAndCalcWeightedProbability(dataVector,weightVector,vectorLength); double entropy = wEntropy(state); freeWeightedProbState(state); return entropy; }/*discAndCalcWeightedEntropy(double *,double *,int)*/ double wJointEntropy(WeightedJointProbState state) { double jointEntropy = 0.0; double tempValue = 0.0; int i; /*H_w(X,Y) = - \sum_x \sum_y w(x,y) p(x,y) \log p(x,y)*/ 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); return jointEntropy; } double calcWeightedJointEntropy(uint *firstVector, uint *secondVector, double *weightVector, int vectorLength) { WeightedJointProbState state = calculateWeightedJointProbability(firstVector,secondVector,weightVector,vectorLength); double jointEntropy = wJointEntropy(state); freeWeightedJointProbState(state); return jointEntropy; }/*calcWeightedJointEntropy(uint *,uint *,double *,int)*/ double discAndCalcWeightedJointEntropy(double *firstVector, double *secondVector, double *weightVector, int vectorLength) { WeightedJointProbState state = discAndCalcWeightedJointProbability(firstVector,secondVector,weightVector,vectorLength); double jointEntropy = wJointEntropy(state); freeWeightedJointProbState(state); return jointEntropy; }/*discAndCalcWeightedJointEntropy(double *,double *,double *,int)*/ double wCondEntropy(WeightedJointProbState state) { double condEntropy = 0.0; double jointValue = 0.0; double condValue = 0.0; int i; /*H_w(X|Y) = - \sum_x \sum_y w(x,y)p(x,y) \log (p(x,y)/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); return condEntropy; } double calcWeightedConditionalEntropy(uint *dataVector, uint *conditionVector, double *weightVector, int vectorLength) { WeightedJointProbState state = calculateWeightedJointProbability(dataVector,conditionVector,weightVector,vectorLength); double condEntropy = wCondEntropy(state); freeWeightedJointProbState(state); return condEntropy; }/*calcWeightedConditionalEntropy(uint *,uint *,double *,int)*/ double discAndCalcWeightedConditionalEntropy(double *dataVector, double *conditionVector, double *weightVector, int vectorLength) { WeightedJointProbState state = discAndCalcWeightedJointProbability(dataVector,conditionVector,weightVector,vectorLength); double condEntropy = wCondEntropy(state); freeWeightedJointProbState(state); return condEntropy; }/*discAndCalcWeightedConditionalEntropy(double *,double *,double *,int)*/ MIToolbox-3.0.0/src/WeightedMutualInformation.c000066400000000000000000000121631303427071400214630ustar00rootroot00000000000000/******************************************************************************* ** 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-2017 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/MIToolbox.h" #include "MIToolbox/ArrayOperations.h" #include "MIToolbox/CalculateProbability.h" #include "MIToolbox/WeightedEntropy.h" #include "MIToolbox/WeightedMutualInformation.h" double wmi(WeightedJointProbState state) { double mutualInformation = 0.0; int firstIndex,secondIndex; int i; /* ** I_w(X;Y) = \sum_x \sum_y w(x,y)p(x,y) * \log (p(x,y)/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)) { mutualInformation += state.jointWeightVector[i] * state.jointProbabilityVector[i] * log(state.jointProbabilityVector[i] / state.firstProbabilityVector[firstIndex] / state.secondProbabilityVector[secondIndex]); } } mutualInformation /= log(LOG_BASE); return mutualInformation; } double calcWeightedMutualInformation(uint *dataVector, uint *targetVector, double *weightVector, int vectorLength) { WeightedJointProbState state = calculateWeightedJointProbability(dataVector,targetVector,weightVector,vectorLength); double mutualInformation = wmi(state); freeWeightedJointProbState(state); return mutualInformation; }/*calcWeightedMutualInformation(uint *,uint *,double *,int)*/ double discAndCalcWeightedMutualInformation(double *dataVector, double *targetVector, double *weightVector, int vectorLength) { WeightedJointProbState state = discAndCalcWeightedJointProbability(dataVector,targetVector,weightVector,vectorLength); double mutualInformation = wmi(state); freeWeightedJointProbState(state); return mutualInformation; }/*discAndCalcWeightedMutualInformation(double *,double *,double *,int)*/ double calcWeightedConditionalMutualInformation(uint *dataVector, uint *targetVector, uint *conditionVector, double *weightVector, int vectorLength) { double mutualInformation = 0.0; double firstCondition, secondCondition; uint *mergedVector = (uint *) checkedCalloc(vectorLength,sizeof(uint)); 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 = calcWeightedConditionalEntropy(dataVector,conditionVector,weightVector,vectorLength); secondCondition = calcWeightedConditionalEntropy(dataVector,mergedVector,weightVector,vectorLength); mutualInformation = firstCondition - secondCondition; FREE_FUNC(mergedVector); mergedVector = NULL; return mutualInformation; }/*calcWeightedConditionalMutualInformation(double *,double *,double *,double *,int)*/ double discAndCalcWeightedConditionalMutualInformation(double *dataVector, double *targetVector, double *conditionVector, double *weightVector, int vectorLength) { double mutualInformation = 0.0; double firstCondition, secondCondition; uint *dataNormVector = (uint *) checkedCalloc(vectorLength,sizeof(uint)); uint *targetNormVector = (uint *) checkedCalloc(vectorLength,sizeof(uint)); uint *conditionNormVector = (uint *) checkedCalloc(vectorLength,sizeof(uint)); uint *mergedVector = (uint *) checkedCalloc(vectorLength,sizeof(uint)); normaliseArray(dataVector,dataNormVector,vectorLength); normaliseArray(targetVector,targetNormVector,vectorLength); normaliseArray(conditionVector,conditionNormVector,vectorLength); mergeArrays(targetNormVector,conditionNormVector,mergedVector,vectorLength); /* I_w(X;Y|Z) = H_w(X|Z) - H_w(X|YZ) */ /* double calculateWeightedConditionalEntropy(double *dataVector, double *conditionVector, double *weightVector, int vectorLength); */ firstCondition = calcWeightedConditionalEntropy(dataNormVector,conditionNormVector,weightVector,vectorLength); secondCondition = calcWeightedConditionalEntropy(dataNormVector,mergedVector,weightVector,vectorLength); mutualInformation = firstCondition - secondCondition; FREE_FUNC(dataNormVector); FREE_FUNC(targetNormVector); FREE_FUNC(conditionNormVector); FREE_FUNC(mergedVector); dataNormVector = NULL; targetNormVector = NULL; conditionNormVector = NULL; mergedVector = NULL; return mutualInformation; }/*discAndCalcWeightedConditionalMutualInformation(double *,double *,double *,double *,int)*/ MIToolbox-3.0.0/test/000077500000000000000000000000001303427071400143465ustar00rootroot00000000000000MIToolbox-3.0.0/test/testMIToolbox.c000077500000000000000000000057451303427071400173040ustar00rootroot00000000000000#include #include #include #include "MIToolbox/ArrayOperations.h" #include "MIToolbox/Entropy.h" #include "MIToolbox/MutualInformation.h" int main(int argc, char *argv[]) { int i; double length, miTarget, entropyTarget, cmiTarget; double firstEntropy, secondEntropy, thirdEntropy, targetEntropy; double firstMItarget, secondMItarget, thirdMItarget, targetMItarget; int *testFirstVector, *testSecondVector, *testThirdVector, *testMergedVector; struct timeval start,end; int *firstVector = (int *) calloc(4,sizeof(int)); int *secondVector = (int *) calloc(4,sizeof(int)); int *thirdVector = (int *) calloc(4,sizeof(int)); int *targetVector = (int *) calloc(4,sizeof(int)); firstVector[0] = 0; firstVector[1] = 0; firstVector[2] = 1; firstVector[3] = 1; secondVector[0] = 0; secondVector[1] = 1; secondVector[2] = 0; secondVector[3] = 1; thirdVector[0] = 0; thirdVector[1] = 1; thirdVector[2] = 1; thirdVector[3] = 1; targetVector[0] = 0; targetVector[1] = 1; targetVector[2] = 1; targetVector[3] = 0; firstEntropy = calcEntropy(firstVector,4); secondEntropy = calcEntropy(secondVector,4); thirdEntropy = calcEntropy(thirdVector,4); targetEntropy = calcEntropy(targetVector,4); printf("Entropies - first: %f, second: %f, third: %f, target %f\n",firstEntropy,secondEntropy,thirdEntropy,targetEntropy); firstMItarget = calcMutualInformation(firstVector,targetVector,4); secondMItarget = calcMutualInformation(secondVector,targetVector,4); thirdMItarget = calcMutualInformation(thirdVector,targetVector,4); targetMItarget = calcMutualInformation(targetVector,targetVector,4); printf("MIs - first: %f, second: %f, third: %f, target %f\n",firstMItarget,secondMItarget,thirdMItarget,targetMItarget); testFirstVector = (int *) calloc(10000,sizeof(int)); testSecondVector = (int *) calloc(10000,sizeof(int)); testThirdVector = (int *) calloc(10000,sizeof(int)); testMergedVector = (int *) calloc(10000,sizeof(int)); for (i = 0; i < 10000; i++) { testFirstVector[i] = i % 2; testSecondVector[i] = i % 4; testThirdVector[i] = i % 3; } /* struct timeval * { * time_t tv_sec seconds * suseconds_t tv_usec microseconds * } */ gettimeofday(&start, NULL); for (i = 0; i < 1000; i++) { miTarget = calcMutualInformation(testFirstVector,testSecondVector,10000); entropyTarget = calcEntropy(testFirstVector,10000); cmiTarget = calcConditionalMutualInformation(testFirstVector,testSecondVector,testThirdVector,10000); mergeArrays(testFirstVector,testSecondVector,testMergedVector,10000); } gettimeofday(&end, NULL); printf("I(X;Y) = %f, H(X) = %f, I(X;Y|Z) = %f\n",miTarget,entropyTarget,cmiTarget); length = end.tv_sec - start.tv_sec; length = length + (end.tv_usec - start.tv_usec) / 1000000.0; printf("Time taken for a thousand I(X;Y), H(X), I(X;Y|Z), merge(X,Y) is %lf seconds\n",length); }/*main(int, char **)*/