splitx <- split(x,grouping)
splitx
x=Matrix
load("Data_Matrix_Tutorial.RData") #Data already loaded with DNCImper package. Can be call like that: DNCImper::Matrix
library(vegan)
library(plyr)
library(ggplot2)
library(DNCImper)
load("Data_Matrix_Tutorial.RData") #Data already loaded with DNCImper package. Can be call like that: DNCImper::Matrix
DNCImper::Matrix
DNCImper::Group
grouping=Group
x=Matrix
group.combinations <- combn(unique(sort(grouping)),2)
ddelta <- NULL
for(i in 1:NCOL(group.combinations)) {
splitx <- split(x,grouping)
#Ici symmetrize:
if(symmetrize == TRUE)
{
Add <- which(c(NROW(splitx[[group.combinations[1,i]]]),
NROW(splitx[[group.combinations[2,i]]])) == max(c(NROW(splitx[[group.combinations[1,i]]]),
NROW(splitx[[group.combinations[2,i]]]))))
if(Add == 1)
{
sampled_lines <- sample(1:length(splitx[[group.combinations[1,i]]][,1]),
length(splitx[[group.combinations[2,i]]][,1]))
splitx[[group.combinations[1,i]]] <- splitx[[group.combinations[1,i]]][sampled_lines,]
}
if(Add == 2)
{
sampled_lines <- sample(1:length(splitx[[group.combinations[2,i]]][,1]),
length(splitx[[group.combinations[1,i]]][,1]))
splitx[[group.combinations[2,i]]] <- splitx[[group.combinations[2,i]]][sampled_lines,]
}
}
paired.x <- rbind(splitx[[group.combinations[1,i]]],
splitx[[group.combinations[2,i]]])
# remove empty species
ifzero <- which(apply(paired.x, 2, sum) == 0)
if(length(ifzero > 0)){
paired.x <- paired.x[,-which(colSums(paired.x)==0)]}
if(length(which(rowSums(paired.x) == 0)) != 0){stop("ERROR : A row/sample is empty")}
group.pair <- c(rep(group.combinations[1,i], NROW(splitx[[group.combinations[1,i]]])),
rep(group.combinations[2,i], NROW(splitx[[group.combinations[2,i]]])))
ddelta <- rbind(ddelta, DNCImper:::DNCI.ses(x=paired.x,grouping=group.pair,id=id, Nperm = Nperm, count = count, plotSIMPER = plotSIMPER, dataTYPE = dataTYPE)) #here is the part that calculates the index based on PERSIMPER
}
splitx <- split(x,grouping)
splitx
Add <- which(c(NROW(splitx[[group.combinations[1,i]]]),
NROW(splitx[[group.combinations[2,i]]])) == max(c(NROW(splitx[[group.combinations[1,i]]]),
NROW(splitx[[group.combinations[2,i]]]))))
splitx[[group.combinations[1,i]]][,1]
i
group.combinations[1,i]
group.combinations
(splitx[[group.combinations[1,i]]]
)
splitx[[group.combinations[1,i]]]
splitx[[1]]
splitx
splitx[[group.combinations[1,i]]][,1]
grouping=groups
grouping
x=comm_df_for_DNCI
x
group.combinations <- combn(unique(sort(grouping)),2)
group.combinations
splitx <- split(x,grouping)
splitx
Matrix
Group
?split()
x
grouping
split(x,grouping)
comm_df_for_DNCI<-read.csv("~/.julia/dev/MetaCommunityMetrics/benchmarks/benchmark_r/data/DNCI_comm.csv")
comm_df_for_DNCI
split(x,grouping)
grouping
x
grouping
split(x,grouping)
x=Matrix
grouping=Group
x
grouping
split(x,grouping)
x=comm_df_for_DNCI
x
x=as.matrix(comm_df_for_DNCI)
x
grouping
grouping=groups
grouping
split(x,grouping)
str(Matrix)
str(comm_df_for_DNCI)
split(comm_df_for_DNCI, groups)
x=comm_df_for_DNCI
str(X)
str(x)
grouping=groups
split(x,grouping)
group.combinations <- combn(unique(sort(grouping)),2)
group.combinations
ddelta <- NULL
for(i in 1:NCOL(group.combinations)) {
splitx <- split(x,grouping)
#Ici symmetrize:
if(symmetrize == TRUE)
{
Add <- which(c(NROW(splitx[[group.combinations[1,i]]]),
NROW(splitx[[group.combinations[2,i]]])) == max(c(NROW(splitx[[group.combinations[1,i]]]),
NROW(splitx[[group.combinations[2,i]]]))))
if(Add == 1)
{
sampled_lines <- sample(1:length(splitx[[group.combinations[1,i]]][,1]),
length(splitx[[group.combinations[2,i]]][,1]))
splitx[[group.combinations[1,i]]] <- splitx[[group.combinations[1,i]]][sampled_lines,]
}
if(Add == 2)
{
sampled_lines <- sample(1:length(splitx[[group.combinations[2,i]]][,1]),
length(splitx[[group.combinations[1,i]]][,1]))
splitx[[group.combinations[2,i]]] <- splitx[[group.combinations[2,i]]][sampled_lines,]
}
}
paired.x <- rbind(splitx[[group.combinations[1,i]]],
splitx[[group.combinations[2,i]]])
# remove empty species
ifzero <- which(apply(paired.x, 2, sum) == 0)
if(length(ifzero > 0)){
paired.x <- paired.x[,-which(colSums(paired.x)==0)]}
if(length(which(rowSums(paired.x) == 0)) != 0){stop("ERROR : A row/sample is empty")}
group.pair <- c(rep(group.combinations[1,i], NROW(splitx[[group.combinations[1,i]]])),
rep(group.combinations[2,i], NROW(splitx[[group.combinations[2,i]]])))
ddelta <- rbind(ddelta, DNCImper:::DNCI.ses(x=paired.x,grouping=group.pair,id=id, Nperm = Nperm, count = count, plotSIMPER = plotSIMPER, dataTYPE = dataTYPE)) #here is the part that calculates the index based on PERSIMPER
}
splitx <- split(x,grouping)
splitx
Add <- which(c(NROW(splitx[[group.combinations[1,i]]]),
NROW(splitx[[group.combinations[2,i]]])) == max(c(NROW(splitx[[group.combinations[1,i]]]),
NROW(splitx[[group.combinations[2,i]]]))))
Add
sampled_lines <- sample(1:length(splitx[[group.combinations[1,i]]][,1]),
length(splitx[[group.combinations[2,i]]][,1]))
splitx[[group.combinations[1,i]]] <- splitx[[group.combinations[1,i]]][sampled_lines,]
paired.x <- rbind(splitx[[group.combinations[1,i]]],
splitx[[group.combinations[2,i]]])
length(which(rowSums(paired.x)
)
)
length(which(rowSums(paired.x) == 0))
if(symmetrize == TRUE)
{
Add <- which(c(NROW(splitx[[group.combinations[1,i]]]),
NROW(splitx[[group.combinations[2,i]]])) == max(c(NROW(splitx[[group.combinations[1,i]]]),
NROW(splitx[[group.combinations[2,i]]]))))
if(Add == 1)
{
sampled_lines <- sample(1:length(splitx[[group.combinations[1,i]]][,1]),
length(splitx[[group.combinations[2,i]]][,1]))
splitx[[group.combinations[1,i]]] <- splitx[[group.combinations[1,i]]][sampled_lines,]
}
if(Add == 2)
{
sampled_lines <- sample(1:length(splitx[[group.combinations[2,i]]][,1]),
length(splitx[[group.combinations[1,i]]][,1]))
splitx[[group.combinations[2,i]]] <- splitx[[group.combinations[2,i]]][sampled_lines,]
}
}
splitx <- split(x,grouping)
#Ici symmetrize:
if(symmetrize == TRUE)
{
Add <- which(c(NROW(splitx[[group.combinations[1,i]]]),
NROW(splitx[[group.combinations[2,i]]])) == max(c(NROW(splitx[[group.combinations[1,i]]]),
NROW(splitx[[group.combinations[2,i]]]))))
if(Add == 1)
{
sampled_lines <- sample(1:length(splitx[[group.combinations[1,i]]][,1]),
length(splitx[[group.combinations[2,i]]][,1]))
splitx[[group.combinations[1,i]]] <- splitx[[group.combinations[1,i]]][sampled_lines,]
}
if(Add == 2)
{
sampled_lines <- sample(1:length(splitx[[group.combinations[2,i]]][,1]),
length(splitx[[group.combinations[1,i]]][,1]))
splitx[[group.combinations[2,i]]] <- splitx[[group.combinations[2,i]]][sampled_lines,]
}
}
i
paired.x <- rbind(splitx[[group.combinations[1,i]]],
splitx[[group.combinations[2,i]]])
ifzero <- which(apply(paired.x, 2, sum) == 0)
if(length(ifzero > 0)){
paired.x <- paired.x[,-which(colSums(paired.x)==0)]}
if(length(which(rowSums(paired.x) == 0)) != 0){stop("ERROR : A row/sample is empty")}
length(which(rowSums(paired.x) == 0))
rowSums(paired.x)
paired.x
rowSums(paired.x)
which(rowSums(paired.x) == 0)
rowSums(paired.x) == 0
which(rowSums(paired.x) == 0)
length(which(rowSums(paired.x) == 0))
apply(paired.x, 2, sum)
paired.x
which(apply(paired.x, 2, sum) == 0)
comm_df_for_DNCI
library(data.table)
library(tidyverse)
library(bench)
## R packages that provide equivalent functions of my package
library(adespatial)
library(DNCImper)
comm_df_for_DNCI<-read.csv("~/.julia/dev/MetaCommunityMetrics/benchmarks/benchmark_r/data/DNCI_comm_t50.csv")
grouping_for_DNCI<-read.csv("~/.julia/dev/MetaCommunityMetrics/benchmarks/benchmark_r/data/cluster_list_t50.csv")
groups <- grouping_for_DNCI$Group
groups
omm_df_for_DNCI
comm_df_for_DNCI
groups
DNCI_multigroup_result <- mark(DNCImper:::DNCI_multigroup(comm_df_for_DNCI,
groups, Nperm = 100,
symmetrize = FALSE,
plotSIMPER = FALSE),
iterations = 1000,
check = TRUE,
time_unit = "ms")
DNCI_multigroup_result
comm_df_for_DNCI
DNCI_multigroup_result <- mark(DNCImper:::DNCI_multigroup(comm_df_for_DNCI,
groups, Nperm = 100,
symmetrize = FALSE,
plotSIMPER = FALSE),
iterations = 100,
check = TRUE,
time_unit = "ms")
DNCI_multigroup_result <- mark(DNCImper:::DNCI_multigroup(comm_df_for_DNCI,
groups, Nperm = 100,
symmetrize = FALSE,
plotSIMPER = FALSE),
iterations = 10,
check = TRUE,
time_unit = "ms")
DNCI_multigroup_result <- mark(DNCImper:::DNCI_multigroup(comm_df_for_DNCI,
groups, Nperm = 10,
symmetrize = FALSE,
plotSIMPER = FALSE),
iterations = 10,
check = TRUE,
time_unit = "ms")
DNCI_multigroup_result <- mark(DNCImper:::DNCI_multigroup(comm_df_for_DNCI,
groups, Nperm = 100,
symmetrize = FALSE,
plotSIMPER = FALSE),
iterations = 1,
check = TRUE,
time_unit = "ms")
DNCI_multigroup_result <- mark(DNCImper:::DNCI_multigroup(comm_df_for_DNCI,
groups, Nperm = 100,
symmetrize = FALSE,
plotSIMPER = FALSE),
iterations = 2,
check = FALSE,
time_unit = "ms")
DNCI_multigroup_result
DNCI_multigroup_result <- mark(DNCImper:::DNCI_multigroup(comm_df_for_DNCI,
groups, Nperm = 10,
symmetrize = FALSE,
plotSIMPER = FALSE),
iterations = 100,
check = FALSE,
time_unit = "ms")
DNCI_multigroup_result <- mark(DNCImper:::DNCI_multigroup(comm_df_for_DNCI,
groups, Nperm = 100,
symmetrize = FALSE,
plotSIMPER = FALSE),
iterations = 10,
check = FALSE,
time_unit = "ms")
DNCI_multigroup_result
execution_time_millisecond <- as.numeric(mean(temporal_beta_div_3$time[[1]])) * 1000
execution_time_millisecond <- as.numeric(mean(DNCI_multigroup_result $time[[1]])) * 1000
memory_usage_mib <- as.numeric(DNCI_multigroup_result $mem_alloc) / 1.048576e+6
cat("Execution Time (Milliseconds):", execution_time_millisecond, "\n")
cat("Memory Usage (MiB):", memory_usage_mib, "\n")
as.numeric(DNCI_multigroup_result $mem_alloc)
DNCI_multigroup_result $mem_alloc
11073560448/1.048576e+6
mean(DNCI_multigroup_result $time[[1]])
41.2*1000
DNCI_multigroup_result_2 <- mark(DNCImper:::DNCI_multigroup(comm_df_for_DNCI,
groups, Nperm = 100,
symmetrize = FALSE,
plotSIMPER = FALSE),
iterations = 10,
check = FALSE,
time_unit = "ms")
DNCI_multigroup_result_2
DNCI_multigroup_result_3 <- mark(DNCImper:::DNCI_multigroup(comm_df_for_DNCI,
groups, Nperm = 100,
symmetrize = FALSE,
plotSIMPER = FALSE),
iterations = 100,
check = FALSE,
time_unit = "ms")
DNCI_multigroup_result_3
execution_time_millisecond <- as.numeric(mean(DNCI_multigroup_result_3$time[[1]])) * 1000
memory_usage_mib <- as.numeric(DNCI_multigroup_result_3$mem_alloc) / 1.048576e+6
# Print the results
cat("Execution Time (Milliseconds):", execution_time_millisecond, "\n")
cat("Memory Usage (MiB):", memory_usage_mib, "\n")
saveRDS(DNCI_multigroup_result_3, "~/.julia/dev/MetaCommunityMetrics/benchmarks/benchmark_r/benchmarking_result_R/DNCI_multigroup_result_100iter.rds")
DNCI_multigroup_result_3$time
df
CV_simple_function(df$Species, df$Sampling_date_order, df$plot, df$Abundance)
CV_simple_function <- function(species, time, plot, abundance){
# Extract unique values for Species, Sampling_date_order, and plot
species_vals <- unique(species)
date_vals <- unique(time)
plot_vals <- unique(plot)
# Create the array with dimensions based on unique Species, Sampling_date_order, and plot
metacomm_tsdata <- array(0, dim = c(length(species_vals), length(date_vals), length(plot_vals)))
# Populate the array with values from the data frame
for(i in 1:nrow(df)){
# Map Species, Sampling_date_order, and plot to their corresponding indices
species_index <- which(species_vals == species[i])
date_index <- which(date_vals == time[i])
plot_index <- which(plot_vals == plot[i])
# Assign the value from the data frame (Abundance) to the array
metacomm_tsdata[species_index, date_index, plot_index] <- abundance[i]
}
result <- var.partition(metacomm_tsdata)
return(result)
}
CV_simple_function(df$Species, df$Sampling_date_order, df$plot, df$Abundance)
var.partition <- function(metacomm_tsdata){
## The function "var.partition" performs the partitioning of variability
## across hierarchical levesl within a metacommunity.
## The input array "metacomm_tsdata" is an N*T*M array. The first dimension represents N species,
## the second represents time-series observations of length T, and the third represents M local communities.
## The output includes four variability and four synchrony metrics as defined in the main text.
## Note that, to be able to handle large metacommunities, this code has avoided calculating all covariance.
ts_metacom <- apply(metacomm_tsdata,2,sum) #summing the total biovolume across all species and patch at every time point
ts_patch <- apply(metacomm_tsdata,c(2,3),sum) #summing the total biovolume across all species in every patch at every time point
ts_species <- apply(metacomm_tsdata,c(1,2),sum)#summing the total biovolume across all patches for every species at every time point
sd_metacom <- sd(ts_metacom)
sd_patch_k <- apply(ts_patch,2,sd)
sd_species_i <- apply(ts_species,1,sd)
sd_species_patch_ik <- apply(metacomm_tsdata,c(1,3),sd)
mean_metacom <- mean(ts_metacom)
CV_S_L <- sum(sd_species_patch_ik)/mean_metacom
CV_C_L <- sum(sd_patch_k)/mean_metacom
CV_S_R <- sum(sd_species_i)/mean_metacom
CV_C_R <- sd_metacom/mean_metacom
phi_S_L2R <- CV_S_R/CV_S_L
phi_C_L2R <- CV_C_R/CV_C_L
phi_S2C_L <- CV_C_L/CV_S_L
phi_S2C_R <- CV_C_R/CV_S_R
partition_3level <- c(CV_S_L=CV_S_L, CV_C_L=CV_C_L, CV_S_R=CV_S_R, CV_C_R=CV_C_R,
phi_S_L2R=phi_S_L2R, phi_C_L2R=phi_C_L2R, phi_S2C_L=phi_S2C_L,
phi_S2C_R=phi_S2C_R)
return(partition_3level)
}
CV_simple_function(df$Species, df$Sampling_date_order, df$plot, df$Abundance)
df
??sd()
ts_metacom
apply(metacomm_tsdata,2,sum)
species= df$Species
time=df$Sampling_date_order
plot=df$plot
species_vals <- unique(species)
date_vals <- unique(time)
plot_vals <- unique(plot)
metacomm_tsdata <- array(0, dim = c(length(species_vals), length(date_vals), length(plot_vals)))
# Populate the array with values from the data frame
for(i in 1:nrow(df)){
# Map Species, Sampling_date_order, and plot to their corresponding indices
species_index <- which(species_vals == species[i])
date_index <- which(date_vals == time[i])
plot_index <- which(plot_vals == plot[i])
# Assign the value from the data frame (Abundance) to the array
metacomm_tsdata[species_index, date_index, plot_index] <- abundance[i]
}
abundance=df$Abundance
# Populate the array with values from the data frame
for(i in 1:nrow(df)){
# Map Species, Sampling_date_order, and plot to their corresponding indices
species_index <- which(species_vals == species[i])
date_index <- which(date_vals == time[i])
plot_index <- which(plot_vals == plot[i])
# Assign the value from the data frame (Abundance) to the array
metacomm_tsdata[species_index, date_index, plot_index] <- abundance[i]
}
metacomm_tsdata
ts_metacom <- apply(metacomm_tsdata,2,sum) #summing the total biovolume across all species and patch at every time point
ts_metacom
df
df <- read.csv("/Users/yc2864/Documents/research/MetaCommunityMetrics.jl/data/metacomm_rodent_df.csv")
df
species=df$Species
time=df$Sampling_date_order
plot=df$plot
abundance=df$Abundance
species_vals <- unique(species)
date_vals <- unique(time)
plot_vals <- unique(plot)
# Create the array with dimensions based on unique Species, Sampling_date_order, and plot
metacomm_tsdata <- array(0, dim = c(length(species_vals), length(date_vals), length(plot_vals)))
# Populate the array with values from the data frame
for(i in 1:nrow(df)){
# Map Species, Sampling_date_order, and plot to their corresponding indices
species_index <- which(species_vals == species[i])
date_index <- which(date_vals == time[i])
plot_index <- which(plot_vals == plot[i])
# Assign the value from the data frame (Abundance) to the array
metacomm_tsdata[species_index, date_index, plot_index] <- abundance[i]
}
ts_metacom <- apply(metacomm_tsdata,2,sum) #summing the total biovolume across all species and patch at every time point
ts_metacom
df
library(tidyverse)
df%>%
filter(Species =="BA")
df%>%
filter(Species =="BA")%>%summarise(total_abundacne=sum(Abundance))
df%>%
filter(Species =="BA")%>%summarise(total_abundance=sum(Abundance))
data <- read.csv("/Users/yc2864/Documents/research/MetaCommunityMetrics.jl/data/metacomm_rodent_df.csv")
df <- data%>%
filter(Species =="BA")%>%summarise(total_abundance=sum(Abundance))
species=df$Species
time=df$Sampling_date_order
plot=df$plot
abundance=df$Abundance
species_vals <- unique(species)
date_vals <- unique(time)
plot_vals <- unique(plot)
# Create the array with dimensions based on unique Species, Sampling_date_order, and plot
metacomm_tsdata <- array(0, dim = c(length(species_vals), length(date_vals), length(plot_vals)))
# Populate the array with values from the data frame
for(i in 1:nrow(df)){
# Map Species, Sampling_date_order, and plot to their corresponding indices
species_index <- which(species_vals == species[i])
date_index <- which(date_vals == time[i])
plot_index <- which(plot_vals == plot[i])
# Assign the value from the data frame (Abundance) to the array
metacomm_tsdata[species_index, date_index, plot_index] <- abundance[i]
}
metacomm_tsdata
df
df <- data%>%
filter(Species =="BA")
species=df$Species
time=df$Sampling_date_order
plot=df$plot
abundance=df$Abundance
species_vals <- unique(species)
date_vals <- unique(time)
plot_vals <- unique(plot)
# Create the array with dimensions based on unique Species, Sampling_date_order, and plot
metacomm_tsdata <- array(0, dim = c(length(species_vals), length(date_vals), length(plot_vals)))
# Populate the array with values from the data frame
for(i in 1:nrow(df)){
# Map Species, Sampling_date_order, and plot to their corresponding indices
species_index <- which(species_vals == species[i])
date_index <- which(date_vals == time[i])
plot_index <- which(plot_vals == plot[i])
# Assign the value from the data frame (Abundance) to the array
metacomm_tsdata[species_index, date_index, plot_index] <- abundance[i]
}
metacomm_tsdata
ts_metacom <- apply(metacomm_tsdata,2,sum) #summing the total biovolume across all species and patch at every time point
ts_metacom
View(data)
unique(data$Species)
data <- read.csv("/Users/yc2864/Documents/research/MetaCommunityMetrics.jl/data/metacomm_rodent_df.csv",na.strings = c(""), stringsAsFactors = FALSE)
unique(data$Species)
df <- data
species=df$Species
time=df$Sampling_date_order
plot=df$plot
abundance=df$Abundance
CV_simple_function(df$Species, df$Sampling_date_order, df$plot, df$Abundance)
unique(df$Species)
# Extract unique values for Species, Sampling_date_order, and plot
species_vals <- unique(species)
date_vals <- unique(time)
plot_vals <- unique(plot)
# Create the array with dimensions based on unique Species, Sampling_date_order, and plot
metacomm_tsdata <- array(0, dim = c(length(species_vals), length(date_vals), length(plot_vals)))
# Populate the array with values from the data frame
for(i in 1:nrow(df)){
# Map Species, Sampling_date_order, and plot to their corresponding indices
species_index <- which(species_vals == species[i])
date_index <- which(date_vals == time[i])
plot_index <- which(plot_vals == plot[i])
# Assign the value from the data frame (Abundance) to the array
metacomm_tsdata[species_index, date_index, plot_index] <- abundance[i]
}
ts_metacom <- apply(metacomm_tsdata,2,sum) #summing the total biovolume across all species and patch at every time point
ts_metacom
ts_species <- apply(metacomm_tsdata,c(1,2),sum)#summing the total biovolume across all patches for every species at every time point
ts_species
sd_species_i <- apply(ts_species,1,sd)
sd_species_i
ts_metacom
CV_simple_function(df$Species, df$Sampling_date_order, df$plot, df$Abundance)
