From 56b3db2c89afad280d5c6d5399236e245d4ea9bc Mon Sep 17 00:00:00 2001
From: Tristan Stenner <ttstenner@gmail.com>
Date: Mon, 17 May 2021 12:39:58 +0200
Subject: [PATCH 1/3] Re-implement stat_surrogate_pvals with an O(n) algorithm
 and support for left-tailed tests

---
 functions/statistics/stat_surrogate_pvals.m | 42 +++++++--------------
 1 file changed, 14 insertions(+), 28 deletions(-)

diff --git a/functions/statistics/stat_surrogate_pvals.m b/functions/statistics/stat_surrogate_pvals.m
index 98ae01dfd..eebff491d 100644
--- a/functions/statistics/stat_surrogate_pvals.m
+++ b/functions/statistics/stat_surrogate_pvals.m
@@ -52,36 +52,22 @@
 % THE POSSIBILITY OF SUCH DAMAGE.
 
 function pvals = stat_surrogate_pvals(distribution,observed,tail)
-
-numDims = myndims(distribution);
-
-% append observed to last dimension of surrogate distribution
-distribution = cat(numDims,distribution,observed);
-
-numDims = myndims(distribution);
-
-% sort along last dimension (replications)
-[tmp idx] = sort( distribution, numDims,'ascend');
-[tmp mx]  = max( idx,[], numDims);
-
-len = size(distribution,  numDims );
-pvals = 1-(mx-0.5)/len;
-if strcmpi(tail, 'both')
-    pvals = min(pvals, 1-pvals);
-    pvals = 2*pvals;
+numDims = ndims(distribution);
+if iscolumn(distribution)
+	numDims = 1;
 end
 
+n = size(distribution, numDims);
+pvals = sum(distribution >= observed, numDims) / n;
 
-% get the number of dimensions in a matrix
-function val = myndims(a)
-if ndims(a) > 2
-    val = ndims(a);
+if strcmpi(tail, 'both')
+	pvals = 2 * min(pvals, 1 - pvals);
+elseif strcmpi(tail, 'left')
+	pvals = 1 - pvals;
+elseif any(strcmpi(tail, {'right', 'one'}))
+	% nothing to be done
 else
-    if size(a,1) == 1,
-        val = 2;
-    elseif size(a,2) == 1,
-        val = 1;
-    else
-        val = 2;
-    end
+	error('invalid value for tail: "%s", should be left, right, one or both', tail);
 end
+
+end
\ No newline at end of file

From 90eec096d090888fc758585f14899b345fabca94 Mon Sep 17 00:00:00 2001
From: Tristan Stenner <ttstenner@gmail.com>
Date: Mon, 17 May 2021 11:47:56 +0200
Subject: [PATCH 2/3] Compatibility with Matlab < R2016b

---
 functions/statistics/stat_surrogate_pvals.m | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/functions/statistics/stat_surrogate_pvals.m b/functions/statistics/stat_surrogate_pvals.m
index eebff491d..0c051b3df 100644
--- a/functions/statistics/stat_surrogate_pvals.m
+++ b/functions/statistics/stat_surrogate_pvals.m
@@ -58,7 +58,9 @@
 end
 
 n = size(distribution, numDims);
-pvals = sum(distribution >= observed, numDims) / n;
+% once support for matlab <= R2016b is dropped:
+% pvals = sum(distribution >= observed, numDims) / n;
+pvals = sum(bsxfun(@ge, distribution, observed), numDims) / n;
 
 if strcmpi(tail, 'both')
 	pvals = 2 * min(pvals, 1 - pvals);

From 91e761a531f8cfc863708f5d370dd054de2ac01f Mon Sep 17 00:00:00 2001
From: Tristan Stenner <ttstenner@gmail.com>
Date: Mon, 17 May 2021 12:36:15 +0200
Subject: [PATCH 3/3] Handle observed==distribution for left- and two-tailed
 tests

---
 functions/statistics/stat_surrogate_pvals.m | 13 +++++++++----
 1 file changed, 9 insertions(+), 4 deletions(-)

diff --git a/functions/statistics/stat_surrogate_pvals.m b/functions/statistics/stat_surrogate_pvals.m
index 0c051b3df..158800968 100644
--- a/functions/statistics/stat_surrogate_pvals.m
+++ b/functions/statistics/stat_surrogate_pvals.m
@@ -62,12 +62,17 @@
 % pvals = sum(distribution >= observed, numDims) / n;
 pvals = sum(bsxfun(@ge, distribution, observed), numDims) / n;
 
+if any(strcmpi(tail, {'right', 'one'}))
+	% nothing to be done
+	return;
+end
+
+p_left = 1 - pvals + sum(bsxfun(@eq, distribution, observed), numDims) / n;
+
 if strcmpi(tail, 'both')
-	pvals = 2 * min(pvals, 1 - pvals);
+	pvals = 2 * min(pvals, p_left);
 elseif strcmpi(tail, 'left')
-	pvals = 1 - pvals;
-elseif any(strcmpi(tail, {'right', 'one'}))
-	% nothing to be done
+	pvals = p_left;
 else
 	error('invalid value for tail: "%s", should be left, right, one or both', tail);
 end