From c8d4e2196af8e491f9df1fd0db6487119995ce99 Mon Sep 17 00:00:00 2001
From: Arthur <arn203@exeter.ac.uk>
Date: Fri, 3 Sep 2021 10:05:50 +0100
Subject: [PATCH] test against R

---
 test/runtests.jl | 60 +++++++++++++++++++++++++++++++++++++++---------
 1 file changed, 49 insertions(+), 11 deletions(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index 86c60d0..78f3673 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,14 +1,52 @@
-using BetaDispersion, Distances, StatsBase,BenchmarkTools
+using BetaDispersion, Distances, StatsBase,BenchmarkTools, RCall
 using Test
-
+vegan = "vegan"
+R"install.packages($vegan)"
+R"library(vegan)"
+# Rcall uses a single thread, so make sure Julia is using only ine thread as well
 @testset "BetaDispersion.jl" begin
-    x = rand(1000,50) 
-    y =rand(1:10,1000)
-    d = dispersion(x,y,Euclidean)
-    dispersion(x,y, metric = true)
-    bench = @benchmark dispersion($x,$y, BrayCurtis)
-    @test mean(bench.times) < 2*10^9
-    bench = @benchmark permutest($d)
-    @test mean(bench.times) < 4*10^9
-    
+    #test against vegan on small data set
+    x = rand(30,10) 
+    y =rand(1:3,30)
+    Dj = pairwise(BrayCurtis(),x,dims = 1)
+    Dr =R"vegdist($x)"
+    j = @benchmark dispersion(Dj,y)
+    r = @benchmark R"betadisper($Dr,as.factor($y))"
+    @test mean(j.times) < mean(r.times)
+    dispj = dispersion(Dj,y)
+    dispr = R"betadisper($Dr,as.factor($y))"
+    j = @benchmark permutest(dispj,999)
+    r = @benchmark R"permutest($dispr,pairwise = TRUE, permutations = 999)"
+    @test mean(j.times) < mean(r.times)
+
+    #try with more variables/species
+    x = rand(30,100) 
+    y =rand(1:3,30)
+    Dj = pairwise(BrayCurtis(),x,dims = 1)
+    Dr =R"vegdist($x)"
+    j = @benchmark dispersion(Dj,y)
+    r = @benchmark R"betadisper($Dr,as.factor($y))"
+    @test mean(j.times) < mean(r.times)
+    dispj = dispersion(Dj,y)
+    dispr = R"betadisper($Dr,as.factor($y))"
+    j = @benchmark permutest(dispj,999) 
+    r = @benchmark R"permutest($dispr,pairwise = TRUE, permutations = 999)" 
+    @test mean(j.times) < mean(r.times)
+
+   #now with a larger data set and more groups  
+   x = rand(1000,100) 
+   y =rand(1:20,1000)
+   Dj = pairwise(BrayCurtis(),x,dims = 1) #  Distances.jl returns a distance matrix insantaneously
+   Dr =R"vegdist($x)"   # vegan takes a very long time.
+   j = @benchmark dispersion(Dj,y) #318.505 ms
+   r = @benchmark R"betadisper($Dr,as.factor($y))" #4.477 s
+   @test mean(j.times) < mean(r.times)
+   dispj = dispersion(Dj,y)
+   dispr = R"betadisper($Dr,as.factor($y))"
+   j = @benchmark permutest(dispj,999) #279.505 ms on my machine (single thread)
+   r = @benchmark R"permutest($dispr,pairwise = TRUE, permutations = 999)" #11.515 s
+   @test mean(j.times) < mean(r.times)
+
+
+  
 end