n.row <- n.col <- 50

##parameters
f<-1 ##fertility advantage of females
n<-50 ##ovule number
s<-.1 ##selfing rate
ib<-0.9 ##inbreeding depression
q<-0.01 ## number of seed going to random cells
h<-0.01 ## number of pollen grains going to random cells
meanp<-10 ## mean number of trips from each cell
N<- 100 ##number of simulations
B<- 100 ##burn in length
I<- 3 ##number of increments
num<-0 ## number of cycles in total

results<-matrix(0, (I*N), 9)
par<-matrix(0, (I*N),7)

fertility<-c(f)
is.vector(fertility)
generation<-c(0)
is.vector(generation)
cycle<-c(0)
is.vector(cycle)
nherm<-c(0)
is.vector(nherm)
sim<-c(0)
is.vector(sim)
e<-c(0)
is.vector(e)


##change parameters

for(var in 1:I)												{

bankf<-matrix(0, n.row, n.col)
bankh<-matrix(0, n.row, n.col)
if(var==2) s=.5
if(var==3) s=.9




for(A in 1:N){
g<-1 ##generation
##matrix of herm
m <- matrix(1, n.row, n.col)

##sum up number of hermaphrodites
summ<-length((m)[m==1])
empty<-length((m)[m==3])
sumf<-length((m)[m==0])

ph<-(summ/(sumf+summ))	

burn=0
while(burn<=B & empty<(n.col*n.row)){
burn=burn+1

fem<-matrix(0,n.row, n.col)
out<-matrix(0,n.row, n.col)
self<-matrix(0,n.row, n.col)

i<-m==1
pollen <- matrix(0,n.row * n.col)
pollen[i] <- meanp
pollen <- matrix(pollen,n.row, n.col)

self[i]=rpois((sum(i)), (s*n*(1-ib)))  
out[i]=rpois((sum(i)), ((1-s)*n))  
						 						 
lpollen<-matrix(0, n.row, n.col)
rpollen<-matrix(0, n.row, n.col)
rpollen<-(h*pollen)
lpollen<-((1-h)*pollen) 
grains<-matrix(0, n.row, n.col)

rec <- matrix(rpois((n.row*n.col)^2, rpollen/(n.row*n.col-1)), n.row*n.col)
diag(rec) <- 0
grains <- matrix(colSums(rec),n.row,n.col)

v<-matrix(0,1,n.col)
u<-matrix(0,(n.row-1),1)
a<-matrix(rpois((n.row*n.col),(lpollen/8)), n.row,n.col)
b<-a[2:n.row, 1:(n.col-1)]
grains<-grains+rbind(cbind(u,b),v)

a<-matrix(rpois((n.row*n.col),(lpollen/8)), n.row,n.col)
b<-a[2:n.row, 2:n.col]
grains<-grains+rbind(cbind(b,u),v)

a<-matrix(rpois((n.row*n.col),(lpollen/8)), n.row,n.col)
b<-a[1:(n.row-1), 1:(n.col-1)]
grains<-grains+rbind(v,(cbind(u,b)))

a<-matrix(rpois((n.row*n.col),(lpollen/8)), n.row,n.col)
b<-a[1:(n.row-1), 2:n.col]
grains<-grains+rbind(v,(cbind(b,u)))

v<-matrix(0,n.row,1) 
a<-matrix(rpois((n.row*n.col),(lpollen/8)), n.row,n.col)
b<-a[1:n.row, 2:n.col]
grains<-grains+cbind(b,v)

a<-matrix(rpois((n.row*n.col),(lpollen/8)), n.row,n.col)
b<-a[1:n.row, 1:(n.col-1)]
grains<-grains+cbind(v,b)

v<-matrix(0,1,n.col)
a<-matrix(rpois((n.row*n.col),(lpollen/8)), n.row,n.col)
b<-a[2:n.row,1:n.col]
grains<-grains+rbind(b,v)

a<-matrix(rpois((n.row*n.col),(lpollen/8)), n.row,n.col)
b<-a[1:(n.row-1),1:n.col]
grains<-grains+rbind(v,b)

herm<-ifelse(m==1, ifelse(grains<1, self, out+self), 0)

herml<-round((1-q)*herm, digits=0) 
hermr<-herm-herml
						 
rec <- matrix(rpois((n.row*n.col)^2, hermr/(n.row*n.col-1)), n.row*n.col)
diag(rec) <- 0
toth <- matrix(colSums(rec),n.row,n.col)

a<-matrix(rpois((n.row*n.col),(herml/9)), n.row,n.col)
toth<-a+toth

v<-matrix(0,1,n.col)
u<-matrix(0,(n.row-1),1)
a<-matrix(rpois((n.row*n.col),(herml/9)), n.row,n.col)
b<-a[2:n.row, 1:(n.col-1)]
toth<-toth+rbind(cbind(u,b),v)

a<-matrix(rpois((n.row*n.col),(herml/9)), n.row,n.col)
b<-a[2:n.row, 2:n.col]
toth<-toth+rbind(cbind(b,u),v)

a<-matrix(rpois((n.row*n.col),(herml/9)), n.row,n.col)
b<-a[1:(n.row-1), 1:(n.col-1)]
toth<-toth+rbind(v,(cbind(u,b)))

a<-matrix(rpois((n.row*n.col),(herml/9)), n.row,n.col)
b<-a[1:(n.row-1), 2:n.col]
toth<-toth+rbind(v,(cbind(b,u)))

v<-matrix(0,n.row,1) ##row 10 long
a<-matrix(rpois((n.row*n.col),(herml/9)), n.row,n.col)
b<-a[1:n.row, 2:n.col]
toth<-toth+cbind(b,v)

a<-matrix(rpois((n.row*n.col),(herml/9)), n.row,n.col)
b<-a[1:n.row, 1:(n.col-1)]
toth<-toth+cbind(v,b)

v<-matrix(0,1,n.col) 
a<-matrix(rpois((n.row*n.col),(herml/9)), n.row,n.col)
b<-a[2:n.row,1:n.col]
toth<-toth+rbind(b,v)

a<-matrix(rpois((n.row*n.col),(herml/9)), n.row,n.col)
b<-a[1:(n.row-1),1:n.col]
toth<-toth+rbind(v,b)
	
seed<-matrix(0, n.row, n.col)
seed<-(toth)
pherm<-(toth/seed)
					
m<-(ifelse(seed > 0, 1, 3))
		
												
##print("the matrix of empty and hermaphrodites is")
##print(m)

summ<-length((m)[m==1])
empty<-length((m)[m==3])
sumf<-length((m)[m==0])

ph<-(summ/(sumf+summ))	



}


##sample one cell 
i <- sample(1:(n.row*n.col), 1)
##and make female
m[i]=0

##sum up number of hermaphrodites
summ<-length((m)[m==1])
empty<-length((m)[m==3])
sumf<-length((m)[m==0])
burne<-empty

ph<-(summ/(sumf+summ))							
g<-1
seedf<-matrix(0,n.row,n.col)
seedh<-matrix(0,n.row,n.col)																																																	
##if number of hermaphrodites is less than total and all of the cells are not empty(i.e. pop extinction)
while(g<= 100 &( empty+summ )<(n.row*n.col)) {


i<-m==1
pollen <- matrix(0,n.row * n.col)
pollen[i] <- meanp
pollen <- matrix(pollen,n.row, n.col)

fem<-matrix(0,n.row, n.col)
out<-matrix(0,n.row, n.col)
self<-matrix(0,n.row, n.col)

i<-m==1
self[i]=rpois((sum(i)), (s*n*(1-ib)))  
out[i]=rpois((sum(i)), ((1-s)*n))  
i<-m==0
fem[i]=rpois((sum(i)), ((1+f)*n)) 
						 
lpollen<-matrix(0, n.row, n.col)
rpollen<-matrix(0, n.row, n.col)
rpollen<-(h*pollen)
lpollen<-((1-h)*pollen) 
grains<-matrix(0, n.row, n.col)

rec <- matrix(rpois((n.row*n.col)^2, rpollen/(n.row*n.col-1)), n.row*n.col)
diag(rec) <- 0
grains <- matrix(colSums(rec),n.row,n.col)


v<-matrix(0,1,n.col)
u<-matrix(0,(n.row-1),1)
a<-matrix(rpois((n.row*n.col),(lpollen/8)), n.row,n.col)
b<-a[2:n.row, 1:(n.col-1)]
grains<-grains+rbind(cbind(u,b),v)

a<-matrix(rpois((n.row*n.col),(lpollen/8)), n.row,n.col)
b<-a[2:n.row, 2:n.col]
grains<-grains+rbind(cbind(b,u),v)

a<-matrix(rpois((n.row*n.col),(lpollen/8)), n.row,n.col)
b<-a[1:(n.row-1), 1:(n.col-1)]
grains<-grains+rbind(v,(cbind(u,b)))

a<-matrix(rpois((n.row*n.col),(lpollen/8)), n.row,n.col)
b<-a[1:(n.row-1), 2:n.col]
grains<-grains+rbind(v,(cbind(b,u)))

v<-matrix(0,n.row,1) 
a<-matrix(rpois((n.row*n.col),(lpollen/8)), n.row,n.col)
b<-a[1:n.row, 2:n.col]
grains<-grains+cbind(b,v)

a<-matrix(rpois((n.row*n.col),(lpollen/8)), n.row,n.col)
b<-a[1:n.row, 1:(n.col-1)]
grains<-grains+cbind(v,b)

v<-matrix(0,1,n.col) 
a<-matrix(rpois((n.row*n.col),(lpollen/8)), n.row,n.col)
b<-a[2:n.row,1:n.col]
grains<-grains+rbind(b,v)

a<-matrix(rpois((n.row*n.col),(lpollen/8)), n.row,n.col)
b<-a[1:(n.row-1),1:n.col]
grains<-grains+rbind(v,b)

		
fem<-ifelse(m==0, ifelse(grains<1, 0, fem), 0)
herm<-ifelse(m==1, ifelse(grains<1, self, out+self), 0)

toth<-matrix(0,n.col,n.row)
totf<-matrix(0,n.col,n.row)
herml<-round((1-q)*herm, digits=0) 
feml<-round((1-q)*fem, digits=0)
hermr<-herm-herml
femr<-fem-feml
						 
rec <- matrix(rpois((n.row*n.col)^2, hermr/(n.row*n.col-1)), n.row*n.col)
diag(rec) <- 0
toth <- matrix(colSums(rec),n.row,n.col)

rec <- matrix(rpois((n.row*n.col)^2, femr/(n.row*n.col-1)), n.row*n.col)
diag(rec) <- 0
totf <- matrix(colSums(rec),n.row,n.col)

a<-matrix(rpois((n.row*n.col),(herml/9)), n.row,n.col)
toth<-toth+a

v<-matrix(0,1,n.col)
u<-matrix(0,(n.row-1),1)
a<-matrix(rpois((n.row*n.col),(herml/9)), n.row,n.col)
b<-a[2:n.row, 1:(n.col-1)]
toth<-toth+rbind(cbind(u,b),v)

a<-matrix(rpois((n.row*n.col),(herml/9)), n.row,n.col)
b<-a[2:n.row, 2:n.col]
toth<-toth+rbind(cbind(b,u),v)

a<-matrix(rpois((n.row*n.col),(herml/9)), n.row,n.col)
b<-a[1:(n.row-1), 1:(n.col-1)]
toth<-toth+rbind(v,(cbind(u,b)))

a<-matrix(rpois((n.row*n.col),(herml/9)), n.row,n.col)
b<-a[1:(n.row-1), 2:n.col]
toth<-toth+rbind(v,(cbind(b,u)))

v<-matrix(0,n.row,1) 
a<-matrix(rpois((n.row*n.col),(herml/9)), n.row,n.col)
b<-a[1:n.row, 2:n.col]
toth<-toth+cbind(b,v)

a<-matrix(rpois((n.row*n.col),(herml/9)), n.row,n.col)
b<-a[1:n.row, 1:(n.col-1)]
toth<-toth+cbind(v,b)

v<-matrix(0,1,n.col)
a<-matrix(rpois((n.row*n.col),(herml/9)), n.row,n.col)
b<-a[2:n.row,1:n.col]
toth<-toth+rbind(b,v)

a<-matrix(rpois((n.row*n.col),(herml/9)), n.row,n.col)
b<-a[1:(n.row-1),1:n.col]
toth<-toth+rbind(v,b)

a<-matrix(rpois((n.row*n.col),(feml/9)), n.row,n.col)
totf<-totf+a

v<-matrix(0,1,n.col)
u<-matrix(0,(n.row-1),1)
a<-matrix(rpois((n.row*n.col),(feml/9)), n.row,n.col)
b<-a[2:n.row, 1:(n.col-1)]
totf<-totf+rbind(cbind(u,b),v)

a<-matrix(rpois((n.row*n.col),(feml/9)), n.row,n.col)
b<-a[2:n.row, 2:n.col]
totf<-totf+rbind(cbind(b,u),v)

a<-matrix(rpois((n.row*n.col),(feml/9)), n.row,n.col)
b<-a[1:(n.row-1), 1:(n.col-1)]
totf<-totf+rbind(v,(cbind(u,b)))

a<-matrix(rpois((n.row*n.col),(feml/9)), n.row,n.col)
b<-a[1:(n.row-1), 2:n.col]
totf<-totf+rbind(v,(cbind(b,u)))

v<-matrix(0,n.row,1) 
a<-matrix(rpois((n.row*n.col),(feml/9)), n.row,n.col)
b<-a[1:n.row, 2:n.col]
totf<-totf+cbind(b,v)

a<-matrix(rpois((n.row*n.col),(feml/9)), n.row,n.col)
b<-a[1:n.row, 1:(n.col-1)]
totf<-totf+cbind(v,b)

v<-matrix(0,1,n.col) 
a<-matrix(rpois((n.row*n.col),(feml/9)), n.row,n.col)
b<-a[2:n.row,1:n.col]
totf<-totf+rbind(b,v)

a<-matrix(rpois((n.row*n.col),(feml/9)), n.row,n.col)
b<-a[1:(n.row-1),1:n.col]
totf<-totf+rbind(v,b)

##add on seed bank females and herm seeds
totf<-(totf+bankf)
toth<-(toth+bankh)

seed<-matrix(0, n.row, n.col)
seed<-(totf+toth)
pherm<-(toth/seed)
						   
for(i in 1:(n.row*n.col)) {
if(seed[i] > 0) m[i]=rbinom(1, 1, pherm[i]) else m[i]=3 

		 }
		 	 
##print(pherm)		 
seed<-matrix(0, n.row, n.col)		 
##save seeds in seed bank for herm and females


##total number of seeds of each type in the seed bank


seedf<-ifelse((m==0), (totf-1)+seedf, totf+seedf)
seedh<-ifelse((m==1), (toth-1)+seedh, toth+seedh) 
						
print("the matrix of females and hermaphrodites is")
print(m)
g=g+1
print("in gerneration")
print(g)
print("in number")
print(num)
##sum up total number of hermaphrodites

summ<-length((m)[m==1])
empty<-length((m)[m==3])
sumf<-length((m)[m==0])

ph<-(summ/(sumf+summ))		

}		
seedfs<-sum(seedf)
seedhs<-sum(seedh)

num<-num+1
results[num, 1]=num
results[num,2]=g
results[num, 3]=var
results[num,4]=summ
results[num, 5]=A
results[num, 6]=empty
results[num,7]=burne
results[num,8]=seedfs
results[num,9]=seedhs

par[num,1]=num
par[num,2]=f
par[num,3]=s
par[num,4]=ib
par[num,5]=h
par[num,6]=q
par[num,7]=meanp

fertility<-append(fertility, f)
generation<-append(generation, g)
cycle<-append(cycle, var)
nherm<-append(nherm, summ)
sim<-append(sim, A)
e<-append(e, empty)

}}
											
print(results)
print(fertility)
print(generation)
print(cycle)
print(nherm)
print(sim)
print(e)

plot(fertility, generation)

write(results, file = "results",
	ncolumns = 6,
	append = FALSE, sep = " ")

write(par, file = "results",
	ncolumns = 6,
	append = TRUE, sep = " ")
		
write(m, file = "results",
	ncolumns = 6,
	append = TRUE, sep = " ")