Nest Migration Script 05-06-19

remove_pika_aipom_05-06-19.R

filename <- "remove_pika_aipom"

# remove taxa
removeTaxa <- c("Pikachu", "Aipom") 



####################################
###################################

currentNesting<-c(

# pre-halloween-2018 nesting list, as estimated by u/Clyphox and u/SnipeThemAll

 ###############
 # Gen 1
 ##############
 "Bulbasaur", "Charmander", "Squirtle",
	# "Caterpie", "Weedle", "Pidgey", "Rattata", "Spearow", "Ekans", #REMOVED - u/Clyphox
 "Pikachu",
	# "Sandshrew", "Nidoran-F", "Nidoran-M", #REMOVED - u/Clyphox
 "Clefairy", "Vulpix",
	# "Jigglypuff", "Zubat", #REMOVED - u/Clyphox
 "Oddish",
	# "Paras", "Venonat", "Diglett", "Meowth", #REMOVED - u/Clyphox
 "Psyduck",
	# "Mankey", #REMOVED - u/Clyphox
 "Growlithe", "Poliwag", "Abra", "Machop",
	# "Bellsprout", #REMOVED - u/Clyphox
 "Tentacool",
	# "Geodude", #REMOVED - u/Clyphox
 "Ponyta",
	# "Slowpoke", #REMOVED - u/Clyphox
 "Magnemite", "Doduo", "Seel", "Shellder", "Gastly", "Onix",
	# "Drowzee", "Krabby", #REMOVED - u/Clyphox
 "Voltorb", "Exeggcute", "Cubone", "Rhyhorn", "Horsea",
	# "Goldeen", "Staryu", #REMOVED - u/Clyphox
 "Scyther", "Jynx", "Electabuzz", "Magmar", "Pinsir",
 "Magikarp", "Eevee", 
 "Omanyte", "Kabuto",
   		 #"Dratini", #CD Pokemon
 #################
 ###############
 # Gen 2
 ##############
 "Chikorita", "Cyndaquil", "Totodile",
	# "Sentret", "Hoothoot", "Ledyba", "Spinarak", #REMOVED - u/Clyphox
 "Chinchou",
	# "Natu", #REMOVED - u/Clyphox
		 #"Mareep", # CD Pokemon
	# "Marill", #REMOVED - u/Clyphox
  "Sudowoodo", #### ADDED - u/Clyphox #####
	# "Hoppip", #REMOVED - u/Clyphox
 "Aipom",
	# "Sunkern", #REMOVED - u/Clyphox
 "Yanma",
	# "Wooper", #REMOVED - u/Clyphox
 "Misdreavus", "Wobbuffet", "Girafarig", "Dunsparce", "Snubbull",
 "Qwilfish", "Shuckle", "Sneasel", "Teddiursa",
	# "Slugma", #REMOVED - u/Clyphox
 "Swinub",
	# "Remoraid",  #REMOVED - u/Clyphox
 "Mantine", "Skarmory", #### ADDED - u/Clyphox #####
	# "Houndour", #REMOVED - u/Clyphox
 "Phanpy", "Stantler", #### ADDED - u/Clyphox #####
		  #"Larvitar", #CD Pokemon
 ###################
 ###############
 # Gen 3
 ##############
 "Treecko", "Torchic", "Mudkip", "Poochyena",
	# "Zigzagoon", "Wurmple", #REMOVED - u/Clyphox
 "Seedot",
	# "Tailow", #REMOVED - u/Clyphox
 "Wingull", "Surskit", "Shroomish",
	# "Whismur", #REMOVED - u/Clyphox
 "Makuhita", "Nosepass", "Skitty", 
	"Sableye", #### ADDED - u/Clyphox #####
	# "Aron", #REMOVED - u/Clyphox
 "Meditite", "Electrike", 
	"Plusle", "Minun", #### ADDED - u/Clyphox #####
	# "Gulpin", #REMOVED - u/Clyphox
 "Carvanha", "Wailmer", "Numel", "Spoink", "Swablu",
	# "Barboach", #REMOVED - u/Clyphox
 "Corphish", "Baltoy",
	 "Lileep", "Anorith", #### ADDED - u/Clyphox #####
 "Shuppet", "Duskull",
	# "Spheal", #REMOVED - u/Clyphox
 "Luvdisc",
		  #"Beldum", #CD Pokemon
 ##################
 ###############
 # Gen 4
 ##############
 "Turtwig", "Chimchar", "Piplup", "Bidoof"
)

