library(Matching)

# before Sc, Sc, ScR, R, RRs1, Rs1
veins <- c(1, 7, 15, 7, 2, 8) # semicanescens at Sc
veins <- c(7, 15, 7, 3, 8) # semicanescens bt R and Rs1
veins <- c(7, 15, 7, 2, 8) # semicanescens removed

b.probs <- c(dbinom(0,5,0.5), dbinom(1,5,0.5), dbinom(2,5,0.5), dbinom(3,5,0.5), dbinom(4,5,0.5), dbinom(5,5,0.5))
u.probs <- rep(1/6,6)

s <- 1000
b <- round(b.probs*s, 0)
u <- round(u.probs*s, 0)

m <- ks.boot(veins,b,1000) # p-value = 0.004958
m <- ks.boot(veins,u,1000) # p-value = 0.004958

m <- ks.test(veins,b,1000) # p-value = 0.004958
m <- ks.test(veins,u,1000) # p-value = 0.004958



# before Sc, Sc, ScR, R, RRs1, Rs1
veins <- c(7, 15, 7, 3, 8) # semicanescens bt R and Rs1

b.probs <- c(dbinom(0,4,0.5), dbinom(1,4,0.5), dbinom(2,4,0.5), dbinom(3,4,0.5), dbinom(4,4,0.5))
u.probs <- rep(1/5,5)

s <- 1000
b <- round(b.probs*s, 0)
u <- round(u.probs*s, 0)

m <- ks.boot(veins,b,1000) # p-value = 0.01348
m <- ks.boot(veins,u,1000) # p-value = 0.01348

m <- ks.test(veins,b,1000) # p-value = 0.01348
m <- ks.test(veins,u,1000) # p-value = 0.01348




# before Sc, Sc, ScR, R, RRs1, Rs1
veins <- c(7, 15, 7, 2, 8) # semicanescens removed

b.probs <- c(dbinom(0,4,0.5), dbinom(1,4,0.5), dbinom(2,4,0.5), dbinom(3,4,0.5), dbinom(4,4,0.5))
u.probs <- rep(1/5,5)

s <- 1000
b <- round(b.probs*s, 0)
u <- round(u.probs*s, 0)

m <- ks.boot(veins,b,1000) # p-value = 0.01348
m <- ks.boot(veins,u,1000) # p-value = 0.01348

m <- ks.test(veins,b,1000) # p-value = 0.01348
m <- ks.test(veins,u,1000) # p-value = 0.01348





library(ggplot2)

new <- ggplot() +
theme_bw() +
theme(axis.title=element_text(size=10), plot.title=element_text(size=10, hjust = 0.5), axis.text.x = element_text(angle = 45, hjust = 1), axis.text.y = element_text(angle = 90, hjust = 0.5), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_rect(color = "black", size=0.5)) +
xlab("Distal edge of CSS along costa") +
ylab("Number of specimens") +
annotate("rect", xmin=seq(0.5,5.5,1), xmax=seq(1.5,6.5,1), ymin=0, ymax=veins, col="black", fill="#BBBBBB") +
scale_x_continuous(limits=c(0.5,6.5),breaks=seq(1,6,1),labels=c("Before Sc","At Sc","B/t Sc & R","At R","B/t R and Rs1","At Rs1")) +
scale_y_continuous(breaks=seq(0,16,2)) +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
theme(axis.text.y = element_text(angle = 0, hjust = 1))

pdf("histogram.pdf", height=3.5, width=3.25)
new
dev.off()

# No assumptions about consistency within species, either within or between sexes


my.unif <- seq(1,6,1)

# returns percentage of trials in which any value is 1 - any character state is represented by only 1 specimen
test.u <- function(pop, tries) {
	m <- rep(my.unif,round(pop/length(my.unif)))
	out <- NULL
	for (i in 1:tries) {
		this <- m[sample(pop,40)]
		this <- as.numeric(table(this))
		ret <- 0
		if (min(this < 3)) {
			ret <- 1
		}
		out <- c(out,ret)
	}
	out <- (length(which(out > 0)) / length(out))
	return(out)
}

u1000 <- test.u(1000,10000) # 0
u10000 <- test.u(10000,10000) # 0


my.norm <- c(5, 27, 68, 68, 27, 5)

# returns percentage of trials in which an intermediate value is 2 or 1 - any intermediate character state is represented by only 1 or 2 specimens
test.n <- function(pop, tries) {
	m <- rep(my.norm,round(pop/length(my.norm)))
	out <- NULL
	for (i in 1:tries) {
		this <- m[sample(pop,40)]
		this <- this[which(this > 1)]
		this <- this[which(this < 6)]
		this <- as.numeric(table(this))
		ret <- 0
		if (min(this < 3)) {
			ret <- 1
		}
		out <- c(out,ret)
	}
	out <- (length(which(out > 0)) / length(out))
	return(out)
}

n1000 <- test.n(1000,10000) # 0
n10000 <- test.n(10000,10000) # 0