From fc8986d0aba89a68576bd06c9f4e8d2f05b57ab1 Mon Sep 17 00:00:00 2001 From: "joseph.lizier" Date: Sat, 12 Jan 2013 12:38:19 +0000 Subject: [PATCH] Made the DetectingInteractionLags demos compatible with Matlab --- .../coupledLogisticMap.m | 58 +++++----- .../transferWithSourceMemory.m | 109 ++++++++++++------ 2 files changed, 100 insertions(+), 67 deletions(-) diff --git a/demos/octave/DetectingInteractionLags/coupledLogisticMap.m b/demos/octave/DetectingInteractionLags/coupledLogisticMap.m index 404cacb..5098ebd 100755 --- a/demos/octave/DetectingInteractionLags/coupledLogisticMap.m +++ b/demos/octave/DetectingInteractionLags/coupledLogisticMap.m @@ -51,7 +51,7 @@ function coupledLogisticMap() cYToX = 0.2; cXToY = 0.5; T = 512; - printf('For 1000 repeats, expect the calculations to take ~5 minutes ...\n'); + fprintf('For 1000 repeats, expect the calculations to take ~5 minutes ...\n'); repeats = 1000; % General results visible for 100 repeats if you want to see them faster (~20 sec) k = 1; % history length @@ -66,7 +66,7 @@ function coupledLogisticMap() % The actual values of the results slightly change with the Kraskov K parameter, % but we note that delay 2 is higher for K=1,2,3,4,10... - KraskovK ="4"; % Use Kraskov parameter K=4 for 4 nearest points + KraskovK ='4'; % Use Kraskov parameter K=4 for 4 nearest points % For K=2 we get: %TE(X->Y,delay=1) = 0.7773 nats (+/- std 0.0959, stderr 0.0096) %TE(X->Y,delay=2) = 1.8759 nats (+/- std 0.0332, stderr 0.0033) @@ -75,7 +75,7 @@ function coupledLogisticMap() tic; % Add utilities to the path - addpath(".."); + addpath('..'); % Assumes the jar is two levels up - change this if this is not the case % Octave is happy to have the path added multiple times; I'm unsure if this is true for matlab @@ -92,8 +92,8 @@ function coupledLogisticMap() if (ts > 1) % We've just been running transients in X and Y, so copy the % last seedSteps up into the first few rows: - X(1:seedSteps,:) = X(rows(X)-seedSteps+1:rows(X),:); - Y(1:seedSteps,:) = Y(rows(Y)-seedSteps+1:rows(Y),:); + X(1:seedSteps,:) = X(size(X,1)-seedSteps+1:size(X,1),:); + Y(1:seedSteps,:) = Y(size(Y,1)-seedSteps+1:size(Y,1),:); end % Now run the process from these initial states for n = 1 : T @@ -117,24 +117,24 @@ function coupledLogisticMap() teCalc=javaObject('infodynamics.measures.continuous.kraskov.ConditionalMutualInfoCalculatorMultiVariateKraskov2'); % Perform calculation for X -> Y (lag 1) teCalc.initialise(1,1,k); % Use history length k (Schreiber k) - teCalc.setProperty("k", KraskovK); - teCalc.setObservations(octaveToJavaDoubleMatrix(X(seedSteps:rows(X)-1,r)), ... - octaveToJavaDoubleMatrix(Y(seedSteps+1:rows(Y),r)), ... - octaveToJavaDoubleMatrix(Y(seedSteps:rows(Y)-1,r))); + teCalc.setProperty('k', KraskovK); + teCalc.setObservations(octaveToJavaDoubleMatrix(X(seedSteps:size(X,1)-1,r)), ... + octaveToJavaDoubleMatrix(Y(seedSteps+1:size(Y,1),r)), ... + octaveToJavaDoubleMatrix(Y(seedSteps:size(Y,1)-1,r))); resultsLag1(r) = teCalc.computeAverageLocalOfObservations(); % Perform calculation for X -> Y (lag 2) teCalc.initialise(1,1,k); % Use history length k (Schreiber k) - teCalc.setProperty("k", KraskovK); - teCalc.setObservations(octaveToJavaDoubleMatrix(X(seedSteps:rows(X)-2,r)), ... - octaveToJavaDoubleMatrix(Y(seedSteps+2:rows(Y),r)), ... - octaveToJavaDoubleMatrix(Y(seedSteps+1:rows(X)-1,r))); + teCalc.setProperty('k', KraskovK); + teCalc.setObservations(octaveToJavaDoubleMatrix(X(seedSteps:size(X,1)-2,r)), ... + octaveToJavaDoubleMatrix(Y(seedSteps+2:size(Y,1),r)), ... + octaveToJavaDoubleMatrix(Y(seedSteps+1:size(X,1)-1,r))); resultsLag2(r) = teCalc.computeAverageLocalOfObservations(); % Perform calculation for X -> Y (lag 2) teCalc.initialise(1,1,k); % Use history length k (Schreiber k) - teCalc.setProperty("k", KraskovK); - teCalc.setObservations(octaveToJavaDoubleMatrix(X(seedSteps:rows(X)-3,r)), ... - octaveToJavaDoubleMatrix(Y(seedSteps+3:rows(Y),r)), ... - octaveToJavaDoubleMatrix(Y(seedSteps+2:rows(X)-1,r))); + teCalc.setProperty('k', KraskovK); + teCalc.setObservations(octaveToJavaDoubleMatrix(X(seedSteps:size(X,1)-3,r)), ... + octaveToJavaDoubleMatrix(Y(seedSteps+3:size(Y,1),r)), ... + octaveToJavaDoubleMatrix(Y(seedSteps+2:size(X,1)-1,r))); resultsLag3(r) = teCalc.computeAverageLocalOfObservations(); %% TransferEntropyCalculatorKraskov uses two MI calculators, so doesn't eliminate all bias. @@ -142,13 +142,13 @@ function coupledLogisticMap() % teCalc=javaObject('infodynamics.measures.continuous.kraskov.TransferEntropyCalculatorKraskov'); %% Perform calculation for X -> Y (lag 1) %teCalc.initialise(k); % Use history length k (Schreiber k) - %teCalc.setProperty("k", KraskovK); - %teCalc.setObservations(X(seedSteps:rows(X),r), Y(seedSteps:rows(Y),r)); + %teCalc.setProperty('k', KraskovK); + %teCalc.setObservations(X(seedSteps:size(X,1),r), Y(seedSteps:size(Y,1),r)); %resultsLag1(r) = teCalc.computeAverageLocalOfObservations(); %% Perform calculation for X -> Y (lag 2) %teCalc.initialise(k); % Use history length k (Schreiber k) - %teCalc.setProperty("k", KraskovK); - %teCalc.setObservations(X(seedSteps-1:rows(X)-1,r), Y(seedSteps:rows(Y),r)); + %teCalc.setProperty('k', KraskovK); + %teCalc.setObservations(X(seedSteps-1:size(X,1)-1,r), Y(seedSteps:size(Y,1),r)); %resultsLag2(r) = teCalc.computeAverageLocalOfObservations(); % Kernel estimator returns the correct ordering of lag 1 and 2 for @@ -158,23 +158,23 @@ function coupledLogisticMap() % we are grouping together, then we start to see the same effect as we did with % the symbolic coding. %teCalc=javaObject('infodynamics.measures.continuous.kernel.TransferEntropyCalculatorKernel'); - %kernelWidth = "0.25"; % normalised units + %kernelWidth = '0.25'; % normalised units %% Perform calculation for X -> Y (lag 1) %teCalc.initialise(k); % Use history length k (Schreiber k) - %teCalc.setProperty("EPSILON", kernelWidth); - %teCalc.setObservations(X(seedSteps:rows(X),r), Y(seedSteps:rows(Y),r)); + %teCalc.setProperty('EPSILON', kernelWidth); + %teCalc.setObservations(X(seedSteps:size(X,1),r), Y(seedSteps:size(Y,1),r)); %resultsLag1(r) = teCalc.computeAverageLocalOfObservations(); %% Perform calculation for X -> Y (lag 2) %teCalc.initialise(k); % Use history length k (Schreiber k) - %teCalc.setProperty("EPSILON", kernelWidth); - %teCalc.setObservations(X(seedSteps-1:rows(X)-1,r), Y(seedSteps:rows(Y),r)); + %teCalc.setProperty('EPSILON', kernelWidth); + %teCalc.setObservations(X(seedSteps-1:size(X,1)-1,r), Y(seedSteps:size(Y,1),r)); %resultsLag2(r) = teCalc.computeAverageLocalOfObservations(); end - printf("TE(X->Y,delay=1) = %.4f nats (+/- std %.4f, stderr %.4f) or %.4f bits\n", ... + fprintf('TE(X->Y,delay=1) = %.4f nats (+/- std %.4f, stderr %.4f) or %.4f bits\n', ... mean(resultsLag1), std(resultsLag1), std(resultsLag1)./sqrt(repeats), mean(resultsLag1)./log(2) ); - printf("TE(X->Y,delay=2) = %.4f nats (+/- std %.4f, stderr %.4f) or %.4f bits\n", ... + fprintf('TE(X->Y,delay=2) = %.4f nats (+/- std %.4f, stderr %.4f) or %.4f bits\n', ... mean(resultsLag2), std(resultsLag2), std(resultsLag2)./sqrt(repeats), mean(resultsLag2) ./log(2)); - printf("If measured:\nTE(X->Y,delay=3) = %.4f nats (+/- std %.4f, stderr %.4f) or %.4f bits\n", ... + fprintf('If measured:\nTE(X->Y,delay=3) = %.4f nats (+/- std %.4f, stderr %.4f) or %.4f bits\n', ... mean(resultsLag3), std(resultsLag3), std(resultsLag3)./sqrt(repeats), mean(resultsLag3) ./log(2)); toc; diff --git a/demos/octave/DetectingInteractionLags/transferWithSourceMemory.m b/demos/octave/DetectingInteractionLags/transferWithSourceMemory.m index f080859..1ba29a2 100755 --- a/demos/octave/DetectingInteractionLags/transferWithSourceMemory.m +++ b/demos/octave/DetectingInteractionLags/transferWithSourceMemory.m @@ -13,6 +13,9 @@ % showing that it is not maximised % for the correct delay even in such simple unidirectional coupling. % +% NOTE: You may need to increase the Java heap space in Matlab for this +% to work (you will get a "java.lang.OutOfMemoryError: Java heap space" +% error if this is a problem). % % Inputs: % - savePlot - true if you want eps files of the plots saved @@ -22,7 +25,7 @@ function transferWithSourceMemory(savePlot) tic; % Add utilities to the path - addpath(".."); + addpath('..'); % Assumes the jar is two levels up - change this if this is not the case % Octave is happy to have the path added multiple times; I'm unsure if this is true for matlab @@ -93,7 +96,7 @@ function transferWithSourceMemory(savePlot) 0, 1, octaveToJavaIntArray([2])); mitXnminus1ToYnplus1(deltaIndex) = compTeCalc.computeAverageLocalOfObservations(); - printf("delta=%.3f, TE(X_{n} -> Y_{n+1})=%.3f, TE(X_{n-1} -> Y_{n+1})=%.3f, MIT(X_{n} -> Y_{n+1})=%.3f, MIT(X_{n-1} -> Y_{n+1})=%.3f\n", ... + fprintf('delta=%.3f, TE(X_{n} -> Y_{n+1})=%.3f, TE(X_{n-1} -> Y_{n+1})=%.3f, MIT(X_{n} -> Y_{n+1})=%.3f, MIT(X_{n-1} -> Y_{n+1})=%.3f\n', ... delta, teXnToYnplus1(deltaIndex), teXnminus1ToYnplus1(deltaIndex), ... mitXnToYnplus1(deltaIndex), mitXnminus1ToYnplus1(deltaIndex)); end @@ -109,47 +112,77 @@ function transferWithSourceMemory(savePlot) % plot types: h - open diamonds, d - closed diamonds, s - closed squares, p - open squares markersize = 15; figure; - if (savePlot) - set(gca, 'fontsize', 32); % do this first to get fontsize right for the key - end - plot(deltas, teXnToYnplus1, "x1;Empirical: TE_{SPO}(X \rightarrow Y, 1);", "markersize", markersize); - hold on; - plot(deltas, teXnToYnplus1_an, "p1;Analytic: TE_{SPO}(X \rightarrow Y, 1);", "markersize", markersize); % h plots open diamonds - plot(deltas, teXnminus1ToYnplus1, "+2;Empirical: TE_{SPO}(X \rightarrow y, 2);", "markersize", markersize); - plot(deltas, teXnminus1ToYnplus1_an, "o2;Analytic: TE_{SPO}(X \rightarrow Y, 2);", "markersize", markersize); - hold off; - % Make figures ok for plotting: - legend('location', 'east'); - if (savePlot) - xlabel('\eta', 'fontsize', 32); - ylabel('Information (bits)', 'fontsize', 32); - print('te.eps', '-deps', '-color'); + if (exist ('OCTAVE_VERSION', 'builtin')) + % Make the full plots only on octave + if (savePlot) + set(gca, 'fontsize', 32); % do this first to get fontsize right for the key + end + plot(deltas, teXnToYnplus1, 'x1;Empirical: TE_{SPO}(X \rightarrow Y, 1);', 'markersize', markersize); + hold on; + plot(deltas, teXnToYnplus1_an, 'p1;Analytic: TE_{SPO}(X \rightarrow Y, 1);', 'markersize', markersize); % h plots open diamonds + plot(deltas, teXnminus1ToYnplus1, '+2;Empirical: TE_{SPO}(X \rightarrow y, 2);', 'markersize', markersize); + plot(deltas, teXnminus1ToYnplus1_an, 'o2;Analytic: TE_{SPO}(X \rightarrow Y, 2);', 'markersize', markersize); + hold off; + % Make figures ok for plotting: + legend('location', 'east'); + if (savePlot) + xlabel('\eta', 'fontsize', 32); + ylabel('Information (bits)', 'fontsize', 32); + print('te.eps', '-deps', '-color'); + else + % do these without fontsize to have more stable displays in octave + xlabel('\eta'); + ylabel('Information (bits)'); + end else - % do these without fontsize to have more stable displays in octave - xlabel('\eta'); - ylabel('Information (bits)'); + % We're on Matlab - just make a quick plot + plot(deltas, teXnToYnplus1, 'rx', 'markersize', markersize); % Empirical: TE_{SPO}(X \rightarrow Y, 1); + hold on; + plot(deltas, teXnToYnplus1_an, 'rp', 'markersize', markersize); % Analytic: TE_{SPO}(X \rightarrow Y, 1) + plot(deltas, teXnminus1ToYnplus1, 'g+', 'markersize', markersize); % Empirical: TE_{SPO}(X \rightarrow y, 2) + plot(deltas, teXnminus1ToYnplus1_an, 'go', 'markersize', markersize); % Analytic: TE_{SPO}(X \rightarrow Y, 2) + hold off; + % Make figures ok for plotting: + legend('Empirical: TE_{SPO}(X \rightarrow Y, 1)', 'Analytic: TE_{SPO}(X \rightarrow Y, 1)', ... + 'Empirical: TE_{SPO}(X \rightarrow y, 2)', 'Analytic: TE_{SPO}(X \rightarrow Y, 2)', ... + 'Location', 'East'); end figure; - if (savePlot) - set(gca, 'fontsize', 32); % do this first to get fontsize right for the key - end - plot(deltas, mitXnToYnplus1, "x1;Empirical: MIT(X \rightarrow Y, 1);", "markersize", markersize); - hold on; - plot(deltas, mitXnToYnplus1_an, "p1;Analytic: MIT(X \rightarrow Y, 1);", "markersize", markersize); - plot(deltas, mitXnminus1ToYnplus1, "+2;Empirical: MIT(X \rightarrow Y, 2);", "markersize", markersize); - plot(deltas, mitXnminus1ToYnplus1_an, "o2;Analytic: MIT(X \rightarrow Y, 2);", "markersize", markersize); - hold off; - % Make figures ok for plotting: - legend('location', 'east') - if (savePlot) - xlabel('\eta', 'fontsize', 32); - ylabel('Information (bits)', 'fontsize', 32); - print('mit.eps', '-deps', '-color'); + if (exist ('OCTAVE_VERSION', 'builtin')) + % Make the full plots only on octave + if (savePlot) + set(gca, 'fontsize', 32); % do this first to get fontsize right for the key + end + plot(deltas, mitXnToYnplus1, 'x1;Empirical: MIT(X \rightarrow Y, 1);', 'markersize', markersize); + hold on; + plot(deltas, mitXnToYnplus1_an, 'p1;Analytic: MIT(X \rightarrow Y, 1);', 'markersize', markersize); + plot(deltas, mitXnminus1ToYnplus1, '+2;Empirical: MIT(X \rightarrow Y, 2);', 'markersize', markersize); + plot(deltas, mitXnminus1ToYnplus1_an, 'o2;Analytic: MIT(X \rightarrow Y, 2);', 'markersize', markersize); + hold off; + % Make figures ok for plotting: + legend('location', 'east') + if (savePlot) + xlabel('\eta', 'fontsize', 32); + ylabel('Information (bits)', 'fontsize', 32); + print('mit.eps', '-deps', '-color'); + else + % do these without fontsize to have more stable displays in octave + xlabel('\eta'); + ylabel('Information (bits)'); + end else - % do these without fontsize to have more stable displays in octave - xlabel('\eta'); - ylabel('Information (bits)'); + % We're on Matlab - just make a quick plot + plot(deltas, mitXnToYnplus1, 'rx', 'markersize', markersize); % Empirical: MIT(X \rightarrow Y, 1) + hold on; + plot(deltas, mitXnToYnplus1_an, 'rp', 'markersize', markersize); % Analytic: MIT(X \rightarrow Y, 1) + plot(deltas, mitXnminus1ToYnplus1, 'g+', 'markersize', markersize); % Empirical: MIT(X \rightarrow Y, 2) + plot(deltas, mitXnminus1ToYnplus1_an, 'go', 'markersize', markersize); % Analytic: MIT(X \rightarrow Y, 2) + hold off; + % Make figures ok for plotting: + legend('Empirical: MIT(X \rightarrow Y, 1)', 'Analytic: MIT(X \rightarrow Y, 1)', ... + 'Empirical: MIT(X \rightarrow Y, 2)', 'Analytic: MIT(X \rightarrow Y, 2)', ... + 'Location', 'East'); end toc;