newNesting <- currentNesting
whichRemove <- sapply(newNesting,
	function(x) all(x != removeTaxa))
newNesting <- newNesting[whichRemove ]

length(currentNesting)
length(newNesting)



#########################################

# for nest migrations - source functions, modify vectors of lists of nesting species, then run
# use kable to produce markdown table for pasting to Silph Road

#approximate ranges for each species in the nesting lists
getRangesForNestingSpecies<-function(nestingSpecies, arbMax=256){
	nNesting<-length(nestingSpecies)
	intLength<-arbMax/nNesting
	lowBound<-((1:nNesting)-1)*intLength
	upperBound<-lowBound+intLength
	res<-cbind(lowBound,upperBound)
	rownames(res)<-nestingSpecies
	return(res)
	}

# find new species for each old nesting seed range
findNew<-function(oldRange,newRanges){
	newRanges<-round(newRanges,5)
	oldLow<-round(oldRange[1],5)
	oldHigh<-round(oldRange[2],5)
	#
	newLowRow<-which(newRanges[,1]<=oldLow & newRanges[,2]>oldLow)
	newHighRow<-which(newRanges[,1]<oldHigh & newRanges[,2]>=oldHigh)
	#
	#print(c(newLowRow,newHighRow))
	newRows<-newLowRow:newHighRow
	newNestNames<-rownames(newRanges)[newRows]
	# calculate proportion of overlap (prob of shift)
	propOverlap<-sapply(newRows,function(x)
		 getPropOverlap(X=newRanges[x,],Y=oldRange))
	# paste them together
	res<-paste0(newNestNames," (",round(propOverlap*100,1),"%)")
	return(res)
	}
	
getPropOverlap<-function(X,Y){
	# calculate how much range X overlaps with range Y
	overlapRange<-c(max(X[1],Y[1]),min(X[2],Y[2]))
	overlap<-overlapRange[2]-overlapRange[1]
	# divide by length of Y
	propOverlap<-overlap/(Y[2]-Y[1])
	return(propOverlap)
	}

# master function
predictNestingShifts<-function(oldNestingList,newNestingList){
	oldRanges<-getRangesForNestingSpecies(oldNestingList)
	newRanges<-getRangesForNestingSpecies(newNestingList)
	#
	nestShifts<-apply(oldRanges,1,findNew,newRanges=newRanges)
	if(length(dim(nestShifts))!=2){
		# need to add NAs when there are uneven number of possible species
			# for different nests to potentially migrate to
		maxNSlength<-max(sapply(nestShifts,length))
		nestShifts<-sapply(nestShifts,function(x){
			if(length(x)<maxNSlength){
				res<-c(x,rep(NA,maxNSlength-length(x)))
			}else{
				res<-x
				}
			return(res)
			})
		}
	nestShifts<-t(nestShifts)
	#print(nestShifts)
	nestShifts<-cbind(oldNestingSpecies=rownames(nestShifts),nestShifts)
	colnames(nestShifts)[2:ncol(nestShifts)]<-paste0("NewNestingSpecies",
		1:(ncol(nestShifts)-1))
	rownames(nestShifts)<-NULL
	return(nestShifts)
	}

##############################
##########################################

nestShifts<-predictNestingShifts(
	oldNestingList=currentNesting,
	newNestingList=newNesting)

nestShifts


# get path
filename <- paste0("nesting_species_predictions_",filename)
filename <- paste0(filename,"_",format(Sys.time(), "%m-%d-%y"))
path<-"D:\\dave\\fun\\Pokemon\\analyses of pokemon stuff\\pokemon nest modelling\\"
filename<-paste0(path,filename)
# use knitr to convert the table of nest shifts to markdown
x<-knitr::kable(nestShifts)
cat(x,file=paste0(filename, ".txt"), sep = "\n")
rmarkdown::render(paste0(filename, ".txt"))
shell.exec(paste0(filename, ".txt"))
shell.exec(paste0(filename, ".html"))