# Analysing Flower shape variation on SPUR and PETAL datasets separate

# Script written by Sara van de Kerke (sara.vandekerke@wur.nl)
# Written December 2018, updated December 2018
# Accompanying van de Kerke et al. 

library(geomorph) # using version ...

### set working directory 
setwd("D:/PhD/OneDrive - WageningenUR/morphometrics_chapter_2-3/Rscripten")



###########################################
############### Data import ###############


load("pelargonium_SPUR_rotated.Rda")
load("pelargonium_PETAL_rotated.Rda")

classifier.variables.PETAL.matching <- read.csv("Variables_PETAL_overlap.csv")
classifier.variables.SPUR.matching <- read.csv("variables_SPUR_overlap.csv", header=TRUE)




###########################################
############### Resize data ###############


### We put the size component back into the coordinate data. For this, we use the Csize extracted in the GPA analysis.

### SPUR

pelargonium.SPUR.rotated.initial <- 
  simplify2array(lapply(1:dim(pelargonium.SPUR.rotated$Coords.raw)[3], 
                        function(j) pelargonium.SPUR.rotated$Coords.GPA$coords[,,j] * 
                          pelargonium.SPUR.rotated$Coords.GPA$Csize[[j]]) )                                        


dimnames(pelargonium.SPUR.rotated.initial)[[3]] <- 
  dimnames(pelargonium.SPUR.rotated$Coords.raw)[[3]]


plotAllSpecimens(pelargonium.SPUR.rotated$Coords.GPA$coords)
plotAllSpecimens(pelargonium.SPUR.rotated.initial)

### PETAL

pelargonium.PETAL.rotated.initial <- 
  simplify2array(lapply(1:dim(pelargonium.PETAL.rotated$Coords.raw)[3], 
                        function(j) pelargonium.PETAL.rotated$Coords.GPA$coords[,,j] * 
                          pelargonium.PETAL.rotated$Coords.GPA$Csize[[j]]) )                                        


dimnames(pelargonium.PETAL.rotated.initial)[[3]] <- 
  dimnames(pelargonium.PETAL.rotated$Coords.raw)[[3]]


plotAllSpecimens(pelargonium.PETAL.rotated$Coords.GPA$coords)
plotAllSpecimens(pelargonium.PETAL.rotated.initial)


#plot(pelargonium.PETAL.rotated.initial[,,100])



###################################################################################
############### Bootstrapping & Combine 2 2D matrices in 1 3D matrix###############
###################################################################################

### First, we select the data that is present in both datasets (complete overlap)

petal<-classifier.variables.PETAL.matching
spur<-classifier.variables.SPUR.matching


### Then, we check whether there are any mismatches between the variable data and the coordinate data

petal$ID[petal$ID %in% dimnames(pelargonium.PETAL.rotated.initial)[[3]]==F]
spur$ID[spur$ID %in% dimnames(pelargonium.SPUR.rotated.initial)[[3]]==F]

### We make empty vectors to store the variable data in. 

individuals<-numeric(0)
nr.spur<-numeric(0)
nr.petal<-numeric(0)
sumnr.spur<-0
sumnr.petal<-0

### nr.spur holds the counted number of individuals in the original variable dataset
### sumnr.spur gives the rownumber each next species starts at in the variable dataset

for (i in 1:68) {
  nr.s<- sum(spur[,3]==i)
  nr.p<- sum(petal[,3]==i)
  
  nr.spur<- c(nr.spur,nr.s)
  nr.petal<-c(nr.petal, nr.p)
  sumnr.spur<-c(sumnr.spur, sum(nr.spur))
  sumnr.petal<-c(sumnr.petal, sum(nr.petal))
}

nr.spur
nr.petal

sumnr.spur
sumnr.petal

### Start of the bootstrapping loop.

### First, we make a matrix to store all the bootstrap replicates in.

#outsboot<-matrix(ncol=3*100,nrow = 6*56*100)
#outsboot <- data.frame(matrix(nrow = 16800, ncol = 4))



