Created
February 4, 2013 21:37
-
-
Save anonymous/4709926 to your computer and use it in GitHub Desktop.
Code for Trap Simulations using Shiny
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #StrategicPestTechnologies: Functions | |
| #~~~~~~~~~~~Calc trap info~~~~~~~~~~~~ | |
| trap.info.func<-function(file){ | |
| #Take file | |
| #Return: traps, n.traps, ylims, xlims, total.area | |
| traps<-read.table(file, header=FALSE) | |
| colnames(traps)<-c("ID","X","Y") | |
| n.traps<-dim(traps)[1] | |
| xlims<-c(min(traps$X)-buffer,max(traps$X)+buffer) | |
| ylims<-c(min(traps$Y)-buffer,max(traps$Y)+buffer) | |
| total.area<-(diff(xlims)/100)*(diff(ylims)/100) | |
| trap.area<-((max(traps$X)-min(traps$X))/100) *((max(traps$Y)-min(traps$Y))/100) | |
| list(traps=traps, n.traps=n.traps, ylims=ylims, xlims=xlims, total.area=total.area, trap.area=trap.area) | |
| } | |
| #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ | |
| #~~~~~~~~~~~Create traps info~~~~~~~~~~~~ | |
| trap.make.func<-function(traps.in.x, traps.in.y, trap.space.x, trap.space.y, buffer, start.xy){ | |
| #Take: traps.in.x, traps.in.y, trap.space, buffer, start.xy | |
| #Return: traps, n.traps, ylims, xlims, total.area | |
| x.t<-seq(from=start.xy[1], by=trap.space.x, length.out=traps.in.x) | |
| y.t<-seq(from=start.xy[2], by=trap.space.y, length.out=traps.in.y) | |
| xy<-expand.grid(x.t,y.t) | |
| ID<-1:length(xy[,1]) | |
| traps<-cbind(ID,xy) | |
| colnames(traps)<-c("ID","X","Y") | |
| n.traps<-dim(traps)[1] | |
| xlims<-c(min(traps$X)-buffer,max(traps$X)+buffer) | |
| ylims<-c(min(traps$Y)-buffer,max(traps$Y)+buffer) | |
| total.area<-(diff(xlims)/100)*(diff(ylims)/100) | |
| trap.area<-((max(traps$X)-min(traps$X))/100) *((max(traps$Y)-min(traps$Y))/100) | |
| list(traps=traps, n.traps=n.traps, ylims=ylims, xlims=xlims, total.area=total.area, trap.area=trap.area) | |
| } | |
| #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ | |
| #~~~~~~~~~~Place animals on grid~~~~~~~~~~~~ | |
| animal.loc.func<-function(xlims, ylims, total.area, density, dist.type="F"){ | |
| #Take: xlims, ylims, total.area, density, dist.type, | |
| #Return: animals, n.animals | |
| n.fixed<-round(total.area*density,0) | |
| n.animals<-n.fixed | |
| if(dist.type=="P"){ | |
| n.pois<-rpois(1,n.fixed) | |
| n.animals<-n.pois | |
| } | |
| xloc<-runif(n.animals,min=xlims[1],max=xlims[2]) | |
| yloc<-runif(n.animals,min=ylims[1],max=ylims[2]) | |
| animals<-as.data.frame(cbind(1:n.animals,xloc,yloc)) | |
| list(animals=animals, n.animals=n.animals, dist.type=dist.type) | |
| } | |
| #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ | |
| #~~~~~~~~Plot distance vs prob from g0 and sigma~~~~~~~~~ | |
| plot.g0.sigma<-function(g0,sigma){ | |
| #Calculate what the probability of capture will be: | |
| #i.e. P(Cap) = dnorm(distance,sd=sigma)*g0/dnorm(0,sigma) | |
| k<-g0/dnorm(0,sd=sigma) | |
| # X11() | |
| print(plot(c(0:(3*sigma)),dnorm(seq(from=0, by=1, to=(3*sigma)), sd=sigma)*k, type="l", ylim=c(0,min(c(1,max(c(g0+.1,.2))))), xlab="Distance (m)", ylab="Prob of capture"), xaxs="r", yaxs="r") | |
| by.y<-ifelse(g0<.3,.1,.2) | |
| abline(h=seq(from=0, to=1, b=by.y), col="red", lty=2) | |
| abline(v=seq(from=0, by=20, to=(3*sigma)), col="red", lty=2) | |
| abline(h=0) | |
| abline(v=0) | |
| mtext(paste("g0 = ",g0,", sigma = ",sigma, sep=""),3,padj=-1) | |
| } | |
| #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ | |
| #This runs better than the original as it only calculates the prob when distance < x*Sigma | |
| #~~~~~Calculate distance and probability of capture~~~~~~~~ | |
| dist.prob.func<-function(n.animals, animals, n.traps, traps, g0, sigma, max.sd=3){ | |
| #For each trap[i], the distance^2 and prob to each animal[j]: | |
| ax2<-(outer(traps$X,animals$xloc,'-'))^2 #Distance in x direction (squared) | |
| ay2<-(outer(traps$Y,animals$yloc,'-'))^2 #Distance in y direction (squared) | |
| dist2.xy<-(ax2+ay2) #Combined | |
| max.dist<-(sigma*max.sd)^2 | |
| dist2.xy.tmp<-dist2.xy[dist2.xy<(max.dist)] | |
| # sum(dist2.xy<(max.dist)) | |
| prob.xy3<-matrix(0,n.traps,n.animals) | |
| prob.xy3[dist2.xy<(max.dist)]<-g0*exp(-(dist2.xy.tmp)/(2*sigma^2)) | |
| prob.xy<-prob.xy3 | |
| } | |
| #~~~~~End distance~~~~~~~~~~~~~~ | |
| #~~~~Resample function which works when length(x)==1 | |
| resamp <- function(x,...){if(length(x)==1) x else sample(x,...)} | |
| #~~~~End resample | |
| #~~~~~Calculate which animals are caught where~~~~~~~~~~ | |
| sim.captures.func<-function(n.traps, n.nights, n.animals, max.catch, prob.xy, type="KILL", check="END"){ | |
| #type can be KILL or CHEW | |
| #check can be "END" or "DAY" i.e. every day or at the end of the study | |
| #Needs: n.traps, n.nights, n.animals, prob.xy | |
| #Returns: trap.catch, trap.animals | |
| trap.catch<-matrix(0,n.traps,n.nights) #Trap.catch stores the captures in each trap each night | |
| trap.animals<-rep(0,n.animals) #Indicator of whether an animal has been caught | |
| trap.sum<-rep(0,n.traps) #Record the number caught in each trap | |
| trap.rem<-rep(T,n.traps) #Record if the trap remains in operation...TRUE, FALSE | |
| trap.pot.tmp<-rep(FALSE,n.traps) | |
| trap.animal.mat<-matrix(0,n.traps, n.animals) #Record which trap each animal is detected in... | |
| trap.animal.night<-vector("list",n.nights) #Record which trap each animal is detected in each night. | |
| for (t in 1:n.nights){ #For each night | |
| trap.animal.night[[t]]<-matrix(0,n.traps, n.animals) | |
| if (check =="DAY"){ #If checking each day, then this gets reset... | |
| trap.rem<-rep(T,n.traps) #Record if the trap remains in operation... | |
| trap.pot.tmp<-rep(FALSE,n.traps) | |
| } | |
| not.caught<-(1:n.animals)[trap.animals==0] #Animals not already caught | |
| if (type=="CHEW"){ | |
| not.caught<-(1:n.animals) #For a chew card, all animals can be redetected... | |
| } | |
| for (j in not.caught){ #For each animal not already caught | |
| prob.tmp<-prob.xy[,j]*(trap.rem) #Adjust the probabilities so previous traps cannot catch anything | |
| prob.gt.0<-prob.tmp>0 #"Probabilities greater than 0" | |
| trap.pot<-trap.pot.tmp | |
| trap.pot[prob.gt.0]<-prob.tmp[prob.gt.0]>runif(sum(prob.gt.0)) #Traps that potentially catch something | |
| #...Only does runif for prob>0 | |
| if(sum(trap.pot)>0){ #If at least one trap can catch it... | |
| trap.tmp<-which(trap.pot) #Which traps potentially catch something? | |
| if (type=="KILL"){ | |
| trap.i<-resamp(trap.tmp,1) #Sample which trap it is. | |
| trap.catch[trap.i,t]<-trap.catch[trap.i,t]+1 #Set trap.catch to current + 1 | |
| trap.sum[trap.i]<-trap.sum[trap.i]+1 #Updates the sum of animals caught in trap (gets rid of rowSums) | |
| trap.animals[j]<-1 #Set trapped animal to 1 | |
| } | |
| if (type=="CHEW"){ | |
| trap.i<-trap.tmp #Every card gets chewed as possums can eat more than one | |
| trap.catch[trap.i,t]<-trap.catch[trap.i,t]+1 #Set trap.catch to current + 1 | |
| trap.sum[trap.i]<-trap.sum[trap.i]+1 #Updates the sum of animals caught in trap (gets rid of rowSums) | |
| trap.animals[j]<-trap.animals[j]+1 #Set trapped animal to 1 | |
| trap.animal.mat[trap.i,j]<-trap.animal.mat[trap.i,j]+1 | |
| trap.animal.night[[t]][trap.i,j]<-trap.animal.night[[t]][trap.i,j]+1 | |
| } | |
| trap.rem<-(trap.sum<max.catch) #Calculate which traps are not full | |
| } | |
| } | |
| } | |
| caught.traps<-(rowSums(trap.catch)>0)*1 | |
| list(trap.animals=trap.animals, trap.catch=trap.catch, caught.traps=caught.traps, trap.animal.mat=trap.animal.mat,trap.animal.night=trap.animal.night) | |
| } | |
| #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ | |
| #~~~~~Calculate which animals are caught where~~~~~~~~~~ | |
| sim.captures.func.month<-function(n.traps, n.months, n.animals, max.catch, prob.xy, type="KILL", check=c(1,2)){ | |
| #type can be KILL or CHEW | |
| #check can be "END" or "DAY" i.e. every day or at the end of the study | |
| #Needs: n.traps, n.nights, n.animals, prob.xy | |
| #Returns: trap.catch, trap.animals | |
| n.nights<-n.months*30 | |
| n.checks<-(check*30)+1 | |
| trap.catch<-matrix(0,n.traps,n.nights) #Trap.catch stores the captures in each trap each night | |
| trap.animals<-rep(0,n.animals) #Indicator of whether an animal has been caught | |
| trap.sum<-rep(0,n.traps) #Record the number caught in each trap | |
| trap.rem<-rep(T,n.traps) #Record if the trap remains in operation...TRUE, FALSE | |
| trap.pot.tmp<-rep(FALSE,n.traps) | |
| trap.animal.mat<-matrix(0,n.traps, n.animals) #Record which trap each animal is detected in... | |
| trap.animal.night<-vector("list",n.nights) #Record which trap each animal is detected in each night. | |
| for (t in 1:n.nights){ #For each night | |
| trap.animal.night[[t]]<-matrix(0,n.traps, n.animals) | |
| if (sum(t==n.checks)>0){ #If checking time, then this gets reset... | |
| trap.rem<-rep(T,n.traps) #Record if the trap remains in operation... | |
| trap.pot.tmp<-rep(FALSE,n.traps) | |
| } | |
| not.caught<-(1:n.animals)[trap.animals==0] #Animals not already caught | |
| for (j in not.caught){ #For each animal not already caught | |
| prob.tmp<-prob.xy[,j]*(trap.rem) #Adjust the probabilities so previous traps cannot catch anything | |
| prob.gt.0<-prob.tmp>0 #"Probabilities greater than 0" | |
| trap.pot<-trap.pot.tmp | |
| trap.pot[prob.gt.0]<-prob.tmp[prob.gt.0]>runif(sum(prob.gt.0)) #Traps that potentially catch something | |
| #...Only does runif for prob>0 | |
| if(sum(trap.pot)>0){ #If at least one trap can catch it... | |
| trap.tmp<-which(trap.pot) #Which traps potentially catch something? | |
| if (type=="KILL"){ | |
| trap.i<-resamp(trap.tmp,1) #Sample which trap it is. | |
| trap.catch[trap.i,t]<-trap.catch[trap.i,t]+1 #Set trap.catch to current + 1 | |
| trap.sum[trap.i]<-trap.sum[trap.i]+1 #Updates the sum of animals caught in trap (gets rid of rowSums) | |
| trap.animals[j]<-1 #Set trapped animal to 1 | |
| } | |
| trap.rem<-(trap.sum<max.catch) #Calculate which traps are not full | |
| } | |
| } | |
| } | |
| caught.traps<-(rowSums(trap.catch)>0)*1 | |
| list(trap.animals=trap.animals, trap.catch=trap.catch, caught.traps=caught.traps, trap.animal.mat=trap.animal.mat,trap.animal.night=trap.animal.night) | |
| } | |
| #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ | |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| # Define server logic required to summarize and view the selected dataset | |
| shinyServer(function(input, output) { | |
| data<-reactive(function(){ | |
| #Read in the values from the ui | |
| trap.space.x<-input$trap.space.x | |
| trap.space.y<-input$trap.space.y | |
| traps.in.x<-input$traps.in.x | |
| traps.in.y<-input$traps.in.y | |
| g0<-input$g0 | |
| sigma<-input$sigma | |
| max.catch<-input$max.catch | |
| n.nights<-input$n.nights | |
| density<-input$density | |
| start.xy<-c(1000,1000) #Don't need to change this | |
| buffer<-50 | |
| dist.type<-input$dist.type | |
| reps<-input$reps | |
| #Call functions | |
| source("pest.tech.functions.R") | |
| dens.vector<-density | |
| d<-1 | |
| #Set up matrices to hold the results | |
| n.caught<-rep(NA,reps) | |
| p.caught<-rep(NA,reps) | |
| p.summ<-matrix(NA,length(dens.vector),3) | |
| t.catch<-rep(0,13) | |
| hist.brks<-c(-.5,0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5,10.5,11.5,12.5) | |
| pb <- winProgressBar(title = "progress bar", min = 0, max = length(dens.vector)*reps, width = 300) | |
| for(s in 1:reps){ | |
| #---------Make the traps from these input parameters above | |
| trap.info<-trap.make.func(traps.in.x, traps.in.y, trap.space.x, trap.space.y, buffer, start.xy) | |
| # print(paste(trap.info$trap.area, "ha")) | |
| xlims<-trap.info$xlims | |
| ylims<-trap.info$ylims | |
| traps<-trap.info$traps | |
| n.traps<-trap.info$n.traps | |
| total.area<-trap.info$total.area | |
| trap.area<-trap.info$trap.area | |
| #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ | |
| #----1. Simulate animal locations. | |
| animal.locs<-animal.loc.func(xlims, ylims, total.area, density, dist.type) | |
| animals<-animal.locs$animals | |
| n.animals<-animal.locs$n.animals | |
| dist.type<-animal.locs$dist.type | |
| #---2. Calculate the distance between each trap/animal and the prob of capture | |
| #Problem here with allocation of memory...fixed with setting sd to 3? | |
| dist.prob<-dist.prob.func(n.animals, animals, n.traps, traps, g0, sigma,3) | |
| prob.xy<-dist.prob #$prob.xy, # dist2.xy<-dist.prob$dist2.xy | |
| #---3. Simulate the captures | |
| #Can also get a problem here with memory allocation | |
| sim.captures<-sim.captures.func(n.traps, n.nights, n.animals,max.catch, prob.xy) | |
| trap.animals<-sim.captures$trap.animals | |
| trap.catch<-sim.captures$trap.catch | |
| caught.traps<-sim.captures$caught.traps | |
| n.caught[s]<-sum(trap.animals) | |
| p.caught[s]<-sum(trap.animals)/n.animals | |
| t.catch<- t.catch+hist(rowSums(trap.catch),hist.brks, plot=FALSE)$counts | |
| #Store histogram of trap catches (after n.nights) | |
| # t.catch[d,]<- t.catch[d,]+hist(rowSums(trap.catch),hist.brks, plot=FALSE)$counts | |
| si<-((d-1)*reps)+s | |
| setWinProgressBar(pb, si, title=paste("Density = ", density, ", Rep =",s," (", round(si/(length(dens.vector)*reps)*100, 0),"% done)")) | |
| } | |
| #The number of animals that each trap catches... | |
| names(t.catch)<-c(0,1,2,3,4,5,6,7,8,9,10,11,12) | |
| t.catch<-t.catch/reps | |
| t.catch<-t.catch[1:(max.catch+1)] | |
| close(pb) | |
| #Summary | |
| data.summary<-as.data.frame(matrix(NA,2,2)) | |
| rownames(data.summary)<-c("Animals caught", "Traps with captures") | |
| colnames(data.summary)<-c("Mean", "SD") | |
| data.summary[1,1]<-mean((p.caught))#Proportion of animals caught (number would be n.caught) | |
| data.summary[1,2]<-sd((p.caught))#Proportion of animals caught (number would be n.caught) | |
| # data.summary[2,1]<-sum(caught.traps)/n.traps#Prop of traps with captures | |
| # trap.catch is traps*nights but number of animals per trap | |
| # caught,traps is whether the trap caught anything during the study... | |
| # data.summary<-n.caught | |
| #Return the values... | |
| return(list(animals=animals,traps=traps,xlims=xlims, ylims=ylims, trap.catch=trap.catch, n.animals=n.animals,n.traps=n.traps, total.area=total.area, trap.area=trap.area, g0=g0,sigma=sigma, max.catch=max.catch, n.nights=n.nights, data.summary=data.summary, n.caught=n.caught, p.caught=p.caught, t.catch=t.catch)) | |
| # return(list(animals=animals,traps=traps,xlims=xlims, ylims=ylims, n.animals=n.animals,n.traps=n.traps, total.area=total.area, trap.area=trap.area, g0=g0,sigma=sigma, max.catch=max.catch, n.nights=n.nights)) | |
| })#End of function.... | |
| #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ | |
| output$main_plot <- reactivePlot(function() { | |
| #~~~~~~~~~~~~~~~~~~~~ Do the plots ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ | |
| par(mar=c(4,5,5,1)) | |
| par(layout(matrix(c(1,1,2,1,1,3,1,1,4),3,3,byrow=TRUE))) | |
| traps<-data()$traps | |
| pch.trap<-ifelse(data()$max.catch==1,3,8) # | |
| animals<-data()$animals | |
| plot(traps[,2:3], pch=pch.trap, col="red", xlim=data()$xlims, ylim=data()$ylims, asp=1, cex=.6) | |
| points(animals$xloc, animals$yloc, col="blue", pch=20, cex=.8) | |
| mtext(paste("Number of animals = ",data()$n.animals),3, padj=-1, adj=0) | |
| mtext(paste("Number of traps = ",data()$n.traps),3, padj=-2.5, adj=0) | |
| mtext(paste("Total area = ",round(data()$total.area,0),"ha",", Trap area = ",round(data()$trap.area,0),"ha"),3, padj=-4, adj=0) | |
| #~~~~~~~~Plot distance vs prob from g0 and sigma~~~~~~~~~ | |
| #Calculate what the probability of capture will be: | |
| #i.e. P(Cap) = dnorm(distance,sd=sigma)*g0/dnorm(0,sigma) | |
| g0<-data()$g0 | |
| sigma<-data()$sigma | |
| plot.g0.sigma(g0,sigma) | |
| #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ | |
| # #---Plot of cumulative traps occupied by time | |
| plot(1:data()$n.nights,cumsum(colSums(data()$trap.catch))/(data()$n.traps*data()$max.catch),pch=20,xlim=c(0,data()$n.nights), ylim=c(0,1), xlab="Nights", ylab="Proportion of traps full traps" ) | |
| mtext("Cumulative traps full up (by time)",3) | |
| # hist((data()$n.caught)) #Number of animals caught | |
| hist((data()$p.caught)) #Propotion of animals caught | |
| }) | |
| # Generate a summary of the dataset | |
| output$summary1 <- reactivePrint(function() { | |
| (round(data()$data.summary,3)) | |
| }) | |
| output$summary2 <- reactivePrint(function() { | |
| data()$t.catch | |
| }) | |
| }) | |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| shinyUI(pageWithSidebar( | |
| # Application title | |
| headerPanel("Trap simulations v0.1"), | |
| # Sidebar with controls to select a dataset and specify the number | |
| # of observations to view | |
| sidebarPanel( | |
| wellPanel( | |
| h5("Trap Parameters"), | |
| numericInput("trap.space.x", "Spacing (x):", 100), | |
| numericInput("trap.space.y", "Spacing (y):", 50), | |
| numericInput("traps.in.x", "Number of columns:", 20), | |
| numericInput("traps.in.y", "Number of rows:", 20), | |
| numericInput("max.catch", "Catch per trap:", 1), | |
| numericInput("n.nights", "Number of nights:", 20) | |
| ), | |
| wellPanel( | |
| h5("Species Parameters"), | |
| numericInput("g0", "g0:", 0.05), | |
| numericInput("sigma", "Home range sigma (m):", 63), | |
| numericInput("density", "Density (per ha):", 1), | |
| radioButtons("dist.type", "Distribution of animals", | |
| list("Fixed" = "F", "Poisson" = "P")) | |
| ), | |
| wellPanel(h5("Simulation parameters"), | |
| numericInput("reps","Iterations:",10)), | |
| submitButton("Update") | |
| ), | |
| mainPanel( | |
| plotOutput(outputId = "main_plot", height = "800px"), | |
| h5("Probability of capture"), | |
| verbatimTextOutput("summary1"), | |
| h5("Animals captured per trap"), | |
| verbatimTextOutput("summary2") | |
| ) | |
| ) | |
| ) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment