From 734b2ac44a9f681404178fc4bce1c78aee88aba3 Mon Sep 17 00:00:00 2001 From: "joseph.lizier" Date: Thu, 1 Nov 2012 12:09:02 +0000 Subject: [PATCH] CA demos: added demos for entropy, entropy rate and excess entropy (even though EE calculator is not uploaded to the toolkit yet). Added capability for a moving frame of reference to the toolkit. Sped up javaMatrixToOctave conversion for a sub-matrix. --- build.xml | 1 + .../plotLocalInfoMeasureForCA.m | 177 +++++++++++++++++- demos/octave/javaMatrixToOctave.m | 38 ++-- 3 files changed, 190 insertions(+), 26 deletions(-) diff --git a/build.xml b/build.xml index 731ab51..27e529e 100755 --- a/build.xml +++ b/build.xml @@ -146,6 +146,7 @@ + diff --git a/demos/octave/CellularAutomata/plotLocalInfoMeasureForCA.m b/demos/octave/CellularAutomata/plotLocalInfoMeasureForCA.m index 68eb4fc..26aa2b4 100755 --- a/demos/octave/CellularAutomata/plotLocalInfoMeasureForCA.m +++ b/demos/octave/CellularAutomata/plotLocalInfoMeasureForCA.m @@ -17,9 +17,13 @@ % - "transfer", 1 - apparent transfer entropy (requires measureParams.k and j) % - "transfercomplete", 2 - complete transfer entropy (requires measureParams.k and j) % - "separable", 3 - separable information (requires measureParams.k) +% - "entropy", 4 - excess entropy (requires no params) +% - "entropyrate", 5 - excess entropy (requires measureParams.k) +% - "excess", 6 - excess entropy (requires measureParams.k) % - "all", -1 - plot all measures % - measureParams - a structure containing options as described for each measure above: % - measureParams.k - history length for information dynamics measures +% (for excess entropy, predictive information formulation, this is the history length and future length) % - measureParams.j - we measure information transfer across j cells to the right per time step % - options - a stucture containing a range of other options, i.e.: % - plotOptions - structure as defined for the plotRawCa function @@ -28,7 +32,7 @@ % We set rand("state", options.seed) if options.seed is supplied, and restore the previous seed afterwards. % - plotRawCa - default true % - saveImages - whether to save the plots or not (default false) -% - movingFrameSpeed - to investigate a moving frame of reference (default 0) (as in Lizier & Mahoney paper) +% - movingFrameSpeed - moving frame of reference's cells/time step, as in Lizier & Mahoney paper (default 0) function plotLocalInfoMeasureForCA(neighbourhood, base, rule, cells, timeSteps, measureId, measureParams, options) @@ -67,6 +71,22 @@ function plotLocalInfoMeasureForCA(neighbourhood, base, rule, cells, timeSteps, end figNum = 2; toc + % The offsets of the parents (see runCA for how this is computed, especially for even neighbourhood): + fullSetOfParents = ceil(-neighbourhood / 2) : ceil(-neighbourhood / 2) + (neighbourhood-1); + + if (isfield(options, "movingFrameSpeed")) + % User has requested us to evaluate information dynamics with a moving frame of reference + % (see Lizier and Mahoney paper). + % The shift for each row of the CA is the negative of the movingFrameSpeed, since + % frame moving at 1 cell/time step is same as next row being shifted backwards by 1 cell: + caStates = accumulateShift(caStates, -options.movingFrameSpeed); + % Also account for the moved frame of reference in the offsets of parents to the destination: + fullSetOfParents = fullSetOfParents - options.movingFrameSpeed; + if (isfield(measureParams, "j")) + % Also account for the moved frame of reference in the j parameter for transfer across j cells + measureParams.j = measureParams.j - options.movingFrameSpeed; + end + end % convert the states to a format usable by java: caStatesJInts = octaveToJavaIntMatrix(caStates); @@ -74,7 +94,10 @@ function plotLocalInfoMeasureForCA(neighbourhood, base, rule, cells, timeSteps, plottedOne = false; - % Make the local information dynamics measurement + % Make the local information dynamics measurement(s) + + %============================ + % Active information storage if ((ischar(measureId) && (strcmpi("active", measureId) || strcmpi("all", measureId))) || \ ((measureId == 0) || (measureId == -1))) % Compute active information storage @@ -84,15 +107,25 @@ function plotLocalInfoMeasureForCA(neighbourhood, base, rule, cells, timeSteps, avActive = activeCalc.computeAverageLocalOfObservations(); printf("Average active information storage = %.4f\n", avActive); javaLocalValues = activeCalc.computeLocalFromPreviousObservations(caStatesJInts); + localValues = javaMatrixToOctave(javaLocalValues); + if (isfield(options, "movingFrameSpeed")) + % User has requested us to evaluate information dynamics with a moving frame of reference + % (see Lizier and Mahoney paper). + % Need to shift the computed info dynamics back (to compensate for earlier shift to CA states: + localValues = accumulateShift(localValues, options.movingFrameSpeed); + end toc figure(figNum) figNum = figNum + 1; - plotLocalInfoValues(javaLocalValues, options.plotOptions); + plotLocalInfoValues(localValues, options.plotOptions); if (options.saveImages) print("figures/active.eps", "-color", "-deps"); end plottedOne = true; end + + %============================ + % Apparent transfer entropy if ((ischar(measureId) && (strcmpi("transfer", measureId) || strcmpi("all", measureId))) || \ ((measureId == 1) || (measureId == -1))) % Compute apparent transfer entropy @@ -105,15 +138,25 @@ function plotLocalInfoMeasureForCA(neighbourhood, base, rule, cells, timeSteps, avTransfer = transferCalc.computeAverageLocalOfObservations(); printf("Average apparent transfer entropy (j=%d) = %.4f\n", measureParams.j, avTransfer); javaLocalValues = transferCalc.computeLocalFromPreviousObservations(caStatesJInts, measureParams.j); + localValues = javaMatrixToOctave(javaLocalValues); + if (isfield(options, "movingFrameSpeed")) + % User has requested us to evaluate information dynamics with a moving frame of reference + % (see Lizier and Mahoney paper). + % Need to shift the computed info dynamics back (to compensate for earlier shift to CA states: + localValues = accumulateShift(localValues, options.movingFrameSpeed); + end toc figure(figNum) figNum = figNum + 1; - plotLocalInfoValues(javaLocalValues, options.plotOptions); + plotLocalInfoValues(localValues, options.plotOptions); if (options.saveImages) print(sprintf("figures/transfer-%d.eps", measureParams.j), "-color", "-deps"); end plottedOne = true; end + + %============================ + % Complete transfer entropy if ((ischar(measureId) && (strcmpi("transfercomplete", measureId) || strcmpi("completetransfer", measureId) || strcmpi("all", measureId))) || \ ((measureId == 2) || (measureId == -1))) % Compute complete transfer entropy @@ -123,47 +166,149 @@ function plotLocalInfoMeasureForCA(neighbourhood, base, rule, cells, timeSteps, transferCalc = javaObject('infodynamics.measures.discrete.CompleteTransferEntropyCalculator', ... base, measureParams.k, neighbourhood - 2); transferCalc.initialise(); - % The offsets of the parents (see runCA for how this is computed, especially for even neighbourhood): - fullSetOfParents = ceil(-neighbourhood / 2) : ceil(-neighbourhood / 2) + (neighbourhood-1); % Offsets of all parents can be included here - even 0 and j, these will be eliminated internally: transferCalc.addObservations(caStatesJInts, measureParams.j, octaveToJavaIntArray(fullSetOfParents)); avTransfer = transferCalc.computeAverageLocalOfObservations(); printf("Average complete transfer entropy (j=%d) = %.4f\n", measureParams.j, avTransfer); javaLocalValues = transferCalc.computeLocalFromPreviousObservations(caStatesJInts, ... measureParams.j, octaveToJavaIntArray(fullSetOfParents)); + localValues = javaMatrixToOctave(javaLocalValues); + if (isfield(options, "movingFrameSpeed")) + % User has requested us to evaluate information dynamics with a moving frame of reference + % (see Lizier and Mahoney paper). + % Need to shift the computed info dynamics back (to compensate for earlier shift to CA states: + localValues = accumulateShift(localValues, options.movingFrameSpeed); + end toc figure(figNum) figNum = figNum + 1; - plotLocalInfoValues(javaLocalValues, options.plotOptions); + plotLocalInfoValues(localValues, options.plotOptions); if (options.saveImages) print(sprintf("figures/transferComp-%d.eps", measureParams.j), "-color", "-deps"); end plottedOne = true; end + + %============================ + % Separable information if ((ischar(measureId) && (strcmpi("separable", measureId) || strcmpi("all", measureId))) || \ ((measureId == 3) || (measureId == -1))) % Compute separable information separableCalc = javaObject('infodynamics.measures.discrete.SeparableInfoCalculator', ... base, measureParams.k, neighbourhood - 1); separableCalc.initialise(); - % The offsets of the parents (see runCA for how this is computed, especially for even neighbourhood): - fullSetOfParents = ceil(-neighbourhood / 2) : ceil(-neighbourhood / 2) + (neighbourhood-1); % Offsets of all parents can be included here - even 0 and j, these will be eliminated internally: separableCalc.addObservations(caStatesJInts, octaveToJavaIntArray(fullSetOfParents)); avSeparable = separableCalc.computeAverageLocalOfObservations(); printf("Average separable information = %.4f\n", avSeparable); javaLocalValues = separableCalc.computeLocalFromPreviousObservations(caStatesJInts, ... octaveToJavaIntArray(fullSetOfParents)); + localValues = javaMatrixToOctave(javaLocalValues); + if (isfield(options, "movingFrameSpeed")) + % User has requested us to evaluate information dynamics with a moving frame of reference + % (see Lizier and Mahoney paper). + % Need to shift the computed info dynamics back (to compensate for earlier shift to CA states: + localValues = accumulateShift(localValues, options.movingFrameSpeed); + end toc figure(figNum) figNum = figNum + 1; - plotLocalInfoValues(javaLocalValues, options.plotOptions); + plotLocalInfoValues(localValues, options.plotOptions); if (options.saveImages) print(sprintf("figures/separable-k%d.eps", measureParams.k), "-color", "-deps"); end plottedOne = true; end + %============================ + % Entropy + if ((ischar(measureId) && (strcmpi("entropy", measureId) || strcmpi("all", measureId))) || \ + ((measureId == 4) || (measureId == -1))) + % Compute entropy + entropyCalc = javaObject('infodynamics.measures.discrete.EntropyCalculator', ... + base); + entropyCalc.initialise(); + entropyCalc.addObservations(caStatesJInts); + avEntropy = entropyCalc.computeAverageLocalOfObservations(); + printf("Average entropy = %.4f\n", avEntropy); + javaLocalValues = entropyCalc.computeLocalFromPreviousObservations(caStatesJInts); + localValues = javaMatrixToOctave(javaLocalValues); + if (isfield(options, "movingFrameSpeed")) + % User has requested us to evaluate information dynamics with a moving frame of reference + % (see Lizier and Mahoney paper). + % Need to shift the computed info dynamics back (to compensate for earlier shift to CA states): + % (Note for entropy, the shifts to and back don't make any difference) + localValues = accumulateShift(localValues, options.movingFrameSpeed); + end + toc + figure(figNum) + figNum = figNum + 1; + plotLocalInfoValues(localValues, options.plotOptions); + if (options.saveImages) + print("figures/entropy.eps", "-color", "-deps"); + end + plottedOne = true; + end + + %============================ + % Entropy rate + if ((ischar(measureId) && (strcmpi("entropyrate", measureId) || strcmpi("all", measureId))) || \ + ((measureId == 5) || (measureId == -1))) + % Compute entropy rate + entRateCalc = javaObject('infodynamics.measures.discrete.EntropyRateCalculator', base, measureParams.k); + entRateCalc.initialise(); + entRateCalc.addObservations(caStatesJInts); + avEntRate = entRateCalc.computeAverageLocalOfObservations(); + printf("Average entropy rate = %.4f\n", avEntRate); + javaLocalValues = entRateCalc.computeLocalFromPreviousObservations(caStatesJInts); + localValues = javaMatrixToOctave(javaLocalValues); + if (isfield(options, "movingFrameSpeed")) + % User has requested us to evaluate information dynamics with a moving frame of reference + % (see Lizier and Mahoney paper). + % Need to shift the computed info dynamics back (to compensate for earlier shift to CA states: + localValues = accumulateShift(localValues, options.movingFrameSpeed); + end + toc + figure(figNum) + figNum = figNum + 1; + plotLocalInfoValues(localValues, options.plotOptions); + if (options.saveImages) + print(sprintf("figures/entrate-k%d.eps", measureParams.k), "-color", "-deps"); + end + plottedOne = true; + end + + %============================ + % Excess entropy + if ((ischar(measureId) && (strcmpi("excess", measureId) || strcmpi("all", measureId))) || \ + ((measureId == 6) || (measureId == -1))) + % Compute excess entropy + + error('infodynamics.measures.discrete.ExcessEntropyCalculator is not yet implemented'); + + excessEntropyCalc = javaObject('infodynamics.measures.discrete.ExcessEntropyCalculator', base, measureParams.k); + excessEntropyCalc.initialise(); + excessEntropyCalc.addObservations(caStatesJInts); + avExcessEnt = excessEntropyCalc.computeAverageLocalOfObservations(); + printf("Average excess entropy = %.4f\n", avExcessEnt); + javaLocalValues = excessEntropyCalc.computeLocalFromPreviousObservations(caStatesJInts); + localValues = javaMatrixToOctave(javaLocalValues); + if (isfield(options, "movingFrameSpeed")) + % User has requested us to evaluate information dynamics with a moving frame of reference + % (see Lizier and Mahoney paper). + % Need to shift the computed info dynamics back (to compensate for earlier shift to CA states: + localValues = accumulateShift(localValues, options.movingFrameSpeed); + end + toc + figure(figNum) + figNum = figNum + 1; + plotLocalInfoValues(localValues, options.plotOptions); + if (options.saveImages) + print(sprintf("figures/excessentropy-k%d.eps", measureParams.k), "-color", "-deps"); + end + plottedOne = true; + end + if (not(plottedOne)) error(sprintf("Supplied measureId %s did not match any measurement types", measureId)); end @@ -171,3 +316,15 @@ function plotLocalInfoMeasureForCA(neighbourhood, base, rule, cells, timeSteps, toc end +% Perform a circular shift of each row of the matrix, shifting the first row by zero, +% the second by shiftPerRow, the third by 2*shiftPerRow, and so on +function shiftedMatrix = accumulateShift(matrix, shiftPerRow) + % Allocate required space up front (makes this much faster): + shiftedMatrix = zeros(rows(matrix), columns(matrix)); + shiftedMatrix(1,:) = matrix(1,:); + for r = 2 : rows(matrix) + % circshift operates on shifting rows, so we transpose the input and output to it: + shiftedMatrix(r,:) = circshift(matrix(r,:)', shiftPerRow*(r-1))'; + end +end + diff --git a/demos/octave/javaMatrixToOctave.m b/demos/octave/javaMatrixToOctave.m index 5a1457f..498186a 100755 --- a/demos/octave/javaMatrixToOctave.m +++ b/demos/octave/javaMatrixToOctave.m @@ -2,9 +2,9 @@ % % Convert a java matrix (1 or 2D, double or int - but not Integer!!) to an octave matrix % -% Unfortunately octave-java doesn't seem to handle the conversion properly, +% Octave-java doesn't seem to handle the conversion natively, +% so we either use org.octave.Matrix (built-in) to do it, or % and we must convert each individual array item ourselves (This can be very slow for large matrices) -% or with org.octave.Matrix (built-in). % function octaveMatrix = javaMatrixToOctave(javaMatrix, startRow, startCol, numRows, numCols) @@ -25,25 +25,31 @@ function octaveMatrix = javaMatrixToOctave(javaMatrix, startRow, startCol, numRo if (exist ('OCTAVE_VERSION', 'builtin')) % We're in octave: % Using 'org.octave.Matrix' is much faster than conversion cell by cell - % (only use if we're converting the whole array though, too many issues - % with type checking etc to bother otherwise) - if (nargin < 2) - tmp = javaObject('org.octave.Matrix', javaMatrix); - % Make sure tmp.ident() is converted to native octave: - oldFlag = java_convert_matrix (1); - unwind_protect - octaveMatrix = tmp.ident(tmp); - unwind_protect_cleanup - % restore to non-default conversion, otherwise we get - % bad errors on other calls - java_convert_matrix(oldFlag); - end_unwind_protect + % Convert whole matrix first: + tmp = javaObject('org.octave.Matrix', javaMatrix); + % Make sure tmp.ident() is converted to native octave: + oldFlag = java_convert_matrix (1); + converted = false; + unwind_protect + octaveMatrix = tmp.ident(tmp); + converted = true; + unwind_protect_cleanup + % restore to non-default conversion, otherwise we get + % bad errors on other calls + java_convert_matrix(oldFlag); + end_unwind_protect + if (converted) + if (nargin >= 2) + % Do some resizing: + octaveMatrix = octaveMatrix(startRow:startRow+numRows-1, startCol:startCol+numCols-1); + end return; end + % else fall through to cell by cell conversion, as per for matlab end - % Else, either we couldn't be bothered doing the resizing for java, + % Else, either we encountered an error in the octave resizing, % or we were in Matlab all along. (If there's a fast way for matlab, tell me) octaveMatrix = zeros(numRows, numCols);