for (j in 1:20) {
  
  ### We make empty vectors to store the bootstrap sampling data in.  
  
  bootspur.number <-numeric(0)
  bootpetal.number<-numeric(0)
  
  ### Bootstrap 6 specimens for each of the 68 species we have in the data.
  
  for (i in 1:68) {
    bootspur.number <-c(bootspur.number,sumnr.spur[i]+sample(nr.spur[i],6, replace = TRUE))
    bootpetal.number <-c(bootpetal.number,sumnr.petal[i]+sample(nr.petal[i],6, replace = TRUE))
  }
  
  ### Check content of bootspur and bootpetal
  
  bootspur.number
  bootpetal.number
  
  
  
  ### Make sure the output of bootspur is in the same order as the variable data
  bootspur<-as.vector(spur$ID)[bootspur.number]
  bootpetal<-as.vector(petal$ID)[bootpetal.number]
  
  ### Check
  head(bootspur)
  head(bootpetal)
  
  ### Combine bootspur and bootpetal (and write to a file)
  
  bbbb<-cbind(bootspur,bootpetal)
  
  #write.csv(bbbb, file="bbbb.txt", append = TRUE ,sep = '\t') 
  
  
  ### Connect numbers of bootspur/bootpetal to the individual numbers of the variable dataset and write to a file
  #boots <- pelargonium.side.initial.GPA$coords[,,bootspur]
  #bootp <- pelargonium.front.initial.GPA$coords[,,bootpetal]
  
  ### Connect numbers of bootspur/bootpetal to the individual numbers of the variable dataset and write to a file
  boots <- pelargonium.SPUR.rotated.initial[,,bootspur]
  bootp <- pelargonium.PETAL.rotated.initial[,,bootpetal]
  
  
  write.csv(boots, file="bootstrap.csv", append = TRUE ,sep = '\t')
  
  
  
  
  ### Combine the coordinates of the 2 2D matrices in 1 3D matrix. For this, we define 2 anchor points and add a third
  ### coordinate to each coordinate system
  
  
  for (i in 1:length(bootspur)){
    #  for (i in 235:264){
    ### First, we extract the individual names
    
    a <- dimnames(boots)[[3]][i]
    b <- dimnames(bootp)[[3]][i]
    
    ### For the spur dataset, the x1-x1[,anchor] should be the new z-coordinate and y1-y1[,anchor] is the 
    ### difference from the reference to y2 (this dataset is 'flipped. We do this for each of the combinations made by the bootstrap 
    ### sampling above.The x coordinate is 0 for all (object is flat)  and is added. In this case, the anchor point is the first landmark.  
    
    ### Extract the coordinate data from the first/spur bootstrap dataset and select the first (for x coordinate) and second (for
    ### y coordinate) column and store in object. 
    
    xy1 <- boots[,,i]
    
    x1 <- xy1[,1]
    head(x1)
    
    y1 <- xy1[,2]
    head(y1)  
    
    
    ### Put new spur x ( = 0), y, and z coordinates in one matrix. Var1 now contains the coördinate system correct for after 
    ### 'gluing' the data.
    
    var1<-cbind(0,(y1-y1[1]),(x1-x1[1]))
    
    
    ### For the petal dataset, the x[anchor] is defined as the average between the landmarks 22 and 23.
    
    ### Extract the coordinate data from the second/petal bootstrap dataset and select the first (for x coordinate) and second (for
    ### y coordinate) column and store in object. We add a z coördinate which is 0 for all (object is flat).
    
    xy2 <- bootp[,,i]
    
    x2 <- xy2[,1]
    head(x2)
    
    y2 <- xy2[,2]
    head(y2)  
    
    z2 <- rep(0,length(y2))
    
    ### We rotate this set of landmarks (the x2 and y2 coordinates using) a transformation matrix before "gluing" the data.
    
    rotmat<- matrix(nrow=2,ncol=2,c(0,1,-1,0))
    
    rotx2y2<-t(rotmat%*%rbind(x2,y2) )
    
    refx2<- (rotx2y2[22,1] + rotx2y2[23,1])/2
    refy2<- (rotx2y2[22,2] + rotx2y2[23,2])/2
    
    #refx2<- (x2[22] + x2[23])/2
    #refy2<- (y2[22] + y2[23])/2
    
    ### Put new petal x, y and z coördinates in one matrix. var2 now contains the coördinate system correct for after 'gluing' 
    ### the data
    
    var2<-cbind(rotx2y2,z2)
    
    #var2<-cbind(x2,y2,z2)
    
    ### Now the anchor point of the petal dataset (defined in refx2) is used to define the new x coordinate of the spur dataset. 
    ### This x was previously defined to be zero, but has to moved to the reference point of the other system.
    ### The same goes for the y coördinate, which is defined as its original value in var1 plus the y-value defined in refy2.
    
    var1[,1]<-refx2
    var1[,2]<-var1[,2]+refy2
    
    
    ### Then we combine var1 and var2 and write is to a file with a make-up that is suitable for further use in Geomorph. 
    
    var<-rbind(var1,var2)
    
    
    write("LM=114", file="68species3.txt", append = TRUE)
    write(t(var),file="68species3.txt", ncolumns = 3, sep = '\t', append = TRUE)
    write(paste("IMAGE="), "68species3.txt", append = TRUE)
    write(paste("ID=", a, b, j, i, sep = "_"), "68species3.txt", append = TRUE)
    
    
    
    
  }
  ### Connect numbers of bootspur to species names of the variable dataset and write to a file
  individuals.raw <- as.vector(spur$Species[bootspur.number])
  individuals <- cbind(j, i, bootspur, bootpetal, individuals.raw)
  write.table(individuals, file="individuals68species3.csv", append=TRUE, sep = ",", row.names = TRUE, col.names = FALSE)
  
  
  #outsboot[i] <- var
}


#write.csv(t(outsboot), file="outsboot.txt", append = TRUE)




### import TPS file with landmark data of the COMBINED dataset
pelargonium.combined.original <- 
  readland.tps("68species3.txt", specID = "ID", readcurves = TRUE, warnmsg = TRUE)



### Import .csv file containing individual species info

classifier.individuals <- read.csv("individuals68species3.csv", header=FALSE)
classifier.individuals$numbers <- paste(classifier.individuals$V1, 
                                        classifier.individuals$V2, sep = "_")
classifier.individuals$species <- classifier.individuals$V6

classifier.individuals$combination <- 
  paste(classifier.individuals$V4,classifier.individuals$V5, sep = "_")


######################################################
############### Combine data in a list ############### 

pelargonium.combined.initial <- list(pelargonium.combined.original, 
                                     classifier.individuals$species,
                                     classifier.individuals$combination)
names(pelargonium.combined.initial) <- c("Coords",   "Species", "Combination")


save(pelargonium.combined.initial, file = "Pelargonium_combined_initial.Rda")