Skip to content

Commit

Permalink
Get the y coordinat for the clade branches, progress #42
Browse files Browse the repository at this point in the history
  • Loading branch information
richelbilderbeek committed Jan 13, 2022
1 parent f99975e commit bc38d94
Show file tree
Hide file tree
Showing 2 changed files with 199 additions and 16 deletions.
50 changes: 37 additions & 13 deletions R/plot_daisie_data.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,22 +9,50 @@
plot_daisie_data <- function(daisie_data) {
t <- DAISIEmainland::daisie_data_to_tables(daisie_data)

p <- ggplot2::ggplot() +
p <- ggplot2::ggplot(t$colonists_general) +
ggplot2::scale_x_continuous(
name = "Time",
limits = c(0, t$header$island_age)
) + ggplot2::theme(
axis.title.y = ggplot2::element_blank(),
axis.text.y = ggplot2::element_blank(),
axis.ticks.y = ggplot2::element_blank()
)
# Sometimes, no colonisations take place, making an island history dull ...
if (nrow(t$colonists_general) == 0) {
return(p)
}

# Remove the redundant first branching time, which equals the island
# age, as it is already the header
are_redundant <- t$colonists_branching_times$branching_times == t$header$island_age
t$colonists_branching_times <- t$colonists_branching_times[!are_redundant, ]
# Remove the redundant first branching time, which equals the island
# age, as it is already the header
are_redundant <- t$colonists_branching_times$branching_times == t$header$island_age
t$colonists_branching_times <- t$colonists_branching_times[!are_redundant, ]

# Draw the branches. Note we only have extant species
n_branches_per_clade_index <- dplyr::summarise(
dplyr::group_by(t$colonists_branching_times, clade_index),
n = dplyr::n()
)
branches_horizontal <- t$colonists_branching_times
branches_horizontal$y <- NA
cur_clade_index <- 0
delta_y <- 0
y <- 0
for (row_index in seq_along(branches_horizontal)) {
if (branches_horizontal$clade_index[row_index] != cur_clade_index) {
cur_clade_index <- branches_horizontal$clade_index[row_index]
delta_y <- 1.0 / n_branches_per_clade_index[n_branches_per_clade_index$clade_index == cur_clade_index]$n
y <- delta_y / 2.0
} else {
y <- y + delta_y
}
branches_horizontal$y[row_index] <- y
}
branches_horizontal

# Only obtain the colonisations, i.e. the first branching time

# Only obtain the colonisations, i.e. the first branching time,
# plot these as points with a shape depending on stac_str
first_branching_times <- dplyr::slice(
dplyr::group_by(
t$colonists_branching_times,
Expand All @@ -34,12 +62,8 @@ plot_daisie_data <- function(daisie_data) {
)
colonisations <- merge(first_branching_times, t$colonists_general)

ggplot2::ggplot() +
ggplot2::geom_point(
p + ggplot2::geom_point(
data = colonisations,
ggplot2::aes(x = branching_times, y = clade_index, shape = stac_str)
) + ggplot2::scale_x_continuous(
name = "Time",
limits = c(0, t$header$island_age)
)
ggplot2::aes(x = branching_times, y = 0.5, shape = stac_str)
) + ggplot2::facet_grid(clade_index ~ .)
}
165 changes: 162 additions & 3 deletions tests/testthat/test-plot_daisie_data.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ test_that("No extant colonists", {
expect_silent(plot_daisie_data(empirical_daisie_data))
})

test_that("One colonist clade", {
test_that("One colonist clade, stac = 4: Non_endemic", {
set.seed(
1,
kind = "Mersenne-Twister",
Expand All @@ -42,6 +42,165 @@ test_that("One colonist clade", {
ideal_daisie_data <- daisie_mainland_data$ideal_multi_daisie_data[[1]]
empirical_daisie_data <- daisie_mainland_data$empirical_multi_daisie_data[[1]]
daisie_data <- ideal_daisie_data
expect_silent(plot_daisie_data(ideal_daisie_data))
expect_silent(plot_daisie_data(empirical_daisie_data))
plot_daisie_data(ideal_daisie_data)
plot_daisie_data(empirical_daisie_data)
})

test_that("stac == 2: Endemic (and Non_endemic)", {
set.seed(
4,
kind = "Mersenne-Twister",
normal.kind = "Inversion",
sample.kind = "Rejection"
)
daisie_mainland_data <- sim_island_with_mainland(
total_time = 1,
m = 10,
island_pars = c(1, 1, 10, 0.1, 1),
mainland_ex = 1,
mainland_sample_prob = 1,
mainland_sample_type = "complete",
replicates = 1,
verbose = FALSE
)
ideal_daisie_data <- daisie_mainland_data$ideal_multi_daisie_data[[1]]
empirical_daisie_data <- daisie_mainland_data$empirical_multi_daisie_data[[1]]
testthat::expect_true(
ideal_daisie_data[[2]]$stac != 4 || empirical_daisie_data[[2]]$stac != 4
)
plot_daisie_data(daisie_data = ideal_daisie_data)
plot_daisie_data(empirical_daisie_data)
})

test_that("stac == 5: Endemic_singleton_MaxAge in empirical_daisie_data", {
set.seed(
7,
kind = "Mersenne-Twister",
normal.kind = "Inversion",
sample.kind = "Rejection"
)
daisie_mainland_data <- sim_island_with_mainland(
total_time = 1,
m = 10,
island_pars = c(1, 1, 10, 0.1, 1),
mainland_ex = 1,
mainland_sample_prob = 1,
mainland_sample_type = "complete",
replicates = 1,
verbose = FALSE
)
ideal_daisie_data <- daisie_mainland_data$ideal_multi_daisie_data[[1]]
empirical_daisie_data <- daisie_mainland_data$empirical_multi_daisie_data[[1]]
plot_daisie_data(daisie_data = ideal_daisie_data)
plot_daisie_data(empirical_daisie_data)
})

test_that("stace == 6: Endemic_clade_MaxAge, in empirical_daisie_data", {
set.seed(
40,
kind = "Mersenne-Twister",
normal.kind = "Inversion",
sample.kind = "Rejection"
)
daisie_mainland_data <- sim_island_with_mainland(
total_time = 1,
m = 10,
island_pars = c(1, 1, 10, 0.1, 1),
mainland_ex = 1,
mainland_sample_prob = 1,
mainland_sample_type = "complete",
replicates = 1,
verbose = FALSE
)
ideal_daisie_data <- daisie_mainland_data$ideal_multi_daisie_data[[1]]
empirical_daisie_data <- daisie_mainland_data$empirical_multi_daisie_data[[1]]
plot_daisie_data(daisie_data = ideal_daisie_data)
plot_daisie_data(empirical_daisie_data)
})

test_that("stac == 3: Endemic&Non_Endemic", {
set.seed(
179,
kind = "Mersenne-Twister",
normal.kind = "Inversion",
sample.kind = "Rejection"
)
daisie_mainland_data <- sim_island_with_mainland(
total_time = 1,
m = 10,
island_pars = c(1, 1, 10, 0.1, 1),
mainland_ex = 1,
mainland_sample_prob = 1,
mainland_sample_type = "complete",
replicates = 1,
verbose = FALSE
)
ideal_daisie_data <- daisie_mainland_data$ideal_multi_daisie_data[[1]]
empirical_daisie_data <- daisie_mainland_data$empirical_multi_daisie_data[[1]]
plot_daisie_data(daisie_data = ideal_daisie_data)
plot_daisie_data(empirical_daisie_data)
})


test_that("much branching", {
set.seed(
1018,
kind = "Mersenne-Twister",
normal.kind = "Inversion",
sample.kind = "Rejection"
)
daisie_mainland_data <- sim_island_with_mainland(
total_time = 1,
m = 10,
island_pars = c(1, 1, 10, 0.1, 1),
mainland_ex = 1,
mainland_sample_prob = 1,
mainland_sample_type = "complete",
replicates = 1,
verbose = FALSE
)
ideal_daisie_data <- daisie_mainland_data$ideal_multi_daisie_data[[1]]
empirical_daisie_data <- daisie_mainland_data$empirical_multi_daisie_data[[1]]
plot_daisie_data(daisie_data = ideal_daisie_data)
plot_daisie_data(empirical_daisie_data)
})


test_that("Search for interesting scenarions", {
skip("Only run locally")
seed <- 0
while (1) {
seed <- seed + 1
message(seed)
set.seed(
seed,
kind = "Mersenne-Twister",
normal.kind = "Inversion",
sample.kind = "Rejection"
)
daisie_mainland_data <- sim_island_with_mainland(
total_time = 1,
m = 10,
island_pars = c(1, 1, 10, 0.1, 1),
mainland_ex = 1,
mainland_sample_prob = 1,
mainland_sample_type = "complete",
replicates = 1,
verbose = FALSE
)
ideal_daisie_data <- daisie_mainland_data$ideal_multi_daisie_data[[1]]
empirical_daisie_data <- daisie_mainland_data$empirical_multi_daisie_data[[1]]
if (length(ideal_daisie_data) == 1) next
if ("look for interesting stac" == "what I want") {
uninteresting_set <- c(2, 4, 3, 5, 6)
if (!ideal_daisie_data[[2]]$stac %in% uninteresting_set || !empirical_daisie_data[[2]]$stac %in% uninteresting_set) {
message(seed)
stop(seed)
}
}
if (length(ideal_daisie_data[[2]]$branching_times) > 3) {
message(seed)
stop(seed)
}
}
})

0 comments on commit bc38d94

Please sign in to comment.